The Dynamics of DNA Methylation in Maize Roots under Pb Stress

Plants adapt to adverse conditions through a series of physiological, cellular, and molecular processes, culminating in stress tolerance. However, little is known about the associated regulatory mechanisms at the epigenetic level in maize under lead (Pb) stress. Therefore, in this study, we aimed to compare DNA methylation profiles during the dynamic development of maize roots following Pb treatment to identify candidate genes involved in the response to Pb stress. Methylated DNA immunoprecipitation-sequencing (MeDIP-seq) was used to investigate the genome-wide DNA methylation patterns in maize roots under normal condition (A1) and 3 mM Pb(NO3)2 stress for 12 h (K2), 24 h (K3) and 48 h (K4). The results showed that the average methylation density was the highest in CpG islands (CGIs), followed by the intergenic regions. Within the gene body, the methylation density of the introns was higher than those of the UTRs and exons. In total, 3857 methylated genes were found in 4 tested samples, including 1805 differentially methylated genes for K2 versus A1, 1508 for K3 versus A1, and 1660 for K4 versus A1. Further analysis showed that 140 genes exhibited altered DNA methylation in all three comparisons, including some well-known stress-responsive transcription factors and proteins, such as MYB, AP2/ERF, bZIP, serine-threonine/tyrosine-proteins, pentatricopeptide repeat proteins, RING zinc finger proteins, F-box proteins, leucine-rich repeat proteins and tetratricopeptide repeat proteins. This study revealed the genome-scale DNA methylation patterns of maize roots in response to Pb exposure and identified candidate genes that potentially regulate root dynamic development under Pb stress at the methylation level.


Introduction
Plants have developed the ability to adapt to environment changes that adversely affect growth, development and reproduction through various biochemical and physiological processes. Lead (Pb), which is one of the most common pollutants in the environment, readily accumulates in soils and is then absorbed by plants, accumulating in different plant tissues, with the highest amounts generally found in root tissues [1,2]. Some plants have evolved detoxification mechanisms that result in a natural tolerance to heavy metals [3]. These species can accumulate an inordinate amount of heavy metals and inhabited heavy metal-enriched or -contaminated soil, extracting large concentrations of heavy metal Pb into their roots and translocating it to above-ground shoots to produce large quantities of plant biomass [4]. In our previous study, Pb concentrations were measured in the roots and above-ground parts of 19 inbred lines of maize seedlings [5]. Among these lines, line 9782 lacked the ability to hyperaccumulate Pb and showed increased tolerance to Pb stress in the roots and above-ground parts following growth in soil contaminated with 750 mg·kg −1 Pb, precluding the threat of Pb entry into the food chain [5]. Moreover, protein catabolic-related genes and transcription factor families, such as bZIP, ERF and GARP, accumulated predominantly in the maize roots during development in response to Pb stress as shown by RNA-seq [6].
Interestingly, maize inbred line 9782 is capable of modulating gene expression in response to Pb stress in a flexible and short-term manner. We speculated that this response occurs via epigenetic modifications, including reversible DNA methylation, histone modifications and chromatin remodeling, and may be associated with hereditable and transgenerational alterations in gene expression [7]. Recently, epigenetic factors, especially DNA methylation, have received considerable attention because of their potential influences on complex traits and responses to adverse environmental conditions, such as drought [8], salt stress [9], and others [10]. Epigenome modifications can occur in stress-related genes through the hypomethylation or hypermethylation of DNA at specific loci or at random loci [11]. In plants, DNA methylation occurs at both symmetric (CpG and CpNpG) and asymmetric (CpNpN) sites through the action of specific de novo and maintenance methyltransferases [12]. In recent years, genome-wide analyses of methylation patterns have aided in the detection of differentially methylated regions, and it is probable that these high-throughput methods will provide valuable information in a wide range of fields in plant biology. The genome-wide profiling of DNA methylation levels in different tissues of rice (Oryza sativa L.) had few differences in DNA methylation among vegetative tissues compared with those observed between endosperm and other tissues [13]. Robert et al. have reported the epigenome-wide inheritance of cytosine methylation variants in a recombinant inbred population in maize [13], and Eichten et al. have found epigenetic and genetic influences on DNA methylation variation in maize populations [14]. Nevertheless, the epigenetic mechanisms underlying the response of maize roots to Pb stress remain poorly understood.
The objective of the present study was to use MeDip-seq to assess genome-wide DNA methylation patterns in maize roots and to identify tissues methylated during dynamic development in response to Pb stress using four treatments, including a mock treatment (A1), Pb treatment for 12 h (K2), Pb treatment for 24 h (K3), and Pb treatment for 48 h (K4). The evaluation of the distribution of DNA methylation in the genome of this plant may reveal a large number of differentially methylated genes among these dynamic, responsive time points and identify genes involved in the regulation of the response of maize roots to Pb stress.

