Cytosine Methylation Alteration in Natural Populations of Leymus chinensis Induced by Multiple Abiotic Stresses

Background Human activity has a profound effect on the global environment and caused frequent occurrence of climatic fluctuations. To survive, plants need to adapt to the changing environmental conditions through altering their morphological and physiological traits. One known mechanism for phenotypic innovation to be achieved is environment-induced rapid yet inheritable epigenetic changes. Therefore, the use of molecular techniques to address the epigenetic mechanisms underpinning stress adaptation in plants is an important and challenging topic in biological research. In this study, we investigated the impact of warming, nitrogen (N) addition, and warming+nitrogen (N) addition stresses on the cytosine methylation status of Leymus chinensis Tzvel. at the population level by using the amplified fragment length polymorphism (AFLP), methylation-sensitive amplified polymorphism (MSAP) and retrotransposon based sequence-specific amplification polymorphism (SSAP) techniques. Methodology/Principal Findings Our results showed that, although the percentages of cytosine methylation changes in SSAP are significantly higher than those in MSAP, all the treatment groups showed similar alteration patterns of hypermethylation and hypomethylation. It meant that the abiotic stresses have induced the alterations in cytosine methylation patterns, and the levels of cytosine methylation changes around the transposable element are higher than the other genomic regions. In addition, the identification and analysis of differentially methylated loci (DML) indicated that the abiotic stresses have also caused targeted methylation changes at specific loci and these DML might have contributed to the capability of plants in adaptation to the abiotic stresses. Conclusions/Significance Our results demonstrated that abiotic stresses related to global warming and nitrogen deposition readily evoke alterations of cytosine methylation, and which may provide a molecular basis for rapid adaptation by the affected plant populations to the changed environments.


Introduction
The influences of human activities on global environments have been studied extensively in the past years. It is well documented that human activity has increased the atmospheric concentrations of greenhouse gases which have successively elevated global surface temperatures over the past decades [1][2][3]. The consequences of this temperature changes are more frequent occurrence of extreme weather and climate events leading to global environmental changes. In addition, the global nitrogen cycle has also been altered by human activities such as excessive use of nitrogen fertilizers, legume crops and fossil fuel combustion [4][5][6][7]. This added nitrogen has a profound effect on the chemistry of the atmosphere, aquatic and terrestrial ecosystems, and ultimately results in changes in global environments and ecosystems. These global environmental changes will result in species distribution shifts, behavioral changes and altered phenol-ogy, and these phenomena have been observed in diverse ecological settings [8]. Taken together, these previous studies have demonstrated that human activities caused alterations in global environments and the increase in extreme events under global change will impose episodic stress upon organisms, such as heat, drought and salt [9].
To survive, plants need to continuously adjust their genomes to external stimuli to adapt to the changing and stressful environments [10,11]. Therefore, the use of molecular techniques to investigate the mechanisms of plant responses to environmental changes has attracted a numerous research efforts. According to the Modern Evolutionary Synthesis, random genetic variations were thought to be the primary sources of heritable adaptation in response to altered environments. Indeed, some studies illustrated that genetic mutations have played a pivotal role in organismal adaptation to various abiotic stress conditions [12][13][14][15][16][17]. In recent years, however, there are increasing numbers of studies addressing how and to what extent the epigenetic alterations might contribute to the plants' ability to cope with various abiotic stresses [18][19][20][21]. These studies demonstrated that epigenetic modifications, including DNA methylation, histone modification and RNA interference, could rapidly alter the gene expression levels and chromatin structure, ultimately lead to heritable changes in biochemical, physiological and morphological traits, and some of which play critical roles in response to a particular stress condition [11,[22][23][24]. Nonetheless, the heritability of such induced methylation alterations in plants remain to be fully addressed. Therefore, several recent studies have investigated the transgenerational methylation changes, and documented that the stress-induced alterations in DNA methylation are common and some of these alterations could be stably transmitted to subsequent generations, and which may help plants to adapt to the same or similar stresses their progenitors once experienced [25][26][27].
Despite accumulated studies on how plants sense and adapt to abiotic stresses, the complexity of environmental variations often makes it difficult to distinguish the real underlying causes under laboratory conditions [28,29]. Therefore, the use of environmentally realistic conditions to address the molecular mechanisms of plants in response and eventual adaptation to abiotic stress is emerging as an important step to evaluate the roles of epigenetic variations. The obvious advantages of such experiments are that the abiotic stresses are much more amenable to be controlled and the interactions between different abiotic stresses can also be assessed.
In this study, we employed the amplified fragment length polymorphism (AFLP), methylation-sensitive amplified polymorphism (MSAP) and retrotransposon based sequence-specific amplification polymorphism (SSAP) techniques to investigate the impact of warming, nitrogen (N) addition and warming+N addition stress on natural populations of a grass species Leymus chinensis Tzvel. under the experimental field conditions. Here, we have asked the following questions: (1) To what extent the abiotic stresses like warming and nitrogen addition would affect the cytosine methylation levels and patterns of L. chinensis? (2) Would the alterations in cytosine methylation show heterogeneity among different genomic regions?

