Predominant intragenic methylation is associated with gene expression characteristics in a bivalve mollusc

Introduction

Epigenetic marks such as DNA methylation are important for genome regulation (Bell & Felsenfeld, 2000; Li, Bestor & Jaenisch, 1992; Hsieh, 1994). DNA methylation has been well-studied in mammals and plants where it has been shown to play important roles in temporal and spatial regulation of gene expression (Okano et al., 1999; Zhang et al., 2006), reduction of transcriptional noise (Bird, 1995), and genome stabilization (Wolffe & Matzke, 1999). However, the distribution and context of DNA methylation varies greatly among phylogenetic groups indicating that these functional roles may not be conserved (Colot & Rossignol, 1999).

In contrast to the heavily methylated vertebrate genomes, invertebrate genomes exhibit a relatively low level of methylation that is primarily confined to gene bodies (Zemach et al., 2010). It is perhaps because of these differences that DNA methylation has remained largely understudied in invertebrates. Recently, however, there has been a renewed interest in invertebrate DNA methylation patterns as it is now being recognized that invertebrates are exceptional models to study functions and evolutionary roles of gene body methylation. Furthermore, it has been shown that DNA methylation mediates phenotypes in response to environmental cues in the invertebrate Apis mellifera (Kucharski et al., 2008; Lyko et al., 2010), indicating an important role in integrating environmental signals. Thus, understanding both the distribution and function of DNA methylation in diverse invertebrate lineages will contribute significantly to our understanding of the evolution of genome regulation and environmental physiology.

The focus of the work presented here is to explore the role of DNA methylation in bivalve molluscs. The presence of DNA methylation has been confirmed in several bivalve species (Wang et al., 2008; Petrovic et al., 2009; Gavery & Roberts, 2010). A majority of the research on DNA methylation in molluscs has focused on the Pacific oyster (Crassostrea gigas), an economically and ecologically important species. Previous studies in the Pacific oyster identified a relationship between gene function and methylation pattern. Specifically, it was shown that genes with housekeeping functions are more methylated than genes involved in inducible functions (i.e., genes involved in response to environment, embryonic development or tissue-specific functions) (Gavery & Roberts, 2010; Roberts & Gavery, 2012). More recently, Riviere et al. (2013) determined that DNA methylation plays a critical role in development as indicated by differential methylation patterns throughout embryogenesis. This was further supported by their observation that 5-aza-cytidine, a potent demethylating agent, significantly disrupts embryonic development (Riviere et al., 2013).

These recent studies on DNA methylation in oysters provide important foundational information on DNA methylation in bivalves. However, previous studies were not able to provide fine scale resolution of DNA methylation patterns, nor examine the relationship with gene expression at the genome-wide level. Here, we provide the first high resolution methylome of a mollusc and examine this in relationship to gene expression data to get a better understanding of the role of DNA methylation in invertebrates.

Methods

Results

DNA methylation mapping

Bisulfite treated DNA sequence reads (139,728,554 total reads; 36 bp) are available in the NCBI Short Read Archive under the accession number SRX32737. A total of 120,734,949 reads (86%) mapped to the C. gigas genome. Fifty-six percent of the 164,873,219 cytosines in the C. gigas genome had at minimum of 1× coverage. Of the 9,978,551 CpG dinucleotides in the genome, 2.6 million (26%) had ≥5× coverage. The distribution of methylation ratios found at CpG dinucleotides ranged between 0.0 and 1.0, but a majority of the loci were either heavily methylated or unmethylated. Specifically, 55% (1,453,752) were methylated (i.e., ≥ 0.50) and another 28% were unmethylated (i.e., = 0.0) (Fig. 1). Genome feature track files (i.e., GFF) representing (1) all CpG dinucleotides and (2) methylated CpG dinucleotides (>50%) for this dataset were developed and are available (Gavery & Roberts, 2013).

Methylation landscape of genomic features

