Genome-wide methylation is modified by caloric restriction in Daphnia magna

Genome-wide percentage of DNA methylation is low in arthropods

DNA methylation, a reversible covalent modification that regulates gene expression, is the best-studied epigenetic mechanism. DNA methylation of cytosines occurs when DNA methyltransferase enzymes (DNMTs) transfer a methyl group onto cytosine [26] to create a 5-methylcytosine. This can occur at a cytosine immediately followed by guanine (CpG site), or as occurs in plants at a CHH or CHG site; where H can be any of A, T or C. In mammals around 70% of CpGs are methylated [27], whereas the genome-wide rate of CpG methylation in invertebrate species sampled to date is much lower, from 0% in Drosophila to 15% in the oyster Crassostrea gigas [18, 28]. The crustaceans Daphnia magna and Daphnia pulex (Arthropoda: Crustacea) have genomic CpG methylation of 0.52–.74 and 0.7% respectively [25, 29].

In mammals, certain insect orders [18], and the crustacean Daphnia, there are three DNMT enzymes. The gene products of DNMT3a and DNMT3b establish methylation de novo, while DNMT1 maintains methylation, and DNMT2 has no known role in DNA methylation. DNMT3 has been lost the most among insects and is present in only four orders [18]. Despite the apparent key role of this gene, the genome of the silkworm (Bombyx mori) which lacks DNMT3 along with other Lepidoptera shows evidence for DNA methylation [8, 18].

CpG methylation can increase or decrease gene expression dependent on the location of the methylation. In mammalian promoter regions, which can be rich in CpGs and are known as CpG islands, it represses expression of the gene. Further to this, many mammalian CpG islands are also enriched for permissive chromatin modification, which condenses the structure of chromatin and further prevents transcription. In contrast, methylation of gene bodies in mammals leads to an increase in expression of the effected gene. Invertebrates have few CpG Islands, and methylation predominantly occurs in gene bodies, and is enriched in exonic sequence to form a mosaic-like pattern genome-wide [30–33]. Daphnia magna is no exception to this [25, 29]. The role of gene-body methylation in invertebrates is not as clear as in mammals, but it is implicated in facilitating stable transcription of housekeeping genes and mediating alternative splicing [33–35]. This link to alternative splicing has been identified in D. magna [22, 25]. In Daphnia magna and pulex there is a also positive correlation between highly methylated genes and gene expression in the germline [25]. More specifically, methylation of exons 2–4 is significantly associated with this effect, as it is in humans [25, 36]. In honeybees, differential methylation of CpGs is potentially associated with differential gene expression through ‘priming’ of genes in response to intruders [20]. Interestingly, in silkworms there is no correlation between methylation in promoters and gene expression [8], suggesting invertebrates and vertebrates differ in their usage of CpG methylation.

Identifying differentially methylated regions in organisms with low-levels of methylation, like most arthropods, lacks a standardised approach at present, but DMRs are detectable as a growing body of research shows. This may reflect the nature of invertebrate methylation, with clusters of methylated CpGs interspersed by long, fallow, regions. In the above honeybee study, the authors identified differentially methylated regions between the brains of aggressive and control bees [20]. In Daphnia, transgenerational inheritance of CpG modifications has been shown in response to salinity, changes in which impact life-history traits, toxic cyanobacteria, and radioactivity via a whole genome bisulphite sequencing approach [22–24].

Caloric restriction and DNA methylation

The relationship between diet and CpG methylation, and subsequent impact on ageing and health, is well established in mammals. Indeed, DNA methylation may be a predictor of biological age [37, 38], and a reduction in expression levels of DNMT enzymes leading to a global loss of genomic methylation is associated with ageing humans [39]. Although promoter regions of certain genes become hypermethylated with age [40]. This is in contrast to honey bees in which an increase in de novo methylation activity through DNMT3 is observed [19]. Inhibition of DNMT activity increased lifespans of treated bees, but in a CR independent manner [19]. Specific examples of a diet by methylation interaction include the expression of DNMTs which have elevated expression in response to CR in cancer cells, counteracting the global hypomethylation [41] observed during ageing. CR also causes a reduction in lipid metabolism gene expression by DNA methylation of gene bodies in mouse livers [42]. As a result, older mice undergoing CR were protected from fatty degeneration, visceral fat accumulation, and hepatic insulin resistance compared to controls. In rats and monkeys short-term CR in older individuals ameliorates the effects of ageing with respect to disease markers, oxidative stress and damage, and increases the expression of longevity related genes [43, 44].

