Epigenetic, transcriptional and phenotypic responses in two generations of Daphnia magna exposed to the DNA methylation inhibitor 5-azacytidine.

The water flea Daphnia magna is a keystone species in freshwater ecosystems and has been widely used as a model organism in environmental ecotoxicology. This aquatic crustacean is sensitive to environmental stressors and displays considerable plasticity in adapting to changing environmental conditions. Part of this plasticity may be due to epigenetic regulation of gene expression, including changes to DNA methylation and histone modifications. Because of the generally hypomethylated genome of this species, we hypothesized that the histone code may have an essential role in the epigenetic control and that histone modifications might be an early marker for stress. This study aims to characterize the epigenetic, transcriptional and phenotypic responses and their causal linkages in directly exposed adult (F0) Daphnia and peritoneal exposed neonates (F1) after a chronic (7-day) exposure to a sublethal concentration (10 mg/l) of 5-azacytidine, a well-studied vertebrate DNA methylation inhibitor. Exposure of the F0 generation significantly reduced the cumulative fecundity, accompanied with differential expression of genes in the one-carbon-cycle metabolic pathway. In the epigenome of the F0 generation, a decrease in global DNA methylation, but no significant changes on H3K4me3 or H3K27me3, were observed. In the F1 offspring generation, changes in gene expression, a significant reduction in global DNA methylation and changes in histone modifications were identified. The results indicate that exposure during adulthood may result in more pronounced effects on early development in the offspring generation, though interpretation of the data should be carefully done since both the exposure regime and developmental period is different in the two generations examined. The obtained results improve our understanding of crustacean epigenetics and the tools developed may promote use of epigenetic markers in hazard assessment of environmental stressors.

Introduction

Cladocerans, (water fleas), are common crustaceans found in freshwater lakes and ponds. Due to their high susceptibility to environmental stressors, Daphnia, in particular, have been extensively used as model organisms in ecotoxicological studies [1]. Moreover, current research indicates that Daphnids is an excellent model for environmental epigenetic research [2, 3]. The programming of each cell’s transcriptome is influenced by alterations in chromatin structures orchestrated by epigenetic modifications. These alterations are grouped into DNA base modification types, including 5-methyl-2′-deoxycytidine (mdC), non-coding RNAs and post-translational modifications (PTMs) of histone proteins, complexed with the DNA in chromatin [4]. To interpret the biological consequence of histone PTMs, a histone code [5] or language [6] has been suggested, where histone PTMs are associated with a transcriptional outcome. For example trimethylation of lysine (K) at the fourth position at the protruding tail of histone H3 (H3K4me3) is associated with open chromatin and active transcription, while trimethylation of H3K27 is associated with condensed chromatin and transcriptional repression. Until now, epigenetic studies in Daphnia have focused mainly on changes in mC content, primarily in the species Daphnia magna [3, 7–13] and Daphnia pulex [13, 14]. Although a small subset of genes have DNA methylation levels >50% [15], the global DNA methylation level is typically <1% [12, 16]. Detection of H3K4me3, H3K14ac, H4K20me2, H3K9me3 and H3K27me3 in Daphnia embryos [17, 18] suggests that D. magna uses histone PTMs as mechanism for transcriptional control. Compared to vertebrates like zebrafish, which exhibit ∼8% DNA methylation [19], invertebrate genomes, including that of D. magna, are hypomethylated. This suggests that histone PTMs might play a more dominant role in transcriptional control in this species, without precluding an active role for DNA methylation in epigenetic control.

Several epigenetic marks are dependent on methyl donors provided by the one-carbon metabolism (OCM) cycle. The OCM (Fig. 1) is a complex and essential pathway for cell biosynthesis and maintenance of the overall methylation state of cells [20]. The OCM consists of two cycles, where the folate cycle is linked to the methionine cycle through the reduction of 5,10-methyleneTHF to 5-methylTHF by Methylenetetrahydrofolate reductase (Mthfr). 5-methylTHF donates a methyl group to homocysteine resulting in the formation of methionine and THF. The methionine adenosyltransferase (Metk) generates S-adenosylmethionine (SAM), the principal donor of methyl groups required for the methylation of histones, non-histone proteins and DNA [21]. A methyl group from SAM is donated via the DNA methyltransferases (Dnmts) and histone methylasestomethylate DNA and histones. During this process, SAM is converted to S-adenosylhomocysteine (SAH) and further converted via the S-adenosylhomocysteine hydrolase (Sahh) to homocysteine, looping the methionine cycle [22]. The ratio of SAM and SAH is central for the balance between the folate and methionine cycles [23]. Changes to methionine concentration affects this ratio and thereby the activity of Dnmts and glycine N-methyltransferase (Gnmt) [24]. Furthermore, SAM inhibits Mthfr [25] and 5-methylTHF inhibits Gnmt [26]. Due to the importance of the OCM for cellular homeostasis, and its contribution to the energy balance by providing molecules of ATP and NADPH [27, 28], changes in OCM are expected to have an impact on other processes, including reproduction [29, 30]. A recent study showed that 5-azacytidine (5-AzaC) reduces DNA methylation in D. magna at several life stages and also affects the OCM [11]. Therefore, OCM genes related to DNA methylation and demethylation, the methionine cycle, the folate cycle and reproductive genes including Dnmts, methylcytosine dioxygenase (Tet2), Gnmt, Metk, Sahh and Mthfr were selected for evaluating epigenetic responses to 5-AzaC exposure.