Analysis of Methylated DNA Immunoprecipitation-Sequencing (MeDIP-seq) Reads
In the present study, four maize roots tissues were used to generate one pooled DNA sample for each group, including a mock-treated group (A1) and those exposed to Pb1000 stress for 12 h (K2), 24 h (K3), and 48 h (K4) (see Experimental Section for details). A range of 6,955,318 to 8,616,236 raw reads was generated for the four groups. More than 80% of the reads were mapped for each group, and approximately 16% of the reads were uniquely mapped to the maize genome. The uniquely mapping reads of A1, K2, K3, and K4 covered 16.13%, 15.95%, 16.27%, and 16.36% of the maize genome, respectively. The proportions of reads uniquely mapped to CpG islands (CGIs) in A1, K2, K3, and K4 were approximately 71.16%, 66.55%, 68.95%, and 63.81%, respectively (Table S1). In addition, analysis of read distributions in different components of the genome showed that the uniquely mapped reads were mainly present in intergenic regions, which contained 80% unique reads, followed by the introns, promoters, and downstream regions. Few reads were mapped to exons, 5' UTRs and 3' UTRs ( Figure 1).

DNA Methylation Profiles of Maize Roots
To decipher the genome-wide DNA methylation profiles of the maize roots, we used the uniquely mapped reads to detect the methylated peaks and further analyzed the peak distribution in the different components of the genome through a comparison of methylation densities. We obtained 48,412, 24,599, 34,380, and 38,318 methylated peaks in A1, K2, K3, and K4, respectively (Table S2).
The majority of peaks were present in intergenic regions followed by introns and promoters. The comparison of the average methylation densities of the different components of the genome showed that the methylation levels significantly differed ( Figure 2). Among all classes, the average methylation density of the intergenic regions was the highest followed by CGIs. The intergenic regions exhibited significantly higher methylation levels than the exon and intron regions (p > 0.01). Within the gene body, the methylation density of the introns was significantly higher than those of the UTRs and exons (p > 0.01).

Distribution of DNA Methylation in CGIs
It has been reported that CGIs were associated with the majority of the annotated gene promoters. The CGIs were classified into two types based on their methylation statuses. Those containing methylated peaks were regarded as methylated CGIs, and the rest were considered as unmethylated ones. In this study, a total of 356,833 CGIs were scanned by CpGPlot software and detected in the maize genome. Of these, 223,321 were found to be methylated, and 161,777 (72.44%), 153,445 (68.71%), 162,476 (72.75%) and 164,296 CGIs (73.56%) were methylated in A1, K2, K3, and K4, respectively. In addition, most of the methylated CGIs were present in intergenic regions. Within the gene body, the exons showed more methylated CGIs than the UTRs and introns (Table S3). Furthermore, we found that methylated CGIs were enriched in intergenic regions compared with other classes (25%).

Validation of Differentially Methylated Genes by Quantitative Real-Time PCR (qRT-PCR)
To confirm the low-Pb-responsive differentially methylated genes detected by MeDIP-seq, we performed qPCR using three replicates to assess these randomly selected, differentially methylated genes associated with functional categories. The qPCR results showed that these genes were significantly differentially expressed and exhibited contrasting expression compared to the MeDIP-seq data (Figure 4), confirming that these genes were induced under conditions of Pb stress.