Caloric restriction and CpG methylation in Daphnia

Our work aims to determine if an experimentally controlled nutritional environment directs changes in methylation status in the crustacean Daphnia magna. Species of Daphnia are well characterised for their response to CR [45], including clear maternal effect phenotypes in D. magna [9, 10, 46–48]. However, increased lifespan is not a universal response among D. magna strains (T. Little, unpublished data) which can exhibit inter-clonal differences in lifespan and other age-related parameters such as age and size at first reproduction when sampled from the same habitat [49]. This is also true for D. pulex where increased lifespans are observed under CR [50–52] but not always [53, 54]. We performed whole genome bisulphite sequencing of CR and normally-fed (NF) replicates of D. magna of the the same Kaimes strain as prior studies demonstrating a maternal effect response to CR. We identify regions of differential methylation using an approach that successful in previous studies focused on differential modifications at CpG in response to environment-derived treatments [22–25] including salinity and exposure to the toxic cyanobacterium Microcystis aureginosa. D. magna show the classic response of longer life under CR and strong maternal effects; the offspring of calorie-restricted mothers being larger and more resistant to pathogens than their counterparts from better fed mothers. Provisioning of offspring, e.g. with carbohydrates, protein or fats, is one explanation for these maternal-effect phenotypes, and epigenetic processes, such as methylation, are also potentially key regulators in these plastic responses to fluctuating environments.

Daphnia have many attributes that make them favourable for epigenetic study. It has long been known that they undergo adult tissue regeneration unlike other invertebrates, akin to humans [54–56]. Furthermore, they reproduce clonally, which permits the study of epigenetic changes in the absence of genetic variation. Clonal replicates are equivalent to identical twin studies, but with an experimentally chosen number of -uplets. CpG-based methylation occurs in Daphnia, its genome encodes all three DNMT enzymes orthologous to mammalian DNMT enzymes, and a growing body of research is showing that global methylation patterns in Daphnia change in response to a variety of environmental factors [21–25, 29, 57, 58].

Methylated sites prediction

We performed whole genome bisulphite sequencing (WGBS) of calorie-restricted (CR) and normally-fed (NF) D. magna to generate data for identification of methylated C residues and exploration of the C methylation response to CR. Trimmed bisulphite-converted reads aligned to the D. magna strain 32 converted genome [48] using Bismark [59] exhibited lower mapping efficiencies than standard short-read alignments, as is typical of WGBS [60], with 20–32% of reads not aligning to the reference genome (Additional files 1 and 2). Of the aligned reads, 29–38% of reads were discarded as PCR duplicates, and 6–10% of the remainder contained predicted CHH or CHG methylated sites. These CHH and CHG site containing reads were removed from analyses, as D. magna contains negligible CHH/G methylation [29] meaning they are indicative of bisulphite unconverted reads. This resulted in average replicate read coverages of 8- to 12-fold. Global CpG methylation was 0.7–0.9% in all samples, and no difference in overall methylation levels was observed between NF and CR replicates paired t-test two-tailed p-value = 0.17). Removal of polymorphic CpG sites using variants predicted from the bisulphite-unconverted data had little effect on the total number of sites (Additional file 1). This was necessary as at a polymorphic C/T site the T allele can be miscalled as an unmethylated C causing an underestimation of the methylation profile at that site. After filtering, 99.2% of sites were retained per replicate, with average total sites across replicates going from 6.9 million to 6.85 million. Hierarchical clustering of replicates by methylation status demonstrated that mother had a stronger effect on global methylation status than nutritional treatment (Fig. 1). Thus, it was necessary to incorporate pairing information into our statistical test. We did this by performing paired t-tests of CpG sites in the Bioconductor R package “Bsseq” [61].

Differential methylation assessment

