Cadmium-induced genome-wide DNA methylation changes in growth and oxidative metabolism in Drosophila melanogaster

Background Cadmium (Cd)-containing chemicals can cause serious damage to biological systems. In animals and plants, Cd exposure can lead to metabolic disorders or death. However, for the most part the effects of Cd on specific biological processes are not known. DNA methylation is an important mechanism for the regulation of gene expression. In this study we examined the effects of Cd exposure on global DNA methylation in a living organism by whole-genome bisulfite sequencing (WGBS) using Drosophila melanogaster as model. Results A total of 71 differentially methylated regions and 63 differentially methylated genes (DMGs) were identified by WGBS. A total of 39 genes were demethylated in the Cd treatment group but not in the control group, whereas 24 showed increased methylation in the former relative to the latter. In most cases, demethylation activated gene expression: genes such as Cdc42 and Mekk1 were upregulated as a result of demethylation. There were 37 DMGs that overlapped with differentially expressed genes from the digital expression library including baz, Act5C, and ss, which are associated with development, reproduction, and energy metabolism. Conclusions DNA methylation actively regulates the physiological response to heavy metal stress in Drosophila in part via activation of apoptosis. Electronic supplementary material The online version of this article (10.1186/s12864-019-5688-z) contains supplementary material, which is available to authorized users.


Background
Cadmium (Cd)-based chemicals are essential in many industries, including plastics and battery manufacturing and non-ferrous metallurgy [1]. As a result of their widespread use, large amounts of Cd have been released into the environment over many decades, causing pollution that threatens global ecosystems as well as human health [2,3]. Through the food chain, these chemicals can accumulate in organisms inhabiting contaminated environments [4], resulting in genetic damage, reduced reproductive capacity, growth inhibition, and even death [5,6].
Given their ubiquitous presence, there is an urgent need to better understand the biochemical impacts of Cd-based chemicals and develop effective detoxification mechanisms [7]. Many studies have addressed not only the repair of genetic damage caused by Cd but also apoptosis and oxidative stress [8,9]. However, there is little known about how Cd affects DNA methylation, a type of epigenetic modification that is important for gene regulation [10][11][12].
Drosophila melanogaster is considered a suitable model species for investigating biological responses to toxic chemicals [13]. Genes in D. melanogaster have many homologs in mammals including humans, with many genes being structurally and functionally conserved; however, Drosophila has the advantage of a simpler genome that makes it more amenable to studies of complex biological mechanisms [14][15][16]. Although global DNA methylation level is lower overall in the genome of Drosophila as compared to mammals, there are also fewer methylases. DNA methylation is an important epigenetic mechanism for the regulation of gene expression in development, reproduction, and stress resistance [17][18][19][20].
Although it is presumed that DNA methylation is involved in the response to Cd stress in Drosophila, there have been no detailed surveys of DNA methylation profiles following exposure to heavy metal stress and many questions remain unanswered, including the number and identity of methylated genes and how methylation affects gene expression. To address these points, in this study we used whole-genome bisulfite sequencing (WGBS) to evaluate genome-wide DNA methylation changes in D. melanogaster subjected to Cd stress. We identified many differentially methylated genes (DMGs) and demonstrated their relationship to gene expression. Our results provide evidence for the broad involvement of DNA methylation in the response to heavy metal stress in animals.