Figure 1:

simplified OCM with enzyme genes investigated marked in red. Folic acid is the substrate in the folate cycle and used to produce SAM. SAM is used as an universal methyl donor and therefore essential for the methylation of DNA and histones (modified from [31])

We evaluate D. magna as a relevant model for the role of epigenetic controls in ecotoxicology through targeted analysis of genes related to methylation with or within the OCM pathway (Fig. 1) and by linking gene expression and histone modification responses across two generations (F0 and F1) following direct 5-AzaC exposure. To assess the possibility of additional interference with normal processes relevant for fecundity, gene expression of the methoprene-tolerant (Met) protein [32] and the two known Vitellogenin genes, i.e. Vtg1 and Vtg2 [33], were also assessed. Met proteins are involved in juvenile hormone signalling related to growth, metamorphosis and reproduction [32, 34] while Vtg encodes the precursors for phospholipo-proteins used by oviparous organisms as nutrients for the developing embryo [33]. These genes are expected to reflect reproductive health status, i.e. fertility and fecundity. In classic toxicology, hazardous chemicals are assessed on the basis of concentration at which effects appear on the individual or at the cohort level. No observed effect concentrations (NOEC) and lowest observed effect concentrations (LOEC) in addition to acute toxicity are measured by lethality or absence of response. NOEC and LOEC thereby serve to indicate the concentration at which no effect, or the first effects may occur in an organism. Molecular responses are expected to be more sensitive to stress and appear earlier than behavioural, organismal or lethal responses. Because of this, the use of NOEC and LOEC in toxicological research has been thoroughly debated [35].

The objective of this study was to determine epigenetic stress responses in Daphnia exposed to the demethylation chemical, 5-AzaC [13, 36], at concentrations that were non-toxic (e.g. lower than the NOEC). Since evolutionary conserved genes seem to be highly methylated in daphnids [13] and other invertebrates [37, 38], conserved genes key to the OCM, development and reproduction in crustaceans were targeted. In order to generate more detailed knowledge about epigenetic stress responses, a series of experiments comparing exposure induced changes in gene expression, histone PTMs and global DNA methylation with biological responses was carried out.

Results

Acute Toxicity

A concentration-dependent increase in D. magna mortality was observed for 5-AzaC after 48 h exposure, where NOEC and LOEC were 10 and 30 mg/l, respectively. The estimated EC50 was 230.1 ± 98.4 mg/l (R2 = 0.945, CI 95%, Fig. 2).

Figure 2:

mortality (mean ± SEM, % of control, n = 6) in cultures of juvenile D. magna exposed to 5-AzaC. Asterisks (*) indicate significant differences compared to the control (P < 0.05). Dashed line is the control value. NOEC and LOEC are indicated

Reproductive Effects

Significant (P = 0.0255) reduction in cumulative fecundity was observed in F0 of D. magna after 7-day exposure to 10 mg/l 5-AzaC (Fig. 3).

Figure 3:

cumulative fecundity in adult female D. magna (mean indicated, n = 6) after 7-day exposure to 10 mg/l 5-AzaC. Asterisk denotes significant difference from the control (P = 0.0255, unpaired t-test)

Global Methylation

Global DNA methylation analysis in both directly 5-AzaC exposed F0 and in utero exposed F1 showed a significant decline in DNA methylation (Fig. 4). We also screened for 5-hydroxymethyl-2′-deoxycytidine (hmdC), but this modification was not detected within the detection limit of 0.001%.

Figure 4:

global DNA methylation level in controls (black) compared to 7-day direct (F0) or perinatal exposure (F1) to 10 mg/l 5-AzaC (grey). Asterisks (*) indicate significant differences compared to the control (P < 0.0001)

Histone PTM Enrichment Profiles in Whole Adult Daphnids