After filtering for coverage, 4.15 million CpG sites were retained for differential methylation testing. Bsseq testing of differential methylation revealed 8764 potential regions using t-statistic cut-offs of − 4.6, 4.6. This was filtered to 453 DMRs with at least three methylated CpGs and a mean difference greater than 0.1. These DMRS were further filtered after combining raw per region p-values into one p-value (false discovery rate correction < 0.05) in comb-p [62]. This resulted in 333 regions (Fig. 2. Volcano plot of FDR corrected p-values) adapting Herb et al. [20]. Of these 117 (35%) were hypermethylated in the CR group versus normal food, and 216 (65%) were hypomethylated in the CR group versus normal food. Only a fraction of these DMRs exhibited a mean difference of greater than 20% in either direction at 4% (5/117) for hyper- and 6% (12/216) for hypomethylated regions respectively. DMRs were 161 base-pairs (bp) and 193 bp long on average for CR hyper- and CR hypo-methylated regions, respectively and ranged from 18 to 602 bp in length. There were from three to 20 CpGs per cluster with an average of six. CR hypermethylated DMRs overlapped 165 gene predictions in the D. magna gene set, while CR hypomethylated DMRs overlapped 284. The majority of these overlapped exonic sequences for both hypo- and hypermethylated regions: 82% (234/384) for CR hypomethylated DMRs and 81% (133/165) for CR hypermethylated DMRs. Only 33 DMRs (10%) did not overlap a predicted gene body. The increase in number of DMR containing genes versus total DMRs reflects overlapping or redundant predictions in the D. magna genome annotation version 2.4 [63]. Because of this we used a non-redundant set of genes from the D. magna gene set as input to GO term enrichment, which resulted in 192 hypo- and 115 hypermethylated genes being considered for functional enrichment. Annotations for these genes are given in Additional files 3 and 4, a greater portion of genes are uncharacterised for hypo- than hypermethylated genes at 28 to 10% respectively. Two genes contain DMRs in both directions (Table 1), one encodes for a latent nuclear antigen and the other a heat shock cognate (Hsc-70) interacting protein. For both genes the exon (exon 2) closer to the predicted start of the gene is hypomethylated, and the 3′ exon is hypermethylated (exon 4 and 6).

Gene	Genomic Position	DMR	Length of DMR	Condition in CR	Number of sites	Corrected P-value	Exon	Strand
Dapma7bEVm000710	scaffold02569: 76199–79,126	76,729–77,106	377	hypo	16	0.0007575	2 of 8	+
Dapma7bEVm000710	scaffold02569: 76199–79,126	77,694–77,926	232	hyper	5	0.001717	6 of 8	+
Dapma7bEVm003762	scaffold01005: 527546–535,133	534,741–534,857	116	hypo	9	0.009821	2 of 6	–
Dapma7bEVm003762	scaffold01005: 527546–535,133	534,106–534,213	106	hyper	4	0.006839	4 of 6	–

Gene: gene name, Genomic Position: gene coordinates, DMR: coordinates of the differentially methylated region, Length of DMR: length in bp, Number of site: number of methylated CpGs in region, Corrected P-value: FDR corrected p-value for that DMR from comb-p analysis, Exon: which exon of the gene DMR occurs in followed by number of exons for that gene, Strand: which genomic strand the gene occurs on

GO term enrichment in methylated genes and DMRs

Global methylation and differentially methylated regions

The global D. magna CpG methylome of ~ 0.7%, consistent across replicates, is similar to that reported previously for different clones of D. magna (0.74 and 0.5% respectively) [25, 29]. Slight differences between studies may be related to D. magna strain differences, or differences in different data filtering methods [29]. In line with previously observed methylation patterns in Daphnia [22–25] and arthropods [30, 32, 64] in general, the majority of DMRs were identified in exonic regions of gene bodies. The magnitude of change in methylation of DMRs is mainly between 10 and 20% (Fig. 2), with some outliers. Indeed, the greatest difference between treatments was only 33% for a DMR hypomethylated under CR overlapping both exonic and intronic sequence in 40S ribosomal protein S26 (Additional files 3 and 4, genes with exons overlapping DMRs of greater than 20% mean difference). This could reflect our whole organism-based study in which large changes in specific tissues are moderated by no change the majority of tissues for that genomic region. These differences in mean methylation are in line with that for D. magna response to Microcystis stress, where differences up to 40% were observed [22]. However, the comparison is not direct as that study compared difference between methylation at gene bodies rather than DMRs. The extent to which changes in magnitude of DMRs or gene bodies is functionally relevant is currently unknown in Daphnia, and arthropods in general. It is possible that the low, but significant, differences in mean methylation > 10% have little functional effect, and as such requires validation through future studies. In Daphnia and other arthropods there is a correlation between DMRs and alternative splicing [22, 24, 25, 30], and is proposed to form an adaptive response to environmental pressures [22, 24, 25]. Increased gene expression is associated with high-levels of germline methylation in Daphnia [25], but the impact of environmentally-induced changes on gene expression is not yet understood beyond alternative splicing [22].