Results
DNA methylation state of the Drosophila genome WGBS yielded 35.5 Gb of raw data from six different samples (three repeats for each of the two groups) comprising about 38.2 billion nucleotides, all with Q20 values above 95% ( Table 1). The raw reads numbered more than 37.6 million among the six samples, and after removing those of low quality (i.e., those with a high number of 'N' , poly-A contamination, and contamination by adaptor sequences), at least 98% of the reads were retained and were taken as the high-quality (HQ) clean reads. Given the number of retained HQ reads, we expected an average genome coverage of about 30×. For all samples, between 63.56 and 74.60% of the HQ reads mapped uniquely to the genome, giving an average genome coverage between 27.28× and 35.67× (Table 1).
The average number of methylated cytosines detected in the Cd treatment and control groups was about 0.1% of all cytosines in the Drosophila genome. There were 12,397 methylated cytosines for CG, 9880 for CHG, and 30,678 for CHH (where H represents A, C, or T) in the treatment group ( Fig. 1a and Table 2), which was significantly lower (P < 0.05, Fisher's exact test) than in the control group (15,854,12,243, and 37,246, respectively, Fig. 1b and Fig.   1c), indicating that Cd treatment reduced global methylation levels.

Preferred sequences flanking the methylation site
We analyzed the relationship between the type of methylation and surrounding sequences by identifying the features of the 9-mer sequence around the methylation site ( Fig. 2a and b). For CHH, the Cd treatment and control groups showed identical sequence enrichment at each genomic region, with "TTG" and "TTT" being the preferred sequences around the methylation site. In the CG and CHG environment, sequences around methylation were slightly different. At the CG locus, with "TTT" and "AAAA" being the preferred sequences of the treatment group and "TTT" and "AATT" being those of the control group. Judging by such pattern, there seem to be almost equal preference for A, T, C, and G around all types of methylation sites. Thus, there doesn't seem to be any significantly enriched motifs in any of the treatments. Methylation occurred at similar sequence environments.

DNA methylation levels in different genomic regions
DNA methylation levels generally show a varied distribution across different functional regions of the genome. We examined the distribution of DNA methylation sites and found that the methylation levels in the promoter, 5′ untranslated region (UTR), exon, intron, and 3′ UTR were similar between Cd treatment and control groups (Fig. 3). The promoter region had the fewest methylation sites (0.03% of all sites), whereas those in introns accounted for over 65% of total sites (Fig. 3a, b). We used the sliding window method to examine DNA methylation levels in these five gene components. Methylation levels were similarly distributed in the treatment and control groups (Fig. 4). Compared to other genetic components, changes in methylation level were observed in the promoter region, but the overall methylation level was high. Methylation levels did not differ significantly across regions, and there was little change in the exon and intron, which showed a stable distribution of methylation marks.

DMRs and related genes
We used swDMR software with stringent parameters to identify DMRs between Cd treatment and control groups.
The methlytion signals, along with QQ-plots of P values associated with the DMRs, and the range of P values are shown in Additional file 1: Figure S1, Additional file 2: Figure S2 and Additional file 3: Figure S3. The QQ-plots shows that all dots represent observed log p-values of CG\CHG\CHG formed an almost straight line that away from the line that represent log p-values under the null expections, which indicate these DMRs are actually siginificate deviated. A total of 71 DMRs were detected throughout the genome (Additonal file 4 Table S1). To identify the methylated genes, we used the genomic localization of each DMR and information on D. melanogaster genome structure annotation to label the DMRs, and determined that they belong to 63 genes (Additional file 1: Table S1).
In the treatment group, 30 DMRs in 24 genes were hypermethylated and 41 DMRs in 39 genes were demethylated relative to the control group. Thus, the rate of demethylation was greater than the rate of hypermethylation. A box plot analysis of the DMRs showed that the methylation level was slightly lower in the treatment as compared to the control group ( Fig. 5), indicating that in addition to the number of de−/hypermethylated sites, demethylation occurred at a higher rate during Cd exposure. Moreover, the DMRs were mainly distributed in introns and exons-i.e., 21 and 29, respectively ( Fig. 6a), with most located on chromosome 3R, followed by chromosomes 3 L and X. With the exception of chromosome 2 L, hypermethylation was less frequently observed than demethylation on all chromosomes, with chromosome 2R having the lowest level of demethylation (Fig. 6b). Additionally, more DMRs were demethylated than were methylated, and methylation sites of the CHH type were mostly demethylated (Fig. 6c).

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses of DMRs
We carried out GO and KEGG enrichment analyses for all DMRs to clarify the functional significance of differential methylation (Additional file 4: Table S1). For all DMRs, the enriched GO terms were related to critical biological processes in D. melanogaster including reproduction, locomotion, development, growth, and response to stimulus,  indicating that Cd exposure affects the methylation of genes related to the basic physiology of Drosophila ( Fig. 7 and Additional file 5: Table S2). This was true for both hyperand demethylation, suggesting that DNA methylation broadly affects gene regulation in intricately connected biological processes. We also found other GO terms that were enriched such as immune system, single-organism process, and biological regulation. In the cell component and molecular function categories, organelle and catalytic activity were significantly enriched. Thus, genes regulated by DNA methylation are not limited to those directly involved in the response to Cd toxicity; instead, epigenetic modifications are associated with overall regulation of gene expression. Cd exposure altered a variety of pathways in the KEGG enrichment analysis; the top pathways are shown in Fig. 8, and included phagosome, phototransduction, and Hippo and Notch signaling pathways (Additional file 5 Table S2). In agreement with the enriched terms identified by GO analysis, these pathways are associated with reproduction and development. Thus, the results of the KEGG analysis demonstrate that multiple cellular mechanisms are activated in the response to Cd exposure and that DNA methylation is actively involved in their regulation.

Association between DMGs and differentially expressed genes (DEGs)
The GO and KEGG enrichment analyses of DMRs provided insight into the processes affected by DNA methylation in D. melanogaster in response to Cd treatment on a large scale. To identify the specific genes involved in these processes, we compared the complete gene sequences of these DMRs-that is, DMGs-with Drosophila digital expression library (DGE) data obtained under the same Cd treatment conditions as those of the present study. In aim to test whether these overlaps are meaningful, same number of genes and genomic regions were randomly picked and counted for the overlap for 25 times. When we are sure about that the random overlap will not likely to affect the results, we finally identified 1971 DEGs associated with heavy metal Cd stress in D. melanogaster, of which 37 were associated with 62 DMRs (Fig. 9 and Additional file 6: Table . This represented only a small proportion of all DEGs (1.87%); on the other hand, the fraction of DMGs was very high (59.6%), indicating that changes in DNA methylation state regulate gene expression but are not the main regulatory process in Drosophila. The observed correspondence between methylation and gene expression levels provide further evidence that DNA methylation regulates gene expression in combination with other mechanisms, and may only occur at specific sites in genes. In most of the 37 DMRs, methylation levels were negatively correlated with gene expression level, that is, methylation repressed gene expression, which in turn activated expression. Exceptions to this trend include Eip75B, a gene related to female gamete production whose expression increased with methylation level.
An analysis of the 37 overlapping genes showed that 27 of these had critical functions (Fig. 10) that were related to development and reproduction according to the GO and KEGG pathway analyses, including cnn, ssh, Act5C, pot, baz, Cdc42, Hem, Eip75B, and cv-c. We also found four genes (ade3, CG6729, Slbp, and CG8878) related to metabolic biosynthesis and 13 involved in resistance to Cd stress (cenG1A, Cyp6u1, AGO3, betaTub60D, alphaTub84B, Act79B, Act88F, CG43102, dx, Ant2, CG6470, Mekk1, CG4020, and Cdc42). These genes have binding or transferase activity and are associated with the immune system or intracellular signaling pathways, with functions in antioxidant and metal ion binding as well as resistance to external stimuli and initiation of apoptosis. Of these genes, Mekk1 has been linked to the response to Cd toxicity through positive regulation of the mitogen-activated protein kinase (MAPK) cascade, whereas Cdc42 is closely related to cell cycle arrest and apoptosis.
We focused on 12 genes for which there was a negative correlation between changes in methylation and expression patterns-namely, dx, Cdc42, CG8878, Ant2, Hem, pot, AGO3, ssh, Mekk1, Slbp, baz, and CG6470 (Fig. 10). The expression levels of these genes were upregulated (except for dx, which was downregulated) in response to Cd stress through DNA methylation. Discussion DNA methylation is an epigenetic regulatory mechanism that controls gene expression through modification of cytosine bases that alters chromatin structure and stability and DNA-protein interactions [21,22]. There is increasing evidence that DNA methylation is a mechanism in animals that allows adaptation to environmental stress or trauma [23] through controlled changes in gene expression levels [24][25][26][27].
In this study, we determined that Cd ion stress altered DNA methylation patterns in the Drosophila genome by WGBS. Although the DNA methylation rate in the genome was very low (~0.1%), it affected genes related to the stress response to Cd exposure. We also identified 71 DMRs encompassing 63 genes. In general, demethylation of genes in these regions was associated with increased gene expression in response to Cd treatment, which is contrary to previous findings that DNA methylation has a strictly inhibitory role in transcription [28,29]. It is worth noting that the demethylated genes had functions associated with essential biological processes such as development, reproduction, cellular defense and repair, antioxidant stress, and apoptosis.
Detoxification proteins are continuously synthesized in cells exposed to toxic elements. The results of our study suggest that this is regulated by DNA demethylation in Drosophila exposed to the heavy metal Cd, leading to the activation of stress resistance genes. Our results provide new evidence for the biological importance of DNA methylation and insight into how gene expression is regulated by epigenetic modifications under conditions of stress. GO and KEGG enrichment analyses can be used to analyze the functions of DEGs [30]. In this study, we carried out a functional enrichment analysis of GO terms and KEGG pathways for all DMRs [31,32] and found that DNA methylation during Cd stress affects genes that are involved in basic physiological functions. Enriched GO terms included reproduction, locomotion, development process, growth, immune system, and response to stimulus. This is in agreement with previous reports that Cd inhibits development and leads to decreased fertility [33] and reduced immunity [34,35]. Similar results were obtained by KEGG pathway analysis, which identified pathways associated with the phagosome, Hippo and Notch signaling, and phototransduction as those affected by changes in DNA methylation as a result of Cd stress; these processes and pathways are implicated in the regulation of immunity, somite development, ocular development, neurogenesis, and embryogenesis. Thus, DNA methylation can directly affect biological mechanisms such as development and immunity to counteract Cd toxicity, which has not yet been demonstrated; most previous studies have suggested that the mechanism of resistance to Cd stress in vivo is indirect, involving free radical scavenging (e.g., glutathione, heat shock protein, and metallothionein) or apoptosis.
We examined genes showing the greatest differences in expression due to changes in DNA methylation [36] and identified 27 including AGO3, Myo81F, and Cdc42 from the set of 37 DMGs overlapping with the DGEs. These 27 genes covered all the biological mechanisms identified by GO and KEGG enrichment analyses; 12 showed changes in methylation that were consistent with the change in their expression level, while 11 were upregulated as a result of demethylation following Cd treatment.
Previous studies have demonstrated that DNA methylation is implicated in development and reproduction [37,38]. Two DMGs in this study-namely, ssh and Act5C-are involved in eye and brain development, respectively. In addition, baz, Eip75B, cv-c, and cnnwhich are involved in oocyte axis specification, female gamete generation, embryonic morphogenesis, and embryo development, respectively-were also differentially methylated, indicating that epigenetic regulation of genes involved in resistance to heavy metal toxicity begins when the fertilized egg begins to form and is passed on to offspring.
We also identified four DMGs related to metabolic biosynthesis, namely ade3, CG6729, Slbp, and CG8878. Although the methylation patterns of these genes was inconsistent, all showed increased expression. These four genes are associated with purine nucleotide metabolism, nuclear-transcribed mRNA catabolism, histone mRNA metabolism, and protein modification, and their upregulation implied that Cd exposure induced base utilization, mRNA recovery, and protein synthesis rate and consequently, gene transcription and protein translation in Drosophila. This expression profile is consistent with the mechanism of stress resistance and demonstrates that it is not possible to predict changes in the regulation of gene expression based solely on changes in DNA methylation levels.
The most important findings of this study are that we identified eight genes related to the immune system and intracellular signaling that were differentially methylated by Cd treatment (Cdc42, cenG1A, CG43102, Mekk1, beta-Tub60D, alphaTub84B, Act79B, and Act88F). These genes are implicated in cell death or apoptosis and their methylation has been linked to a variety of malignancies and metabolic diseases [39,40]. Apoptosis can be activated by stressful stimuli such as Cd exposure and determines cell fate in organisms [41,42]. On the other hand, alphaTub84B and Act79B are involved in regulation of the actin cytoskeleton and mitotic spindle. Mekk1 is a previously reported Cd stress resistance gene encoding a zinc finger protein that binds Cd ions and induces mitotic arrest, triggering apoptosis. Cdc42 (also known as cenG1A or DG43102) is a well-known gene associated with cancer cell proliferation [43,44] belonging to the Rho family of small GTPases that regulate mitosis, establishment of cell polarity [45], cell migration [46], and MAPK signaling [47]. Upregulation of Cdc42 is a marker for cellular responses to external stimuli; in this study, we found that Cdc42 was demethylated, which corresponded to an increase in gene expression.
In conclusion, the results of this study demonstrate that changes in DNA methylation in Drosophila caused by exposure to Cd activate genes involved in apoptosis and other basic cellular processes. These findings provide novel insight into physiological response to heavy metal stress in a multicellular organism and a basis for the development of measures to alleviate the effects of these toxic compounds in humans and other animals.

Experimental Drosophila
The D. melanogaster line used in this study was maintained at the Institute of Genetics, School of Life Sciences, Shaanxi Normal University. The strain was  from the University of Cambridge (Cambridge, UK) and is guaranteed to have a consistent genetic background.

Cd treatment
Adult female flies were exposed to Cd at a concentration of 52 mg l − 1 and maintained on standard Drosophila cornmeal-sucrose-agar-yeast medium at 25°C ± 1°C on a 12:12-h light/dark cycle for 10 days.

Genomic DNA extraction
Genomic DNA was extracted from samples comprising approximately 50 Drosophila whole muscle tissues using the DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Genomic DNA contamination and degradation were verified by agarose gel electrophoresis [48,49]. DNA purity (OD260/280 ratio) was measured using a Nano-Photometer spectrophotometer (Implen, München, Germany) and DNA concentration was measured using a Qubit 2.0 fluorometer (Life Technologies, Carlsbad, CA, USA) [50,51].

Database preparation and quantification
The isolated genomic DNA (5 μg) was used to construct a library for the Cd treatment and control groups. The DNA was sonicated to 200-to 300-bp fragments using a S220 sonicator (Covaris, Woburn, MA). After purification, the DNA fragments were end-repaired, with an adenine added to the two 3′ termini. The cytosine methylation sequencing linker was then ligated to both ends of the DNA fragments, which were treated with the EZ DNA Methylation-Gold kit (Zymo Research, Irvine, CA, USA) and bisulfite and PCR amplified. The length of the insert was assessed using an Agilent Bioanalyzer 2100 system with a Qubit 2.0 fluorometer and library concentration was determined by quantitative PCR. The library was subjected to two-terminal sequencing of 125-bp fragments using the HiSeq 2500 platform (Illumina, San Diego, CA, USA). The image data were converted to the original sequence (sequencing reads) by base calling and were stored as a FASTQ file. The sequencing error rate and substrate content distribution along each read were analyzed with internal Perl scripts. The original reads in FastQ format were processed using Trimmomatic software by (1) removing the joint; (2) rejecting > 10% of reads containing N (unknown substrate); and (3) removing low-quality reads (low quality score with Phred score ≤ 20). At the same time, the Q20, Q30, and GC contents of the data were calculated using internal scripts.

Read mapping to the reference genome
After filtering low-quality reads, bisulfite-treated reads were aligned to the downloaded D. melanogaster genome sequence using Bismark v.0.12.5 software with default parameters (NCBI, ftp://ftp.ncbi.nih.gov/genomes/drosophila melanogaster). Prior to mapping, the D. melanogaster genome sequence and post-treatment reads were converted to a bisulfite-transformed version (C-to-T and G-to-A conversion), and the transformed genome sequences were indexed using Bowtie2 software (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml). Reads that were perfectly matched from the forward and reverse sequencing data were retained for further mapping against the reference Drosophila genome.
The methylation status of all cytosine positions in the reads was inferred from the mapping results. Identical reads that were aligned to the same location in the Drosophila genome were considered as replicate reads; these were used to summarize the sequencing depth and Fig. 10 Methylation and expression levels of 27 candidate genes that are critical for the response to Cd stress. Changes in methylation (orange bars) and expression (blue bars) are shown. Values were calculated for the Cd treatment vs. the control group. The depth of methylation and transcripts per million were used in the calculation and log2 analysis was used to standardize the data. Values less and greater than zero represent down-and upregulated genes, respectively. P values are shown for each fold change; those in white and black fonts were calculated based on methylation and expression fold change, respectively coverage of each sample. The results of the methylated extract were converted to bigWig format so that they could be viewed using the IGV browser. The unconverted sodium bisulfite was determined as the percentage of sequenced cytosines that were sequenced at the reference position in the phage genome.

DMR analysis
Based on the methylation information for each site, DMRs were confirmed using swDMR software [7,8]. The genomes of the treatment and control groups were first scanned using the sliding window method with a window size of 1000 bp and step size of 100 bp. Only windows containing more than 10 cytosine sites were retained and used to calculate the average methylation level. Those with a fold change (i.e., a difference in average methylation level) between the two samples of > 2 and > 0.1 and a P value < 0.05 (reflecting a significant change) with Fisher's exact probabilistic test were considered as potential DMRs. The above procedure was repeated until the genomes of all potential DMRs were confirmed, and their P values were corrected by the false discovery rate method (corrected P < 0.05). Thereafter, overlapping potential DMRs were subjected to one-time merging and statistical analysis, and the final merged tends were considered as alternative DMRs (Additional file 4: Table S1).
We compared the chromosomal location information of the DMR with the standard gene file of the D. melanogaster reference genome. When a DMR overlapped with a gene or functional component of a gene (such as the 5′ or 3′ UTR, exon, or intron), it was assigned to that gene and its components. The 5′ and 3′ UTR, exon, and intron positional information for each gene was obtained from the standard gene file. The promoter region contained a 2-kb region upstream of the transcription start site.

GO and KEGG enrichment analysis of DMGs
We used GOseq R package to perform GO and KEGG enrichment analyses of DMGs [52,53]. Using the hypergeometric distribution test, a corrected P value below 0.05 was set to identify DMRs significantly enriched for GO terms. KEGG pathway enrichment analysis was carried out using the whole genome as background to calculate the significance; according to the background, enriched metabolic pathways were identified based on the number of DMGs. The P-values of these GO and KEGG terms were shown in Additional file 5: Table S2.

Funding
This work was financially supported by the Excellent Doctor Innovation Project of Shaanxi Normal University (no. S2015YB03); Natural Science Foundation of Shaanxi Province, China (program no. 2017JM3032); and Fundamental Research Funds for Central Universities (program no. GK201703040). These funding bodies had no role in the design of the study; collection, analysis, and interpretation of data; or in the writing of the manuscript.

Availability of data and materials
The data that support the findings of this study are available from the repository of NCBI Sequence Read Archive (SRA) with the GenBank accession No: SRR7825323, SRR7825324, SRR7825321, SRR7825322, SRR7825319, and SRR7825320.
Authors' contributions DLG and MZ designed the experiments. RRD and XYH performed the experiments. DLG and XRY analyzed the data and prepared the figures. DLG drafted the manuscript. WG, SQX, and MZ edited and revised the manuscript. All authors read and approved the final manuscript.
Ethics approval and consent to participate The D. melanogaster line used in this study was maintained at the Institute of Genetics, School of Life Sciences, Shaanxi Normal University. The specimen collection did not ask for any ethical or institutional approvals, because D. melanogaster is not endangered or protected by any current law. The care and treatment of animals in this study were performed according to the Guideline for the Care and Use of Laboratory Animals in China.

Consent for publication
Not applicable.