Skip to main content

Advertisement

High-resolution DNA methylome reveals that demethylation enhances adaptability to continuous cropping comprehensive stress in soybean

Abstract

Background

Continuous cropping stress involves such factors as biological barriers, allelopathic autotoxicity, deterioration of soil physicochemical properties, and soil fertility imbalance and is regarded as a kind of comprehensive stress limiting soybean yield and quality. Genomic DNA methylation is an important regulatory mechanism for plants to resist various environmental stresses. Therefore, it is especially worthwhile to reveal genomic methylation characteristics under stress and clarify the relationship between DNA methylation status and continuous cropping stress adaptability in soybean.

Results

We generated a genome-wide map of cytosine methylation induced by this kind of comprehensive stress in a tolerant soybean variety (Kang Xian 2, KX2) and a sensitive variety (He Feng, HF55) using whole-genome bisulfite sequencing (WGBS) technology. The expression of DNA demethylase genes was detected using real-time quantitative PCR (qRT-PCR). The functions of differentially methylated genes (DMGs) involved in stress response in biochemical metabolism and genetic information transmission were further assessed based on Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The results showed that genomic DNA demethylation was closely related to continuous cropping comprehensive stress adaptability in soybean, which was further verified by the increasing expression of DNA demethylases ROS1 and DML. The demethylation of mCpG and mCpHpG (mCpApG preferred) contexts was more critical, which mainly occurred in gene-regulatory regions at the whole-chromosome scale. Moreover, this kind of stress adaptability may be related to various stress responders generated through strengthened glucose catabolism and amino acid and fatty acid anabolism, as well as fidelity transmission of genetic information.

Conclusions

Genomic DNA demethylation was closely associated with continuous cropping comprehensive stress adaptability, highlighting the promising potential of screening continuous cropping-tolerant cultivars by DNA methylation index and further exploring the application of DNA demethylases in soybean breeding.

Background

Soybean (Glycine max [L.] Merr.), an agricultural product used for grain, cooking oil, fodder, and important industrial raw materials, is a continuous global staple crop [1, 2]. Soybean plants are also important for soil fertility because they can fix atmospheric nitrogen through symbiosis with microbes in the rhizosphere [3]. However, due to salinization, desertification, the growing population, and other reasons, the area of arable land has decreased considerably over the last few decades [4, 5]. The increasing demand for soy products and reduced cultivated land acreage have resulted in large areas of soybean coming under continuous cropping stress, especially in China [6,7,8]. For instance, the acreage devoted to soybean cultivation under continuous cropping accounted for more than 40% of the whole soybean planting area in 2012 in Heilongjiang Province, Northeast China [9]. After long-term continuous cropping, the crop may have poor growth due to continuous cropping obstacles including biological barriers, allelopathic autotoxicity of plants, the deterioration of soil physicochemical properties, and soil fertility imbalance, leading to low yields and poor quality [10,11,12,13]. Therefore, the obstacle of continuous cropping, a kind of comprehensive adversity, has been one of the bottlenecks restricting soybean yield increases and quality improvement.

When crops are exposed to stressful conditions, they will resort to various strategies to minimize the effects of stress, such as tolerance, resistance and avoidance. These strategies usually arise from changes in related gene expression [14, 15]. DNA methylation is an indispensable epigenetic mechanism for normal plant development under adverse conditions that can result in stable alterations in gene expression without changes in the underlying DNA sequence [16,17,18,19].

In plants, DNA methylation commonly occurs at cytosine sites (where a methyl group is added at the 5′ position to form 5-methylcytosine) in either symmetrical CpG and CpHpG sequence contexts or asymmetrical CpHpH (where H is A, C, or T) contexts [20, 21]. Cytosine methylation in all sequence contexts is established through the RNA-directed DNA methylation (RdDM) pathway guided by 24-nt small interfering RNAs (siRNAs), in which the DNA methyltransferase Domains Rearranged Methyltransferase 2 (DRM2) is recruited to mediate de novo methylation. CpG and CpHpG methylations are maintained during subsequent rounds of DNA replication because of their symmetrical nature by DNA methyltransferase 1 (MET1) and chromomethylase 3 (CMT3), respectively, whereas asymmetrical CpHpH methylation is maintained by the RdDM pathway and chromomethylase 2 (CMT2) [22]. DNA methylation can also be removed by either passive DNA demethylation (failure to maintain methylation after replication) or active DNA demethylation (active removal by some enzymes). In plants, active DNA demethylation is mediated by members of the bifunctional DNA glycosidase subfamily, including Demeter (DME), Repressor of Silencing 1 (ROS1) and Demeter-like (DML). These enzymes can not only catalyse the hydrolysis of a glycosylic bond between the methylcytosine base and deoxyribose but also cleave the DNA backbone at abasic sites to form a single-nucleotide gap that will be filled with an unmethylated cytosine nucleotide by polymerase and ligase [21, 23,24,25,26,27,28].

By these DNA methylation or demethylation processes, DNA methylation can be dynamically regulated and maintained at a proper level in plants. When crop plants encounter environmental stresses, genomic DNA methylation levels will be changed to adapt to the challenge [29]. Diverse environmental stresses, such as cold, aluminium, herbicide, salt, drought stress, and nutrient stress, can induce heritable alteration in DNA methylation in plants [16, 30,31,32,33,34]. For example, cold stress causes strong genome-wide DNA demethylation in maize seedlings, and the transcription of some demethylated genes increases in response to cold [29, 35]. Salinity stress also induces DNA demethylation events in a tolerant Setaria italica L. cultivar, and the expression of stress-responsive genes is modulated [36]. However, in Medicago truncatula, salinity stress increases DNA methylation level as a stress-adaptive response [30]. Consequently, it is reasonable for us to assume that methylation alteration might be important for soybean plants to adapt to the comprehensive stress of continuous cropping. Up to now, there has been no report discussing the relationship between soybean continuous cropping adaptation and genomic DNA methylation.

In this study, we generated DNA methylomic maps of soybean varieties sensitive (He Feng 55, HF55) and tolerant (Kang Xian 2, KX2) to continuous cropping and examined their methylation changes induced by continuous cropping comprehensive stress using whole-genome bisulfite sequencing (WGBS) technology, which can accurately quantify whole-genome methylation at single-base resolution. The methylation levels of total C, CpG, CpHpG, and CpHpH contexts; their distribution characteristics in the genome, on every chromosome and in specific gene regions; and their changes under continuous cropping comprehensive stress were revealed. Based on demethylation obtained in this work, the expression levels of DNA demethylase genes, including ROS1, DME and DML, were detected by real-time quantitative PCR (qRT-PCR). Moreover, analysis of Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of differentially methylated genes (DMGs) were performed to identify some genes involved in metabolic processes and fidelity transmission of genetic information associated with response to continuous cropping comprehensive stress. This work made it clear that genomic DNA demethylation was closely associated with continuous cropping comprehensive stress adaptability, reinforcing the understanding of DNA demethylation knowledge applied to continuous cropping-tolerant cultivars breeding in soybean. To the best of our knowledge, this report is the first detailing the soybean genomic DNA methylation characteristics induced by continuous cropping comprehensive stress, corresponding demethylase activity, and possible metabolic mechanism and gene regulation caused by demethylation related to this kind of comprehensive stress adaptation.

Results

Bisulfite sequencing reveals that demethylation is closely related to continuous cropping resistance in soybean

Before analysing genomic methylation, we verified the different adaptabilities of sensitive HF55 and tolerant KX2 by investigating some morphological indexes including plant height, leaf area, stem and leaf dry weight, nodule number, nodule dry weight and chlorophyll content under continuous cropping comprehensive stress (Fig. 1a, Additional file 1: Table S1). On the basis of those data, bisulfite-seq analysis was performed to obtain base-pair-resolution DNA methylomes using an Illumina HiSeq 2000 and the WGBS method. The WGBS library is the most comprehensive, highest-resolution method for detecting cytosine methylation (5mC) to reveal DNA methylation patterns and variation on a genome-wide scale, especially in the CpHpG and CpHpH contexts [37, 38]. Cytosine methylation of non-CpG sites is extensive (more than 30% of the total 5mC) in plant genomes [39, 40]. Two cultivars were used to construct the WGBS libraries under continuous and non-continuous cropping conditions, and 33 Gb of sequence data per sample was generated, which covered more than 80% of sequences of each chromosome and gene region (Additional file 1: Figure S1, Table S2, Table S3).

Fig. 1
figure1