There are two genes that contain DMRs methylated in different directions (Table 1), in both cases the first exon (exon two in both genes), has the greater number of methylated CpGs versus the later exon. This is consistent with Kvist et al. [25] as exons two and three in D. magna are the most methylated exons genome-wide, with a decline in methylation in subsequent exons. It remains for future work to dissect if the differential methylation in these genes results in any change in gene expression or variation in alternative splicing at these loci.

Enriched GO terms hypomethylated under CR

The Biological Process GO term ‘methylation’ was enriched in or DMRs and associated with two ribosomal RNA methyltransferase genes and a Cap-specific mRNA (Nucleoside-2′-O-)-methyltransferase gene. The biological roles of these genes are poorly understood, but one hypothesis is that RNA methyltransferases are regulators of global protein translation [65]. One rRNA methyltransferase, NSUN5, is part of a conserved mechanism that modulates ageing [66], as reduced levels of this gene increases lifespan and stress resistance in yeast, worms and flies [66]. Another mechanism by which RNA methyltransferases may act is in protein synthesis that determine cell size, which is hypothesised to explain the small stature of organisms lacking rRNA methyltransferases [65]. This is concordant with the observed phenotype of adult CR Daphnia which are smaller than NF individuals, although functional validation would be required to link RNA methyltransferases causally to size difference phenotypes in Daphnia. Cap-specific mRNA (Nucleoside-2′-O-)-methyltransferase adds a methyl group to at the ribose 2′-O position of the first transcribed nucleotide, which is important for discriminating self RNAs from foreign RNAs produced by viruses [67, 68]. There is no clear link to caloric restriction with this gene.

The enriched molecular function GO term, Acyl CoA dehydrogenase activity, is associated with two medium-chain specific acyl-CoA dehydrogenase genes. These genes are involved in the breakdown of fatty acid molecules through β-oxidation. Reduced activity of this class of genes in humans impairs the β-oxidation pathway which leads to an intolerance to fasting [69]. There are no studies linking these genes to diet in worms or flies as yet, however they are implicated in heat adaptation through the stabilization of lipid membranes [70]. Broadly interpreted, lipid metabolism is strongly linked to caloric restriction-derived lifespan increase in flies [71], but the specific role of acyl CoA dehydrogenases remains to be elucidated. DNA methylation may act in concert with other mechanisms in response to CR as has been shown with miRNAs [48], but most-likely also through histone modifications [72]. A focus for future research will be the interaction of these different mechanism and their impact, or lack thereof, on global gene expression changes that result in the observed Daphnia CR phenotypes.

This experiment was performed with clonally reproducing Daphnia. This has the large advantage of a genetically uniform background on which to base results [72], but caution is required when extrapolating our results to sexually reproducing organisms. As discussed by Kvist et al. [24], it is not known if differential methylation is maintained through both parthenogenesis (mitotic) and sexual reproduction (meiosis) as only parthenogenetically transmitted phenotypic changes have been studied. This is also the case for maternal effects resulting from CR in the strain of D. magna used here [9, 10, 46, 48]. Jeremias et al. [23] speculate that parthenogenetic species may use epigenetic modifications to introduce phenotypic variation rapidly in response to a stimulus versus obligately sexual reproducers. They note that Daphnia species are well known for their phenotypic plasticity, especially in defensive modifications to the exoskeleton. The ‘resetting’ of CpG methylation in mammals for example could prevent such a mechanism of phenotypic plasticity propagating [73]. However, the prevalence of this meiotic reset does vary among vertebrate taxa [73–75], suggesting the potential exists for propagation of DNA methylation changes through sexual reproduction. The mode of methylation in the Daphnia genome: low global levels of methylation interspersed with highly methylated gene bodies, is in keeping with sexually reproducing invertebrates [23, 25].