Promoter DNA Methylation and Gene Expression Level
We found that most of the promoter regions were associated with CpG islands and were highly methylated. It is well known that promoter DNA methylation is a repressive signal for gene transcription. We obtained gene differential expression profiles for K2, K3, K4 compared with A1 respectively, using RNA-seq [6]. In the present study, we defined the genomic regions 2 kb upstream and downstream of the gene body as the proximal promoters, and the p value of the methylation peaks was used for the methylation level measurements to detect the differentially methylated genes in K2, K3, and K4 compared with A1. We observed that gene expression levels were negatively correlated with DNA methylation in the proximal promoter regions in K2, K3 and K4, whereas there was a relatively lower level of methylation in K3 compared with K2 and K4 ( Figure 5). and K4 (C). Genes were assessed according to differential expression levels. DNA methylation level was measured by the log ratio of the p value of the methylation peaks, with each point representing the mean expression level and the relative methylation level.

DNA Methylation Profiles
The characterization of genome-wide patterns of methylation in plant systems has largely been carried out using the model organism Arabidopsis. Eichten et al. [14] have recently performed methylated DNA immunoprecipitation (ChIP) analysis to locate differentially methylated regions (DMRs) in 51 maize and teosinte inbred genotypes. However, the present study is the first to systematically compare genome-wide maize root methylation profiles in response to Pb stress. Considering that Pb and phosphorus (Pi) would be interact and precipitate in the plant roots, we appropriately increased the concentration of Pi (0.5 mM) by time-course observations of the phenotype and SOD, POD enzyme activity under Pb stress with 3 mM (Figures S2 and S3), which could avoid a Pi deficiency stress. We aimed to identify methylated genes affecting maize roots growth under only heavy metal stress. We used the MeDIP-seq method to investigate genome-wide methylation during dynamic root development in response to Pb treatment. Read distribution analysis found that uniquely mapped reads were enriched in the intergenic regions. In addition, to investigate global methylation patterns, we used Model-based Analysis for ChIP-Seq (MACS) to scan the methylation-enriched regions (called peaks) detected by MeDIP-seq. Peak distribution analysis demonstrated that the promoters were high methylated, whereas the methylation levels in gene body regions were relatively low. Methylation upstream or downstream of genes had repressive effects on gene expression [15][16][17]. Our results support this notion because the promoter-methylated genes had lower expression than those that were not found to be methylated at any component ( Figure 5). In addition, the patterns of methylation within and around protein-coding genes were consistent with those observed in previous studies [18][19][20][21]. The 5' and 3' UTRs contained high levels of methylation. Within the transcribed region, methylation was lowest near the transcription start and stop sites and increased away from these sites within the gene body ( Figure S4).