DNA methylomes of continuous cropping sensitive and tolerant soybean varieties and expression analysis of some demethylase genes. a Soybean plant morphology of sensitive HF55 and tolerant KX2 under continuous cropping stress. The tolerant variety was verifiable by its phenotype 60 days after sowing. b Whole-genome methylation levels (mC/C ratios) in both varieties under continuous cropping comprehensive stress and the non-continuous cropping condition. c Expression analysis of the demethylase genes DME, DML and ROS1 under continuous cropping stress. (NCC: non-continuous cropping; CC: continuous cropping; HF55: He Feng 55; KX2: Kang Xian 2; same as below)

In this work, the whole-genome methylation levels (the ratio of methylcytosine to cytosine) of each sample were first calculated. The results showed that the methylation ratio of tolerant KX2 and sensitive HF55 was 16.78 and 18.57% under non-continuous cropping conditions, respectively. In contrast, both cultivars displayed demethylation when they were exposed to continuous cropping stress. The degree of demethylation in KX2 was 18.77%, which was higher than that in HF55 (8.35%) (Fig. 1b). These results indicate that the tolerant soybean variety KX2 had a higher ability than the sensitive variety HF55 to adjust its DNA methylation levels upon exposure to continuous cropping stress, suggesting a link between the plasticity of DNA methylation and plant performance under continuous cropping stress.

To confirm the reliability of the above results, we further detected the expression of some important DNA demethylase genes, including DME, DML and ROS1, by qRT-PCR. The results revealed that these genes were all up-regulated in both KX2 and HF55 under continuous cropping conditions (Fig. 1c). In particular, compared to the DML and ROS1 gene expression levels in HF55, those in KX2 were increased by 85.9 and 60.1%, respectively (both P < 0.05). Therefore, genomic DNA demethylation was enhanced by increasing DML and ROS1 expression, further indicating that demethylation was closely related to continuous cropping resistance of soybean.

Demethylation of mCpG and mCpHpG contexts is more important than that of the mCpHpH context for resisting continuous cropping comprehensive stress in soybean

The total DNA methylation levels and characteristics were further estimated at CpG and non-CpG (CpHpG and CpHpH) sites. The level of mCs in CpG dinucleotides was 62.86 and 60.08% in HF55 and KX2, respectively, which was higher than that of CpHpG sequences (40.02% in HF55 and 37.41% in KX2). CpHpH contexts only had a low methylation rate (6.92% in HF55 and 5.65% in KX2) (Fig. 2a). However, interestingly, the percentages of total methylcytosine events that occurred in the three sequence contexts were not consistent with the DNA methylation levels. Among the contexts, mCpHpH accounted for the highest proportion (53.99 and 50.83% in HF55 and KX2, respectively) (Fig. 2b), while the proportion for mCpG was only 23.37 and 24.92% in HF55 and KX2, respectively, similar to that in CpHpG (22.64 and 24.45% in HF55 and KX2). These contrasting results indicated that the methylation ratio in the CpG context was highest in the soybean genome, but in terms of the number of mCs, the mCpG and mCpHpG contexts had much lower rates than the mCpHpH context. Moreover, mCpG and mCpHpG were mainly located in high methylation level (Fig. 2c).

Fig. 2
figure2

DNA methylation status of each sequence context in two soybean cultivars under different conditions. a Percentage of mCs at CpG, CpHpG, and CpHpH sites (H = C, T, or A). b The relative fraction of mCs identified for each sequence context (CpG, CpHpG and CpHpH). c Distribution of mC methylation level in each sequence context. Only mCs covered by at least 5 reads were used to calculate methylation level

Under continuous cropping stress, when genomic DNA demethylation occurred in both varieties, all the CpG, CpHpG and CpHpH sites in tolerant KX2 were more demethylated than those in sensitive HF55 (Fig. 2a). Specifically, the methylation levels of CpG, CpHpG and CpHpH sites decreased by 6.03, 8.71 and 0.95%, respectively, in tolerant KX2, all of which were higher than those of sensitive HF55 (2.07, 1.90 and 0.77% at CpG, CpHpG and CpHpH sites, respectively). Therefore, the demethylation action of mCpG and mCpHpG contexts may be a more important factor than that of the mCpHpH context in resisting continuous cropping comprehensive stress in soybean.

Further, we analysed the relationship between CpHpG and CpHpH sequence contexts and methylation preference. We calculated the numbers of 9-mer sequences in which the mC was at the fourth position (Fig. 3). In non-methylation sequence contexts, CpTpG and CpApT were the most frequent, and there was no difference between HF55 and KX2. However, in the methylated CpHpG context, mCpApG was the most frequent in tolerant KX2. Under continuous cropping stress, the mCpApG ratio decreased, while the mCpTpG ratio became dominant, showing that mCpApG demethylation was more important than demethylation of other contexts for the continuous cropping adaptation of tolerant soybean. The methylated CpHpH context was very different from non-methylated sequences. The CpT dinucleotide was the most enriched upstream of mC, followed by ApA, indicating that surrounding sequences may also be important in determining CpHpH methylation.

Fig. 3
figure3

Sequence preferences of methylated and non-methylated CpHpG and CpHpH contexts in two soybean cultivars under different conditions. The methylated cytosine is in the fourth position

Soybean chromosomal DNA demethylation exhibits whole-scale regular characteristics under continuous cropping comprehensive stress

We generated a set of methylation data for each of the twenty chromosomes in the four samples, which revealed the methylation characteristics at the chromosome level. As shown in Table 1, the total methylation levels of C, CpG, CpHpG and CpHpH sites were lowest on chromosome 13 and highest on chromosome 18 in both KX2 and HF55. Moreover, the methylation levels of all chromosomes in tolerant KX2 were lower than those in sensitive HF55. Under continuous cropping stress, cytosine (C, CpG, CpHpG and CpHpH) demethylation occurred on all chromosomes of both tolerant KX2 and sensitive HF55, and tolerant KX2 was more demethylated on all chromosomes than sensitive HF55. Specifically, total C demethylation ranged from 2.68 to 3.56 in KX2 but only from 1.15 to 1.8 in HF55, and chromosome 6 contained the most total C demethylation. In addition, in KX2, the demethylation of CpG and CpHpG contexts (5.22–7.61 and 6.80–10.08, respectively) on all chromosomes was much higher than that of the CpHpH context (0.85–1.02) under continuous cropping stress, and the CpHpG demethylation rate was the highest.

Table 1 Methylation characteristics at the chromosome level in sensitive and tolerant soybean varieties under different conditions

To characterize the distribution of all kinds of C methylation at the chromosome scale in greater detail, we generated dot plots of average C methylation levels in sliding 100-kb windows along each chromosome (Fig. 4, Additional file 1: Figures S2-1–S2-4). Moreover, we acquired the information for gene models and centromere location from the Soybean Genome Browser at SoyBase.org http://soybase.org/gbrowse/cgi-bin/gbrowse/gmax1.01/. [41], and then matched the concrete methylation characteristics in the CpG, CpHpG and CpHpH contexts on all chromosomes. When comparing all kinds of cytosine methylation levels based on the density of genes, we found that DNA methylation was roughly negatively correlated with the density of genes. The density peaks of DNA methylation were most likely to be located in the regions containing few genes, which were mainly located in the pericentromeric region. In contrast, chromosome end regions containing more genes showed the opposite pattern. Further, most chromosomes had a higher methylation level in the telomeric region in the CpG context than in the other two contexts (Additional file 1: Figure S2), suggesting that the methylcytosine distribution at the chromosome level had regional characteristics. However, we did not find obvious, large fluctuations when the methylation distribution was compared between continuous and non-continuous cropping conditions within the same variety, indicating that the demethylation caused by continuous cropping likely occurred at the whole-chromosome scale.

Fig. 4
figure4

Methylcytosine density throughout chromosome 6 in sensitive HF55 and tolerant KX2 under different conditions. Normalized methylated cytosine over total cytosine positions in 10-kb windows (blue dots, left axis) and normalized methylated CpG, CpHpG, and CpHpH contexts in 100-kb windows (smoothed lines, right axis)

DNA demethylation induced by continuous cropping comprehensive stress mainly occurs in gene-regulatory regions

