- Research article
- Open Access
Genome-wide methylation is modified by caloric restriction in Daphnia magna
BMC Genomics volume 20, Article number: 197 (2019)
The degradation of epigenetic control with age is associated with progressive diseases of ageing, including cancers, immunodeficiency and diabetes. Reduced caloric intake slows the effects of ageing and age-related disease in vertebrates and invertebrates, a process potentially mediated by the impact of caloric restriction on epigenetic factors such as DNA methylation. We used whole genome bisulphite sequencing to study how DNA methylation patterns change with diet in a small invertebrate, the crustacean Daphnia magna. Daphnia show the classic response of longer life under caloric restriction (CR), and they reproduce clonally, which permits the study of epigenetic changes in the absence of genetic variation.
Global cytosine followed by guanine (CpG) methylation was 0.7–0.9%, and there was no difference in overall methylation levels between normal and calorie restricted replicates. However, 333 differentially methylated regions (DMRs) were evident between the normally fed and CR replicates post-filtering. Of these 65% were hypomethylated in the CR group, and 35% were hypermethylated in the CR group.
Our results demonstrate an effect of CR on the genome-wide methylation profile. This adds to a growing body of research in Daphnia magna that demonstrate an epigenomic response to environmental stimuli. Specifically, gene Ontology (GO) term enrichment of genes associated with hyper and hypo-methylated DMRs showed significant enrichment for methylation and acyl-CoA dehydrogenase activity, which are linked to current understanding of their roles in CR in invertebrate model organisms.
Epigenetic modifications play a key role in patterning gene expression and organismal development. This is particularly evident when epigenetic control degrades, resulting in progressive diseases in humans, including cancers, immunodeficiency and diabetes . The degradation of epigenetic control with age is proposed to occur in a drift-like process. One mechanism that may rescue age-related epigenetic dysregulation is caloric restriction (CR) – reduced caloric intake without malnutrition or loss of micro-nutrients. CR slows the effects of aging and postpones the development of age-related diseases [2,3,4]. Extension of lifespan is observed in invertebrates and vertebrates, including in yeast, spiders, copepods, worms, silkworms, fish, non-human primates and Daphnia [4,5,6,7,8,9,10]. The marine copepod crustacean Paracartia grani is an example of this , and the observed lifespan increase is coupled to lower cumulative egg production in comparison with cohorts fed a high-food diet. This is interpreted as an evolutionary trade-off that allows survival during food shortages that increases the probability of recruitment for some offspring [7, 11]. In rhesus monkeys and mice, a caloric restriction of 30 and 40% respectively appears to reduce epigenetic drift in methylation and increases lifespan. In rodents lifespan can be extended by up to 50% . CR may also delay onset of a spectrum of diseases including cancer, kidney disease, autoimmune disease and diabetes [13,14,15], as well as neurodegenerative diseases [16, 17].
Mammals have been the focus of studies linking methylation, CR, and ageing due to the high rates of methylation in vertebrates, and the regulatory effect of methylation on disease progression in humans. Conversely, the lower genome-wide rates of methylation in arthropods  has hindered study of methylation-derived responses to the environment in this group, but this is now changing [8, 18,19,20,21,22,23,24,25]. Especially as the model organism Drosophila melanogaster, along with other surveyed Dipterans, exhibits negligible CpG methylation .
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  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 , 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 , 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 . 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,31,32,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,34,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 . 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 . Interestingly, in silkworms there is no correlation between methylation in promoters and gene expression , 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 . 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,23,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 . Although promoter regions of certain genes become hypermethylated with age . This is in contrast to honey bees in which an increase in de novo methylation activity through DNMT3 is observed . Inhibition of DNMT activity increased lifespans of treated bees, but in a CR independent manner . 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  observed during ageing. CR also causes a reduction in lipid metabolism gene expression by DNA methylation of gene bodies in mouse livers . 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 , including clear maternal effect phenotypes in D. magna [9, 10, 46,47,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 . This is also true for D. pulex where increased lifespans are observed under CR [50,51,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,23,24,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,55,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,22,23,24,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  using Bismark  exhibited lower mapping efficiencies than standard short-read alignments, as is typical of WGBS , 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  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” .
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 . This resulted in 333 regions (Fig. 2. Volcano plot of FDR corrected p-values) adapting Herb et al. . 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 . 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).
GO term enrichment in methylated genes and DMRs
GO term enrichment was explored using the “weight01” algorithm in topGO (downloaded from: https://bioconductor.org/packages/release/bioc/html/topGO.html) and molecular function (MF) and biological process (BP) terms are reported. The enrichment analysis showed significant over representation in one molecular function “acyl-CoA dehydrogenase activity” (GO:0003995) GO term and one biological process “methylation” (GO:0032259) GO term for genes associated with the DMRs (Table 2, Additional file 5), both of which were hypomethylated under CR in all categories of gene tested (all, > 1 CpG, > 3 CpGs). This is not altogether surprising as the number of gene containing methylated CpGs was a high fraction of the total annotated genes (75% for > 1 CpG, and 68% for > 3 CpGs) and is keeping with previous observations . There were two genes associated with acyl-CoA dehydrogenase activity and three with methylation, none of which overlap between terms. The methylation genes consist of two ribosomal RNA methyltransferases and a Cap-specific mRNA (Nucleoside-2′-O-)-methyltransferase, none of which are a member of the DNMT gene family. Both of the acyl-CoA dehydrogenase genes are medium-chain specific acyl-CoA dehydrogenases.
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 . In line with previously observed methylation patterns in Daphnia [22,23,24,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 . 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 , but the impact of environmentally-induced changes on gene expression is not yet understood beyond alternative splicing .
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.  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 . One rRNA methyltransferase, NSUN5, is part of a conserved mechanism that modulates ageing , as reduced levels of this gene increases lifespan and stress resistance in yeast, worms and flies . 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 . 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 . 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 . Broadly interpreted, lipid metabolism is strongly linked to caloric restriction-derived lifespan increase in flies , 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 , but most-likely also through histone modifications . 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 , but caution is required when extrapolating our results to sexually reproducing organisms. As discussed by Kvist et al. , 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.  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 . However, the prevalence of this meiotic reset does vary among vertebrate taxa [73,74,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].
We have shown that caloric restriction effects the methylation status of a subset of genes in Daphnia magna, despite the low overall CpG methylation found in this species. This is in line with recent work into other environmental stressors of Daphnia. Functional enrichment analysis of the annotation of loci with changed methylation status identified hypomethylation in the CR treatment of genes involved in methylation and acyl-coenzyme A dehydrogenases as the functionally enriched mode of the response. While we have focused on the effect of caloric restriction on DNA methylation status, there are alternative potential epigenetic responses to CR, including differential expression of small RNAs (sRNAs) and histone modifications. Previously, we established that CR induces differential micro RNA (miRNA) expression in D. magna under an equivalent experimental design , but other sRNAs, for example piwi interacting RNAs (piRNAs) and transfer RNA-derived small RNA (tsRNA), could also have a role in CR-dependent gene regulation [76, 77]. Histone modification in response to CR or protein restriction are proposed to increase longevity  by delaying and repressing ageing-related processes and diseases. Future studies using the tractable D. magna model could vary a range of dietary components (overall calories, proteins or fatty acids), and examine the joint effects of a range of epigenetic mechanisms.
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 . 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  at 20 °C and on a 12 h:12 h light:dark cycle and fed 2.5 × 106 cells of the single-celled green algae Chlorella vulgaris daily. Following three generations of acclimatisation (detailed in ), 40 offspring from each mother were isolated and split to form a replicate. Twenty were fed a normal diet of 5 × 106 algal cells/day and the remaining 20 that were fed a CR diet of 1 × 106 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 . 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)  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 . 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 . 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  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 .
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  and pvclust . 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  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 . 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”  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. . We took raw p-values per methylated CpG calculated in bsseq, as opposed to bumphunter in , and overlapped them with all of the bsseq predicted DMR regions, not the filtered set. These p-values were then combined using comb-p  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.
ATP-binding cassette sub-family D member 1
Calcium/calmodulin-dependent protein kinase kinases
Cytosine followed by nucleotides where H is adenine, thymine or cytosine and G is guanine
Cytosine followed by guanine
Differentially methylated region
European Nucleotide Archive
False discovery rate
Intracisternal A particle
Minichromosome maintenance complex component 2
Nicotinamide adenine dinucleotide
Polymerase chain reaction
piwi interacting RNA
Structural maintenance of chromosomes protein 2
Tris (hydroxymethyl) aminomethane ethylenediaminetetraacetic acid
transfer RNA-derived small RNA
Whole genome bisulphite sequencing
Egger G, Liang G, Aparicio A, Jones PA. Epigenetics in human disease and prospects for epigenetic therapy. Nature. Nature Publishing Group. 2004;429:457.
Cooper TM, Mockett RJ, Sohal BH, Sohal RS, Orr WC. Effect of caloric restriction on life span of the housefly, Musca domestica. FASEB J FASEB. 2004;18:1591–3.
Forster MJ, Morris P, Sohal RS. Genotype and age influence the effect of caloric intake on mortality in mice. FASEB J FASEB. 2003;17:690–2.
Colman RJ, Anderson RM, Johnson SC, Kastman EK, Kosmatka KJ, Beasley TM, et al. Caloric restriction delays disease onset and mortality in rhesus monkeys. Science (80-. ). American Association for the Advancement of Science. 2009;325:201–4.
Weindruch R. The retardation of aging by caloric restriction: studies in rodents and primates. Toxicol. Pathol. Sage publications sage CA: thousand oaks, CA. 1996;24:742–5.
Weindruch R, Walford RL. Retardation of aging and disease by dietary restriction. Springfield: CC Thomas; 1988.
Saiz E, Calbet A, Griffell K, Bersano JGF, Isari S, Solé M, et al. Ageing and caloric restriction in a marine planktonic copepod. Sci Rep. 2015;5:14962.
Xiang H, Zhu J, Chen Q, Dai F, Li X, Li M, et al. Single base-resolution methylome of the silkworm reveals a sparse epigenomic map. Nat Biotechnol. 2010;28:516–20.
Clark J, Garbutt JS, McNally L, Little TJ. Disease spread in age structured populations with maternal age effects. Ecol Lett Wiley Online Library. 2017;20:445–51.
Garbutt JS, Little TJ. Maternal food quantity affects offspring feeding rate in Daphnia magna. Biol. Lett. The Royal Society. 2014;10:20140356.
Kirkwood TBL, Shanley DP. Food restriction, evolution and ageing. Mech. Ageing Dev. Elsevier. 2005;126:1011–6.
Sohal RS, Weindruch R. Oxidative Stress, Caloric Restriction, and Aging. Science (80-. ). 1996;273:59–63.
Fernandes G, Yunis EJ, Good RA. Suppression of adenocarcinoma by the immunological consequences of calorie restriction. Nature. Nature Publishing Group. 1976;263:504.
Sarkar NH, Fernandes G, Telang NT, Kourides IA, Good RA. Low-calorie diet prevents the development of mammary tumors in C3H mice and reduces circulating prolactin level, murine mammary tumor virus expression, and proliferation of mammary alveolar cells. Proc. Natl. Acad. Sci. National Acad Sciences. 1982;79:7758–62.
Kubo C, Johnson BC, Good RA. A crucial influence of total calorie intake on autoimmune-prone mice: influence of diets of grossly different composition on immunologic functions. Fed Proc. 1984;3:2249.
Duan W, Mattson MP. Dietary restriction and 2-deoxyglucose administration improve behavioral outcome and reduce degeneration of dopaminergic neurons in models of Parkinson’s disease. J Neurosci Res Wiley Online Library. 1999;57:195–206.
Zhu H, Guo Q, Mattson MP. Dietary restriction protects hippocampal neurons against the death-promoting action of a presenilin-1 mutation. Brain Res Elsevier. 1999;842:224–9.
Bewick AJ, Vogel KJ, Moore AJ, Schmitz RJ. Evolution of DNA methylation across insects. Mol Biol Evol. 2017;34:654–65.
Cardoso-Júnior CAM, Guidugli-Lazzarini KR, Hartfelder K. DNA methylation affects the lifespan of honey bee (Apis mellifera L.) workers-evidence for a regulatory module that involves vitellogenin expression but is independent of juvenile hormone function. Insect Biochem. Mol. Biol. Elsevier. 2018;92:21–9.
Herb BR, Shook MS, Fields CJ, Robinson GE. Defense against territorial intrusion is associated with DNA methylation changes in the honey bee brain. BMC Genomics. BioMed Central. 2018;19:216.
Asselman J, De Coninck DIM, Vandegehuchte MB, Jansen M, Decaestecker E, De Meester L, et al. Global cytosine methylation in Daphnia magna depends on genotype, environment, and their interaction. Environ Toxicol Chem Wiley Online Library. 2015;34:1056–61.
Asselman J, De Coninck DIM, Beert E, Janssen CR, Orsini L, Pfrender ME, et al. Bisulfite sequencing with Daphnia highlights a role for epigenetics in regulating stress response to Microcystis through preferential differential methylation of serine and threonine amino acids. Environ. Sci. Technol. 2017;51:924–31.
Jeremias G, Barbosa J, Marques SM, De Schamphelaere KAC, Van Nieuwerburgh F, Deforce D, et al. Transgenerational inheritance of DNA Hypomethylation in Daphnia magna in response to salinity stress. Environ. Sci. Technol. 2018;52:10114–23.
Trijau M, Asselman J, Armant O, Adam-Guillermin C, De Schamphelaere KAC, Alonzo F. Transgenerational DNA methylation changes in Daphnia magna exposed to chronic γ irradiation. Environ Sci Technol. 2018;52:4331–9.
Kvist J, Gonçalves Athanàsio C, Shams Solari O, Brown JB, Colbourne JK, Pfrender ME, et al. Pattern of DNA methylation in Daphnia: evolutionary perspective. Genome Biol. Evol. 2018;10:1988–2007.
Bestor TH. Activation of mammalian DNA methyltransferase by cleavage of a Zn binding regulatory domain. EMBO J Wiley Online Library. 1992;11:2611–7.
Ehrlich M, Gama-Sosa MA, Huang L-H, Midgett RM, Kuo KC, McCune RA, et al. Amount and distribution of 5-methylcytosine in human DNA from different types of tissues or cells. Nucleic Acids Res. Oxford University Press. 1982;10:2709–21.
Olson CE, Roberts SB. Genome-wide profiling of DNA methylation and gene expression in Crassostrea gigas male gametes. Front Physiol. 2014;5:224.
Asselman J, De Coninck DIM, Pfrender ME, De Schamphelaere KAC. Gene body methylation patterns in Daphnia are associated with gene family size. Genome Biol Evol. 2016;8:1185–96.
Zemach A, IE MD, Silva P, Zilberman D. Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science (80-. ). American Association for the Advancement of Science. 2010;328:916–9.
Sarda S, Zeng J, Hunt BG, Yi SV. The evolution of invertebrate gene body methylation. Mol. Biol. Evol. Oxford University Press. 2012;29:1907–16.
Bonasio R, Li Q, Lian J, Mutti NS, Jin L, Zhao H, et al. Genome-wide and caste-specific DNA methylomes of the ants Camponotus floridanus and Harpegnathos saltator. Curr Biol Elsevier. 2012;22:1755–64.
Suzuki MM, Kerr ARW, De Sousa D, Bird A. CpG methylation is targeted to transcription units in an invertebrate genome. Genome Res Cold Spring Harbor Lab. 2007;17:625–31.
Flores K, Wolschin F, Corneveaux JJ, Allen AN, Huentelman MJ, Amdam GV. Genome-wide association between DNA methylation and alternative splicing in an invertebrate. BMC Genomics. 2012;13:480.
Li-Byarlay H, Li Y, Stroud H, Feng S, Newman TC, Kaneda M, et al. RNA interference knockdown of DNA methyl-transferase 3 affects gene alternative splicing in the honey bee. Proc Natl Acad Sci National Acad Sciences. 2013;110:12750–5.
Li S, Zhang J, Huang S, He X. Genome-wide analysis reveals that exon methylation facilitates its selective usage in the human transcriptome. Brief Bioinform. 2017;19:754–64.
Horvath S. DNA methylation age of human tissues and cell types. Genome Biol. BioMed Central. 2013;14:3156.
Horvath S. Erratum to: DNA methylation age of human tissues and cell types. Genome Biol. BioMed Central. 2015;16:96.
Pal S, Tyler JK. Epigenetics and aging. Sci Adv. 2016;2:e1600584.
Jung M, Pfeifer GP. Aging and DNA methylation. BMC Biol. BioMed Central. 2015;13:–7.
Li Y, Liu L, Tollefsbol TO. Glucose restriction can extend normal cell lifespan and impair precancerous cell growth through epigenetic control of hTERT and p16 expression. FASEB J. FASEB. 2010;24:1442–53.
Hahn O, Grönke S, Stubbs TM, Ficz G, Hendrich O, Krueger F, et al. Dietary restriction protects from age-associated DNA methylation and induces epigenetic reprogramming of lipid metabolism. Genome Biol. 2017;18:56.
Lane MA, Tilmont EM, De Angelis H, Handy A, Ingram DK, Kemnitz JW, et al. Short-term calorie restriction improves disease-related markers in older male rhesus monkeys (Macaca mulatta). Mech Ageing Dev Elsevier. 2000;112:185–96.
Kim CH, Lee EK, Choi YJ, An HJ, Jeong HO, Park D, et al. Short-term calorie restriction ameliorates genomewide, age-related alterations in DNA methylation. Aging Cell. 2016;15:1074–81.
Latta LC, Frederick S, Pfrender ME. Diet restriction and life history trade-offs in short-and long-lived species of Daphnia. J. Exp. Zool. Part a Ecol. Genet. Physiol. Wiley Online Library. 2011;315:610–7.
Garbutt JS, Little TJ. Bigger is better: changes in body size explain a maternal effect of food on offspring disease resistance. Ecol Evol Wiley Online Library. 2017;7:1403–9.
Coakley CM, Nestoros E, Little TJ. Testing hypotheses for maternal effects in Daphnia magna. J Evol Biol Wiley Online Library. 2018;31:211–6.
Hearn J, Chow FW-N, Barton H, Tung M, Wilson P, Blaxter M, et al. Daphnia magna microRNAs respond to nutritional stress and ageing but are not transgenerational. Mol Ecol. 2018;27:1402–12.
Pietrzak B. Interclonal differences in age-specific performance in Daphnia magna. J Limnol. 2011;70:345–52.
Lynch M. The life history consequences of resource depression in Daphnia pulex. Ecology Wiley Online Library. 1989;70:246–56.
Dudycha JL. A multi-environment comparison of senescence between sister species of Daphnia. Oecologia Springer. 2003;135:555–63.
McCauley E, Murdoch WW, Nisbet RM. Growth, reproduction, and mortality of Daphnia pulex Leydig: life at low food. Funct Ecol JSTOR. 1990;4:505–14.
Schwartz TS, Pearson P, Dawson J, Allison DB, Gohlke JM. Effects of fluctuating temperature and food availability on reproduction and lifespan. Exp Gerontol Elsevier. 2016;86:62–72.
Kim E, Ansell CM, Dudycha JL. Resveratrol and food effects on lifespan and reproduction in the model crustacean Daphnia. J. Exp. Zool. A. Ecol. Genet. Physiol. 2013/10/16. 2014;321:48–56.
Anderson BG. The number of pre-adult instars, growth, relative growth, and variation in Daphnia magna. Biol Bull Marine Biological Laboratory. 1932;63:81–98.
Agar WE. A statistical study of regeneration in two species of Crustacea. J Exp Biol The Company of Biologists Ltd. 1930;7:349–69.
Vandegehuchte MB, De Coninck D, Vandenbrouck T, De Coen WM, Janssen CR. Gene transcription profiles, global DNA methylation and potential transgenerational epigenetic effects related to Zn exposure history in Daphnia magna. Environ Pollut Elsevier. 2010;158:3323–9.
Vandegehuchte MB, Kyndt T, Vanholme B, Haegeman A, Gheysen G, Janssen CR. Occurrence of DNA methylation in Daphnia magna and influence of multigeneration cd exposure. Environ Int Elsevier. 2009;35:700–6.
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. Oxford University Press. 2011;27:1571–2.
Porter J, Sun MAN, Xie H, Zhang L. Investigating bisulfite short-read mapping failure with hairpin bisulfite sequencing data. BMC Genomics. BioMed Central Ltd. 2015;16:S2.
Hansen KD, Langmead B, Irizarry RA. BSmooth: from whole genome bisulfite sequencing reads to differentially methylated regions. Genome Biol. BioMed Central. 2012;13:R83.
Storey JD, Tibshirani R. Statistical significance for genomewide studies. Proc. Natl. Acad. Sci. 2003;100:9440–5.
Orsini L, Gilbert D, Podicheti R, Jansen M, Brown JB, Solari OS, et al. Daphnia magna transcriptome by RNA-Seq across 12 environmental stressors. Sci. Data. The Author(s). 2016;3:160030.
Hunt BG, Glastad KM, Yi SV, Goodisman MAD. The function of intragenic DNA methylation: insights from insect epigenomes. Integr. Comp. Biol. 2013;53:319–28.
Blanco S, Frye M. Role of RNA methyltransferases in tissue renewal and pathology. Curr Opin Cell Biol Elsevier. 2014;31:1–7.
Schosserer M, Minois N, Angerer TB, Amring M, Dellago H, Harreither E, et al. Methylation of ribosomal RNA by NSUN5 is a conserved mechanism modulating organismal lifespan. Nat. Commun. Nature Publishing Group. 2015;6:6158.
Züst R, Cervantes-Barragan L, Habjan M, Maier R, Neuman BW, Ziebuhr J, et al. Ribose 2′-O-methylation provides a molecular signature for the distinction of self and non-self mRNA dependent on the RNA sensor Mda5. Nat. Immunol. Nature Publishing Group. 2011;12:137.
Daffis S, Szretter KJ, Schriewer J, Li J, Youn S, Errett J, et al. 2′-O methylation of the viral mRNA cap evades host restriction by IFIT family members. Nature. Nature Publishing Group. 2010;468:452.
Matern D, Rinaldo P. In: Adam M, Ardinger H, Pagon R, editors. GeneReviews Medium-Chain Acyl-Coenzyme A Dehydrogenase Deficiency. Seattle (WA): University of Washington, Seattle; 1993.
Ma DK, Li Z, Lu AY, Sun F, Chen S, Rothe M, et al. Acyl-CoA dehydrogenase drives heat adaptation by sequestering fatty acids. Cell. Elsevier. 2015;161:1152–63.
Katewa SD, Demontis F, Kolipinski M, Hubbard A, Gill MS, Perrimon N, et al. Intramyocellular fatty-acid metabolism plays a critical role in mediating responses to dietary restriction in Drosophila melanogaster. Cell Metab Elsevier. 2012;16:97–103.
Harris KDM, Bartlett NJ, Lloyd VK. Daphnia as an emerging epigenetic model organism. Genet Res Int Hindawi. 2012;2012:147892:1–147892:8.
Seisenberger S, Peat JR, Hore TA, Santos F, Dean W, Reik W. Reprogramming DNA methylation in the mammalian life cycle: building and breaking epigenetic barriers. Phil. Trans. R. Soc. B. The Royal Society. 2013;368:20110330.
Potok ME, Nix DA, Parnell TJ, Cairns BR. Reprogramming the maternal zebrafish genome after fertilization to match the paternal methylation pattern. Cell. Elsevier. 2013;153:759–72.
Wang L, Zhang J, Duan J, Gao X, Zhu W, Lu X, et al. Programming and inheritance of parental DNA methylomes in mammals. Cell Elsevier. 2014;157:979–91.
Chen Q, Shi J, Peng H, Zhang X, Zhang Y, Qian J, et al. Sperm tsRNAs contribute to intergenerational inheritance of an acquired metabolic disorder. Science (80-. ). 2016;351:397–400.
Brennecke J, Malone CD, Aravin AA, Sachidanandam R, Stark A, Hannon GJ. An epigenetic role for maternally inherited piRNAs in transposon silencing. Science (80-. ). American Association for the Advancement of Science. 2008;322:1387–92.
Li Y, Daniel M, Tollefsbol TO. Epigenetic regulation of caloric restriction in aging. BMC Med. 2011;9:98.
Auld SKJR, Hall SR, Housley Ochs J, Sebastian M, Duffy MA. Predators and patterns of within-host growth can mediate both among-host competition and evolution of transmission potential of parasites. Am. Nat. University of Chicago Press Chicago, IL. 2014;184:S77–90.
Klüttgen B, Dülmer U, Engels M, Ratte HT. ADaM, an artificial freshwater for the culture of zooplankton. Water Res Elsevier. 1994;28:743–6.
Krueger F. Trim galore. A wrapper tool around Cutadapt FastQC to consistently apply Qual. Adapt. trimming to FastQ files. 2015.
Van Der Auwera GA, Carneiro MO, Hartl C, Poplin R, Levy-Moonshine A, Jordan T, et al. From FastQ data to high confidence varant calls: the genome analysis toolkit best practices pipeline. Curr Protoc Bioinforma. 2014;11:11.10.1–11.10.33.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, et al. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13:R87.
Suzuki R, Shimodaira H. Pvclust: an R package for assessing the uncertainty in hierarchical clustering. Bioinformatics. Oxford University Press. 2006;22:1540–2.
Pedersen BS, Schwartz DA, Yang IV, Kechris KJ. Comb-p : software for combining, analyzing, grouping and correcting spatially correlated P -values. Bioinformatics. 2012;28:2986–8.
We are grateful to the reviewers of the manuscript; their constructive input improved our study substantially. We thank Peter Fields, University of Basel and Don Gilbert, Indiana University for providing scripts, advice, and data on Daphnia magna genome annotation. This research was funded by Leverhulme Trust grant RPG-2015-406 and Wellcome Trust Institutional Strategic Support Fund (round 2, University of Edinburgh).
This research was funded by Leverhulme Trust grant RPG-2015-406 and Wellcome Trust Institutional Strategic Support Fund (round 2, University of Edinburgh). The funders were not involved in study design, and collection, analysis, and interpretation of data and in writing the manuscript.
Availability of data and materials
The raw read data generated during the current study are available in the European Nucleotide Archive under accession PRJEB24784, https://www.ebi.ac.uk/ena/data/view/PRJEB24784.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Reads sequenced per replicate for both converted and unconverted libraries. Alignments analysed is number of read-pairs per bisulphite converted library aligned by Bismark; deduplicated is number of read pairs after removal of PCR duplicates; non CpG filtered is read-pairs remaining after removal of CHH and CHG containing reads; bases remaining are number of bases left for CpG methylation prediction in Bismark; average coverage is for filtered bases at a D. magna genome size of 240 megabases; CpG sites remaining before and after filtering of polymorphic sites, polymorphic sites were identified from unconverted libraries against a reference genome converted to Clone 32; Non CpG % methylated are CHH and CHG sites methylated (100 minus this figure is equivalent to bisulphite conversion efficiency), CpG % methylated is the genome-wide percentage of CpG methylation, and CpG % methylated filtered is the percentage remaining after reads containing CHH or CHG methylation are removed from the analysis. Diet, H: normal food, L: caloric restriction, code: combination of diet and mother. (XLSX 48 kb)
MultiQC report in html format reporting general statistics from the Bismark alignment process for all replicates. Including alignment rates, deduplication effect, overall cytosine methylation and m-bias plot. This plot shows average methylation level per position across reads, demonstrating minimal bias at 5′ and 3′ reads after trimming of first and last nine base pairs of each read. (HTML 1089 kb)
Gene annotations for genes with overlapping significantly hypomethylated regions under caloric restriction. Genes overlapped by a DMR of greater than 20% mean difference in methylation level are indicated in column 3. (XLSX 14 kb)
Gene annotations for genes with overlapping significantly hypermethylated regions under caloric restriction. Genes overlapped by a DMR of greater than 20% mean difference in methylation level are indicated in column 3. (XLSX 12 kb)
GO term hierarchies for the two enriched GO terms: A) Biological Process and B) Molecular Function. Ovals and rectangles contain GO identification, GO name, p-value and number of genes tested and over total number of genes from that category. Red rectangles represent significant terms under the topGO criteria applied. (PS 91 kb)
Raw read files used in this study deposited in the European Nucleotide Archive under accession PRJEB24784. Treatment, H: normal food and L: caloric restricted. (XLSX 34 kb)
R code used in methylKit, pvclust, bsseq, and comb-p sections of this study. (TXT 8 kb)