Daphnia magna preparation and experiment

We used a single clone (known to us as Clone 32) of D. magna, collected from the Kaimes pond near Leitholm in the Scottish Borders [79]. Six replicates of control (i.e. well-fed, or normal food (NF)) D. magna were compared to six replicates of caloric restricted (CR) D. magna to identify differentially methylated regions. Maternal lines were first acclimatized for three generations. For this, individuals were kept in artificial pond medium [80] at 20 °C and on a 12 h:12 h light:dark cycle and fed 2.5 × 10⁶ cells of the single-celled green algae Chlorella vulgaris daily. Following three generations of acclimatisation (detailed in [48]), 40 offspring from each mother were isolated and split to form a replicate. Twenty were fed a normal diet of 5 × 10⁶ algal cells/day and the remaining 20 that were fed a CR diet of 1 × 10⁶ algal cells/day. Each replicate was split and reared in four sub-replicate jars of five animals. These were pooled for DNA extraction. Hence, NF and CR replicates were paired by mother and each consisted of 20 individuals in total. The experiment was ended after the birth of second clutch (approximately day 12 of the treatment generation). D. magna were ground by motorized pestle in Digsol and proteinase K and incubated overnight at 37 °C and stored at − 70 °C until DNA extraction.

DNA extraction and sequencing

DNA was extracted from pooled D. magna per replicate by phenol-chloroform followed by a Riboshredder RNA digestion step and repeat of the phenol-chloroform step. DNA was eluted into 100 μl of TE buffer and quantified by Qubit fluorimeter Sample purity was checked by 260:280 ratio on Nanodrop, and DNA integrity was examined by running approximately 35 ng DNA on a 0.8% agarose gel stained with ethidium bromide. Each DNA extraction was split in two for creation of a bisulphite converted library and corresponding bisulphite unconverted library (all steps the same except bisulphite conversion). Twenty-four libraries were created: 12 bisulphite-converted and 12 corresponding unconverted samples. This was done to identify per replicate mismatches from the reference and remove false positive methylation calls. All libraries were created by Edinburgh Genomics using the Zymogen EZ DNA Methylation-Lightning Kit and Methylseq Library prep Illumina TruSeq DNA Methylation Kit and 125 base pair paired-end sequenced on Illumina HiSeq. Raw read data have been deposited in the European Nucleotide Archive under accession PRJEB24784 (file names and conversion status, Additional file 6).

Quality assessment and mapping