To characterize the abundant C methylation in the regions of the soybean genome, we used heat maps to represent the methylation distribution; CpG, CpHpG and CpHpH density distribution; and the relationship between the two in the regions including the 3′-UTR, 5′-UTR, CDS, and upstream, downstream, intron and genome (Fig. 5). We found that in all regions, the densities of CpG, CpHpG and CpHpH contexts increased in that order in the four samples. However, the number of high-methylation-level windows in upstream, downstream, 5′-UTR, 3′-UTR and genome decreased in the order of CpG, CpHpG and CpHpH contexts. Moreover, continuous cropping comprehensive stress also decreased the number of high-methylation-level windows of CpG and CpHpG in these regions in KX2. Therefore, we compared the average methylation density of each context in specific regions. As shown in Table 2, higher mC densities were located upstream (5.75 in KX2 and 6.89 in HF55) and downstream (5.15 in KX2 and 6.09 in HF55) than in other regions, including the genebody, coding region (CDS), 3′-UTR and 5′-UTR. Under continuous cropping stress, total C methylation levels decreased by different degrees in all regions in both varieties, especially in KX2. The levels of upstream, downstream, 5′-UTR and 3′-UTR were reduced more obviously, with reductions of 22.78, 24.08, 23.42, and 21.15%, respectively.

Fig. 5
figure5

Heat maps of density patterns and methylation distribution of CpG, CpHpG and CpHpH contexts in different genomic regions of two soybean cultivars under different conditions. Each panel represents a separate characteristic, and n represents the number of analysed CpGs, CpHpGs or CpHpHs (per-strand depth ≥ 4) within that feature. X-axis (cytosine density) indicates the number of CpGs, CpHpGs or CpHpHs in 200 bp windows. Y-axis (methylation level) displays the mean methylation level of cytosines in the specific CpGs, CpHpGs or CpHpHs. The black line represents the methylated median value for the specific density of CpGs, CpHpGs or CpHpHs. The red zone reflects the abundance of CpGs, CpHpGs or CpHpHs that fall into bins of given methylation levels and their densities. The top blue bar chart shows the distribution of CpG, CpHpG or CpHpH densities projected onto the horizontal axis of the heat maps. The right green bar chart indicates the distribution of methylation levels, projected onto the vertical axis of the heat maps

Table 2 Relative methylation density (mC/C ratios) in HF55 and KX2 throughout different gene-associated regions under different conditions

Some vital metabolism processes and related DMGs in resisting continuous cropping comprehensive stress are revealed by DMG analysis

To investigate differential methylation under continuous cropping stress, we identified differentially methylated regions (DMRs), which denote genomic regions of adjacent CpG, CpHpG or CpHpH sites that are differentially methylated. DMRs were identified in windows containing at least five CpG (or CpHpG or CpHpH) sites at the same position in two sample genomes. A total of 13,199 DMRs in tolerant KX2 and 4018 DMRs in sensitive HF55 were identified, and the hypomethylated proportion in KX2 (95.33%) was much higher than that in HF55 (65.26%) (Fig. 6a). A number of DMGs were identified in KX2 (4475) and HF55 (1951), which consisted of hypomethylated genes and hypermethylated genes. In KX2, the number of hypomethylated genes was significantly higher (5.67 times) than that of hypermethylated genes under continuous cropping stress, whereas the numbers of hypo- and hypermethylated genes were not very different in HF55 (only 1.61 times) (Fig. 6b). GO annotation revealed that these DMGs are involved in diverse biological processes, such as metabolic processes, response to stimulus, signal, transcription, macromolecular complex and biological regulation (Fig. 6c).

Fig. 6
figure6

Statistical analysis of differential methylation in sensitive HF55 and tolerant KX2 under different conditions. a Number of differentially (hypo- and hyper-) methylated regions (DMRs); b Number of differentially (hypo-, hyper-, hypo- and hyper-) methylated genes (DMGs); c Gene Ontology (GO) categories significantly enriched in the DMRs related to stress

Detailed pathway-based analyses showed that DMGs induced by continuous cropping comprehensive stress in tolerant KX2 are mainly involved in the processes of glucose, amino acid and fatty acid metabolism. As shown in Fig. 7 and Table 3, the HK, aceE, CS, IDH, and OGDH genes encode hexokinase, pyruvate dehydrogenase E1 component, citrate synthase, isocitrate dehydrogenase and 2-oxoglutarate dehydrogenase E1 component, respectively, which are the key enzymes for glucose catabolism. Glucose catabolism can supply energy, reductants and some materials for amino acid and fatty acid synthesis. The amino acid metabolism-related genes ALT, glyA, cysK, metE, DNMT1, and OTC encode alanine transaminase, glycine hydroxymethyltransferase, cysteine synthase A, 5-methyltetrahydropteroyltriglutamate-homocysteine methyltransferase, DNA (cytosine-5)-methyltransferase 1, and ornithine carbamoyltransferase, respectively. These enzymes may catalyse glutamate, glycine, cysteine, methionine, S-adenosylhomo-cysteine (SAH) and citrulline synthesis, which will contribute to the synthesis of polyamines (putrescine, spermidine, and spermine), GSH and phytochelatins (PCs). The GSR and FNR genes encode glutathione reductase (NADPH) and ferredoxin-NADP(+) reductase, which can enhance GSH regeneration. GSH is actively resistant to stress because it is an important antioxidant, which can protect tissues from peroxide damage. Further, the genes glyA, GLDC, and AMT can catalyse one-carbon-unit N5,N10-methylene-THF synthesis from Ser and Gly, which are involved in nucleotide synthesis. The ACACA and fadB genes encode acetyl-CoA carboxylase and [acyl-carrier-protein] S-malonyltransferase, which will enforce the synthesis of fatty acids and their derivatives, such as jasmonic acid (JA), methyl JA (MeJA) and 12-oxo-10,15(Z) phytodienoic acid (OPDA). These compounds are all associated with plant resistance to biotic and abiotic stresses. All of these genes were differentially demethylated in tolerant KX2 but not in sensitive HF55 under continuous cropping stress.

Fig. 7
figure7

Primary metabolic pathways that respond to continuous cropping comprehensive stress. Solid lines represent one-step reactions and corresponding differentially demethylated enzyme genes. Dotted lines indicate that enzyme genes are not differentially demethylated and/or multi-step reactions. HK: hexokinase; GPI: glucose-6-phosphate isomerase; PFK: 6-phosphofructokinase 1; aceE: pyruvate dehydrogenase E1 component; CS: citrate synthase; IDH: isocitrate dehydrogenase; OGDH: 2-oxoglutarate dehydrogenase E1 component; ALT: alanine transaminase; cysK: cysteine synthase A; glyA: glycine hydroxymethyltransferase; GLDC: glycine dehydrogenase; AMT: aminomethyltransferase; metE: 5-methyltetrahydropteroyltriglutamate-homocysteine methyltransferase; DNMT1: DNA (cytosine-5)-methyltransferase 1; OTC: ornithine carbamoyltransferase; ACACA: acetyl-CoA carboxylase; fabD: [acyl-carrier-protein] S-malonyltransferase; GSR: glutathione reductase (NADPH); FNR: ferredoxin-NADP(+) reductase; PCs: Phytochelatins; JA: jasmonic acid; MeJA: methyl JA; OPDA: 12-oxo-10,15(Z) phytodienoic acid

Table 3 Stress-induced DMGs involved in metabolism processes and their demethylated types in soybean

In addition, based on GO annotation, some stress-response genes are involved in the processes of DNA repair (MutLα, XPD, and TFIIH4), RNA surveillance (PP2A and RNGTT), spliceosome formation (U5 snRNA) and protein processing in the endoplasmic reticulum (P97) (Fig. 8 and Table 4). These processes are closely related to plant adaptability to stress. Moreover, these DMGs involved in metabolism processes and fidelity transmission of genetic information were characterized with demethylation of CpG and CpHpG contexts (Table 3 and Table 4).

Fig. 8
figure8

Schematic showing the points at which soybean copes with continuous cropping comprehensive stress by metabolic and genetic regulation resulting from DNA demethylation. The continuous cropping stress signal resulted in DNA demethylation. Some metabolic and genetic processes were modified, causing the plants to adapt to the stress in the continuous cropping environment

Table 4 DMGs involved in fidelity transmission of genetic information within GO annotation and their demethylated types in soybean

Discussion