Methylated CpG dinucleotides, defined as having a methylation ratio of 0.5 or greater, were located predominantly in intragenic regions (exons and introns), but were also present in putative promoters (defined as 1 kb upstream of TSS), transposable elements and unannotated intergenic regions. The distribution of methylated CpG across various genomic regions is significantly different than what would be expected if the methylation were distributed randomly throughout the CpG dinucleotides in the genome (X² = 513, 194.1, df = 4, p < 0.0001). Specifically, DNA methylation appears to be overrepresented in intragenic regions (64% of methylated CpG in combined exons and introns) when compared to the proportion of all CpG in the genome (38%) (Fig. 2). When methylation was examined on a per gene basis a strong positive correlation (R²: 0.86) was observed between exonic and intronic methylation within a gene. Additionally, a strong correlation was observed between the gene methylation measured via high-throughput bisulfite sequencing and the predicted methylation ratio based on the CpG observed to expected ratio (CpGo/e) (Spearman rho: −0.616, p-value: <1 × 10⁻⁴).

Gene expression & DNA methylation

After quality trimming, 45,751,574 RNA-seq reads mapped to the genome. Raw reads are available in the NCBI Short Read Archive under the accession number SRX367081.

The relationship between the proportion of methylation in a gene and overall transcript abundance in gill tissue was characterized (Fig. 3). In general, transcription abundance increases significantly with increasing DNA methylation until the 40th percentile after which it remains relatively stable until the 100th percentile when methylation significantly decreases.

The first two principal components (PC) of the PCA of gene attributes were significant and explained 76.4% of the variation among the genes. This variation was being driven by multiple factors, including DNA methylation (Fig. 4 and Table 1). The only variable that did not contribute significantly to the first two principal components was mean transcript abundance (correlation loading 0.2). The first PC, which explained 50.2% of the variation was loaded heavily by number of CpG dinucleotides, the length of the mRNA and the number of exons. The second PC, which explained 26.1% of the variation was loaded heavily by the %CV of gene expression among tissues and the methylation ratio. DNA methylation is negatively correlated with transcript variance between tissues (%CV) and relatively uncorrelated with attributes such as gene length.

Table 1:

Summary of PCA for gene attributes.

Principal component	%variance explained	Significance value	Significant variable loadings
PC1	50.2	<0.001	Number of CpG	0.9
			Length mRNA	0.9
			Number of exons	0.8
PC2	26.1	<0.0001	Expression %CV	−0.6
			Methylation ratio	0.6

DOI: 10.7717/peerj.215/table-1

Discussion

Here we have used methylation enriched high-throughput bisulfite sequencing in conjunction with genomic feature annotation and transcriptomic data to gain a better understanding of the role of DNA methylation in oysters. This work not only provides new information on DNA methylation in invertebrates but also provides a framework for characterizing DNA methylation in other taxa.

The reduced representation approach was selected to obtain a higher coverage of methylated regions. In addition, since methylation was likely to occur in gene bodies (Zemach et al., 2010), and because transcriptomic data was the primary genomic resource for C. gigas at the time of sequencing (the genome was released soon after), it was expected that methylation enrichment would significantly limit the proportion of unmappable reads. Quantitative methylation data were obtained for both methylated CpG as well as unmethylated CpG that were either interspersed with or flanking these more heavily methylated regions. Therefore, methylation enriched bisulfite-sequencing was effective in generating a comprehensive invertebrate methylome.