Before aligning reads to the D. magna reference genome (version 2.4 downloaded from: http://arthropods.eugenes.org/EvidentialGene/daphnia/daphnia_magna/Genes/function/cddrpsdapmaevg14.gotab2), the reference was edited to match the haplotypes present in clone 32 [48]. This conversion increased read mapping efficiency by reducing polymorphism between the reference (assembled from a different clone) and our data.

Reads from both bisulphite converted and unconverted libraries were trimmed of the first and last nine bases of every read using TrimGalore! (version 0.4.1) [81] after initial inspection of Bismark m-bias plots. TrimGalore! was also used to remove base calls with a Phred score of 20 or lower, adapter sequences, and sequences shorter than 20 bases. FastQC 0.11.4 (downloaded from: http://www.bioinformatics.babraham.ac.uk/projects/fastqc) was used to inspect the data. Bisulphite calls were made with Bismark 0.16.3 [59]. Bismark alignments were performed with option “–score_min L,0,-0.6”. PCR duplicates were removed using the deduplicate_bismark script. Bismark reports indicated that libraries were not fully bisulphite converted and raw methylation calls were approximately 3% for CpG, CHH and CHG sites. Previous research has shown that CHH and CHG methylation is negligible in the D. magna genome [29]. As a result, we used the filter_non_conversion script to remove reads that contained CHH and/or CHG methylation sites as indicative of a non-bisulphite converted read. Finally, methylated sites were identified using bismark_methylation_extractor and reports created with bismark2report.

Further variants were predicted per replicate using the unconverted library reads by following the GATK pipeline [82] and converted strain 32 reference sequence. One round of base quality score recalibration was sufficient using variants previously identified from strain 32. Variants were hard-filtered in GATK using recommendations for single nucleotide polymorphisms (SNPs): quality by depth (QD) > 2, fisher strand (FS) > 60, root mean square of mapping quality (MQ) < 40, mapping quality rank sum test (MQRankSum) < − 12.5 and read position rank sum (ReadPosRankSum) < − 8.0. Sites with SNPs at methylated positions were removed from the analysis using BEDtools [83].

Differential methylation analysis and region validation

All analyses were performed using methylation calls from the bisulphite converted libraries only. Hierarchical sample clustering of genome-wide methylation patterns across replicates was generated using methylKit [84] and pvclust [85]. A matrix was prepared by adapting R code of the “clusterSamples” command in methyKit and used as input to Pvclust (Additional file 7). The dendrogram was then inferred by hierarchical clustering of correlation-based distances using the Ward.D2 method, confidence in relationships between replicates was provided by 10,000 bootstraps. The 6 NF and six CR replicates were then compared using the bsseq [61] Bioconductor package in R to identify regions of differential methylation, the input files were the same as for hierarchical clustering (Additional file 7). Bsseq uses the Bsmooth algorithm to identify local methylation estimates from low-coverage data, it exploits highly-correlated methylation levels across genomes. D. magna is no exception despite low-levels of CpG methylation. It has a “mosaic-like” distribution of methylation as shown by [25]. Regions of differential methylation are then assessed using the biological variation provided by the replicates. We ran bsseq with a paired t-test to incorporate differences between pairs due to mother on methylation. CpG sites were filtered so that at least two replicates in each treatment had six or more reads. We converted locally corrected t-test statistics into p-values (with five degrees of freedom, as six replicate pairs minus one) and a significance threshold of 0.05 was applied. P-values were converted to q-values using the Bioconductor package “qvalue” [62] and a false discovery rate threshold of 0.05 was applied to each CpG site (p- and q-value R scripts: Additional file 7). Bsseq identified DMRs were selected using a t-statistic cutoff of − 4.6 and 4.6, a greater than 0.1 average difference in methylation between groups, and at least three methylated CpGs. Genes overlapping DMRs were identified using the D. magna assembly version 2.4 genome annotation file and BEDTools. This list of overlapping genes was used as a basis for functional enrichment analysis.

To validate bsseq filtered DMR regions we adapted the approach of Herb et al. [20]. We took raw p-values per methylated CpG calculated in bsseq, as opposed to bumphunter in [20], and overlapped them with all of the bsseq predicted DMR regions, not the filtered set. These p-values were then combined using comb-p [86] to identify shared peaks between the bsseq DMR and comb-p approach. Comb-p is a five-step process that first (1) calculates the correlation between proximal p-values of varying distances (the autocorrelation), (2) combines adjacent p-values using the Stouffer-Liptak-Kechris correction which incorporates autocorreclation, (3) adjusts for false discovery by the Benjamini-Hochberg correction, (4) identifies enrichment regions as series of low p-values, and (5) p-value per region is assigned using the Stouffer-Liptak correction. The results of comb-p and bsseq filtered DMRs were intersected and those regions with corrected p-values from comb-p were further considered.

Functional enrichment analysis

Enriched GO (gene ontology) terms were identified using topGO. The weight01 algorithm was applied as it considers the GO hierarchy, a Fisher’s exact test with an α of 0.01 was used to identify enriched genes, and no p-value correction was applied as per topGO author recommendation (R code, Additional file 7). Hypo- and hypermethylated genes of interest were tested against databases of GO terms for all D. magna genes, those with at least one methylated CpG, and those with at least three CpGs to match the bsseq filtering of DMRs. The per-CpG raw methylation (pre-smoothing) rates were extracted from the Bsmooth object using getMeth in bsseq. Molecular function (MF) and biological process (BP) GO terms are reported. D. magna GO terms were downloaded from http://arthropods.eugenes.org/EvidentialGene/daphnia/daphnia_magna/Genes/function/cddrpsdapmaevg14.gotab2.