Functions of Genes Potentially Methylated in Response to Pb Stress
Plants respond to adverse conditions via a series of physiological, cellular, and molecular processes culminating in stress tolerance. Previous studies have indicated that plants have evolved a range of gene regulatory mechanisms to adapt to different stress responses that act together in various response and defense systems [7]. Transcription factors, transport proteins and some other critical genes are involved in certain signal transduction and secondary metabolite pathways and are considered to be the common stress-related transcripts activated under both biotic and abiotic stresses. In the current study, we found that a total of 140 differentially methylated genes that were identified in all three comparisons (K2 versus A1, K3 versus A1, and K4 versus A1) might contribute to the regulation of the response of maize roots to Pb stress. Among these genes, transcription factors play important regulatory roles in stress responses by regulating their target genes via binding to the cognate cis-acting elements [22]. Members of the APETELA2 (AP2), bZIP, NAC, and MYB families have been shown to play regulatory roles in stress responses and have been verified to play significant roles in controlling the expression of specific stress-related genes. In our study, in addition to AP2/ERF, bHLH, MYB, bZIP transcription factors, we also detected the differential methylation of GRAS transcription factor. It has been reported that the salt-and drought-inducible poplar GRAS protein SCL7 confers salt and drought tolerance in Arabidopsis thaliana [23]. Knight and Knight (2001) found that the transcription of bZIP, Myb, and zinc finger transcription factor are induced by Pb [24]. In our study, GRMZM2G406099 (AP2/ERF), GRMZM2G048883 (zinc finger, C2H2-like), GRMZM2G482657 (zinc finger, RING-type) and GRMZM2G062499 (leucine-rich repeat) were validated by qRT-PCR to have decreased methylation levels and thereby increased gene expression levels. It has been demonstrated that plant responses to environmental stresses, including heavy metals, may be regulated by multiple transcription factors. We also found that most of the differentially expressed transcripts were involved in signal transduction and the regulation of gene expression under Pb stress. The first group of genes included those encoding kinases, phosphatases, calcium-binding proteins and proteases that were involved in stress signal transduction. Among them, protein phosphatase (PP) participates in a type of phosphoprotein cascade, resulting in the inactivation of the phosphoprotein. PP has four subunits with different interaction partners for each subunit, whereas PP2B and PP2C are Ca 2+ -dependent. Our results showed that one well-known ABA signal transduction component, GRMZM2G401075 (protein phosphatase 2C (PP2C)-like), was commonly up-regulated [25,26]. F-box protein and U-box protein, which functions as a negative regulator of phytochrome A (phyA)-specific light signaling [27,28], are ubiquitin-related proteins that play important roles in signal transduction in maize during abiotic stress [29]. One gene encoding GRMZM2G406074 (F-box protein) was found and validated in our study. In addition, many tetratricopeptide repeat (TPR) proteins (GRMZM2G082487 and GRMZM2G061876) [30,31] and a ubiquitin-fold modifier 1 (Ufm1) [32] were also identified as commonly up-regulated and directly function as cofactors in Pb stress tolerance without transducing signals. The plant cell wall, which acts as a barrier, plays an important role in regulating heavy metal defense and detoxification by limiting metal uptake and penetration into the protoplast. Many genes involved in cell wall metabolism were found to be repressed in our study. Interestingly, a leucine-rich repeat protein, which has been implicated in cell wall synthesis [33], was present at decreased levels in our study. Two members of the glycoside hydrolase family, GRMZM5G865576 (glycoside hydrolase, family 28) and GRMZM2G145458 (glycosyl transferase, family 48), which decreased in abundance, were identified. These proteins, which are considered to be plasmodesmata-associated [34], are linked to cell elongation and cell wall formation, and the findings herein indicate their involvement in cell wall modification, cell division and growth, enabling a rapid response to Pb stress.

Seed Sterilization and Experiment Design
The seeds of maize (Zea mays) inbred line 9782 were sown on filter paper saturated with distilled water and incubated at 26 °C in the dark. Three days later, seedlings selected for uniform growth were transplanted into an aerated complete nutrient solution (see Table S5 for details) and maintained for 3 days in a growth chamber with a photoperiod of 14 h light/10 h dark at 26 °C and a relative humidity of 70%. Then, the seedlings were randomly divided to two groups as follows: CK-grown seedlings, which were grown only in half-strength Hoagland solution and Pb1000-grown seedlings, which were grown in CK + Pb1000 (3 mM Pb(NO3)2) for Pb stress.