Ethics Statement
No specific permits were required for this study, because the performance of this study was in accordance with guidelines set by the Northeast Normal University, China. No specific permits were required for the described field studies, because the field is owned by Northeast Normal University and the Songnen Grassland Ecological Research Station performs the management. No specific permits were required for these locations/activities, because the location is not privately-owned or protected in any way and the field studies did not involve endangered or protected.

Plant Materials and Experimental Design
The species L. chinensis Tzvel., which belongs to the genus Leymus of family Gramineae, is a dominant species of grassland ecosystem and widely distributed in the Eurasian Steppes [30]. Previous studies demonstrated that this grass species has highly tolerance to several abiotic stresses, and exhibited remarkable plasticity in the morphological and physiological characteristics [31]. These attributes render L. chinensis as an ideal plant to study the molecular mechanisms of adaptation to abiotic stress conditions. In the present study, a field-cultivated experiment of L. chinensis was conducted at the Songnen Grassland Ecological Research Station of Northeast Normal University (44u409 N, 123u449 E). It has a semi-arid, continental climate with mean annual temperature 4.6-6.4uC and annual precipitation 280-400 mm, and soils are mixed saline and alkaline. The species L. chinensis is a clonal perennial grass with large below ground bud bank and could rapid propagate through asexual reproduction. These attributes made it possible to have all the plants used in this study being clonally propagated from a single mother plant.
In detail, twelve 364 meters plots were selected from the same field and each of them were 3 meters apart. Then, these plots were assigned to be as control (CK) and treatment with warming, N addition and warming+N addition stress, respectively. For warming treatment, these plants were warmed continuously from 2006 to 2009 with 200 cm615 cm MSR-2420 infrared radiators (Kalglo Electronics Inc., Bethlehem, PA, USA). To avoid the influence of infrared radiators on nearby plots, the specially designed reflectors behind the heating elements were used to ensure both a nearly uniform irradiation of heated plots and confinement of the heating flux to such plots. In each unwarmed control subplot, there was one 'dummy' heater with the same shape and size as the infrared radiator, suspended with the same height, to simulate the shading effects of the heater. In manipulated warming experiments, infrared radiator is now widely used to increase temperature of warming plots, and there have been a large number of papers published in top journals [32][33][34][35]. With infrared radiators, the wavelength of heater radiation is in the range 800-1100 nm, hence the heaters produce negligible visible light or photosynthetically active radiation [36]. The N addition stress was sprayed with aqueous NH 4 NO 3 (10 g N m 22 ) in mid-May every year. Within each plot, three individuals were randomly selected for further analysis and all of these samples were gathered in 2009. Then, a total of 36 individuals were collected in this study.

DNA Extraction and Molecular Marker Analyses
Fresh leaf tissue was collected from each study individual and dried by silica gel. Total DNA was extracted using the hexadecyltrimethyl-ammonium-bromide procedure [37]. For each individual, the same DNA sample was used as the starting material for the AFLP, MSAP and SSAP analyses described later.
The AFLP fingerprinting was exactly as described in Vos et al. [38], with minor modifications [39]. The MSAP protocol was performed according to Dong et al. [40]. For the SSAP procedure, genomic DNA was completely restricted with Hpa-II and Msp-I restriction enzyme (New England Biolabs, Massachusetts, USA), and simultaneously ligated to Hpa-II/Msp-I adapters (59-GAT-CATGAGTCCTGCT-39 and 59-CGAGCAGGACTCATGA-39). The restricted-ligated products were then pre-amplified with the Hpa-II/Msp-I primer (59-ATCATGAGTCCTGCTCGG-39) and retrotransponson primer (Bare-1:59-CTAGGGCATAATTC-CAACAA-39). For selective amplification, the pre-amplified products were diluted 20-fold and amplified with Bare-1 and Hpa-II/Msp-I selective primers (Table S1). All of the selective primer pairs of AFLP, MSAP and SSAP were listed in the Table  S1.
For all markers, the selective amplified products were separated on 6% denaturing polyacrylamide gels with silver staining, and only clear and completely reproducible bands were scored.