One of our primary findings was the overall level of genome methylation in the oyster. Here we found that 15% of CpG dinucleotides (2% of total cytosines) are methylated in gill tissue. This degree of methylation is much lower than the global methylation patterns seen in mammals where 70%–80% of CpGs are methylated (Bird & Taggart, 1980), but still higher than what has been reported in other invertebrates. For instance, only 0.8% of the CpGs are methylated in the brain of A. mellifera (Lyko et al., 2010) and between 0 and 8% of CpGs are methylated in the nematode, Trichinella spiralis, depending on the life stage (Gao et al., 2012). Although methylation in C. gigas is relatively high for an invertebrate, it is not outside the range of what has been reported in other species by liquid chromatography-mass spectrometry analysis. For example, similar to the oyster, 2% of total cytosines are methylated in the mollusc Biomphalaria glabrata (Fneich et al., 2013). It should be noted that methylation in oysters does likely vary in both a temporal and possibly tissue specific manner, as clearly indicated by Riviere et al. (2013) by characterizing differences in total methylation during development. In addition, because the sample represents a pool of multiple individuals, it cannot be determined whether the variation in methylation at a particular locus represents hemimethylation or differential methylation between individuals. In general, the bimodal pattern observed (Fig. 1) indicates that a CpG locus is either heavily methylated or unmethylated, but future work sequencing individual oysters would provide valuable information regarding individual epigenetic variation in oysters.

This work also provided the first direct evidence in oysters that DNA methylation is prominent in gene bodies (see Fig. 2) and these data are well correlated with previous investigations using an in silico approach (i.e., CpGo/e) to predict methylation in C. gigas (Gavery & Roberts, 2010).The predominance of gene body methylation is consistent with what has been described in other invertebrates (e.g., Suzuki et al., 2007; Zemach et al., 2010; Lyko et al., 2010) and there is increasing evidence that gene body methylation is the ancestral pattern (Lechner et al., 2013). The function of gene body methylation remains unclear, but studies indicate possible active roles in preventing spurious transcription (Bird, 1995; Huh et al., 2013) and regulating alternative splicing (Maunakea et al., 2010; Shukla et al., 2011; Foret et al., 2012), as well as a more passive role for methylation as a byproduct of an open chromatin state (Jjingo et al., 2013). Given the nature of the study design, we are not able to directly test the hypothesis that DNA methylation contributes to spurious transcription or the regulation of alternative isoforms in C.gigas. However, genomic feature tracks have been developed and published (Gavery & Roberts, 2013) so that genome wide methylation can be easily visualized with respect to gene expression patterns (exon-specific RPKM).

Exons are the preferential target for gene body methylation for most species (Feng et al., 2010), and methylation is enriched in exons of the oyster. However, there is also a relatively large amount of intronic methylation in oysters when compared to other invertebrate species. For example, DNA methylation occurs almost exclusively in exons in the honey bee A. mellifera (Lyko et al., 2010). Genome-wide methylation studies in other invertebrate species also report very low levels of intronic methylation relative to other genomic regions (e.g., Gao et al., 2012; Bonasio et al., 2012). Similarly, in plants, methylation is preferably targeted to exons; however, it has been reported that in globally methylated mammalian genomes gene body methylation is not biased toward exons (Feng et al., 2010), although exon/intron boundaries can be marked by differences in DNA methylation (Sati et al., 2012). It appears that bivalves may be unique among the invertebrates examined in terms of the degree of methylation in introns. Intronic methylation has been implicated to be involved in gene regulation through the expression of alternative isoforms of genes in other species (e.g., Maunakea et al., 2010; Foret et al., 2012). Variation in methylation patterns between taxa may indicate that additional model invertebrates are needed to study the function of these epigenetic marks.

The distribution of DNA methylation in the C. gigas genome is consistent with the fractionated or ‘mosaic’ pattern of methylation previously described in invertebrates (Tweedie et al., 1999; Simmen et al., 1999). In oysters, as in other invertebrates, the methylated fraction tends to consist of gene bodies, while other genomic regions exhibit less methylation (Fig. 2). Interestingly, transposable elements (TE) show little methylation in oyster gill tissue. This is in contrast to vertebrate genomes where TE are heavily methylated and function to suppress their activity (Yoder, Walsh & Bestor, 1997). While there is no general consensus regarding the extent of TE methylation across invertebrate taxa, the pattern of sparse TE methylation observed in oysters is similar to what has been described in other invertebrate species (Simmen et al., 1999; Feng et al., 2010; Zemach et al., 2010).