DNA methylation changes in plant genomes can generate novel and heritable phenotypic variation, which will improve tolerance, resistance and adaptation to poor environments by influencing gene expression [42, 43]. The striking phenotypic difference between tolerant KX2 and sensitive HF55 prompted us to profile their methylomes in response to continuous cropping stress. The high-resolution analysis provided unique insight into the plasticity of DNA methylation in response to the comprehensive stress of continuous cropping. We observed that the total methylation ratios of tolerant KX2 and sensitive HF55 were 16.78 and 18.57%, respectively, which are lower than those of both maize (20%) [35] and rice (24.3%) [44] but more than three times that of Arabidopsis (5.26%) [45], indicating specific methylation levels among diverse plant species [29]. Lukens and Zhan reported that a moderate cytosine methylation level played an important role in maintaining genome stability, contributing to the silencing of transposable elements and controlling the transcription of some genes in plants [42]. In this study, we found that continuous cropping comprehensive stress induced genomic DNA demethylation in soybean. The tolerant variety KX2 was able to rapidly reduce its DNA methylation level (by 18.77%), while the sensitive variety HF55 showed a low ability to adjust its DNA methylation level (decrease of 8.35%) upon exposure to continuous cropping stress, suggesting a link between the plasticity of DNA methylation and plant tolerance performance. A decrease in DNA methylation level will cause the activation of some genes, which will enhance the expression of stress-resistant proteins [32, 46,47,48]. Therefore, we infer that the tolerant variety KX2 can better adapt to the comprehensive stress of continuous cropping than the sensitive variety HF55 through further demethylation. This finding is in agreement with earlier results on DNA demethylation in response to adverse conditions, such as biotic stress [32, 49, 50], salt stress [51] and Fe deficiency [52]. The reason is that chromatin demethylation and the relaxation of its structure could serve as a transcriptional switch for many stress-regulated genes [15]. However, hypermethylation of cytosine residues has also been uncovered in pea root tips in response to water deficit stress [53] and in tobacco [54] and potato [55] cell cultures in response to osmotic stress. Therefore, epigenetic change at the DNA methylation level, whether hypomethylation or hypermethylation, plays an important role in plant adaptation to environmental stress and growth and development. Importantly, though, distinct epigenetic responses occur between different plants and in response to different stress stimuli.

DNA methylation in plants occurs in three cytosine contexts including CpG, CpHpG and CpHpH (H replacing A, C or T) [56, 57], and it appears to have various functions, including regulating the expression of some genes, reprogramming and imprinting [58]. The cytosine methylation patterns result not only from the establishment and maintenance of mCs but also from demethylation [59]. The removal of methylcytosine can be accomplished via passive or active processes. In passive demethylation, the mCs are replaced with unmethylated cytosines during DNA replication, while in active demethylation, the methyl mark is removed by 5mC glycosylases such as DME, DML and ROS1. These enzymes, possessing both glycosylase (base excision) and AP lyase activity that are directed towards mCs, are involved in base excision repair [60]. In our work, the relative expression levels of DML and ROS1 in KX2 were significantly higher than those in HF55 under continuous cropping stress. The high expression of these enzymes should be one of the factors causing extensive demethylation in genomic DNA [25]. Active DNA demethylation is important in maintaining epigenomic plasticity to enable efficient response to environmental stresses in a timely manner [17, 59]. These results further confirm that the DNA demethylation is related to the resistance of soybean to the continuous cropping comprehensive stress. Therefore, how to apply demethylase in plant resistance to this kind of stress should be researched in the future.

Cytosine bases of the nuclear genome in higher plants are often extensively methylated [40, 56]. In this work, we described the profile of DNA methylation density for each of the 20 chromosomes in 100-kb windows. Interestingly, chromosome end regions showed lower DNA methylation density, where gene density was high, suggesting that cytosine methylation occurred in the intergenic regions. Figure 4 and Additional file 1: Figure S2(-1,-2,-3,-4) also show low methylation density in the telomeric regions. These results are in agreement with previous reports on Arabidopsis telomere methylation. Vaquero-Sedas et al. found that Arabidopsis telomeres have lower levels of DNA methylation than internal transcribed spacers or subtelomeres [61]. Later, Vega-Vaquero et al. confirmed that DNA methylation is indeed absent in Arabidopsis telomeres based on experiments on high telomeric C-rich strand production efficiencies and methylation-dependent restriction enzyme analyses [62]. In addition, they pointed out that the degree of telomeric DNA denaturation during the process of sequencing or the formation of telomeric C-rich strand secondary structures such as the i-motif might cause the overestimation of telomeric methylation [63, 64]. Schoeftner and Blasco also considered that conserved telomeric repeats remained unmethylated because the asymmetric target units (CCCTAAA) n in plants lack methylable cytosines [65].

In mammals, DNA methylation predominantly occurs at cytosines in CpG sequences, while in plants, methylation of CpH sequences (CpHpG and CpHpH, where H can be A, C or T) is also present and involved in epigenetic regulation and gene expression [66, 67]. Therefore, detecting the context proximal to sites of CpH methylation is essential for determining whether stresses caused some enrichments of particular local sequences and basal changes, as previously reported in Arabidopsis DNA methylomes [45, 68]. No local sequence enrichment was observed upstream of mCpHpG sites, while the base following the methylcytosine tended to be adenine or thymine (mCpA/TpG) (Fig. 3). This trend is consistent with a nearest-neighbour analysis of wheat germ DNA that found a higher level in mCpApG and mCpTpG sites than in mCpCpG sites [69]. Further, we found that the base preference following methylcytosine in mCpHpG changed under continuous cropping comprehensive stress in tolerant KX2, and mCpTpG became more common than mCpApG. Therefore, we infer that mCpApG sequences are more easily demethylated than other sequences under continuous cropping comprehensive stress in tolerant KX2, which resulted in the decrease in the percentage of mCpApG and thus the dominance of mCpTpG. This process may regulate gene expression, but further research is necessary to clarify the role of CHG demethylation, as reviewed by Pelizzola and Ecker [70]. In addition, similar to the mCpHpH context, mC also tended to be followed by an A or a T. This result is consistent with data from mammalian and Arabidopsis genomes indicating that mCpT and mCpA sequences are more frequent than mCpC sequences [45, 68, 71].

DMRs, GO annotation and KEGG pathway analysis indicated that some DMGs in response to continuous cropping comprehensive stress are involved in some vital metabolic processes (Fig. 6). HK, aceE, CS, IDH, and OGDH are key genes for glucose catabolism, which is the central metabolic pathway in all organisms. This pathway produces ATP, reductants and some materials for amino acid and fatty acid synthesis. Additionally, glucose catabolism is directly involved in adaptation of plants to environmental stresses, such as nutrient limitation, drought, low-temperature and osmotic stress [72]. The ALT, glyA, cysK, metE, DNMT1, and OTC gene products catalyse the synthesis of some important amino acids such as SAH, methionine, cysteine, glycine, glutamate and citrulline. These amino acids are necessary for the biosynthesis of some defence compounds, such as GSH, polyamines and PCs. The GSR and FNR gene products also enhance GSH regeneration. GSH has many distinct functions in plant stress defence, including removing harmful H2O2 to protect tissues from peroxide damage, controlling gene expression linked to the redox state of cells and being an important reducing cofactor of many enzymes related to reactive oxygen species (ROS) detoxification [73,74,75]. Some researchers reported that plants exposed to salicylic acid (SA) (which plays a key role in plant stress tolerance) and abscisic acid (ABA) (which is related to environmental stress adaptation) exhibited higher GSH concentrations and glutathione reductase (GR) activity than those without such exposure, which further certifies the relationship between GSH content and stress defence [76,77,78]. Moreover, GSH is also related to the synthesis of PCs, which is involved in stress response [78]. In addition, SAH is the product catalysed by DNMT1 in the S-adenosylmethionine/homocysteine cycle. Fuso and Lu et al. reported that SAH is a strong DNA methyltransferase inhibitor that will reinforce genomic DNA hypomethylation [79, 80]. Polyamine synthesis is an important nitrogen-metabolizing pathway regulating ammonia metabolism and organic nitrogen balance in plant cells [81,82,83,84,85]. During the process of polyamine generation, some intermediates, such as nitric oxide (NO) and γ-aminobutyric acid (GABA), can also be produced [86, 87]. Polyamines, NO and GABA all play important roles in the regulation of plant development and as signal molecules mediating some responses to biotic and abiotic stressors, including pathogens, heavy metal, drought and salt [85, 88,89,90,91]. The ACACA and fadB genes contribute to fatty acid synthesis. Fatty acids not only are crucial components of cellular membranes, suberin, and cutin waxes, but also are important for the remodelling of membrane fluidity, which will increase plant resistance to drought, physical injury and infection by pathogenic microorganisms [92, 93]. Moreover, fatty acids are the materials for the synthesis of substances related to environmental stress adaptation, including JA, MeJA and OPDA, which are also associated with plant basal immunity and responses to pathogens [94, 95]. Therefore, based on the analysis above and our results, we think these substances are involved in the response to continuous cropping comprehensive stress, including biological barriers, allelopathic autotoxicity, the deterioration of soil physicochemical properties and soil fertility imbalance, and DNA demethylation may be the source inducing these complex resistance mechanisms.