DNA Extraction and Preparation for MeDIP-seq
All maize inbred line 9782 root samples were cleaned and immediately frozen in liquid nitrogen for further study. Four libraries were constructed using DNA extracted from the CK-grown (A1) and Pb1000-grown maize roots at 12 h (K2), 24 h (K3) and 48 h (K4) according to the results from POD and SOD assays, respectively [6]. Based on the manufacturer's protocol, genomic DNA was isolated using a TaKaRa Universal Genomic DNA Extraction Kit Ver. 3.0 (DV811A) (TaKaRa, Osaka, Japan), and then DNA quality was evaluated by agarose gel electrophoresis. DNA samples from three randomly replicated roots within each group were mixed in equal amounts to generate a pooled sample using a Quant-iT dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA). Subsequently, the four pooled-samples were sonicated to produce DNA fragments ranging from 100-500 bp in size. After end repair, phosphorylating and A-tailing with a Paired-End DNA Sample Prep Kit (Illumina, San Diego, CA, USA), the DNA was ligated to an Illumina sequencing primer adaptor. Following the manufacturer's recommendations, the fragments were used for MeDIP enrichment with a Magnetic Methylated DNA Immunoprecipitation Kit (Diagenod, Liège, Belgium), and the qualifying DNA was used for PCR amplification. Bands between 220 and 320 bp in size were excised from the gel and then purified and with quantified QIAquick Gel Extraction Kit (Qiagen, Valencia, CA, USA) and Quant-iTTM dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA), respectively, on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Following qPCR, DNA libraries were sequenced with an Illumina Hiseq 2000 (Illumina, San Diego, CA, USA) to generate paired-end 50-bp reads by the Beijing Genomics Institute (Bioss, Beijing, China).

Bioinformatic Analysis
First, adapter sequences were removed, and low-quality tags and contamination due to adapter-adapter ligation were filtered out. Next, sequence reads for each tissue were mapped to v2 of the B73 reference pseudomolecules [35] using Bowtie version 0.12.7 [36]. The uniquely mapped data were retained for read distribution analysis, including assessments of the distribution among maize chromosomes and among the different components of the genome. 5b gene annotation information was downloaded from the maize sequence [35] and the region from transcript start site to transcript end site was defined as gene body region. CpG islands (CGIs) were scanned by CpGPlot [37] with the following criteria: length of 500 bp, GC content of 55%, and observed-to-expected CpG ratio of 0.65. Then, genome-wide methylation peak scanning was conducted using the MACS V 1.4.2 [38,39]. The number of peaks in different components of the maize genome (such as promoters, 5' UTR, 3' UTR, exon, intron, intergenic regions, CGIs, and downstream regions) was analyzed in our study. Moreover, we also analyzed the total peak number in each sample. Overlapping peaks for the different components of the genome were counted as a single peak. The methylation densities of the different components of the genome were compared by calculating the ratio of methylated peaks in a particular component to the total area of that region.

Reverse Transcription, Standard and Real-Time Reverse Transcription PCR
To validate the common differentially methylated genes (DMEs) in the roots as determined by MeDIP-seq, 8 DMEs were subjected to quantitative real-time PCR using an ABI7500 system. 18s rRNA was used as an endogenous control, and cDNA synthesis was carried out using 1 μg of total RNA. The corresponding primers were designed by Primer5 software and are listed in Table S6. According to the standard ABI7500 system protocol, amplification was performed as follows: 40 cycles of 95 °C for 30 s; 95 °C for 5 s, and 60 °C for 30 s, particularly for the verification of amplification specificity, followed by a thermal denaturing step to generate melting curves. All reactions were run in triplicate, including non-template controls. The threshold cycles (Ct) of each tested gene were averaged for triplicate reactions, and the values were normalized according to the Ct of the control products of the 18s rRNA gene. Statistical analysis was performed using the 2 −ΔΔCt method.

GO Annotation of All Genes with Peaks
Differentially methylated genes with peaks were used for the subsequent gene ontology (GO) analysis. Genes exhibiting more than 2-fold changes in methylation levels in the different samples were annotated, with a p < 0.005 and Benjiamini-adjusted p < 0.05, GO functional analysis of the putative target genes was performed by Web Gene Ontology Annotation Plot (WEGO) [40].

Conclusions
In summary, this study provided a comprehensive analysis of DNA methylation profiles of maize roots and revealed 140 differentially methylated genes that might be involved in response to Pb stress. Our observations provide new clues for elucidating the epigenetic mechanisms of the response of maize to Pb stress, and also provide a foundation for the studies of other types of heavy metal stress in plants.