Parthenogenetic proliferation, through clonal reproduction (i.e. eggs), is expected to contribute to increased specificity and reduced noise, an advantage when performing multigenerational sequence-specific epigenetic analysis. A chromatin immunoprecipitation (ChIP) assay for Daphnia was set up with pools of eight individual non-exposed adults. This assay was proven specific and reproducible after specific histone PTM patterns were obtained on different genomic locations. Histone PTM profiles upstream of the actively transcribed Actin gene showed high enrichment of H3K4me3 and H3K9ac together with low levels of H3K27me3 and H3K9me3 (Fig. 5A) in accordance with an actively transcribed gene. In addition, the gene body of 18SrRNA is enriched with the repressive H3K9me3 mark (Fig. 5B). For the Dnmt loci investigated, the upstream sequence regions showed unique histone PTM profiles for all three genes (Dnmt1a, Dnmt3a.1 and Dnmt3a.2) (Fig. 5C–E). From this, we concluded that the ChIP data obtained was reliable and that the upstream sequences analysed harbour regulatory features that can be used to assess enrichment modifications by environmental stressors. For the exposure study, we chose H3K4me3 and H3K27me3 as the most interesting candidate targets for analysis.

Figure 5:

enrichment pattern from pools of eight unexposed whole adult D. magna. ChIP analysis of H3K4me3, H3K27me3, H3K9me3 and H3K9ac on sequences flanking the TSS, ‘upstream’ and ‘gene body’ on indicated genes (A) Actin, (B) 18S rRNA, (C) Dnmt1a), (D) Dnmt3a.1 and (E) Dnmt3a.2. Enrichment is visualized as per cent of input, error margins reflects SEM

Histone PTM Enrichment on 5-AzaC Exposed F0 and F1 Generation

Direct 5-AzaC exposure (F0) did not cause any significant enrichment changes relative to controls forH3K4me3 and H3K27me3 (Fig. 6). However, differences in histone PTM enrichment were observed in perinatal exposed D. magna (F1), where 5-AzaC induced increases in H3K27me3 enrichment on the majority of the tested genes, while H3K4me3 level was evenly distributed with enrichment levels comparable to the control (Fig. 6 and Table 1). In addition, a clear, albeit non-significant, change in histone PTM enrichment was observed in the F1 generation for Dnmt1a, Dnmt3a.1, Dnmt3a.2, Tet2, Gnmt, Metk, Sahh and Mthfr (Table 1).

Figure 6:

relative differences in upstream sequence enrichment of H3K4me3 and H3K27me3 on indicated genes in adult D. magna. Control (grey) and 7-day 10 mg/l 5-AzaC exposed D. magna (black) are presented as direct (F0) and perinatal (F1) exposure

Table 1:

overview of the indicative direction of gene expression changes based on histone PTM enrichment in F0 and F1 generation in OCM cycle genes in D. magna exposed to 10 mg/l 5-AzaC (NOEC) for 7 days

Arrows indicate upregulation (‘↑’) or downregulation (‘↓’) of gene expression relative to the controls, asterisks within the arrows indicate a significant difference compared to the control. Zero (‘0’) indicates no apparent change and ‘n/a’ indicates no data. Colour coding: green represents anticipated changes in gene expression based on histone PTM enrichment, red reflects unanticipated changes based on PTM enrichment. Orange indicates inconsistent change. White arrows indicate only gene expression up or downregulation.

Gene Expression

We assessed four reference genes for downstream normalization of gene expression. There were no significant differences observed between control and exposed groups of the housekeeping gene cyclophilin (Cyp) in F0 (P = 0.565) and in F1 (P = 0.305). In contrast, Actin (F0 P = 0.326, F1 P = 0.003), Ef1a (F0 P = 0.044, F1 P = 0.118) and glyceraldehyde 3-phosphate dehydrogenase (Gapdh) (F0 P = 0.002, F1 P = 0.899) displayed significant differences between control and exposed groups in F0 or F1 generation (data not shown). Consequently, Cyp was used for normalization of the gene expression results. From the panel of genes investigated, significant differences on gene expression level in F0 generation were measured on the transcripts of the OCM cycle-related genes Dnmt1a and Dnmt3a.2, the methionine cycle-related genes Gnmt and Metk and the fecundity-related genes, Vtg1 and Vtg2. Perinatal exposure caused changes on Dnmt3a.1, otherwise no significant changes were observed in the F1 generation (Fig. 7). An apparent, albeit non-significant change in gene expression was observed for Dnmt1a, Dnmt3a.2, Tet2, Gnmt, Metk, Sahh, Mthfr, Vtg1, Vtg2 and Met for F1 daphnids (Table 1).

Figure 7:

relative gene expression in controls, direct (F0) and perinatal exposure (F1) after 7-day exposure to 10 mg/l 5-AzaC in D. magna. The transcripts of Dnmts, Gnmt, Metk, Mthfr, Sahh and Tet2 are one carbon cycle enzymes, while Met, Vtg1 and Vtg2 reflect fertility and fecundity. Error bars indicate SEM and asterisks (*) significance at the P < 0.05 level

Coherence between ChIP and Gene Expression after 5-AzaC Exposure