The stable inheritance of DNA genetic information, RNA transcription and correct protein synthesis are vital molecular processes for ensuring cell-cycle progression and various biofunctions, all of which are involved in plant growth and stress response [96,97,98,99]. In this work, we observed that some DMGs in tolerant KX2 were involved in continuous cropping comprehensive stress responses (GO annotation). These genes took part in the processes of DNA repair, transcription, RNA splicing, RNA surveillance, and protein processing in the endoplasmic reticulum (Fig. 8, Table 4). These molecular processes are capable of adjusting the levels of mRNA and functional proteins that will alter some metabolic processes and thus participate in the stress response [99,100,101]. Consequently, resistant plants can tolerate complicated environmental stresses more efficiently and effectively through metabolic networks and fidelity transmission of genetic information caused by DNA methylation changes, enhancing their adaptability under continuous cropping comprehensive stress (both biotic and abiotic stresses).

Conclusion

Genomic DNA demethylation was closely related to soybean adaptability to the continuous cropping comprehensive stress, which was further verified by the increased expression of DNA demethylases ROS1 and DML. The demethylation of mCpG and mCpHpG (mCpApG preferred) contexts was more important, which mainly occurred in gene-regulatory regions at whole-chromosome scale. Among the DMGs, GO annotation and KEGG pathway analysis further demonstrated that various stress responders generated through strengthened glucose catabolism, amino acid and fatty acid anabolism, as well as fidelity transmission of genetic information, played important roles in soybean adaptability to this kind of adversity.

Methods

Plant materials and growing conditions

The soybean cultivars used in this study included the continuous cropping-tolerant variety KX2 (breeding and provision by Daqing Branch of Heilongjiang Academy of Agriculture Science) and continuous cropping-sensitive variety HF55 (provision by Genetic Breeding Laboratory of Agricultural College of Heilongjiang Bayi Agricultural University; breeding by Hejiang Institute of Agricultural Sciences, Heilongjiang Academy of Agricultural Sciences). Continuous cropping soil was collected from Lindian County, Heilongjiang Province, China, where soybean has been cultivated continuously for 6 years. The control included neighbouring non-continuous cropping soil with physicochemical properties and fertility similar to those of the continuous cropping soil. Soybean seeds were planted in pots (diameter 25 cm, depth 20 cm) filled with continuous cropping or control soil and exposed to 25/18 °C (day/night) and 70% relative humidity for 60 days in a greenhouse. Light, temperature and humidity conditions remained constant throughout the experimental periods. Soybean plants were collected carefully after 60 days of sowing. Four plants per replicate were used for measuring plant height, total leaf area and nodule number. The roots, stems, leaves and nodules of those four soybean plants were packaged in draft paper, oven-dried at 105 °C for 30 min, and then dried at 80 °C to a constant weight. The total leaf area was determined through the dry weight ratio method. Chlorophyll content from leaf tissues was measured from different positions on two of the uppermost, youngest, fully expanded leaves using a portable chlorophyll content meter, CCM-200 PLUS (Opti-sciences, USA). Three repetitions per treatment were used.

Whole-genome bisulfite sequencing library construction and sequencing

Leaf material was excised from the uppermost third of the functional leaves of 60-d-old plants, frozen in liquid nitrogen and stored at − 80 °C until DNA isolation. There were three repetitions per treatment. Genomic DNA was extracted using a QIAamp DNA Micro Kit (Qiagen, USA) according to the manufacturer’s instructions. DNA concentration was quantified using a Qubit Fluorometer (Invitrogen, USA), and DNA integrity was detected by 1% agarose gel electrophoresis. After that, DNA of three repetitions was mixed in equal amounts for library construction.

For normal WGBS library construction, the DNA was fragmented to 100–300 bp by sonication using a Bioruptor (Diagenode, Belgium), following blunt-ending, addition of dA to the 3′-end, and ligation of adaptors to protect against bisulfite conversion. DNA fragments were treated with sodium bisulfite using a Zymo EZ DNA Methylation-Gold kit (Zymo Research, USA). After sodium bisulfite conversion, unmethylated cytosine residues are converted to uracil, whereas 5-methylcytosine (5mC) remains unchanged. Fragments with different insert sizes were excised from the same lane of a 2% TAE agarose gel. Products were purified by using a QIAquick Gel Extraction kit (Qiagen, USA) and amplified by PCR. After PCR amplification, uracil residues were converted to thymine. Finally, the qualified library was sequenced using a HiSeq 2000 platform.

Data analysis

After sequencing, the raw reads were filtered by removing adaptor sequences and low-quality reads using the Illumina analysis pipeline. During this process, due to bisulfite conversion of cytosine to uracil, cytosines on the coding strand were changed to thymidines, and guanines on the template strand were changed to adenosines. Then, the clean and high-quality reads were aligned to the reference genome Phytozome v9.0 (https://genome.jgi.doe.gov/portal/pages/dynamicOrganismDownload.jsf?organism=Phytozome) using SOAP 2.20 software. Only perfectly matched reads were used for methylation analysis. Methylation level was determined by dividing the number of reads covering methylated cytosine (mC) by the total reads covering cytosine (mC/C) [102, 103]. The ratio of mCpG, mCpHpG, or mCpHpH to total mC was used to calculate the proportion of mCpG, mCpHpG, or mCpHpH at all mC sites, respectively. DMRs were identified by comparing the methylation level difference in the same sliding window of genomes between continuous cropping samples and noncontinuous cropping samples. Only windows that contained at least five CpG (or CpHpG or CpHpH) sequences were used, and changes in methylation level after continuous cropping stress had to be at least 2-fold. P-values associated with DMRs were calculated by Fisher’s exact test, and P values < 0.05 were considered significant. All adjacent differentially methylated windows were collapsed into a single DMR. Genes containing DMRs in the genebody and/or 2-kb flanking sequences were considered DMGs. GO enrichment analysis was performed using the BiNGO tool to analyse the molecular functions of DMGs. P-values ≤0.05 after family-wise error rate correction were considered significantly enriched. KEGG pathways were used to analyse the biochemical metabolic processes that DMGs were involved in.

Expression analysis of demethylase genes

Total RNA was extracted from leaves of tolerant KX2 and sensitive HF55 soybean using an RNA extraction kit (UNIQ-10, SK1321, Sangon Biotech, China). The RNA concentration was calculated by measuring the absorbance at 260 nm, and the purity was evaluated by the ratios of 260/280 nm and 260/230 nm using a NanoDrop spectrophotometer (NanoDrop Technologies, USA). cDNA was synthesized using RevertAid Premium Reverse Transcriptase (Thermo Scientific™ EP0733) with 6-mer random primers as recommended. The synthesized cDNA was subjected to qPCR. Specific primers were designed using Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, USA) (Additional file 1: Table S4). The ACTIN gene was used as a reference gene to normalize the amount of RNA in each sample. The qPCR was performed in an ABI StepOne Plus instrument in a 20 μL reaction mixture containing 2 μL cDNA, 4 μmol forward and reverse primers, and 10 μL SybrGreen qPCR Master Mix (Takara Biotech, China). Thermal cycling conditions were as follows: pre-incubation at 95 °C for 3 min, followed by 45 cycles of denaturation at 95 °C for 7 s, annealing at 57 °C for 10 s and elongation at 72 °C for 15 s. A melting curve was analysed at the end of the qPCR to verify specific amplification. The relative expression quantity was determined by the 2–ΔΔCt method, where ΔCt = (Ct target gene − Ct actin gene) and ΔΔCt = (ΔCt treatment–ΔCt control). The experimental data were statistically analysed with the t test using SPSS software (version 25.0, SPSS Inc., USA). P values lower than 0.05 were considered statistically significant.

Abbreviations

5mC:

Cytosine methylation

ABA:

Abscisic acid

CMT:

Chromomethylase 3

DME:

Demeter

DMGs:

Differentially methylated genes

DML:

Demeter-like

DMRs:

Differentially methylated regions

DRM2:

Domains Rearranged Methyltransferase 2

GABA:

γ-aminobutyric acid

GO:

Gene Ontology

GR:

Glutathione reductase

HF55:

He Feng 55

JA:

Jasmonic acid

KEGG:

Kyoto Encyclopedia of Genes and Genomes

KX2:

Kang Xian 2

MeJA:

Methyl JA

MET1:

DNA methyltransferase 1

NO:

Nitric oxide

OPDA:

12-oxo-10,15(Z) phytodienoic acid

PCs:

Phytochelatins

qRT-PCR:

Real-time quantitative PCR

RdDM:

RNA-directed DNA methylation

ROS:

Reactive oxygen species

ROS1:

Repressor of Silencing 1

SA:

Salicylic acid

SAH:

S-adenosylhomo-cysteine

siRNAs:

Small interfering RNAs

WGBS:

Whole-genome bisulfite sequencing

References

  1. 1.

    Van Dam J, Faaij AP, Hilbert J, Petruzzi H, Turkenburg W. Large-scale bioenergy production from soybeans and switchgrass in Argentina: part B. Environmental and socio-economic impacts on a regional level. Renew Sust Energ Rev. 2009;13(8):1679–709.

  2. 2.

    Lu S, Zhao X, Hu Y, Liu S, Nan H, Li X, Fang C, Cao D, Shi X, Kong L, et al. Natural variation at the soybean J locus improves adaptation to the tropics and enhances yield. Nat Genet. 2017;49(5):773–9.

  3. 3.

    Kim E, Hwang S, Lee I. SoyNet: a database of co-functional networks for soybean Glycine max. Nucleic Acids Res. 2016;45(D1):D1082–9.

  4. 4.

    Sojneková M, Chytrý M. From arable land to species-rich semi-natural grasslands: succession in abandoned fields in a dry region of Central Europe. Ecol Eng. 2015;77:373–81.

  5. 5.

    Xie H, Zou J, Jiang H, Zhang N, Choi Y. Spatiotemporal pattern and driving forces of arable land-use intensity in China: toward sustainable land management using emergy analysis. Sustainability. 2014;6(6):3504–20.

  6. 6.

    Zhang Q, Li Z, Han B. Immediate responses of cyst nematode, soil-borne pathogens and soybean yield to one-season crop disturbance after continuous soybean in Northeast China. Int J Plant Prod. 2013;7:341–53.

  7. 7.

    Liu X, Li Y, Han B, Zhang Q. Yield response of continuous soybean to one-season crop disturbance in a previous continuous soybean field in Northest China. Field Crop Res. 2012;138:52–6.

  8. 8.

    Chen X, Wang Y, Li W, Wang Y, Wei D, Wang X, Han X. Impact of long-term continuous soybean cropping on ammonia oxidizing bacteria communities in the rhizosphere of soybean in Northeast China. Acta Agric Scand B—Soil Plant Sci. 2015;65(5):470–8.

  9. 9.

    Chen X. Characterizion of microorganism community in the rhizosphere of continuous cropping soybean in black soil. Univ Chin Acad Sci. 2015;5:12.

  10. 10.

    Li C, Li X, Kong W, Wu Y, Wang J. Effect of monoculture soybean on soil microbial community in the Northeast China. Plant Soil. 2010;330(1–2):423–33.

  11. 11.

    Huang LF, Song LX, Xia XJ, Mao WH, Shi K, Zhou YH, Yu JQ. Plant-soil feedbacks and soil sickness: from mechanisms to application in agriculture. J Chem Ecol. 2013;39(2):232–42.

  12. 12.

    Ruan W, Zhu X, Li H, Zhang X, Guo S, Wang J, Zhang F, Gao Y. Soybean autotoxicity: effects of m-hydroxy-phenylacetic acid on cell ultrastructural changes and gene expression in soybean roots. Allelopath J. 2009;24(2):271–82.

  13. 13.

    Cui J, Wang Y, Han J, Cai B. Analyses of the community compositions of root rot pathogenic fungi in the soybean rhizosphere soil. Chile J Agric Res. 2016;76(2):179–87.

  14. 14.

    Leister D, Wang L, Kleine T. Organellar gene expression and acclimation of plants to environmental stress. Front Plant Sci. 2017;8:387.

  15. 15.

    Boyko A, Kovalchuk I. Epigenetic control of plant stress response. Environ Mol Mutagen. 2008;49(1):61–72.

  16. 16.

    Verhoeven KJ, Jansen JJ, van Dijk PJ, Biere A. Stress-induced DNA methylation changes and their heritability in asexual dandelions. New Phytol. 2010;185(4):1108–18.

  17. 17.

    Zhang M, Kimatu JN, Xu K, Liu B. DNA cytosine methylation in plant development. J Genet Genom. 2010;37(1):1–12.

  18. 18.

    Hewezi T. Editorial: epigenetic regulation of plant development and stress responses. Plant Cell Rep. 2017;37(1):1–2.

  19. 19.

    Ganguly DR, Crisp PA, Eichten SR, Pogson BJ. The Arabidopsis DNA methylome is mtable under transgenerational drought stress. Plant Physiol. 2017;175(4):1893–912.

  20. 20.

    Lang Z, Lei M, Wang X, Tang K, Miki D, Zhang H, Mangrauthia SK, Liu W, Nie W, Ma G, et al. The methyl-CpG-binding protein MBD7 facilitates active DNA demethylation to limit DNA hyper-methylation and transcriptional gene silencing. Mol Cell. 2015;57(6):971–83.

  21. 21.

    He XJ, Chen T, Zhu JK. Regulation and function of DNA methylation in plants and animals. Cell Res. 2011;21(3):442–65.

  22. 22.

    Matzke MA, Mosher RA. RNA-directed DNA methylation: an epigenetic pathway of increasing complexity. Nat Rev Genet. 2014;15(6):394–408.

  23. 23.

    Zhu JK. Active DNA demethylation mediated by DNA glycosylases. Annu Rev Genet. 2009;43:143–66.

  24. 24.

    Park J, Frost J, Park K, Ohr H, Park G, Kim S, Eom H, Lee I, Brooks J, Fischer R, et al. Control of DEMETER DNA demethylase gene transcription in male and female gamete companion cells in Arabidopsis thaliana. Proc Natl Acad Sci. 2017;114(8):2078–83.

  25. 25.

    Agius F, Kapoor A, Zhu J. Role of the Arabidopsis DNA glycosylase/lyase ROS1 in active DNA demethylation. Proc Natl Acad Sci. 2006;103(31):11796–801.

  26. 26.

    Morales-Ruiz T, Ortega-Galisteo A, Ponferrada-Marín M, Martínez-Macías M, Ariza R, Roldán-Arjona T. Demeter and repressor of silencing 1 encode 5-methylcytosine DNA glycosylases. Proc Natl Acad Sci. 2006;103(18):6853–8.

  27. 27.

    Gehring M, Huh JH, Hsieh TF, Penterman J, Choi Y, Harada JJ, Goldberg RB, Fischer RL. DEMETER DNA glycosylase establishes MEDEA polycomb gene self-imprinting by allele-specific demethylation. Cell. 2006;124(3):495–506.

  28. 28.

    Lang Z, Wang Y, Tang K, Tang D, Datsenka T, Cheng J, Zhang Y, Handa AK, Zhu JK. Critical roles of DNA demethylation in the activation of ripening-induced genes and inhibition of ripening-repressed genes in tomato fruit. Proc Natl Acad Sci. 2017;114(22):E4511–9.

  29. 29.

    Viggiano L, de Pinto M: Dynamic DNA methylation patterns in stress response. Plant Epigenetics Springer, Cham 2017:281–302.

  30. 30.

    Yaish MW, Al-Lawati A, Al-Harrasi I, Patankar HV. Genome-wide DNA methylation analysis in response to salinity in the model plant caliph medic (Medicago truncatula). BMC Genomics. 2018;19(1):78.

  31. 31.

    Secco D, Whelan J, Rouached H, Lister R. Nutrient stress-induced chromatin changes in plants. Curr Opin Plant Biol. 2017;39:1–7.

  32. 32.

    Dowen RH, Pelizzola M, Schmitz RJ, Lister R, Dowen JM, Nery JR, Dixon JE, Ecker JR. Widespread dynamic DNA methylation in response to biotic stress. Proc Natl Acad Sci. 2012;109(32):E2183–91.

  33. 33.

    Choi CS, Sano H. Abiotic-stress induces demethylation and transcriptional activation of a gene encoding a glycerophosphodiesterase-like protein in tobacco plants. Mol Gen Genomics. 2007;277(5):589–600.

  34. 34.

    Wang WS, Pan YJ, Zhao XQ, Dwivedi D, Zhu LH, Ali J, Fu BY, Li ZK. Drought-induced site-specific DNA methylation and its association with drought tolerance in rice (Oryza sativa L.). J Exp Bot. 2011;62(6):1951–60.

  35. 35.

    Steward N, Ito M, Yamaguchi Y, Koizumi N, Sano H. Periodic DNA methylation in maize nucleosomes and demethylation by environmental stress. J Biol Chem. 2002;277(40):37741–6.

  36. 36.

    Pandey G, Yadav CB, Sahu PP, Muthamilarasan M, Prasad M. Salinity induced differential methylation patterns in contrasting cultivars of foxtail millet (Setaria italica L.). Plant Cell Rep. 2017;36(5):759–72.

  37. 37.

    Harris RA, Wang T, Coarfa C, Nagarajan RP, Hong C, Downey SL, Johnson BE, Fouse SD, Delaney A, Zhao Y, et al. Comparison of sequencing-based methods to profile DNA methylation and identification of monoallelic epigenetic modifications. Nat Biotechnol. 2010;28(10):1097–105.

  38. 38.

    Adey A, Shendure J. Ultra-low-input, tagmentation-based whole-genome bisulfite sequencing. Genome Res. 2012;22(6):1139–43.

  39. 39.

    Henderson IR, Jacobsen SE. Epigenetic inheritance in plants. Nature. 2007;447(7143):418–24.

  40. 40.

    Vanyushin B, Ashapkin V. DNA methylation in higher plants: past, present and future. Biochi Biophys Acta Gene Regul Mech. 2011;1809(8):360–8.

  41. 41.

    McClean PE, Mamidi S, McConnell M, Chikara S, Lee R. Synteny mapping between common bean and soybean reveals extensive blocks of shared loci. BMC Genomics. 2010;11(1):184.

  42. 42.

    Lukens LN, Zhan S. The plant genome's methylation status and response to stress: implications for plant improvement. Curr Opin Plant Biol. 2007;10(3):317–22.

  43. 43.

    Chinnusamy V, Zhu J-K. Epigenetic regulation of stress responses in plants. Curr Opin Plant Biol. 2009;12(2):133–9.

  44. 44.

    Li X, Zhu JD, Hu FY, Ge S, Ye M. Single-base resolution maps of cultivated and wild rice methylomes and regulatory roles of DNA methylation in plant gene expression. BMC Genomics. 2012;13(1):300.

  45. 45.

    Lister R, O'Malley RC, Tonti-Filippini J, Gregory BD, Berry CC, Millar AH, Ecker JR. Highly integrated single-base resolution maps of the epigenome in Arabidopsis. Cell. 2008;133:523–36.

  46. 46.

    Wada Y, Miyamoto K, Kusano T, Sano H. Association between up-regulation of stress-responsive genes and hypomethylation of genomic DNA in tobacco plants. Mol Gen Genomics. 2004;271(6):658–66.

  47. 47.

    Lippman Z, Gendrel A-V, Black M, Vaughn MW, Dedhia N, McCombie WR, Lavine K, Mittal V, May B, Kasschau KD. Role of transposable elements in heterochromatin and epigenetic control. Nature. 2004;430(6998):471–6.

  48. 48.

    Jaenisch R, Bird A. Epigenetic regulation of gene expression: how the genome integrates intrinsic and environmental signals. Nat Genet. 2003;33(Suppl):245–54.

  49. 49.

    Deleris A, Halter T, Navarro L. DNA methylation and demethylation in plant immunity. Annu Rev Phytopathol. 2016;54:579–603.

  50. 50.

    Boyko A, Kathiria P, Zemp FJ, Yao Y, Pogribny I, Kovalchuk I. Transgenerational changes in the genome stability and methylation in pathogen-infected plants: (virus-induced plant genome instability). Nucleic Acids Res. 2007;35(5):1714–25.

  51. 51.

    Ferreira LJ, Azevedo V, Maroco J, Oliveira MM, Santos AP. Salt tolerant and sensitive rice varieties display differential methylome flexibility under salt stress. PLoS One. 2015;10(5):e0124060.

  52. 52.

    Bocchini M, Bartucca ML, Ciancaleoni S, Mimmo T, Cesco S, Pii Y, Albertini E, Del Buono D. Iron deficiency in barley plants: phytosiderophore release, iron translocation, and DNA methylation. Front Plant Sci. 2015;6:514.

  53. 53.

    Labra M, Ghiani A, Citterio S, Sgorbati S, Sala F, Vannini C, Ruffini-Castiglione M, Bracale M. Analysis of cytosine methylation pattern in response to water deficit in pea root tips. Plant Biol. 2002;4(6):694–9.

  54. 54.

    Kovar A, Koukalova B, Bezde M, Opatrn Z. Hypermethylation of tobacco heterochromatic loci in response to osmotic stress. Theor Appl Genet. 1997;95(1–2):301–6.

  55. 55.

    Sabbah S, Raise M, Tal M. Methylation of DNA in NaCl-adapted cells of potato. Plant Cell Rep. 1995;14(7):467–70.

  56. 56.

    Varriale A. DNA methylation in plants and its implications in development, hybrid vigour, and evolution. Plant Epigenetics. Cham: Springer; 2017. p. 263–80.

  57. 57.

    Takuno S, Ran J, Gaut B. Evolutionary patterns of genic DNA methylation vary across land plants. Nat Plants. 2016;2:15222.

  58. 58.

    Patil V, Ward RL, Hesson LB. The evidence for functional non-CpG methylation in mammalian cells. Epigenetics. 2014;9(6):823–8.

  59. 59.

    Zhu J, Kapoor A, Sridhar VV, Agius F, Zhu JK. The DNA glycosylase/lyase ROS1 functions in pruning DNA methylation patterns in Arabidopsis. Curr Biol. 2007;17(1):54–9.

  60. 60.

    Drohat AC, Coey CT. Role of base excision “repair” enzymes in erasing epigenetic marks from DNA. Chem Rev. 2016;116(20):12711–29.

  61. 61.

    Vaquero-Sedas MI, Gamez-Arjona FM, Vega-Palas MA. Arabidopsis thaliana telomeres exhibit euchromatic features. Nucleic Acids Res. 2011;39(6):2007–17.

  62. 62.

    Vega-Vaquero A, Bonora G, Morselli M, Vaquero-Sedas MI, Rubbi L, Pellegrini M, Vega-Palas M. Novel features of telomere biology revealed by the absence of telomeric DNA methylation. Genome Res. 2016;26(8):1047–56.

  63. 63.

    Guo W, Fiziev P, Yan W, Cokus S, Sun X, Zhang MQ, Chen P-Y, Pellegrini M. BS-Seeker2: a versatile aligning pipeline for bisulfite sequencing data. BMC Genomics. 2013;14(1):774.

  64. 64.

    Fernández S, Eritja R, Aviñó A, Jaumot J, Gargallo R. Influence of pH, temperature and the cationic porphyrin TMPyP4 on the stability of the i-motif formed by the 5′-(C 3 TA 2) 4-3′ sequence of the human telomere. Int J Biol Macromol. 2011;49(4):729–36.

  65. 65.

    Schoeftner S, Blasco MA. A ‘higher order’ of telomere regulation: telomere heterochromatin and telomeric RNAs. EMBO J. 2009;28(16):2323–36.

  66. 66.

    Du J, Johnson LM, Jacobsen SE, Patel DJ. DNA methylation pathways and their crosstalk with histone methylation. Nat Rev Mol Cell Biol. 2015;16(9):519–32.

  67. 67.

    Guo JU, Su Y, Shin JH, Shin J, Li H, Xie B, Zhong C, Hu S, Le T, Fan G, et al. Distribution, recognition and regulation of non-CpG methylation in the adult mammalian brain. Nat Neurosci. 2014;17(2):215–22.

  68. 68.

    Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, Pradhan S, Nelson SF, Pellegrini M, Jacobsen SE. Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature. 2008;452(7184):215–9.

  69. 69.

    Gruenbaum Y, Naveh-Many T, Cedar H, Razin A. Sequence specificity of methylation in higher plant DNA. Nature. 1981;292(5826):860–2.

  70. 70.

    Pelizzola M, Ecker JR. The DNA methylome. FEBS Lett. 2011;585(13):1994–2000.

  71. 71.

    Ramsahoye BH, Biniszkiewicz D, Lyko F, Clark V, Bird AP, Jaenisch R. Non-CpG methylation is prevalent in embryonic stem cells and may be mediated by DNA methyltransferase 3a. Proc Natl Acad Sci. 2000;97(10):5237–42.

  72. 72.

    Plaxton WC. The organization and regulation of plant glycolysis. Annu Rev Plant Biol. 1996;47:185–214.

  73. 73.

    Seth CS, Remans T, Keunen E, Jozefczak M, Gielen H, Opdenakker K, Weyens N, Vangronsveld J, Cuypers A. Phytoextraction of toxic metals: a central role for glutathione. Plant Cell Environ. 2012;35(2):334–46.

  74. 74.

    Noctor G, Mhamdi A, Chaouch S, Han Y, Neukermans J, Marquez-Garcia B, Queval G, Foyer CH. Glutathione in plants: an integrated overview. Plant Cell Environ. 2012;35(2):454–84.

  75. 75.

    Foyer CH, Shigeoka S. Understanding oxidative stress and antioxidant functions to enhance photosynthesis. Plant Physiol. 2011;155(1):93–100.

  76. 76.

    Nazar R, Iqbal N, Syeed S, Khan N. Salicylic acid alleviates decreases in photosynthesis under salt stress by enhancing nitrogen and sulfur assimilation and antioxidant metabolism differentially in two mungbean cultivars. J Plant Physiol. 2011;168(8):807–15.

  77. 77.

    Pál M, Kovács V, Szalai G, Soós V, Ma X, Liu H, Mei H, Janda T. Salicylic acid and abiotic stress responses in rice. J Agron Crop Sci. 2014;200(1):1–11.

  78. 78.

    Capaldi FR, Gratão PL, Reis AR, Lima LW, Azevedo RA. Sulfur metabolism and stress defense responses in plants. Trop Plant Biol. 2015;8(3–4):60–73.

  79. 79.

    Lu SC. S-adenosylmethionine. Int J Biochem Cell Biol. 2000;32(4):391–5.

  80. 80.

    Fuso A, Seminara L, Cavallaro RA, D'Anselmi F, Scarpa S. S-adenosylmethionine/homocysteine cycle alterations modify DNA methylation status with consequent deregulation of PS1 and BACE and beta-amyloid production. Mol Cell Neurosci. 2005;28(1):195–204.

  81. 81.

    Marco F, Alcazar R, Tiburcio AF, Carrasco P. Interactions between polyamines and abiotic stress pathway responses unraveled by transcriptome analysis of polyamine overproducers. OMICS. 2011;15(11):775–81.

  82. 82.

    Tiburcio AF, Altabella T, Bitrián M, Alcázar R. The roles of polyamines during the lifespan of plants: from development to stress. Planta. 2014;240(1):1–18.

  83. 83.

    Pathak MR, Teixeira da Silva JA, Wani SH. Polyamines in response to abiotic stress tolerance through transgenic approaches. GM Crops Food. 2014;5(2):87–96.

  84. 84.

    Minocha R, Majumdar R, Minocha SC. Polyamines and abiotic stress in plants: a complex relationship. Front Plant Sci. 2014;5:175.

  85. 85.

    Gupta K, Dey A, Gupta B. Plant polyamines in abiotic stress responses. Acta Physiol Plant. 2013;35(7):2015–36.

  86. 86.

    Majumdar R, Barchi B, Turlapati SA, Gagne M, Minocha R, Long S, Minocha SC. Glutamate, ornithine, arginine, proline, and polyamine metabolic interactions: the pathway is regulated at the post-transcriptional level. Front Plant Sci. 2016;7:78.

  87. 87.

    Akcay N, Bor M, Karabudak T, Ozdemir F, Turkan I. Contribution of gamma amino butyric acid (GABA) to salt stress responses of Nicotiana sylvestris CMSII mutant and wild type plants. J Plant Physiol. 2012;169(5):452–8.

  88. 88.

    Hussain SS, Ali M, Ahmad M, Siddique K. Polyamines: natural and engineered abiotic and biotic stress tolerance in plants. Biotechnol Adv. 2011;29(3):300–11.

  89. 89.

    Bitrian M, Zarza X, Altabella T, Tiburcio AF, Alcazar R. Polyamines under abiotic stress: metabolic crossroads and hormonal crosstalks in plants. Metabolites. 2012;2(3):516–28.

  90. 90.

    Winter G, Todd CD, Trovato M, Forlani G, Funck D. Physiological implications of arginine metabolism in plants. Front Plant Sci. 2015;6:534.

  91. 91.

    Domingos P, Prado AM, Wong A, Gehring C, Feijo JA. Nitric oxide: a multitasked signaling gas in plants. Mol Plant. 2015;8(4):506–20.

  92. 92.

    Beisson F, Li Y, Bonaventure G, Pollard M, Ohlrogge JB. The acyltransferase GPAT5 is required for the synthesis of suberin in seed coat and root of Arabidopsis. Plant Cell. 2007;19(1):351–68.

  93. 93.

    Upchurch RG. Fatty acid unsaturation, mobilization, and regulation in the response of plants to stress. Biotechnol Lett. 2008;30(6):967–77.

  94. 94.

    Kachroo A, Kachroo P. Fatty acid-derived signals in plant defense. Annu Rev Phytopathol. 2009;47:153–76.

  95. 95.

    Creelman RA, Mulpuri R. The oxylipin pathway in Arabidopsis. Arabidopsis Book. 2002;1:e0012.

  96. 96.

    Tuteja N, Ahmad P, Panda BB, Tuteja R. Genotoxic stress in plants: shedding light on DNA damage, repair and DNA repair helicases. Mutat Res. 2009;681(2–3):134–49.

  97. 97.

    Balestrazzi A, Confalonieri M, Macovei A, Dona M, Carbonera D. Genotoxic stress, DNA repair, and crop productivity. In: Tuteja N, Gill SS, editors. Crop improvement under adverse conditions. Berlin: Springer; 2012. p. 153–70.

  98. 98.

    Filichkin S, Priest HD, Megraw M, Mockler TC. Alternative splicing in plants: directing traffic at the crossroads of adaptation and environmental stress. Curr Opin Plant Biol. 2015;24:125–35.

  99. 99.

    Liu JX, Howell SH. Endoplasmic reticulum protein quality control and its relationship to environmental stress responses in plants. Plant Cell. 2010;22(9):2930–42.

  100. 100.

    Nakaminami K, Matsui A, Shinozaki K, Seki M. RNA regulation in plant abiotic stress responses. Biochim Biophys Acta Gene Regul Mech. 2012;1819(2):149–53.

  101. 101.

    Howell SH. Endoplasmic reticulum stress responses in plants. Annu Rev Plant Biol. 2013;64:477–99.

  102. 102.

    Xiang H, Zhu J, Chen Q, Dai F, Li X, Li M, Zhang H, Zhang G, Li D, Dong Y. Single base-resolution methylome of the silkworm reveals a sparse epigenomic map. Nat Biotechnol. 2010;28(5):516–20.

  103. 103.

    Lister R, Pelizzola M, Dowen RH, Hawkins RD, Hon G, Tonti-Filippini J, Nery JR, Lee L, Ye Z, Ngo QM. Human DNA methylomes at base resolution show widespread epigenomic differences. Nature. 2009;462(7271):315–22.