Statistical Analyses
The presence/absence of AFLP bands were determined by visual inspection and scored as 1/0. Similarly, the bands of MSAP and SSAP were first scored for presence (1) or absence (0) of EcoR-I/Msp-I and EcoR-I/Hpa-II fragments. To confirm the genetic similarity of the CK and other three treatment groups, we calculated the percentage of polymorphic loci based on the AFLP matrix. In addition, the changes of cytosine methylation levels and patterns were also estimated according to the MSAP and SSAP datasets. Then, the cytosine methylation alterations, which change their methylation status in treatment groups contrasting to CK, were divided into four types, CG hypo, CHG hypo, CG hyper and CHG hyper as described by Qi et al. [41]. In detail, the Hpa-II and Msp-I are two isoschizomers which recognize the same restriction site (59-CCGG) but with differential sensitivity to methylation modifications of either C: Hpa-II can not cleave if either of the C is methylated, whereas Msp-I can not cleave if the external C is methylated. Therefore, for a given genotype, methylation of the internal or external C at the 59-CCGG sites can be distinguished in the EcoR-I+Hpa-II/Msp-I-based MSAP fingerprinting profiles. So we hereby define these two major types of cytosine methylation at the 59-CCGG sites as CG methylation (a band present in Msp-I-digest but absent from Hpa-II-digest) and CHG methylation (a band present in Hpa-II-digest but absent from Msp-I-digest), respectively. Accordingly, four patterns of methylation types, namely, CG hypo, CHG hypo, CG hyper and CHG hyper, could be defined. Then, the percentages of the four cytosine methylation types in each of the treatment group were scored. Meanwhile, in order to distinguish the specific loci methylation alterations from the random methylation variations, we also analyzed the percentages of differentially methylated loci (DML) that are polymorphic in CK but show monomorphism in treatment groups. Because if a treatment causing targeted methylation changes at specific loci, there should occur consistently in different sibling plants within the same treatment group [26]. To avoid the bias in parameter estimation, any loci present/ absent in all but one individual were removed from the datasets. Thereafter, these DML of MSAP were excised from the gel and sequenced on an ABI 3730 automatic DNA analyzer (Applied Biosystems, California, USA). In addition, the bands absent in all of individuals in CK but showed monomorphic in treatment groups were also scored and sequenced. The procedure of isolated bands was performed as described by Sha et al. [42], and sequences homology searches were performed using the BLAST program (http://blast.ncbi.nlm.nih.gov/Blast.cgi). Furthermore, analysis of molecular variance (AMOVA) was employed to apportion the variation both within and between the four groups.

Genetic Background of Control and Treatment Groups
Although the L. chinensis plant populations we used in this study were clonal and all derived from a single mother plant, both original genetic heterogeneity and newly arisen genetic changes may occur and may interwoven with the stress-induced epigenetic alterations that we aimed to stress. Thus, we first used the genetic marker AFLP to assess the genetic similarity of control and treatment groups. According to the results of AFLP, 14 primer combinations produced a total of 896 bands, with an average of 59.7 bands per primer pair (Figure 1 and Table S1). As expected, only 11 loci (represent 1.2% of the total bands) showed polymorphic bands across the entire dataset. The AFLP results indicated that there was no obvious genetic differentiation among the four groups. This result enables us to presume that all of these materials used in this study lacked appreciable genetic heterogeneity and the levels and patterns of methylation alterations to be detected should represent pure epigenetic variations.

Levels and Patterns of Methylation Alteration in the Treatment Groups
To assess the levels and patterns of methylation alteration in the treatment groups, the MSAP and SSAP techniques were employed in this study. The data showed that, the 15 primer combinations assayed in the MSAP generated a total of 814 bands, with an average of 54.3 bands per primer pair ( Figure 2 and Table S1). For SSAP analysis, a total of 642 bands were obtained from the 10 primer combinations, with an average of 64.2 bands per primer pair ( Figure 3 and Table S1). These markers enabled us to assess the stress-associated epigenetic changes in the four types of cytosine methylation alterations (CG hypo, CHG hypo, CG hyper and CHG hyper). The results (Figure 4) indicated that, although the three treatment groups had similar frequencies of methylation variation patterns, the types of CG hyper and CHG hyper are higher than the CG hypo and CHG hypo in most of these comparisons. This indicated that the treatment groups have hypermethylated genomes than the untreated control. In addition, the percentages of all four patterns (CG hypo, CHG hypo, CG hyper and CHG hyper) in SSAP analysis were obviously higher than those in MSAP. This result indicated that the regions around the transposable elements had higher frequencies of cytosine methylation changes than other regions of the genome. Furthermore, comparing the percentages of the differentially methylated loci (DML) revealed that the warming+N addition group had the highest proportion of DML and the warming group the least ( Figure 5). Similarly, we also found that the percentages of DML in SSAP apparently exceeded those in MSAP, suggesting higher rates of cytosine methylation alterations adjacent to the transposable elements.
To evaluate whether the treatments caused epigenetic differentiation, AMOVA was performed in this study based on both the MSAP and SSAP datasets ( Table 1). As the MSAP marker revealed, 19.1% of the total variance presented among the four groups, it implies that the four groups had obvious differences in cytosine methylation patterns (F ST = 0.19). In contrast, only a small amount of variation (6.76%) occurred among the four groups in SSAP analysis, suggesting that the regions around the transposable elements of the four groups have diverged only slightly in cytosine methylation variation patterns (F ST = 0.07).

Sequence Analysis of Differently Methylated Loci
Thirty-six differentially methylated DNA fragments of MSAP were recovered, cloned and sequenced, of which, 7 fragments showed significant homology to known sequences or genes (Table 2). Specifically, fragments F1, F21 and F22 were homologous to genomic sequences of Oryza sativa. F5, F16 and F20 showed similarity to Hordeum vulgare subsp. vulgare and F19 shared homology with sequences of Triticum aestivum. Additionally, fragments F1, F19 and F20 were mapped to the upstream and downstream of beta-expansin 1a precursor, glycosyltransferase and ethylene responsive transcription factor, respectively. The other fragments located on the gene body regions of tubby-like Fbox protein, UNR-interacting protein, gag-pol polyprotein and DUF295 family protein, respectively.

Discussion
Investigations of the effects of climate change on the terrestrial ecosystem have been discussed extensively during the past decades. A larger number of studies have illustrated that the changes in global environments have already not only caused distribution shift in plants and animals, but also have profound influences on the behavioral and morphological traits of species [8,[43][44][45]. For example, Pounds et al. [46] have demonstrated that 67% of the 110 or so species of Atelopus are likely to extinction and the large scale warming is a key factor in the disappearances. This conclusion is also supposed by Morin et al. [47] where have illustrated that the climate change has caused 16 North American tree species shift their distribution ranges at a continental scale. These previous studies have increased our understanding of how climate changes influence the terrestrial ecosystem.
In recent years, the use of genomic techniques to address how species adjust their genetic constitution to adapt to the changing environments constitutes an emerging fundamental issue. For instance, Manel et al. [48][49][50] illustrated that changes in temperature and precipitation conditions can caused the adaptive genetic variation in alpine plants at broad scale. Meanwhile, an array of studies have also implicated that heritable epigenetic variation could directly or indirectly contribute to the abilities of species to cope with changing environmental conditions [20,28,[51][52][53][54][55][56][57]. These studies not only revealed the molecular mechanisms of plants in response and ultimate adaptation to changing environmental conditions, but also supplied a wealth of information towards our understanding of the evolutionary roles of epigenetic mechanisms.
In the present study, we have used multiple molecular markers including AFLP, MSAP and SSAP to investigate how natural populations of the perennial grass species L. chinensis alters its cytosine methylation status in response to warming, N addition and warming+N addition stresses. We presume that the warming factor used in this study has similar effects on plants as global warming which caused by human activities. Also, the influences of N addition on plants are similar to these of nitrogen deposition which caused by agricultural production. The purpose of this study is to demonstrate how and to what extent these factors induce the alterations in cytosine methylation modification. Given the theoretical expectations and empirical results, epigenetic variations may usually be affected by the genetic composition. Therefore, it is difficult to distinguish the epigenetic from genetic variations in genetically diverse individuals and populations [26,58]. In considering the effects of genetic variations, we used plants all propagated asexually from a single mother plant. In addition, we used the AFLP marker to assess the genetic similarity of both the control and treatment groups. As revealed in our results, although the materials used in this study were treated with abiotic stresses for four consecutive generations, only 1.2% of a total of 896 bands showed genetic polymorphism, indicating that all the samples used in this study shared the same genetic heritage. As such, the alterations of cytosine methylation status of each individual retrieved from MSAP should be independent of genetic variations.
The observation described above enables us to evaluate the warming, N addition and warming+N addition stress-induced changes in the species L. chinensis. By analyzing the hyper-and hypomethylation of CG and CHG, we found that all the three   treatment groups showed alterations in cytosine methylation pattern. Similar results were also found in some previous studies which illustrated that abiotic stresses could alter the cytosine methylation status through hypo-and hypermethylation of DNA [59][60][61][62][63]. Together, our findings suggested that the changes of cytosine methylation states of each treatment group might be induced by the warming and N addition. Interestingly, in our comparisons we noted that the levels of cytosine methylation alteration around the transposable elements were obviously higher than the other genomic regions. Activation of transposable elements in response and subsequent adaptation to stress has been documented previously [64,65]. Also, several studies have illustrated that the abiotic stresses, including drought, temperature and nutrient-deficiency, could activate the transposable elements as a result of altered cytosine methylation states [66][67][68][69][70][71][72][73]. However, transposable elements are highly unstable and stress treatments may induce both genetic and epigenetic variations. To distinguish the epigenetic from genetic variations, we applied the AFLP technique to analyses the genetic structures of both the CK and treatment groups, and our results illustrated that all the samples used in this study shared a similar genetic background. Therefore, we assume that most of these variations detected by SSAP were retrieved from epigenetic variation. Accordingly, our results have not only verified the previous findings, but also further demonstrated that the different genomic regions of L. chinensis showed differential propensity for cytosine methylation alterations.
As illustrated by our experiment and some previous studies, that the status of cytosine methylation could be modified by abiotic stresses, and some of these induced changes are faithfully transmitted to offspring [26]. More recently, however, several studies in Arabidopsis demonstrated that the alterations in cytosine methylation modification could occur spontaneously without the pressures of stress, and some of these transgenerational cytosine methylation changes could also generate new allelic states that underlie changes in gene transcription, and hence, providing a novel mechanism for phenotypic variation [74,75]. Therefore, in this study, to confirm whether the cytosine methylation alterations were triggered by abiotic stress and distinguish the targeted methylation changes at specific loci from random methylation alterations, the percentages of DML in SSAP and MSAP of each treatment group were scored. As revealed by our analysis, all of these treatment groups, including DML and the percentages of SSAP, were obviously higher than those in MSAP (Fig. 4). These results are concordant with the hyper-and hypomethylation analyses in this study. In general, our findings are in agreement with the study of Taraxacum officinale that the cytosine methylation patterns showed high specificity between the control and treatment groups [26]. Furthermore, according to the gene function analysis, some DML showed significant homologies to known sequences or genes ( Table 2). Take the case of fragment F5, for example, it shared homology with the tubby-like F-box family protein which may participate in the abscisic acid (ABA) signaling pathway [76]. In addition, the fragment F20 was homologous to the 39 regulatory region of the ethylene-responsive transcription factor. Both of these genes could contribute to the abilities of plants to cope with abiotic stresses [77]. Taken together, our findings implied that these DML might play important roles in adaptation by the stressed plants to the particular abiotic stress conditions.
The aforesaid analyses illustrated that the warming, N addition and warming+N addition stresses had a clear impact on the cytosine methylation status of the stressed L. chinensis plants.  However, whether the abiotic stresses could trigger specific methylation changes and then lead to epigenetic differentiation between the control and treatment populations remained to be addressed. However, it is notable that in comparing the cytosine methylation patterns, we found that the control and treatment groups indeed showed epigenetic differentiation (F ST = 0.19 for MSAP and F ST = 0.07 for SSAP). Variation in methylation levels and patterns among natural populations has been reported in several previous studies [20,28,78,79]. These studies revealed that the cytosine methylation status could be triggered by exposure to different environmental conditions, and then lead to rapid epigenetic differentiation. Nonetheless, natural populations inherently harbor genetic variations, and hence, making the single out of epigenetic variations difficult. As illustrated in our analysis, all plants used in this study, being derived from a single mother plant via asexual propagation, shared nearly the same genetic heritage, and hence, the observed epigenetic variations were not associated with genetic variations. It is therefore possible that these pure epigenetic differentiations were occurred de novo as a result of the abiotic stresses.
In conclusion, our findings have clearly shown that warming, N addition and warming+N addition stresses readily induced cytosine methylation changes and the extent of which was variable across different genomic regions. It is reasonable to assume that this flexibility of stress-induced epigenetic modifications are consequential to enhanced adaptation by the stressed plants, and which may bear relevance to global warming and nitrogen deposition. Table S1 The selective amplification primer pairs of the amplified fragment length polymorphism (AFLP), methylation-sensitive amplified polymorphism (MSAP) and retrotransposon based sequence-specific amplification polymorphism (SSAP) used in this study. NB represents the number of bands. (DOCX)