In general, it is expected that the chromatin environment will affect the transcriptional outcome. To interpret our findings, we summarized histone PTM and gene expression results (Table 1). In F0, we observed no concordance between 5-AzaC induced changes on H3K4me3 and H3K27me3 and gene expression. In F1, chromatin responses were observed. If the gene was upregulated or downregulated, H3K4me3 and H3K27me3 tended to be more highly enriched, with the exception of H3K4me3 enrichment on Dnmt1a and Gnmt (Figs 6 and 7 and Table 1).

Discussion

Daphnia magna is an ecologically key species that is widely used in ecological and environmental research because of its well-known biology, sequenced genome, parthenogenic life cycle, easy culturing and short generation time. Taken together this makes daphnids well suited as an in vivo model for laboratory experiments. They are known for their phenotypic plasticity exemplified with the kairomone induced transgenerational ‘helmet’ morphants [39], as well as the shift from parthenogenesis to sexual life cycles [40]. This morphological plasticity has been studied thoroughly, but the mechanistic understanding of epigenetic regulation in daphnids is still limited. Interestingly, crustaceans with low global amount of DNA methylation have been subject mainly to DNA methylation studies, while other epigenetic mechanisms like the histone PTMs have been poorly studied. We therefore hypothesized that histone PTMs could be important additional regulators of gene expression because of the relatively low DNA methylation content in the Daphnia genome [9, 14]. The 5-AzaC effective concentration (EC50) used in this study was comparable with that reported in a previous D. magna study [16]. This suggests, that our Daphnia strain is not over-susceptible to 5-AzaC. The study of molecular changes was thus conducted at concentrations not expected to cause acute toxicity (i.e. below the NOEC), and comparable to other studies [13]. The exposure of D. magna to sublethal concentrations of the DNA methylation inhibiting agent, 5-AzaC, significantly inhibited reproduction, reduced the amount of DNA methylation and changed gene expression in the F0 generation in our study. However, the response of histone PTMs was limited to the F1 generation only, suggesting that DNA methylation is an important player in epigenetic control of gene expression.

Whole Daphnids Reveal Gene-specific Histone PTM Enrichment

To obtain knowledge about the average chromatin landscape in whole daphnids, the enrichment of a library of histone PTMs (H3K4me3, H3K27me3, H3K9me3 and H3K9ac) were measured. To see if the enrichment is specific for upstream regulatory regions, we also included gene body regions, where the enrichment pattern is expected to be different compared to the transcriptional regulatory regions as promoters. At the Actin locus upstream region, we found a typical pattern often observed at promoters of transcribed genes, namely an enrichment of H3K4me3 and H3K9ac in combination with low levels of H3K27me3 and H3K9me3. In contrast, an H3K9me3 enrichment profile was observed for 18S rRNA and was, despite being transcribed by RNA polymerase I, included as an example of specific histone PTM enrichment (Fig. 5). The analysed locus at 18S rRNA is localized downstream of transcription start site (TSS) and not in the upstream sequence, possibly explaining the low H3K4me3 and H3K9ac enrichment in addition to its being an RNA pol I transcribed gene. The different histone PTM patterns at Dnmt1a, Dnmt3a.1 and Dnmt3a.2 indicate differential epigenetic regulation of these genes, with the maintenance methylation enzyme gene Dnmt1 as more enriched, for instance, with H3K4me3 at the promoter compared to the de novo methyltransferase Dnmt3a.1. This is directly reflected in the raw gene expression data, where Dnmt1 has a significant (P = 0.004) lower Cq value, and thus a higher expression level compared to Dnmt3a.1 (results not shown). We conclude that the ChIP results demonstrate methodological robustness to generate reproducible and differentiated results for whole and pooled D. magna, as observed for whole body ChIP analyses on zebrafish (Danio rerio) larvae [41] and the nematode Caenorhabditis elegans [42].

5-AzaC Reduces Fecundity

Some of the most established chronic tests for chemical toxicity is measuring changes to reproduction and development, and studies that link such changes to methylation events are thus highly warranted. The downregulation of the vitellogenin genes in the F0-generation of this study may indicate reduction of reproductive success. Both Vtg1 and Vtg2 genes are required for normal vitellogenin production [43], the key step in enrichment of nutrition to developing embryos, albeit the latter responds better to chemical exposures [33]. As expected based on the findings of Hannas et al. [33], the most pronounced effect of 5-AzaC was measured on Vtg2 expression in F0, whereas Vtg1 was also significantly downregulated. Furthermore, differential gene expression of the two Vtg genes was also observed in the form of zinc induced upregulation of Vtg2 while Vtg1 was not differentially expressed [44]. The downregulation of vitellogenin gene expression by 5-AzaC seems to be associated with reduction in offspring produced, where the F0 generation produced significantly fewer offspring (F1 clones) compared to age-matched control samples (Fig. 3). Another master regulator of reproduction and sex determination in Daphnia, the Met gene, showed a higher mRNA expression in the perinatal exposed F1 after 5-AzaC exposure. The gene expression therefore reflects the current status by indicating compensatory actions in response to damage, e.g. upregulation of Vtg because of vitellogenin protein depletion and more vitellogenin synthesis is needed.