Intragenic DNA methylation is positively correlated with gene expression in C. gigas with moderately and highly expressed genes showing the highest degree of methylation (Fig. 3). This relationship is similar to what has been reported for other invertebrate species (Zemach et al., 2010). Interestingly, Riviere et al. (2013) reported a negative relationship between DNA methylation and expression of certain homeobox (hox) genes during embryonic development in C. gigas. The authors hypothesized that the apparent suppression of hox expression by DNA methylation may be due to repression by DNA methylation proximal to the transcription start site in these genes. Although the results reported here and those of Riviere et al. may appear contradictory, it is possible that depending on the context of the methylation (i.e., whether gene body or promoter methylation) it may play either a repressive or expressive role. This is referred to as the DNA methylation paradox (Jones, 1999) and is observed in a wide range of taxa.

We used an ordination approach to explore genomic attributes or groups of attributes that predictably co-occur with methylated genes in the C. gigas genome. Because multiple factors may be linked with methylation (either through causative or correlative associations), this approach allowed us to identify relationships between multiple variables. The most interesting finding from the PCA analysis is that the amount of methylation in a gene is related to the variance in expression between tissues. Genes that show the least variation in expression between tissues have higher DNA methylation levels than those exhibiting a tissue-specific expression profile (i.e., high %CV between tissues). This observation provides corroboration for previous reports based on in silico analyses in oysters showing that housekeeping genes have the highest amount of methylation in C. gigas (Gavery & Roberts, 2010). Housekeeping genes perform functions required by all cell types, therefore it’s expected that their expression patterns would show low variation across tissues. The results of this study are consistent with the expectation that genes with low expression variation across tissues show a high degree of methylation relative to genes with a more tissue-specific expression pattern. Again, this study supports previous findings (Gavery & Roberts, 2010; Roberts & Gavery, 2012) that heavily methylated genes are enriched in housekeeping functions, which are essential for cellular function. One theory is that the lack of methylation in genes with tissue-specific expression can contribute to phenotypic plasticity by allowing more transcriptional opportunities through processes such as allowing access to alternative TSS, facilitating exon skipping or other alternative splicing mechanisms and allowing for increased sequencing variation (Roberts & Gavery, 2012).

Conclusions

Through the current effort, quantitative methylation data were obtained for over 2.5 million CpG dinucleotides throughout the genome of Crassostrea gigas. These data represent the first high resolution methylome in any mollusc and the analytical approaches provide a framework for DNA methylation characterization in other species. In addition, the dataset developed here will be beneficial for phylogenetic analysis of DNA methylation in invertebrates, which will be more robust with the addition of a lophotrochozoan species. The results of this study highlight similarities in epigenetic profiles of invertebrates such as a predominance of gene body methylation and a positive relationship between intragenic methylation and gene expression. In addition, they highlight interesting differences between invertebrate epigenomes including a higher level of intronic methylation in bivalves than what has been reported, for example, in insects. Although the functional role of DNA methylation in bivalves remains elusive, two scenarios could explain our findings. One possibility is DNA methylation in gene bodies is a byproduct of transcription resulting from an open chromatin state, as proposed by Jjingo et al. (2013). Thus the methylation patterns are influenced by transcriptional activity. The second scenario is DNA methylation is involved in regulating gene activity in bivalves. If in fact DNA methylation does influence transcription, the regulatory role is likely very complex. For instance, DNA methylation could have both a have direct regulatory effect on certain genes as proposed by Riviere et al. (2013), as well as facilitating expanded transcriptional opportunities in other cases. Future studies will certainly be challenging given the dynamic nature of DNA methylation, but will hopefully help better delineate if DNA methylation plays a functional role in regulating genome activity in bivalves and what that role might be.