Download references

Acknowledgements

Not applicable.

Fundings

This work was supported by the National Natural Science Foundation of China [grant number 31201163]; the Science and Technology Development Program of Heilongjiang Agricultural Reclamation Department of China [grant number HNK135–02–06-06]; and PhD research startup foundation of Heilongjiang Bayi Agricultural University [XDB-2016-23]. The funding agents only provided the financial support but did not participate in the design and implementation of the experimental, analysis of data and in the writing of the manuscript.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Author information

XL, SF and ZD conceived and designed this study. XL, SF and YH wrote the paper. YZ and NF analyzed the data. JL supplied materials for experiments. XH, JD and WZ prepared samples for bisulfite sequencing and performed the real-time PCR test. All authors have read and approved the manuscript for publication.

Correspondence to Dianfeng Zheng or Shumei Fang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional file

Additional file 1:

Table S1. Morphological indexes of soybean plants under different conditions. Table S2. Effective coverage of each chromosome in each sample. Table S3. Effective coverage of various regions in each sample. Table S4. The primer sets used in qRT-PCR. Figure S1. Cumulative distribution of effective sequencing depth in total cytosine and three sequence contexts. Figure S2-1. Methylcytosine density throughout chromosome one to five in sensitive HF55 and tolerant KX2 under different conditions. Figure S2-2. Methylcytosine density throughout chromosome six to ten in sensitive HF55 and tolerant KX2 under different conditions. Figure S2-3. Methylcytosine density throughout chromosome eleven to fifteen in sensitive HF55 and tolerant KX2 under different conditions. Figure S2-4. Methylcytosine density throughout chromosome sixteen to twenty in sensitive HF55 and tolerant KX2 under different conditions. (DOCX 10182 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • Continuous cropping comprehensive stress
  • DNA demethylation
  • Demethylase
  • Differentially methylated genes
  • Metabolism analysis
  • Soybean