5-AzaC Reduces DNA Methylation in Both Generations

In a previous study, 5-AzaC was used together with other environmental contaminants in a multigenerational set-up to assess changes in global methylation patterns inherited by subsequent generations [16]. We report similar effects of global demethylation in exposed parents and in the perinatal exposed F1 juveniles. Interestingly, the basal global DNA methylation levels in F1 juveniles were lower compared to F0 adults. Although the biological interpretation of global DNA methylation levels should be made with care, we suggest that the observed differences (Fig. 4), in combination with histone PTM data, can support a model where adult daphnids rely on DNA methylation, whereas histone PTM marks play an active additional role in F1 juveniles. The cause of this discrepancy might imply reprogramming of DNA methylation during the life cycle as commonly observed in vertebrates [45]. Alternatively, the different level of DNA methylation in adults versus neonates could be a function of age, analogous to findings in human blood cells, where global DNA methylation increases a function of age until a decrease in late adulthood [46]. This can be a consequence of the epigenetic clock and epigenetic drift [47]. The current findings illustrate the importance of considering both methylation and histone PTM regulation in studies with rapidly developing species such as daphnids.

5-AzaC Induces Histone PTM Changes in F1 But Not F0

Based on histone PTM enrichment patterns of unexposed adults, the H3K4me3 and H3K27me3 showed the highest enrichment in the amplified regions (Fig. 5) and were chosen for the further experiments. We observed no or little effect on the chosen histone modifications in the F0 generation even though 5-AzaC induced differential gene expression, reduced global levels of DNA methylation and also reduced numbers of F1 offspring. However, in the F1 generation, histone PTMs on the genes of the OCM were affected (Table 1). From our data, it is inconclusive whether these differences are directly a result of exposure or if they are an indirect result arising from different exposure regimes or different developmental stages. Furthermore, exposure regime and developmental period can affect the transcriptional directional changes that are not explained by H3K4me3 and H3K27me3 enrichment levels. The simultaneous enrichment of H3K4me3 and H3K27me3 on several loci has been described as bivalent states in embryonic stem cells [48] and is also a consequence of studying a whole animal with different cell types. Furthermore, we also observed generally higher enrichment of H3K27me3 relative to controls in the F1 generation, as compared to F0, after exposure. Increased H3K27 trimethylation enrichment, which is mediated by polycomb repressive complex 2, suggests a chromatin reorganization or chromatin compactness induced by 5-AzaC. The lack of stress response on H3K4me3 and H3K27me3 enrichment in the directly exposed F0 generation may be specific for the chosen epitopes, where other histone modifications might be more involved in gene regulation caused by stress, as well as loci located at other genes might be more responsive. Future assessments, using ChIP sequencing would help shed light on genome-wide alterations of selected PTMs in Daphnia. In general, genes used for normalization purposes should be expressed at a constant level throughout the experiment. Several of the reference gene candidates were therefore excluded. Caution should be used when evaluating the observed differences of 5-AzaC induced responses between F0 and F1 because of the different exposure regimes and developmental stages. In contrast to the exposure of the adults, the F1 was most likely exposed during sensitive developmental stages and with shorter exposure times. We measured several reference genes to normalize our gene expression analysis. Significant changes in expression of several well-characterized reference genes by 5-AzaC shows its potential to alter even stably expressed genes at concentrations considered non-toxic (i.e. below the NOEC). One of the causes could lie in reduced DNA methylation induced by 5-AzaC, as was shown by our global methylation analysis.

Methylation Cycle Genes Respond to 5-AzaC

In general, Dnmt1a is responsible for the maintenance of DNA methylation, while Dnmt3 catalyse de novo DNA methylation [49], suggesting the significant higher transcription level of Dnmt3a.2 and significant lower expression of Dnmt1a following 5-AzaC exposure. This may be a compensatory response for DNA methylation inhibition (Fig. 7, F0). This change in gene expression did not correlate with anticipated changes in enrichment of H3K4me3 or H3K27me3 in the F0 generation. The Dnmt1a and Dnmt3a.1 of the F1 juveniles displayed higher expression than the F0 adults, implicating developmental-specific gene expression as DNA methylation plays important roles in development [50]. Significant transcriptional changes of the two Dnmt3 homologues indicate that they are involved in methylation regulation in D. magna and may be considered as potential environmental biomarker candidates for stressors causing changes in epigenetic methylation patterns. Higher expression of Tet2 in F1 suggests further oxidation of mdC to hmdC, although global hmdC was not detected in any of the life stages and was possibly below the detection threshold of the assay. Compared to the results in F0, the downregulation of Dnmt3a.2 in F1 is of particular interest, where it coheres with less H3K4me3 and higher H3K27me3 enrichment after exposure to 5-AzaC. This suggests that this gene is under epigenetic control via histone modifications in the perinatal exposed F1 generation.

5-AzaC Induce Higher Expression of Methionine Cycle Genes

Even though the Daphnia genome is hypomethylated, methyl donors are ubiquitously present. The F0 gene expression showed significant Metk upregulation and non-significant upregulation of Gnmt and Mthfr. In F1, Metk is downregulated while Gnmt and Sahh are both upregulated in F0 and F1. Hence, the 5-AzaC induced changes to gene expression in the two generations may reflect different stress responses on the methionine cycle-related genes. The upregulation of the genes for Gnmt and Sahh, which are responsible for conversion of SAM to homocysteine (HCY) (Fig. 1), indicates a high turnover of SAM in both generations. In F0, the formation of SAM from methionine and adenosine triphosphate catalysed by Metk, is significantly increased and thereby works in the same direction as Gnmt and Sahh. However, in the perinatal exposed organisms, downregulation of Metk can lead to less SAM. Since the level of Mthfr product is controlled by the level of SAM [25, 51], the observed difference in expression direction in F0 and F1 Mthfr can be explained by an age-specific response to SAM levels or different susceptibilities to the exposure compound (Table 1). As a universal methyl donor, SAM provides methyl groups to both DNA and histone substrates. Assuming substrate competition between histone methyltransferases and Dnmts, one would expect a larger degree of demethylation of DNA than methylated histones upon 5-AzaC exposure, which works on DNA itself. Interestingly, the upregulation of Gnmt, Metk and Sahh suggests excess SAM caused by DNA methylation inhibition by 5-AzaC in F0. The downregulation of Metk in F1, though non-significant, indicates lower turnover of methionine to SAM. The conversion of SAM to SAH, catalysed by the upregulated Gnmt leads to the upregulation of Mthfr, as also shown by our data (Table 1), due to less inhibition from SAM [25]. In total, this indicates that the methionine cycle genes in F1 are aberrantly expressed as a response to perinatal 5-AzaC exposure. The methionine cycle, represented through Mthfr expression only, appears to affect both F0 and F1. The slightly reduced expression levels of Mthfr in F0 and the small expression increase in F1 reflects the different nature of the exposure in F0 versus F1.

Conclusion

This study used a suite of molecular tools for quantifying changes in methylated histone PTMs, global levels of DNA methylation, gene expression and measuring fecundity to characterize responses to 5-AzaC over two generations. We observed changes on global DNA methylation and both significant and not significant changes on gene expression in F0 and F1 generations. The 5-AzaC induced changes on gene expression were anticipated, although the interpretation of directional changes was not easily explained by H3K4me3 and H3K27me3 enrichment levels upstream of TSS on the selected loci. However, we revealed that effects on histone PTMs were detected in the F1 generation even when there is no change in the parental (F0) generation. Future work on other epigenetic factors in combination with genome-wide sequencing is needed to gain further understanding of how alterations in chromatin as part of stress responses in D. magna correlate with observed phenotypic changes.

Methods

Daphnia Culture

Daphnia magna (Strain from DHI Water Environmental Health, Hørsholm, Denmark) were cultured in M7 medium [52] in a 20 ± 1°C climate room with a photoperiod of 16:8 h. The animals were fed daily with concentrated microalgae Pseudokirchneriella subcapitata (equivalent to 0.1 mg carbon per Daphnia per day), according to the OECD test guideline 211 (D. magna Reproduction Test). The pH in the culture medium was 7.8 ± 0.2 and dissolved oxygen >7 mg/l. The culture medium was renewed twice a week. Under these conditions, adult D. magna reproduce female neonates every 3 days.

Test Compound

5-AzaC was obtained from Sigma (Sigma-Aldrich, St. Louis, MO, USA, CAS 320-67-2, purity ≥98%). To minimize degradation of 5-AzaC, the test solutions were freshly made, shortly before the exposure by dissolving the chemical in the culture medium.

Exposure and Sampling

An acute (48 h) toxicity test was performed according to the OECD test guideline 202 (Daphnia sp. Immobilization Test), to determine the lethal concentrations of 5-AzaC to D. magna. Five concentrations (10, 30, 100, 300 and 1000 mg/l) of 5-AzaC were tested. The daphnids were not fed during the test. The median effective concentration (EC50), LOEC and NOEC were estimated based on a sigmoidal concentration–response curve (variable slope). A 7-day short-term screening assay [53] using synchronized adult (14–15 days old) female D. magna was then modified to investigate the sublethal effects of 5-AzaC on reproduction and molecular endpoints. Two sets of exposure studies were run simultaneously with Daphnia individuals from the same culture batch. In set 1 (reproduction test), D. magna (n = 6) with eggs visible in their brood chambers were individually exposed to 10 mg/l 5-AzaC (∼1/30 of EC50 and 1/3 of NOEC (Fig. 2)), or the culture medium (control) in 50 ml volumes (100 ml beaker, animal load: 50 ml medium/daphnid) for 7 days. The number of neonates produced was recorded every day. The neonates were removed from the test unit to avoid stress from high population density and food deficiency. In set 2 (molecular responses), 6 batches of 20 D. magna with visible eggs were exposed to 10 mg/l 5-AzaC in 1000 ml volumes of test medium (2 l beaker, animal load: 50 ml medium/daphnid) for 7 days. The neonates were removed every day for the first 6 days. After 7 days of exposure, both adults (F0) and juveniles recently released from the brood pouch (F1) were separately sampled for ChIP (10 pooled F0 and 50 pooled F1 from each replicate), global DNA methylation (5 pooled F0 and 20 pooled F1 from each replicate) and gene expression (5 pooled F0 and 10 pooled F1 from each replicate) analyses. Samples for ChIP analysis were immediately cross-linked in 1 ml 1% formaldehyde. The cross-linking reaction was quenched after 10 min by adding glycine to a final concentration of 0.125 M, and subsequently the animals were washed three times in 1 ml M7 medium and frozen in liquid nitrogen. Samples for DNA methylation analysis were snap-frozen in liquid nitrogen and stored at −80°C until use. Samples for gene expression analysis were placed in RNALater (Qiagen, Hilden, Germany) at 4°C overnight, and stored at −80°C until analysis.

Global Methylcytidine and Hydroxymethylcytidine Analysis

Standards for global methylation analysis were procured from different sources. 5-Methyl-2′-deoxycytidine (mdC) was purchased from MP Biomedicals (Santa Ana, CA, USA). 2′-Deoxycytidine, 2′-deoxyguanosine, thymidine and 2′-deoxyadenosine were obtained from Sigma-Aldrich. hmdC was purchased from Carbosynth (Compton, Berkshire, UK). 2′-deoxyguanosine-13C10, 15N5 (13C15N-G) and 5-methyl-2′-deoxycytidine-d3 (DmdC) and 5-hydroxymethyl-2′-deoxycytidine-d3 (DhmdC) were purchased from Toronto research chemicals (North York, Canada). To form six replicates, DNA was purified from pooled F0 (n = 6) and F1 (n = 6) D. magna individuals using the Gentra tissue DNA purification kit according to the manufacturers recommendations (Qiagen). Global DNA methylation analysis was performed as described by Kamstra et al. [54]. In short, 200 ng of DNA in 10 µl TE buffer was digested by adding 10 µl of a mixture of venom phosphodiesterase (60 mU/ml), alkaline phosphatase (40 U/ml) and benzonase (50 U/ml) in buffer (20 mM TRIS, 100 mM NaCl, 20 mM MgCl2, pH 7.9) (all purchased from Sigma-Aldrich). The mixture was incubated for 6 h at 37°C, followed by the addition of a mixture of internal standards (final concentrations of 345 nM 13C10-15N5-G, 20 nM DmdC), in a final volume of 200 µl. A standard curve was made within the expected mdC range for daphnids (0–1% mdC, relative to G), using 200 µl of T and dA (590 nM), dC and dG (345 nM) and the mixture of internal standards (345 nM 13C10-15N5-G, 20 nM DmdC and 0.69 nM of DhmdC).

Primer Design

The mRNA sequences of D. magna DNA (cytosine-5)-methyltransferase 1 (Dnmt1a), DNA (cytosine-5)-methyltransferase 3A homologue 1 (Dnmt3a.1), DNA (cytosine-5)-methyltransferase 3A homologue 2 (Dnmt3a.2), Gnmt, S-adenosylmethionine synthase (Metk), Mthfr, S-adenosylhomocysteine hydrolase (Sahh) and Tet2 were obtained from wfleabase [http://wfleabase.org/ (15 September 2016, date last accessed)]. Gene models used [55] were constructed by the EvidentialGene method [56] based on RNA-sequencing data and chromosome assembly data. The primer sequences of the reference genes Actin, Cyp and Gapdh were obtained from studies reporting these as good reference genes in D. magna [33, 57, 58]. ChIP primers were designed based on sequences around1000 bp upstream of TSSs, with the exception of Actin, where the amplicon was localized at the TSS. Primers were designed using Primer3 v4.0.0 [59, 60]. Primers were purchased from ThermoFisher Scientific (Carlsbad, CA, USA) and after RT-PCR the amplified fragments were tested for specificity with Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). In addition, Vtg1, Vtg2 [33] and Met [61] were included for analysis of gene expression to provide indications of the reproductive disruption potential. These genes were not included in the ChIP analyses due to lack of information on the promoter and coding sequences. All primers are listed in Supplementary Table S1.

Gene Expression Analysis

Total RNA was isolated from pooled F0 (n = 6 individuals) and F1 (n = 6 individuals) D. magna using ZR Tissue & Insect RNA MicroPrep™ kit in combination with on-column DNase I treatment, according to the manufacturer’s protocol (Zymo Research Corp., Irvine, CA, USA), as previously described [62]. The yield (>50 ng/μl) and purity (260/280 > 1.8) of RNA were determined using NanoDropTM ND-1000 spectrometer (NanoDropTechnologies, Wilmington, DE, USA). The integrity (clear peak and flat baseline) of RNA was assessed using RNA 6000 Nano chips on the Bioanalyzer platform (Agilent technologies). Two hundred nanograms of total RNA was reverse transcribed to cDNA using qScript™ cDNA SuperMix (Quanta BioSciencesTM, Gaithersburg, MD, USA) according to the manufacturer’s protocol. Quantitative real-time reverse-transcription polymerase chain reaction (qPCR) was performed using 400 nM concentrations of the forward and reverse primers and 5 ng template cDNA per reaction. Thermocycling conditions consisted of an initial 95°C heating followed by 40 cycles of 95°C for 20 s, 60°C for 20 s and 72°C for 30 s. As target control, an additional melting curve from 65 to 95°C with 0.5°C incremental steps was added to the thermocycling protocol. The relative expression of each gene was calculated using the Pfaffl method [63] and the target gene expression was normalized to the geometric mean expression of the housekeeping genes using the ΔΔCq method [64].

ChIP Assay

The ChIP assay was using a modified protocol originally developed for zebrafish (D. rerio) [65]. Briefly, a 1:100 dilution of protease inhibition mix (Sigma-Aldrich), 1 mM phenylmethylsulfonyl fluoride (Sigma-Aldrich) and 20 mM sodium acetate (Sigma-Aldrich) were added shortly to all solutions before use. F0 (n = 6) and F1 (n = 6) individuals were thawed in 100 µl lysis buffer [1% SDS, 50 mM Tris-HCl, 10 mM EDTA (Sigma-Aldrich)], ground in a mortar and incubated on ice for 30 min. Sonication conditions were optimized to 30 cycles on a Bioruptor Pico (Diagenode, Liege, Belgium), yielding fragmented chromatin of 200–300 bp. The sonicated chromatin was transferred to a QIAshredder column (Qiagen) and centrifuged at 5600 g. The low molecular weight chromatin was subsequently separated from the high molecular weight fraction by centrifugation at 12 000 g at 4°C for 10 min. A 160 µl of the low molecular weight fraction (supernatant) was transferred to a new vial while avoiding transfer of the pellet. DNA concentrations were measured by Qubit (ThermoFisher, Waltham, MA, USA) and diluted in radioimmunoprecipitation assay buffer [10 mM Tris-HCl pH 7.5, 140 mM NaCl, 1 mM EDTA, 0.5 mM EGTA (Sigma-Aldrich), 1% (vol/vol) Triton X-100 (Sigma-Aldrich), 0.1% (wt/vol) SDS (Sigma-Aldrich), 0.1% (wt/vol) Na-deoxycholate (Sigma-Aldrich)] to a DNA concentration equivalent to 3 ng/µl and frozen at −20°C until analysis. For each ChIP reaction, 100 µl chromatin was added to 2.4 µg antibody coupled to 10 µl protein A coated Dynabeads (Life Technologies AS, Oslo, Norway). The enrichment of each ChIP was analysed with qPCR using 2.5 µl of the immunoprecipitated DNA as template and quantified against known standards. The enriched DNA was then visualized as per cent of input DNA or normalized to histone H3 enrichment and visualized as relative effect compared to controls. The ChIP assay was first performed on unexposed adult daphnids with antibodies targeting H3K4me3 (Diagenode, prod. no. C15410003), H3K27me3 (Diagenode, prod. no. C15410069), H3K9me3 (Diagenode, prod. no. C15410056) and H3K9ac (Diagenode, prod. no. C1541009) to determine the dynamics of the histone PTM enrichment pattern, which was subsequently used to choose the histone PTM for evaluating 5-AzaC exposure. For the transgenerational 5-AzaC experiment, we performed immunoprecipitation against H3K4me3 and H3K27me3.

Statistical Analyses

All statistical analyses were carried out in GraphPad Prism (version 6.0.1, GraphPad Software, Inc., La Jolla, CA, USA). Outliers were identified and removed using the ROUT test [66]. In the reproductive test, EC50 and NOEC were calculated using the 4-parameters nonlinear regression and one-way analysis of variance (ANOVA, >2 groups) followed by Tukey’s post-hoc test, respectively. DNA methylation data was analysed using a one-way ANOVA followed by Tukey’s post-hoc test. After normalization to housekeeping genes for gene expression and H3 for ChIP, data that met the prerequisites of normality and equal variance were analysed using the Student t-test (pairwise comparison) or otherwise, the nonparametric Mann–Whitney test. A significance level of P < 0.05 was used in evaluating all statistical tests.

Click here for additional data file.