Skip to main content

Integrated methylome and transcriptome analysis unravel the cold tolerance mechanism in winter rapeseed(Brassica napus L.)

Abstract

Background

Cytosine methylation, the main type of DNA methylation, regulates gene expression in plant response to environmental stress. The winter rapeseed has high economic and ecological value in China's Northwest, but the DNA methylation pattern of winter rapeseed during freezing stress remains unclear.

Result

This study integrated the methylome and transcriptome to explore the genome-scale DNA methylation pattern and its regulated pathway of winter rapeseed, using freezing-sensitive (NF) and freezing-resistant (NS) cultivars.The average methylation level decreased under freezing stress, and the decline in NF was stronger than NS after freezing stress. The CG methylation level was the highest among the three contexts of CG, CHG, and CHH. At the same time, the CHH proportion was high, and the methylation levels were highest 2 kb up/downstream, followed by the intron region. The C sub-genomes methylation level was higher than the A sub-genomes. The methylation levels of chloroplast and mitochondrial DNA were much lower than the B. napus nuclear DNA, the SINE methylation level was highest among four types of transposable elements (TEs), and the preferred sequence of DNA methylation did not change after freezing stress. A total of 1732 differentially expressed genes associated with differentially methylated genes (DMEGs) were identified in two cultivars under 12 h and 24 h in three contexts by combining whole-genome bisulfite sequencing( and RNA-Seq data. Function enrichment analysis showed that most DMEGs participated in linoleic acid metabolism, alpha-linolenic acid metabolism, carbon fixation in photosynthetic organisms, flavonoid biosynthesis, and plant hormone signal transduction pathways. Meanwhile, some DMEGs encode core transcription factors in plant response to stress.

Conclusion

Based on the findings of DNA methylation, the freezing tolerance of winter rapeseed is achieved by enhanced signal transduction, lower lipid peroxidation, stronger cell stability, increased osmolytes, and greater reactive oxygen species (ROS) scavenging. These results provide novel insights into better knowledge of the methylation regulation of tolerance mechanism in winter rapeseed under freezing stress.

Peer Review reports

Introduction

Cold stress mainly includes chilling stress (above 0 °C) and freezing stress (below 0 °C). Freezing stress is one of the most unfavorable abiotic stresses in Northwestern China from late autumn to early spring, seriously affecting the growth, development, quality, and productivity of crops, and even restricting the geographical distribution of plant species, making this region relatively single species. Freezing stress damage includes multiple aspects, for example, cold, dehydration, osmotic, and mechanical stresses, which are more severe than chilling stress damage [1, 2].To cope with freezing injury, plants have evolved various adaptive mechanisms to respond and survive to low freezing conditions, including changes in physiological, biochemical, cellular, and molecular changes [3, 4]. Freezing tolerance is a necessary trait for crops to withstand low temperatures, especially in overwintering crops [5]. Some research showed that the freezing tolerance of plants is associated with cold acclimation, while DNA methylation plays a vital role [6]. At the same time, these studies showed that epigenetic mechanisms form stress memory, which offspring can inherit [7].

DNA methylation, one of the hallmarks of epigenetic mechanisms, plays an important role in plant adaptation to abiotic stresses by regulating transcriptional activities of the stress-response genes [8, 9]. A previous study showed that cold- induced demethylation in maize is involved in many processes, such as hormone regulation, cold response, photosynthesis, and transposon activation [10]. In rice, DNA methylation level was down-regulated in response to chilling stress, and differentially methylated genes were involved in regulating transcription factors families' [11]. In addition, the DNA methylation patterns changed in Celtis bungeana after chilling and freezing stresses [5]. In Brassica rapa, cold acclimation altered DNA methylation patterns and confers tolerance to heat and increases growth rate [12]. In tea plants significant loss of DNA methylation is found at their promoter and gene body regions after chilling stress; DNA methylation in the promoter positively regulates the gene expression, and negatively regulates the gene expression of gene body regions [13].

The rapeseed (Brassica napus L.) is an important oil-yielding seed crops all over the worldwide, with the largest planting area and relatively high yield and oil content. The major area of winter B.napus is the Yangtze River basin, which belongs to semi-winter rapeseed and more than 80% of china’s total area of rapeseed [14]. In a decade year, we have successfully cultivated the freezing-resistant cultivar NS, which can survive normally in winter at the -26℃ in most northwestern areas and contain abundant freezing-resistant genes, originating from the distant hybridization and cold acclimation between B.napus and B. Rapa [15, 16]. It is of great significance to continuously increase the planting areas of winter rapeseed in northern China for exploiting its economic value and ecological value, whereas, to expand the planting areas of winter rapeseed, we must understand molecular genetic mechanism of cold resistance and breed more winter rapeseed cultivars with strong cold tolerance. In recent years, the DNA methylation was the research hotspot on plant response to abiotic stresses, especially in cold stress. However, the epigenetic mechanism of how DNA methylation regulates genes to elevate cold tolerance of winter rapeseed are still little known.

The study first investigated the changes in methylation patterns of a pair of resistant and susceptible winter rapeseed cultivars after freezing stress. Meanwhile, the methylome and transcriptome were used to identify differentially methylated regions (DMRs), differentially methylated genes (DMGs), and differentially expressed genes associated with differentially methylated genes (DMEGs) between the two cultivars. Subsequently, the main metabolic pathways enriched in DMEGs involved in freezing stress were explored. These results will enable us to understand better the methylation-participated molecular regulatory network in winter rapeseed under freezing stress.

Results

Genome-wide DNA methylation profiles of NS and NF

Six genomic DNA libraries were constructed with the leaves of cultivars NS and NF to study the genome-wide DNA methylation pattern of winter rapeseed under freezing treatment.More than 192 million raw reads per DNA library samples were generated by whole-genome bisulfite sequencing (WGBS). After removal of related adapters, low-quality reads, poly-A and Ns-containing reads, 219 million, 190 million, 199 million, 221 million, 243 million, and 288 million were collected in six libraries (Table 1), of which 85.82% (NSt0), 86.38% (NSt1), 85.46% (NSt2), 86.96% (NFt0), 88.92% (NFt1), 86.87% (NFt2) were uniquely mapped to the Brassica napus L reference genome (GCA_000686985.2), respectively. The average sequence depths of all libraries were more than 25.35 fold, and the effective coverage rates of all chromosomes were between 50–60% (Supplementary Table S1). However, practically cytosines conversion rate for all libraries were close to 99%, indicating a high bisulfite sequencing conversion rate in this study.

Table 1 Whole genome DNA bisulfite sequencing data in two rapeseed cultivars

The average methylation level showed a downward trend in winter rapeseed after freezing stress. The trend was similar in CG, CHG, CHH sequence sites(where H = A, T, or C) and more pronounced in NF than NS (Fig. 1). However, comparing the two cultivars showed that methylation levels were hardly different before freezing treatment, but the methylation level of NF was lower than that of NS after freezing treatment, especially in the t2 treatment. It is worth noting that among three sequence contexts, CG has the highest methylation level, followed by CHG and CHH. Furthermore, almost half of 5-methylcytosines (5-mCs) occurred at CHH sequence site in all samples, with a lower proportion occurring at CG and CHG sites (Fig. 2).

Fig. 1
figure 1

The average DNA methylation leves of two rapeseed cultivars in genome-wide methylated cytosine(C), and in the CG,CHG and CHH contexts

Fig. 2
figure 2

The percentages of methylated cytosines under CG, CHG, and CHH contexts. A the percentages of methylated cytosines in NSt0; B the percentages of methylated cytosines in NSt1; C the percentages of methylated cytosines in NSt2; D the percentages of methylated cytosines in NFt0; E the percentages of methylated cytosines in NFt1; F the percentages of methylated cytosines in NFt2

The methylation levels of all 19 chromosomes in B.napus leaves showed that the C sub-genomes methylation level was higher than the A sub-genomes (Fig. 3). After freezing stress, all 19 rapeseed chromosomes showed a downward methylation level trend (Fig. 3; Supplementary Table S1). Interestingly, the average methylation levels of mitochondrion (mt) and chloroplast (cp) DNA were lower than in B.napus nuclear chromosomes; however, the methylation levels of mitDNA and cpDNA increased with treatment. At the same time, the chloroplastt's methylation level was higher than the mitochondrion's under the same conditions.

Fig. 3
figure 3

Distribution of methylation levels in different chromosemes of two rapeseed cultivars

The methylation patterns among gene features in response to freezing stress were compared in different genome regions of six libraries to further identify the methylation patterns. Eight functional regions of the genome, gene body, exons, introns, protein-coding sequences (CDS), 5' untranslated regions (5'UTR), 3' untranslated region (3’UTR), 2 kb upstream, and 2 kb downstream, showed a decline in methylation after freezing stress (Fig. 4). Similarly, the CHH sites showed the lowest methylation levels, while CG had the highest of these functional regions (Supplementary Figure S1). A comparison of different functional regions demonstrated that the methylation levels were highest in the 2 kb upstream and 2 kb downstream, followed by an introns, with the lower the methylation levels in gene body, CDS, Exons, 5'UTR, 3'UTR region. Based on these results, it is speculated that methylation may play a role in the regulatory region of genes to activate gene expression.

Fig. 4
figure 4

Distribution of methylation levels in different genomic functional regions of two rapeseed cultivars

Because different transposable elements (TEs) have different functions, the methylation levels between different TEs can affect their function. Therefore, three types of TEs, long terminal repeats (LTRs), long interspersed nuclear elements (LINEs), short interspersed nuclear elements (SINEs), and DNA methylation levels were analyzed. The result showed that the SINEs had the highest methylation level among the rapeseed TEs in the rapeseed genome, and SINEs had the highest methylation levels in the CHG and CHHsequences but had the lowest methylation level in the CG context (Fig. 5; Supplementary Figure S2). All TEs types in all samples had a similar downward trend after freezing stress, which was closely related to the activation of more stress genes' under freezing stress.

Fig. 5
figure 5

Distribution of methylation levels in different transposable elements of two rapeseed cultivars

Because DNA methylation sites were highly correlated with DNA sequence, the 9 bp sequence information near the methylation sites in three contexts was analyzed to obtain the sequence preference around all mCs in the rapeseed genome. The result suggested that the methylation always occurred at the TCGA sequence in the CG sites (Fig. 6). Moreover, CTG and CAA sequences were the most common sequence motifs at the CHG and CHH sites. Unexpectedly, the sequence preference of cytosine DNA methylation did not differ between two cultivars and between freezing stress and control, and TCGA, CTG, and CAA sequencees dominated in all samples under CG, CHG, and CHH contexts, indicating that the new cytosine preferentially occurs in the above dominate sequences, and the sequence preference of cytosine DNA methylation sits in the rapeseed genome did not alter under freezing stress.

Fig. 6
figure 6

Sequence preference round all methylated cytosine of rapeseed genome under CG, CHG and CHH contexts. A The sequence preference of NF in t0, t1, and t2 treatments under CG, CHG and CHH contexts. B The sequence preference of NS in t0, t1, and t2 treatments under CG, CHG and CHH contexts

Analysis of differentially methylated regions (DMRs) and differentially methylated genes (DMGs) between two cultivars

To understand the difference in DNA methylation in two cultivars with varying tolerance to freezing stress, we compared the DMRs of of NF vs. NS at different treatments at CG, CHG, and CHH sites. In the CG context, 37,898, 47,952, and 60,246 DMRs showed hypermethylation, and 21,258, 25,915, and 23,734 DMRs showed hypomethylation in NFt0 vs. NSt0, NFt1 vs. NSt1, and NFt2 vs. NSt2 respectively (Supplementary Table.S1, S2). For the CHG context, 58,354, 88,277, and 86,366 DMRs were identified in the groups of NFt0 vs. NSt0, NFt1 vs. NSt1, and NFt2 vs. NSt2, including 32,014, 28,595, and 35,333 hypermethylated DMRs and 26,340, 59,682, and 51,033 hypomethylated DMRs, respectively. Similarly, 71,983 (54,417 hyper and 17,566 hypo), 59,691 (29,182 hyper and 30,509 hypos), and 58,658 (42,665 hyper and 15,993 hypos) DMRs were identified in NFt0 vs. NSt0, NFt1 vs. NSt1, and NFt2 vs. NSt2 under CHH context, respectively. Genes overlapping with DMRs for at least 1 bp in the functional region were defined as DMGs. In total of 18,698, 13,847, and 20,981 DMGs were identified in NFt0 vs. NSt0 under the CG, CHG, and CHH context (Fig. 7; Supplementary Table S3). In NFt1 vs. NSt1, 22,261, 16,915, and 19,864 DMGs were obtained in the CG, CHG, and CHH context. Similarly, under three contexts, 24,426,17,088, and 18,452 DMGs were identified in NFt2 vs. NSt2. The number of hypermethylated genes in CG and CHG increased after freezing stress treatment (t1 and t2) compared to the control (t0). However, the number of hypermethylated genes decreased in the CHH context, indicating that hypermethylation mainly occurred in CG and CHG sites during freezing stress.

Fig. 7
figure 7

Number of diferentially methylated genes identified between two rapeseed cultivars under CG, CHG and CHH contexts

RNA sequencing analysis and defining differentially expressed genes (DEGs) between two cultivars

RNA sequencing (ribonucleic acid -seq) was also performed on eighteen libraries from six samples (three biological replicates for per sample) of two cultivars to investigate the association between DNA methylation and gene expression. The high-quality clean reads per library were more than 39 million, the Q30 for all libraries was greater than 94%, and the GC content per library was approximately 48%(Supplementary Table S4). For all libraries, over 77% of reads, were uniquely mapped to the reference genome and the mapping ratio was close to 80%. The correlation analysis showed that the correlation coefficient between the three replicates of any particular sample was extremely high, indicating that the experiment was reliable and reasonable(Supplementary Table S4). A total of 16,441 DEGs were identified in the NFt0 vs NSt0 sample, which including 9204 up-regulated genes and 7237 down-regulated genes(Fig. 8; Supplementary Table S5). Compared to controls, 17,537 (up 11,345 and down 6192) and 19,638 (up 10,222 and down 9416) DEGs were identified in the NFt1 vs. NSt1 and NFt2 vs. NSt2 groups(Fig. 8; Supplementary Table S5). The number of DEGs increased after freezing treatment, indicating that some genes were aroused during freezing stress.

Fig. 8
figure 8

Number of DMEGs identified between two rapeseed cultivars. A The number of DMEGs in NFt0 vs. NSt0; B the number of DMEGs in NFt1 vs. NSt1; C the number of DMEGs in NFt2 vs. NSt2

Comparative analysis of DEGs associated with DMGs between two cultivars

To further explore the possible association between DNA methylation and gene expression, firstly, the overlap between the DMGs and DEGs of two cultivars was examined at CG, CHG, and CHH sites. In the NFt0 vs. NSt0 of 3195, 2347, and 4027 DMEGs were identified in under CG, CHG, and CHH contexts (Fig. 8; Supplementary S6). After merging DMEGs from CG, CHG, and CHH context, 6217 DMEGs were obtained in the NFt0 vs. NSt0 (Fig. 8). After freezing stress for 12 h, 4064, 2933, and 4147 DMEGs were determined in the NFt1 vs. NSt1 sample at the CG, CHG, and CHH sites, and 7157 DMEGs were obtained after combination at CG, CHG, and CHH contexts (Fig. 8; Supplementary S6). After 24 h of freezing stress, 5026, 3324, and 4348 DMEGs were identified in the NFt2 vs. NSt2 under the CG, CHG, and CHH contexts, and 8197 DMEGs were obtained after combination at CG, CHG, and CHH contexts(Fig. 8; Supplementary S6). The DMEGs of NFt0 vs. NSt0 were removed from the overlapping DMEGs of NFt1 vs. NSt1 and NFt2 vs. NSt2 to obtain 1732 DMEGs, which will be used for subsequent functional analysis (Fig. 9).

Fig. 9
figure 9

Number of DMEGs identified between two rapeseed cultivars at three treatments

GO annotation and KEGG pathway enrichment analysis of DMEGs

To understand the potential role of DNA methylation patterns in freezing tolerance, 1732 DMEGs were subjected to gene ontology (GO) classification annotation (Fig. 10A; Supplementary Table S7). Those DMEGs were annotated to 1647 functional categories, including 1143 biological processes, 161 cellular components, and 343 molecular functions. They were found to be mainly involved in biological processes, such as response to external stimulus, response to stimulus, response to osmotic stress, response to abiotic stimulus, and response to stress. In addition, cytoskeletal part, supramolecular fiber, and microtubule were the terms that significantly dominated the cellular components. The endoplasmic reticulum (ER) retention sequence binding, phosphotransferase activity, DNA-3-methyl bases glycosylase activity, and C-methyltransferase activity were the most represented GO categories in the molecular functions.

Fig. 10
figure 10

The functional annotation of DMEGs. A GO annotation of DMEGs; B KEGG pathway enrichment of DMEGs

The 1732 DMEGs identified in winter rapeseed were submitted to KEGG pathway enrichment analysis to determine their complex biological functions (Fig. 10B; Supplementary Table S7). A total of 343 DMEGs were mapped successfully to 113 KEGG pathways. The top 20 KEGG enrichment pathways are shown that linoleic acid metabolism, alpha-linolenic acid metabolism, carbon fixation in photosynthetic organisms, and flavonoid biosynthesis were significantly enriched. Besides, plant hormone signal transduction and some of the amino acid metabolism pathways were also enriched.

Identification of differentially expressed transcription factors for DMEGs

One-hundred and eighty-three 183 transcription factors (TFs) were identified from 1732 overlapped DMEGs in NFt1 vs. NSt1 and NFt2 vs. NSt2 groups.Most TFs were up-regulated, including ethylene response factor (AP2/ERF), transcription factor bHLH (bHLH), WRKY transcription factor (WRKY), zinc finger CCCH domain-containing protein (C3H), transcription factor MYB (MYB), basic leucine zipper (bZIP), homeobox-leucine zipper protein (HB), protein TIFY (TIFY), auxin-responsive protein (AUX/IAA), and GATA transcription factor (C2C2-GATA) families. Most of the down-regulated TFs belonged to the auxin response factor (B3-ARF), heat stress transcription factor (HSF), zinc finger protein (C2H2) (Fig. 11; Supplementary Table S8). Correlation analysis showed that the overall methylation level was significantly negatively correlated with the gene expression level (p < 0.05), but the correlation between different sites and gene regions was inconsistent (Supplementary Table S8). These results suggested that methylation played an important role in the process of activating TFs under freezing stress, especially the up-regulated TFs, thereby regulating the expression of downstream genes in response to freezing stress.

Fig. 11
figure 11

Classification of differentially expressed TFs between two winter rapeseed under freezing treatment. A The unmber of differentially expressed TFs in NSt1 vs. NFt1; B the unmber of differentially expressed TFs in NSt2 vs. NFt2

RNA-Seq validation of DMEGs by quantitative Real-Time PCR

Based on KEGG enrichment analysis, eleven genes related to plant hormone signal transduction and linoleic acid metabolism pathways in rapeseed under freezing stress were selected to validate the reliability of RNA-seq data by quantitative Real-Time PCR (qRT-PCR). Most genes were found to be changed consistently at both the qRT-PCR and RNA-seq data under freezing stress in rapeseed cultivar NS compared to cultivar NF, despite differences in expression levels of some genes. (Fig. 12; Supplementary Table S9). These results suggested that the RNA-seq results were true and reliable, which can be used in the present study.

Fig. 12
figure 12

Comparative analysis of mRNA and RNA-seq levels between two rapeseed cultivars under freezing stress. A The column chart between mRNA and RNA-seq level; B Concordance changes between mRNA and RNA-seq levels. The data is the log2 value of the ratio of NS gene expression value to NF expression value under freezing treatment (t1 and t2).Values are means ± SD (from three biological replicates) of Q-PCR data under freezing treatment

Discussion

Rapeseed is China's primary oilseed crop, including, including the winter rapeseed (over 90%) and little of spring rapeseed [17]. Previous research indicated that the 35° north latitude (Tianshui in the Gansu Province) was the northernmost limit of B. napus winter rapeseed of in China [18]. The cultivated spring rapeseed and B. rapa can survive north of 35°N [19, 20]. Over the past decade, our research group has successfully selected and bred rapeseed cultivars of B. napus that can overwinter safely in most of the northwestern areas from the 38° to the 42° north latitude, where the lowest temperature in winter is -26℃ [15, 16]. Moreover, the cold tolerance characteristics of these rapeseed cultivars at the transcription, protein, physiological, biochemical, and subcellular levels were studied [15, 21, 22]. The strong cold-resistant cultivar had stronger antioxidant enzyme activity and accumulation of cryoprotective molecules, associated with the metabolic pathways, including alpha-linolenic acid metabolism, plant hormone signal transduction, microbial metabolism in diverse environments, linoleic acid metabolism, phenylpropanoid biosynthesis, glutathione metabolism, flavonoid biosynthesis metabolism, and amino acid metabolism. DNA methylation, the main epigenetic modification, was involved in regulating plant adaptation to abiotic stresses and was inherited by offspring [8, 9, 23]. However, the current understanding of the methylation pattern amid freezing stress and the complex mechanisms of methylation regulating freezing tolerance in winter rapeseed is still limited. Therefore, in the current study, the whole-genome bisulfite sequencing (WGBS) was used to analyze the methylation pattern of NS and NF cultivars under freezing stress. The integrative analysis of methylome and transcriptome will pave the pathway to understanding the cold tolerance mechanism of winter rapeseed (Fig. 13).

Fig. 13
figure 13

A possible model of the freezing tolerance mechanism regulated by methylation in winter rapeseed. The solid line arrow represent metabolic pathways and cellular processes. Dashed arrows indicate possible regulation. Heat maps indicate the log2 value of the ratio of NS gene expression value to NF expression value under freezing treatment (t1 and t2)

Methylation level changes of winter rapeseed under freezing stress

Enzymes regulate methylation and demethylation in plants, which are the part of stress reactions to respond to environmental conditions [24, 25]. It is proposed that DNA methylation is related to the stress memory of plants so that their offspring can respond more quickly and accurately to the same stress [7]. This is consistent with our theory of enhancing the cold resistance of rapeseed through cold acclimation, and it also exists in other crops. Current methods have been to detect the methylation dynamics, such as methylation-sensitive amplification polymorphism (MSAP), methylation-specific PCR (MSP), and WGBS. The WGBS method is widely used in plant methylation research due to its high accuracy, wide range, high reliability, and cost-effectiveness [26–, 27,28,29]. In winter rapeseed, previous studies using the MSAP method had shown that the methylation level decreased after low-temperature stress. In this study, the methylation levels of the two cultivars showed a downward trend in all contexts following the freezing stress, and the effect of freezing-sensitive cultivars was stronger than freezing-resistance cultivars. Notably, the methylation level of CG was highest among the three contexts, while the proportion of CHH was abundant, which was reported by other [30]. Previous studies revealed that the DNA methylation level in the CG context was positively correlated with plants' genome size, 30.5% in Arabidopsis and 92.5% in Beta vulgaris [31, 32]. In this study, we found that the methylation level in leaves of Brassica napus was 64.24–71.70% under the CG context (Fig. 1). The methylation level is also related to the tissue parts of the plant [30]. Examination of rapeseed's 19 chromosomes indicated that the C sub-genomes methylation level was higher than in A sub-genomes (Fig. 3). A similar result was obtained in the Genic Male Sterile Line, the Restorer Line, and cultured microspores of rapeseed [30, 33]. Furthermore, we found that the methylation levels of cpDNA and mtDNA were much lower than those of nuclear, and the mtDNA methylation levels were lower than that of cpDNA. These results are consistent with previous report on nuclear DNA, cpDNA and mtDNA methylation [34]. In addition, their research indicated the methylation level first increased and then decreased under freezing stress, which further suggest that the methylation role of chloroplast and mitochondrial DNA control gene expression in plant response to abiotic stress. The comparison between different functional regions depicted that the methylation level was highest in the 2 kb upstream and 2 kb downstream regions, followed by intron, while that in gene body, CDS, exons, 5'UTR, and 3'UTR regions were lower, indicating that the methylation plays a role in the genes' regulatory region (Fig. 4). Simultaneously, the SINEs had the highest average methylation levels, and LTRs had the lowest methylation level among TEs (Fig. 5). Furthermore, the prefered sequence of DNA methylation for winter rapeseed did not change after freezing stress (Fig. 6). In the integrated analysis between DMGs and DEGs, 1732 DMEGs with significant change were yielded under the t1 and t2 treatments for the freezing-resistant (NS) cultivar relative to freezing-sensitive (NF) cultivar (Fig. 8– 9). Function enrichment analysis showed that DMEGs were mainly involved in linoleic acid metabolism, alpha-linolenic acid metabolism, carbon fixation in photosynthetic organisms, nitrogen metabolism, flavonoid biosynthesis, plant hormone signal transduction, and some amino acid metabolism pathways (Fig. 10; Supplementary Table S7).

Methylation-regulated transcription factors and protein kinases

Many transcription factors have been reported participate to convey cold stress signals. Among them, the C-repeat binding transcription factor (CBF) is a central player in freezing tolerance. The inducer of C-repeat-binding factor expression 1 (ICE1) activate CBF and with it the expression of downstream cold-regulated genes, thereby enhancing freezing tolerance [35, 36]. The ICE1 belongs to the transcription factor bHLH family, and it positively regulates the CBF expression to contribute to freezing stress. In addition, many transcription factors have been determined to positively participate in the cold resistance of plants, including AP2/ERF, MYB88 /MYB124, WRKY, bZIP, C3H, and TIFY [15–, 37,38,39,40,41]. In this study, most AP2/ERF, WRKY, MYB, bZIP, C3H, TIFY, AUX, and HB transcription factors were up-regulated in cultivar NS compared to cultivar NF, but some member of TFs were down-regulated(Fig. 12; Supplementary Table S8). Hypermethylation of gene body in CG, CHG, and CHH contexts caused up- or downregulated expression of transcription factors under Phosphate Starvation Stress in Rice [42]. Also, the overall methylation level was significantly negatively correlated with the gene expression level (p < 0.05), but the correlation of different loci and gene positions was different, and the methylation level was also positively correlated with the expression level.(Supplementary Table S8). We find that some members of up-regulated TFs, such as MYBs, WRKY, AP2/ERF, TIFY, and bZIPs, were hypomethylated in cultivar NS compared to cultivar NF in different gene regions, while the some members of up-regulated TFs were hypermethylated. Meanwhile, most of TFs have different methylation level in different gene regions and different contexts, so the methylation effect of TFs also need further research in future works. It has been reported that several protein kinases regulated the cold stress response of plants by a cold-sensing calcium channel [43]. The mitogen-activated protein kinases (MAPKs), Ca2+-dependent protein kinases (CDPKs), calcium/calmodulin-regulated receptor-like kinases (CRLKs), calcineurin-B-like interacting protein kinases (CIPKs), dehydration-responsive element-binding protein kinases (DREB), and receptor-like kinases (RLKs) have also been identified as key protein in the Ca2+ signaling pathway [43,44,45]. Similar results were obtained in this study; more than two-thirds of protein kinases were up-regulated in NS compared to NF after freezing stress (Fig. 13; Supplementary Table S8). The result of methylation showed the methylation regulation patterns of homologous genes or genes' different regions were inconsistent. These results showed that the protein kinases and TFs might be activated by the Ca2+ signaling in rapeseed after cold treatment, and triggering the expressions of downstream freezing responsive genes; the methylation might regulate the Ca2+ signaling by adjusting protein kinases and TFs.

Methylation-regulated lipid metabolism

Under cold stress, the plant cell membrane alters the fatty acid profile and loses fluidity; however, increasing the unsaturation in the membrane lipids can maintain the normal fluidity of membranes and enhance the tolerance to freezing stress [46, 47]. The secretory phospholipase A2 (PLA2G), a key enzyme that converts phosphatidylcholine into linoleic acid and linolenic acid, could alter the lipid composition and improve the membrane lipid unsaturation. In the present study, the expression of PLA2G in NS was higher than that in NF after freezing stress, and it was hypomethylated in the 2 kb downstream region (Fig. 13; Supplementary Table S10). Plant lipoxygenases are a kind of fatty acid dioxygenases with diverse functions, which not only participate in the peroxidation of linolenic and linoleic acids, and start the synthesis of jasmonic acid [48]. It was reported that jasmonic acid (JA) positively modulates the CBF pathway, leading to regulating stomatal closure and maintaining photosynthesis under cold stress by accumulating cryoprotective compounds and interacting with other plant phytohormones [49, 50]. In the study, eight DMEGs, including three lipoxygenase (LOX2S), two allene oxide cyclase (AOC), one acyl-CoA oxidase (AOCX), one 12-oxophytodienoic acid reductase(OPR), and one enoyl-CoA hydratase/3-hydroxyacyl-CoA dehydrogenase (MFP2), involved in linoleic acid, alpha-linolenic acid metabolism and the JA synthesis, and all genes were up-regulated in NS compared to NF after freezing stress. However, 2 kb the downstream region methylation levels were down-regulated (Fig. 13; Supplementary Table S10).In addition, lipid metabolism generates the lipid signal that activates the downstream cold-responsive genes [51]. This study showed that six genes involved the lipid signal transmission, half the number of genes were up-regulated in NS compared to NF after freezing stress, and it was hypermethylated (Fig. 13; Supplementary Table S10). Our results demonstrated that the fatty acid β-oxidation and JA biosynthesis played a central role in plant freezing tolerance, and DNA methylation regulated the expression of the related gene.

Methylation-regulated carbon fixation and sugar metabolism

Photosynthesis is a primary plant energy source and maintains the normal life activities of cells [52]. The fixation capacity of CO2 can reflect the photosynthetic rate of plants. Photosynthesis was easily affected under cold stress, which leads to excessive energy generation and even triggers photoinhibition [53]. The ribulose-bisphosphate carboxylase (RBCs) is the elite enzyme of the Calvin cycle, which converts airborne CO2 into sugars. In addition, the first reaction of the gluconeogenesis pathway can also fix CO2 to sugar, phosphoenolpyruvate carboxylase (PPC) is the core enzyme in this pathway, and phosphoenolpyruvate carboxykinase (PCA) has the opposite function to supplementary amount of pyruvate. The research indicated that the RBCs and PPC were up-regulated, the PCA was down-regulated in NS compared to NF, and they were hypermethylated in 2 kb upstream position (Fig. 13; Supplementary Table S10). The paper's research indicated that the soluble sugars play an important role in maintaining a cellular osmotic balance of plants under abiotic stress [54, 55]. Soluble sugars mainly include monosaccharides and oligosaccharides, which arose from photosynthetic fixation and hydrolysis of polysaccharides and glycoconjugate. The study indicated that many genes encoding amino sugar, nucleotide sugar, starch, and sucrose metabolism, most of the genes were up-regulated expression after freezing stress and were hypermethylated (Fig. 13; Supplementary Table S10). The beta-fructofuranosidase (INV) is a core enzyme from sucrose/sucrose-6P to β-D-fructose or α-D-glucose 6-phosphate, and the α-D-glucose 6-phosphate could enter another metabolism, facilitating the conversion of polysaccharides into soluble sugars. Simultaneously, beta-glucosidase (bgIX) and beta-amylase (amyB), are key enzymes in the conversion of cellulose to α/β glucose or maltose [56]. Our study showed that the INV, amy B, and four out of five bgIX genes were up-regulated in NS compared to NF after freezing stress. The INV genes was hypermethylated 2 kb upstream, the amyB was hypermethylated 2 kb downstream and four genes bgIX genes were hypermethylated in 2 kb up/downstream or at the gene body.

On the other hand, some sucrose synthases (SS) and alpha-amylases (amyA) boost the synthesis of polysaccharides were down-regulated in NS compared to NF after freezing stress, and the SS was hypomethylated in NS compared to NF under t1 and was hypermethylated in t2 at the 2 kb downstream, but amyA was hypermethylated in NS compared to NF at the 2 kb downstream position after freezing stress. These results suggested that the ability of carbon fixation and polysaccharides hydrolysis of NS was higher than that of NF under freezing stress, and DNA methylation also participated in photosynthesis and sugar metabolism.

Methylation-regulated the ROS scavenging-related metabolites

Although moderate reactive oxygen species (ROS), as the signaling role, can be involved in regulating plant stress responses and improving the ability to tolerate stress, excessive ROS can cause peroxidative injury to organisms [52, 57]. Under adverse biotic and abiotic stresses, the plants have ROS-scavenging enzymes and non-enzymatic scavengers to maintain normal cellular redox homeostasis [53]. Peroxidase and peroxisome, as important enzymes, have the function of scavenging ROS in cold stress [14]. In this study, there were ten and eight identified genes involved in the peroxidase and peroxisome pathways, respectively, and most of them were significantly enriched in NS compared to NF after freezing stress (Fig. 13; Supplementary Table S10). However, their methylation was inconsistent in different regions. In addition, flavonoids are an important member of secondary metabolites that contain various important bioactive substances, including flavonols, anthocyanidins, and some sugars; they function as an antioxidant and protect the plant from a wide range of abiotic stresses, such as heavy metal stress, ultraviolet-B, and salt stress, high-temperature damage, and other environmental stresses [54, 55]. The flavonol synthase (FLS) is an enzyme catalyzing the biosynthesis of quercetin, kaempferol, myricetin, and galanin, which were highly active non-enzymatic scavengers of ROS [58, 59]. The dihydroflavonol 4-reductase (DFR) is an enzyme for the biosynthesis of leucocyanidin, leucodelphinidin, leucopelargonidin, luteoforol, and apiforoletc, which also has an antioxidant role. In this study, two DMEGs (FLS, DFR) expressions were increased in NS compared to NF; among them the methylation levels of FLS were up-regulated under the gene body, and DFR was down-regulated under the 2 kb downstream and was up-regulated under 2 kb upstream (Fig. 13; Supplementary Table S10). The naringenin 3-dioxygenase (F3H) is a core enzyme of the flavonol biosynthesis pathway, which was up-regulated expression in NFt2 vs. NSt2, and it was hypermethylated in downstream2kb (Fig. 13; Supplementary Table S10). Before transcriptome studies, the FLS, DFR, and F3H as the key downstream gene of flavonoid biosynthesis were significantly induced under cold stress in Tetrastigma hemsleyanum [60]. Taken together, these results suggested that rapeseed cultivar NS has a stronger ability to scavenge ROS under freezing stress than the cultivar NF, and the methylation regulates the expression of the ROS scavenging-related gene.

Methylation-regulated plant hormone signal transduction related freezing resistance

Phytohormones, small endogenous signaling molecules, interact with ROS to orchestrates the plant response to abiotic stress and drive changes in transcriptomic, metabolic, and proteomic network changes that lead to plant acclimatization and survival [49]. Previous research has shown that the plant produces large amounts of phytohormones under environmental stress, such as IAA, JA, ABA, BR, and ET, which combine endogenous substances with environmental signals to improve the defense ability against stress [61–, 62,63,64]. Similar results were obtained in this research; 28 DMEGs, encoding some of the responsive proteins or receptors related to hormone signals, were enriched in the plant hormone signal transduction pathway. (Fig. 13; Supplementary Table S10). Studies in Arabidopsis and Rice have shown that high levels of auxin expression can improve their frost resistance, and similar results were also obtained in the present study [61, 65]. The 12DMEGs encode auxin-responsive protein (IAA), SAUR family protein (SAUR), auxin influx carrier (AUX1), and auxin response factor (ARF). At the same time, IAA, AUX1, and two-thirds of SAUR genes were up-regulated expression, and ARFwas down-regulated expression in NS compared to NF under freezing stress(Fig. 13; Supplementary Table S10). The methylation data showed that IAA and AUX1 were hypermethylated in NFt1 vs. NSt1 and NFt2 vs. NSt2 under gene body or at the 2 kb up/downstream, and ARF was hypermethylated in NFt1 vs. NSt1 and NFt2 vs. NSt2 under CHH contexts of upstream2kb (Supplementary Table S10). JA regulates the cold tolerance of plants by inducing jasmonate ZIM domain-containing protein (JAZ), a negative regulator of the JA signaling pathway [62]. Two genes encoding JAZ showed decreased expression in the NFt2 vs. NSt2, whereas methylation level increased in the 2 kb upstream and gene body regions. The abscisic acid receptor PYR/PYL family (PYL) is the receptor of abscisic acid (ABA), and protein phosphatase 2C (PP2C) is an ABA co-receptors, which forms the PYL-ABA-PP2C complex to release the SnRK2 from inhibition by the PP2C [36]. When absence of ABA, the PP2C binding SnRK2 kinases, are inactive in SnRK2 kinases [66]. Recently, it was reported that cold stress activates SnRK2.6/OST1, SnRK2.6 interacts with ICE1 and phosphorylates it to activate the CBF-COR gene-expression cascade, and improving the freezing tolerance of plants [67]. In this study, PP2C and PYL have a down-regulated expression in NS compared to NF after freezing stress, and the methylation levels were inconsistent in different regions. Therefore, we hypothesized that PP2C and PYL inhibition would release more ABA to activate SnRK2 (Fig. 13; Supplementary Table S10). The ABA-responsive element binding factor (ABF) is the main ABA signaling pathway element and positively regulates ABA metabolism. We found that the ABF expression level in NS was higher than in NF after freezing stress, and it was hypomethylated in most regions (Fig. 13; Supplementary Table S10).

The paper showed that brassinolides (BRs) increase plant tolerance to protect plants from damage [63]. BRs are perceived by brassinosteroid insensitive 1 (BRI1), a positive regulator of BR signaling [49]. The cyclin D3 (CYCD3) is an element of the promotive effect of BR on cell division. The result showed that BRI1 and CYCD3 were up-regulated expression and hypermethylated in rapeseed cultivar NS compared to NF after freezing stress (Fig. 12; Supplementary Table S10). Recently, it was shown that the ethylene (ET) signaling pathway transcriptionally inhibits CBF/DREB1 to regulate cold responses of soybean via the action of ethylene insensitive3 (EIN3), contributing to poor cold tolerance [64]. Ethylene receptors (ETRs) as the receptor of ET and negatively regulate the ET signaling pathway [68]. The study showed that the expression levels of EIN3 and ETR2 were lower in NS than in NF, and the methylation level was positively correlated with expression (Fig. 13; Supplementary Table S10). Overall, the phytohormones exert important function in plant freezing tolerance, and DNA methylation is involved to regulates related gene expression.

Conclusion

Overall, we have provided the first systematic exploration of winter rapeseed's global DNA methylation patterns under freezing stress and investigated several novels and important DMEGs and pathways related to its freezing tolerance mechanism of winter rapeseed. DNA methylation under freezing stress regulated many metabolisms associated with winter rapeseed tolerance.

Methods

Plant materials and growth conditions

In this experiment, two winter rapeseed cultivars with different freezing tolerance 17NS57 (NS, freezing-resistant cultivar, with more than 90% overwinter survival rate at − 26 °C) and NQF24 (NF, freezing-sensitive cultivar, with 0% overwinter survival rate below − 10 °C), were provided by the Gansu Agricultural University (Lanzhou, China) [15, 16]. The pot experiment was carried out in a greenhouse. Seeds were germinated in a culture dish on filter papers wetted with deionized water and maintained at 22 ± 1 ℃ until the emergence of the cotyledon. Seedlings of two cultivars were transferred to pots (5L) filled with a 3:1 mixture of nutritional soil and vermiculite and were grown in an illumination incubator under normal conditions (22/20 ℃, day/night temperature, 16 h/8 h light/dark cycle). When the seedlings reached the five leaves stage, they were transferred to the chamber offering freezing treatment at -4 ℃ for 12 h (Treatment 1, t1), and 24 h (Treatment 2, t2), while the control group (control, t0) was maintained in normal condition. After treatment, the third fully expanded leaf from the top of the plant was collected, frozen in liquid nitrogen, and stored at − 80 °C for subsequent analysis.

DNA and RNA extraction and library construction

DNA from the leaves of six samples was separately extracted using the QIAGEN genomic DNA extraction kit, following the manufacturer’s protocol. Then using Illumina’s standard DNA methylation analysis protocol and a DNA Methylation-Gold Kit constructed libraries. Total RNA was extracted from the leaves of six samples containing three biological replicates using the Spectrum™ Plant Total RNA Kit, mRNA was enriched by Oligo (dT) beads, the enriched mRNA was reverse transcripted into cDNA, and sequencing libraries were generated using QiaQuick PCR extraction kit following the manufacturer’s instructions. For whole-genome bisulfite sequencing (WGBS) and RNA sequencing, the libraries were sequenced on the Illumina Hiseq TM 2500 platform by Gene Denovo BiotechnologyCo (Guangzhou, China).

Mapping reads to the reference genome

After sequencing, the raw reads were filtered to remove reads containing more than 10% of unknown nucleotides and low-quality reads containing more than 40% of low-quality (Q-value ≤ 20) bases. The obtained clean reads were mapped to the rapeseed reference genome using BSMAP software (version 2.90) by default [69]. Then a custom Perl script was used to call methylated cytosines, and the methylated cytosines were tested with the correction algorithm described in Lister [70]. The methylation level was calculated based on methylated cytosine percentage in the whole genome, in each chromosome, and in different regions of the genome for each sequence context (CG, CHG, and CHH). The methylation profile at flanking 2 kb regions and gene body (or transposable elements) was plotted based on the average methylation levels for each window to assess different methylation patterns in different genomic regions.

Identification of differentially methylated regions

Differential DNA methylation between the two cultivars at each locus was determined using Pearson's chi-square test (χ2) in methyl Kit (version 1.7.10) [71]. To identify DMRs between two cultivars, the minimum read coverage to call a methylation status for a base was set to 4. A sliding-window approach with a 200-bp window sliding at 100-bp intervals was used to identify DMRs. DMRs for each sequence context (CG, CHG, and CHH) according to different criteria: 1) For CG, numbers of GC in each window ≥ 5, the absolute value of the difference in methylation ratio ≥ 0.25, and q ≤ 0.05; 2); For CHG, numbers in a window ≥ 5, the absolute value of the difference in methylation ratio ≥ 0.25, and q ≤ 0.05; 3); For CHH, numbers in a window ≥ 15, the absolute value of the difference in methylation ratio ≥ 0.15, and q ≤ 0.05; 4) For all C, numbers in a window ≥ 20, the absolute value of the difference in methylation ratio ≥ 0.2, and q ≤ 0.05. Genes overlapping with significant DMRs for at least 1 bp in the flanking 2 kb regions and gene body were defined as differentially methylated genes (DMGs), which be used in subsequent analyses.

RNA sequencing and data analysis

Transcriptome sequencing was performed from the same materials as methylation sequencing containing three biological replicates. The reads containing adapters, reads containing poly-N, and low-quality reads from the raw data were removed using the Trimmomatic software (Bolger et al., 2014). The short reads alignment tool Bowtie2 was used for mapping reads to the ribosome RNA (rRNA) database. The rRNA mapped reads were removed. The remaining reads or unmapped reads were aligned to the rapeseed reference genome (GCA_000686985.2) using TopHat2 (version 2.0.9) [72]. The expression level was normalized by calculating the fragments per kilobase of exon model per million mapped fragments (FPKM) value. To identify differentially expressed genes across samples or groups, the edgeR package (http://www.r-project.org/) was used. All significantly altered genes were identified using the criteria of a p-value < 0.001 and a value of |log2foldchange|≥ 2.

Functional annotation and analysis of DMEGs

To analyze the correlation between altered DNA methylation patterns and gene expression, we examined the overlap between the DMGs and DEGs of two cultivars, named differentially methylated regions related to differentially expressed genes (DMEGs). Then Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted for DMEGs.

Validation of DMEGs by quantitative qRT-PCR

Total RNA was extracted from the leaves of each sample using the Steadypure Plant RNA Extraction Kit. The first-strand cDNA was synthesized using RNA, according to the manufacturer’s instructions of Evo M-MLV RT Premix for Qrt-PCR, and qRT-PCR amplification reactions were performed using an SYBR® Green Premix Pro Taq HS qPCR Kit. The 2−ΔΔCT method was used to calculate the gene expression, the actin gene (Bra028615) was used as the reference gene for normalization of gene expression, and each sample was replicated three times. The gene names and specific primer sequences were detailed in Supplementary Table S1.

Availability of data and materials

The sequenced transcriptome raw data and Whole-genome bisulfite sequencing raw data have been deposited to the SRA at NCBI with the accession number PRJNA685002.

Abbreviations

DMRs:

Differentially methylated regions

DMGs:

Differentially methylated genes

DEGs:

Differentially expressed genes

DMEGs:

DEGs associated with DMGs

RNA-Seq:

Ribonucleic acid sequencing

qRT-PCR:

Quantitative Real-Time PCR

WGBS:

Whole-genome bisulfite sequencing

MSAP:

Methylation sensitive amplification polymorphism

MSP:

Methylation specific PCR

TE:

Transposable elements

rRNA:

Ribosome RNA

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

References

  1. Takahashi D, Uemura M, Kawamura Y. Freezing tolerance of plant cells: from the aspect of plasma membrane and microdomain. Adv Exp Med Biol. 2018;1081:61–79.

    Article  CAS  PubMed  Google Scholar 

  2. Dong X, Liu Z, Mi W, Xu C, Xu M, Zhou Y, Zhen G, Cao X, Fang X, Mi C. Overexpression of BrAFP1 gene from winter rapeseed (Brassica rapa) confers cold tolerance in Arabidopsis. Plant Physiol Biochem. 2020;155:338–45.

    Article  CAS  PubMed  Google Scholar 

  3. Ramanjulu S, Bartels D. Drought- and desiccation-induced modulation of gene expression in plants. Plant, Cell Environ. 2002;25(2):141–51.

    Article  CAS  Google Scholar 

  4. Li H, Jiang X, Lv X, Ahammed G, Guo Z, Qi Z, Yu J, Zhou Y. Tomato GLR3.3 and GLR3.5 mediate cold acclimation-induced chilling tolerance by regulating apoplastic H O production and redox homeostasis. Plant, cell & environment. 2019;42(12):3326–39.

    Article  CAS  Google Scholar 

  5. Song Y, Liu L, Feng Y, Wei Y, Yue X, He W, Zhang H, An L. Chilling- and freezing-induced alterations in cytosine methylation and its association with the cold tolerance of an alpine subnival plant, chorispora bungeana. PLoS ONE. 2015;10(8):e0135485.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  6. Mayer B, Ali-Benali M, Demone J, Bertrand A, Charron J. Cold acclimation induces distinctive changes in the chromatin state and transcript levels of COR genes in Cannabis sativa varieties with contrasting cold acclimation capacities. Physiol Plant. 2015;155(3):281–95.

    Article  CAS  PubMed  Google Scholar 

  7. Friedrich T, Faivre L, Bäurle I, Schubert D. Chromatin-based mechanisms of temperature memory in plants. Plant, Cell Environ. 2019;42(3):762–70.

    Article  CAS  Google Scholar 

  8. Chang Y, Zhu C, Jiang J, Zhang H, Zhu J, Duan C. Epigenetic regulation in plant abiotic stress responses. J Integr Plant Biol. 2020;62(5):563–80.

    Article  CAS  PubMed  Google Scholar 

  9. Zhang M, Zhang X, Guo L, Qi T, Liu G, Feng J, Shahzad K, Zhang B, Li X, Wang H, et al. Single-base resolution methylome of cotton cytoplasmic male sterility system reveals epigenomic changes in response to high-temperature stress during anther development. J Exp Bot. 2020;71(3):951–69.

    CAS  PubMed  Google Scholar 

  10. Shan X, Wang X, Yang G, Ying W, Su S, Li S, Liu H, Yuan Y. Analysis of the DNA methylation of maize ( Zea mays L.) in response to cold stress based on methylation-sensitive amplified polymorphisms. J Plant Biol. 2013;56(1):32–8.

    Article  CAS  Google Scholar 

  11. Ahmad F, Farman K, Waseem M, Rana R, Nawaz M, Rehman H, Abbas T, Baloch F, Akrem A, Huang J, et al. Genome-wide identification, classification, expression profiling and DNA methylation (5mC) analysis of stress-responsive ZFP transcription factors in rice (Oryza sativa L.). Gene. 2019;718:144018.

    Article  CAS  PubMed  Google Scholar 

  12. Liu T, Li Y, Duan W, Huang F, Hou X. Cold acclimation alters DNA methylation patterns and confers tolerance to heat and increases growth rate in Brassica rapa. J Exp Bot. 2017;68(5):1213–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Tong W, Li R, Huang J, Zhao H, Ge R, Wu Q, Mallano A, Wang Y, Li F, Deng W, et al. Divergent DNA methylation contributes to duplicated gene evolution and chilling response in tea plants. The Plant journal : for cell and molecular biology. 2021;106(5):1312–27.

    Article  CAS  Google Scholar 

  14. Qi W, Wang F, Ma L, Qi Z, Liu S, Chen C, Wu J, Wang P, Yang C, Wu Y, et al. Physiological and Biochemical Mechanisms and Cytology of Cold Tolerance in Brassica napus. Front Plant Sci. 2020;11:1241.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Wei J, Zheng G, Yu X, Liu S, Dong X, Cao X, Fang X, Li H, Jin J, Mi W, et al. Comparative Transcriptomics and Proteomics Analyses of Leaves Reveals a Freezing Stress-Responsive Molecular Network in Winter Rapeseed (Brassica rapa L.). Front Plant Sci. 2021;12:664311.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Cao X, Liu Z, Mi W, Xu C, Zou Y, Xu M, Zheng G, Fang X, Cui X, Dong X, et al. Analysis on the Adaptability of Northward Planting of Brassica napus. Scientia Agricultura Sinica. 2020;53(20):4164–76.

    Google Scholar 

  17. Lei Y. Physiological and molecular responses to cold stress in rapeseed (Brassica napus L.). J Integr Agric. 2019;18(12):2742–52.

    Article  Google Scholar 

  18. Pu Y, Zhao Y, Wu J, Liu L, Bai J, Ma L, Niu Z, Jin J, Fang Y, Li X, et al. Comprehensive assessment on cold tolerance of the strong winter brassica napus L. cultivated in Northern China. Scientia Agricultura Sinica. 2019;52(19):3291–308.

    Google Scholar 

  19. Ma L, Coulter J, Liu L, Zhao Y, Chang Y, Pu Y, Zeng X, Xu Y, Wu J, Fang Y, et al. Brassica rapatranscriptome analysis reveals key cold-stress-responsive genes in winter rapeseed (Brassica rapa L.). IntJ Mol Sci. 2019;20(5):1071.

    Article  CAS  Google Scholar 

  20. Pu Y, Liu L, Wu J, Zhao Y, Bai J, Ma L, Yue J, Jin J, Niu Z, Fang Y, et al. Brassica napustranscriptome profile analysis of winter rapeseed ( L.) in response to freezing stress, reveal potentially connected events to freezing stress. Int J Mol Sci. 2019;20(11):2271.

    Article  CAS  Google Scholar 

  21. Mi W, Liu Z, Jin J, Dong X, Xu C, Zou Y, Xu M, Zheng G, Cao X, Fang X, et al. Comparative proteomics analysis reveals the molecular mechanism of enhanced cold tolerance through ROS scavenging in winter rapeseed (Brassica napus L.). PloS one. 2021;16(1):e0243292.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Wei J, Zheng G, Dong X, Li H, Liu Z. Integration of Transcriptome and proteome analysis reveals the mechanism of freezing tolerance in winter rapeseed. Plant Growth Regul. 2021;96(1):103–18.

    Article  CAS  Google Scholar 

  23. Xu J, Zhou S, Gong X, Song Y, van Nocker S, Ma F, Guan Q. Single-base methylome analysis reveals dynamic epigenomic differences associated with water deficit in apple. Plant Biotechnol J. 2018;16(2):672–87.

    Article  CAS  PubMed  Google Scholar 

  24. Yaish M, Colasanti J, Rothstein S. The role of epigenetic processes in controlling flowering time in plants exposed to stress. J Exp Bot. 2011;62(11):3727–35.

    Article  CAS  PubMed  Google Scholar 

  25. Du J, Zhong X, Bernatavichute Y, Stroud H, Feng S, Caro E, Vashisht A, Terragni J, Chin H, Tu A, et al. Dual binding of chromomethylase domains to H3K9me2-containing nucleosomes directs DNA methylation in plants. Cell. 2012;151(1):167–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Xing MQ, Zhang YJ, Zho SR, Hu WY, Wu XT. Global analysis reveals the crucial roles of DNA Methylation during rice seed development. Plant Physiol. 2015;168(4):1417–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Kang Y, Bae A, Shim S, Lee T, Lee J, Satyawan D, Kim M, Lee S. Genome-wide DNA methylation profile in mungbean. Sci Rep. 2017;7:40503.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Li Y, Ding X, Wang X, He T, Zhang H, Yang L, Wang T, Chen L, Gai J, Yang S. Genome-wide comparative analysis of DNA methylation between soybean cytoplasmic male-sterile line NJCMS5A and its maintainer NJCMS5B. BMC Genomics. 2017;18(1):596.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  29. Yan H, Kikuchi S, Neumann P, Zhang W, Wu Y, Chen F, Jiang J. Genome-wide mapping of cytosine methylation revealed dynamic DNA methylation patterns associated with genes and centromeres in rice. Plant J. 2010;63(3):353–65.

    Article  CAS  PubMed  Google Scholar 

  30. Wang Z, Wu X, Wu Z, An H, Yi B, Wen J, Ma C, Shen J, Fu T, Tu J. Genome-Wide DNA methylation comparison between brassica napus genic male sterile line and restorer line. Int J Mol Sci. 2018;19(9):2689.

    Article  PubMed Central  CAS  Google Scholar 

  31. Niederhuth C, Bewick A, Ji L, Alabady M, Kim K, Li Q, Rohr N, Rambani A, Burke J, Udall J, et al. Widespread natural variation of DNA methylation within angiosperms. Genome Biol. 2016;17(1):194.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  32. Ausin I, Feng S, Yu C, Liu W, Kuo H, Jacobsen E, Zhai J, Gallego-Bartolome J, Wang L, Egertsdotter U, et al. DNA methylome of the 20-gigabase Norway spruce genome. Proc Natl Acad Sci USA. 2016;113(50):E8106–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Li J, Huang Q, Sun M, Zhang T, Li H, Chen B, Xu K, Gao G, Li F, Yan G, et al. Global DNA methylation variations after short-term heat shock treatment in cultured microspores of Brassica napus cv. Topas Scientific reports. 2016;6:38401.

    Article  CAS  PubMed  Google Scholar 

  34. Mazin A, Boĭko L, Ogarkova O, Vaniushin B. Loss of CpG dinucleotides from DNA. VI. Methylation of mitochondrial and chloroplast genes. Mol biol. 1988;22(6):1688–96.

    CAS  Google Scholar 

  35. Chinnusamy V, Zhu J, Zhu JK. Cold stress regulation of gene expression in plants. Trends Plant Sci. 2007;12(10):444–51.

    Article  CAS  PubMed  Google Scholar 

  36. Zhu JK. Abiotic stress signaling and responses in plants. Cell. 2016;167(2):313–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Park S, Lee C, Doherty C, Gilmour S, Kim Y, Thomashow M. Regulation of the Arabidopsis CBF regulon by a complex low-temperature regulatory network. Plant J Cell Mol Biol. 2015;82(2):193–207.

    Article  CAS  Google Scholar 

  38. Xie Y, Chen P, Yan Y, Bao C, Li X, Wang L, Shen X, Li H, Liu X, Niu C, et al. An atypical R2R3 MYB transcription factor increases cold hardiness by CBF-dependent and CBF-independent pathways in apple. New Phytol. 2018;218(1):201–18.

    Article  CAS  PubMed  Google Scholar 

  39. Pan F, Wu M, Hu W, Liu R, Yan H, Xiang Y. Phyllostachys edulisGenome-Wide Identification and expression analyses of the bZIP transcription factor genes in moso bamboo (Phyllostachys edulis). Int J Mol Sci. 2019;20(9):2203.

    Article  CAS  PubMed Central  Google Scholar 

  40. Pradhan S, Pandit E, Nayak D, Behera L, Mohapatra T. Genes, pathways and transcription factors involved in seedling stage chilling stress tolerance in indica rice through RNA-Seq analysis. BMC Plant Biol. 2019;19(1):352.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  41. Waititu J, Cai Q, Sun Y, Sun Y, Li C, Zhang C, Liu J, Wang H. Transcriptome profiling of maize ( Zea mays L.) leaves reveals key cold-responsive genes, transcription factors, and metabolic pathways regulating cold stress tolerance at the seedling stage. Genes. 2021;12(10):1638.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Kumar S, Seem K, Kumar S, Vinod K, Chinnusamy V, Mohapatra T. Pup1 QTL regulates gene expression through epigenetic modification of dna under phosphate starvation stress in rice. Front Plant Sci. 2022;13:871890.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Yuan P, Yang T, Poovaiah WB. Calcium Signaling-Mediated Plant Response to Cold Stress. Int J Mol Sci. 2018;19(12):3896.

    Article  PubMed Central  CAS  Google Scholar 

  44. Liu Z, Jia Y, Ding Y, Shi Y, Li Z, Guo Y, Gong Z, Yang S. Plasma membrane CRPK1-mediated phosphorylation of 14-3-3 proteins Induces their nuclear Import to fine-tune cbf signaling during cold response. Mol Cell. 2017;66(1):117-128.e115.

    Article  CAS  PubMed  Google Scholar 

  45. Zhao C, Wang P, Si T, Hsu CC, Wang L, Zayed O, Yu Z, Zhu Y, Dong J, Tao WA. MAP Kinase cascades regulate the cold response by modulating ICE1 protein stability. Dev Cell. 2017;43(5):618–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Degenkolbe T, Giavalisco P, Zuther E, Seiwert B, Hincha D, Willmitzer L. Differential remodeling of the lipidome during cold acclimation in natural accessions of Arabidopsis thaliana. Plant J Cell Mol Biol. 2012;72(6):972–82.

    Article  CAS  Google Scholar 

  47. de Mendoza D. Temperature sensing by membranes. Annu Rev Microbiol. 2014;68:101–16.

    Article  PubMed  CAS  Google Scholar 

  48. Zhao Y, Dong W, Zhang N, Ai X, Wang M, Huang Z, Xiao L, Xia G. A wheat allene oxide cyclase gene enhances salinity tolerance via jasmonate signaling. Plant Physiol. 2014;164(2):1068–76.

    Article  CAS  PubMed  Google Scholar 

  49. Verma V, Ravindran P, Kumar P. Plant hormone-mediated regulation of stress responses. BMC Plant Biol. 2016;16:86.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  50. Hu Y, Jiang Y, Han X, Wang H, Pan J, Yu D. Jasmonate regulates leaf senescence and tolerance to cold stress: crosstalk with other phytohormones. J Exp Bot. 2017;68(6):1361–9.

    Article  CAS  PubMed  Google Scholar 

  51. Hou Q, Ufer G, Bartels D. Lipid signalling in plant responses to abiotic stress. Plant, Cell Environ. 2016;39(5):1029–48.

    Article  CAS  Google Scholar 

  52. Mittler R. ROS Are Good. Trends Plant Sci. 2017;22(1):11–9.

    Article  CAS  PubMed  Google Scholar 

  53. Ray P, Huang B, Tsuji Y. Reactive oxygen species (ROS) homeostasis and redox regulation in cellular signaling. Cell Signal. 2012;24(5):981–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Zhang H, Wu Z, Suo Y, Wang J, Zheng L. Gene expression and flavonol biosynthesis are induced by ultraviolet-B and salt stresses in Reaumuria trigyna. Biol Plant. 2017;61(2):246–54.

    Article  CAS  Google Scholar 

  55. Muhlemann J, Younts T, Muday G. Flavonols control pollen tube growth and integrity by regulating ROS homeostasis during high-temperature stress. Proc Natl Acad Sci USA. 2018;115(47):E11188–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Bauer S, Vasu P, Persson S, Mort A, Somerville C. Development and application of a suite of polysaccharide-degrading enzymes for analyzing plant cell walls. Proc Natl Acad Sci USA. 2006;103(30):11417–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Baxter A, Suzuki N, Mittler R. ROS as key players in plant stress signalling. (Special Issue: Oxidative stress and cell death.). J Exp Bot. 2014;65(5):1229–40.

    Article  CAS  PubMed  Google Scholar 

  58. Noratto GD, Kim Y, Talcott ST, Mertens-Talcott SU. Flavonol-rich fractions of yaupon holly leaves (Ilex vomitoria, Aquifoliaceae) induce microRNA-146a and have anti-inflammatory and chemopreventive effects in intestinal myofribroblast CCD-18Co cells. Fitoterapia. 2011;82(4):557–69.

    Article  CAS  PubMed  Google Scholar 

  59. Xu D, Hu M, Wang Y, Cui Y. Antioxidant activities of quercetin and its complexes for medicinal application. Mol (Basel, Switzerland). 2019;24(6):1123.

    Article  CAS  Google Scholar 

  60. Peng X, Wu H, Chen H, Zhang Y, Qiu D, Zhang Z. Transcriptome profiling reveals candidate flavonol-related genes of Tetrastigma hemsleyanum under cold stress. BMC Genomics. 2019;20(1):687.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  61. Jain M, Khurana J. Transcript profiling reveals diverse roles of auxin-responsive genes during reproductive development and abiotic stress in rice. FEBS J. 2009;276(11):3148–62.

    Article  CAS  PubMed  Google Scholar 

  62. Hu Y, Jiang L, Wang F, Yu D. Jasmonate regulates the inducer of cbf expression-C-repeat binding factor/DRE binding factor1 cascade and freezing tolerance in Arabidopsis. Plant Cell. 2013;25(8):2907–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Xia X, Fang P, Guo X, Qian X, Zhou J, Shi K, Zhou Y, Yu J. Brassinosteroid-mediated apoplastic H O -glutaredoxin 12/14 cascade regulates antioxidant capacity in response to chilling in tomato. Plant, Cell Environ. 2018;41(5):1052–64.

    Article  CAS  Google Scholar 

  64. Robison J, Yamasaki Y, Randall S. Glycine maxThe ethylene signaling pathway negatively impacts CBF/DREB-regulated cold response in Soybean (Glycine max). Front Plant Sci. 2019;10:121.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Hannah M, Heyer A, Hincha D. A global survey of gene regulation during cold acclimation in Arabidopsis thaliana. PLoS Genet. 2005;1(2):e26.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  66. Soon F-F, Ng L-M, Zhou XE, West GM, Kovach A, Tan MHE, Suino-Powell KM, He Y, Xu Y, Chalmers MJ, et al. Molecular mimicry regulates ABA Signaling by SnRK2 Kinases and PP2C Phosphatases. Science. 2012;335(6064):85–8.

    Article  CAS  PubMed  Google Scholar 

  67. Ding Y, Li H, Zhang X, Xie Q, Gong Z, Yang S. OST1 kinase modulates freezing tolerance by enhancing ICE1 stability in Arabidopsis. Dev Cell. 2015;32(3):278–89.

    Article  CAS  PubMed  Google Scholar 

  68. Shi Y, Tian S, Hou L, Huang X, Zhang X, Guo H, Yang S. Ethylene signaling negatively regulates freezing tolerance by repressing expression of CBF and type-A ARR genes in Arabidopsis. Plant Cell. 2012;24(6):2578–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Xi Y, Li W. BSMAP: whole genome bisulfite sequence MAPping program. BMC Bioinformatics. 2009;10:232.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

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

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Akalin A, Kormaksson M, Li S, Garrett-Bakelman F, Figueroa M, Melnick A, Mason C. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13(10):R87.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg S. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013;14(4):R36.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

Download references

Acknowledgements

We are grateful to Guangzhou Genedenovo Biotechnology Co, Ltd, for assisting in the sequencing and bioinformatics analysis.

Funding

The research was supported by the State Key Laboratory of Aridland Crop Science, Gansu Agricultural University (GSCS-2021-Z01), Corresponding projects of the overall planning fund of the Institute of Agricultural Sciences of the Tibetan Academy of Agricultural Sciences(2021NKS-TCJJ-001), Industrial Support Plan Project of Gansu (2021CYZC-46), Young Doctoral Fund of Gansu (2021QB-035), Special funds from the central government to guide local scientific and technological development of China (ZCYD-2020–2-3), Ministry of Science and Technology of China (2018YFD0100500), National Natural Science Foundation of China (31660404).

Author information

Authors and Affiliations

Authors

Contributions

The work presented here was carried out in collaboration among all authors. Zigang Liu and Guoqiang Zheng designed the experiment. Guoqiang Zheng executed the omics data analysis and wrote the manuscript. Jiaping Wei, Ali Aslam, and Junmei Cui polished and revised this paper. Xiaoyun Dong, Hui Li, Ying Wang, Haiyan Tian, and Xiaodong Cao participated in the qRT-PCR of DMEGs to validate the reliability of RNA-seq data. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Zigang Liu.

Ethics declarations

Ethics approval and consent to participate

All experimental studies on plants were complied with relevant institutional, national, and international guidelines and legislation.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing financial and non-finanical interests.

Additional information

Publisher’s Note

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

Supplementary Information

Additional file 1. Table S1.1

The qPCR rimers used in this study. Table S1.2 The effective coverage rates of all chromosomes of all samples. Table S1.3 Number of differentially methylated regions identified in CG, CHG, and CHH contexts between two cultivars of rapeseed. Table S1.4 Distribution of methylation levels in different chromosomes.

Additional file 2. Table S2.1

The DMRs between two cultivars under CG context in control. Table S2.2 The DMRs between two cultivars under CHG context in control. Table S2.3 The DMRs between two cultivars under CHH context in control. Table S2.4 The DMRs between two cultivars under CG context in treatment1. Table S2.5 The DMRs between two cultivars under CHG context in treatment1. Table S2.6 The DMRs between two cultivars under CHH context in treatment1. Table S2.7 The DMRs between two cultivars under CG context in treatment2. Table S2.8 The DMRs between two cultivars under CHG context in treatment2. Table S2.9 The DMRs between two cultivars under CHH context in treatment2.

Additional file 3. Table S3.1

The DMGs between two cultivars under CG context in control. Table S3.2 The DMGs between two cultivars under CHG context in control. Table S3.3 The DMGs between two cultivars under CHH context in control. Table S3.4 The DMGs between two cultivars under CG context in treatment1. Table S3.5 The DMGs between two cultivars under CHG context in treatment1. Table S3.6 The DMGs between two cultivars under CHH context in treatment1. Table S3.7 The DMGs between two cultivars under CG context in treatment2. Table S3.8 The DMGs between two cultivars under CHG context in treatment2. Table S3.9 The DMGs between two cultivars under CHH context in treatment2.

Additional file 4. 

Table S4.1 The data of RNA-seq in two cultivars of rapeseed. Table S4.2 Number of differentially expression genes between two cultivars of rapeseed. Table S4.3 Correlation analysis between the three replicates of RNA-seq all libraries.

Additional file 5. Table S5.1

The DEGs between two cultivars in control. Table S5.2 The DEGs between two cultivars in treatment1. Table S5.3 The DEGs between two cultivars in treatment2.

Additional file 6. Table S6.1

The DMEGs between two cultivars under CG context in control. Table S6.2 The DMEGs between two cultivars under CHG context in control. Table S6.3 The DMEGs between two cultivars under CHH context in control. Table S6.4 The DMEGs between two cultivars under CG context in treatment1. Table S6.5 The DMEGs between two cultivars under CHG context in treatment1. Table S6.6 The DMEGs between two cultivars under CHH context in treatment1. Table S6.7 The DMEGs between two cultivars under CG context in treatment2. Table S6.8 The DMEGs between two cultivars under CHG context in treatment2. Table S6.9 The DMEGs between two cultivars under CHH context in treatment2. 

Additional file 7. Table S7.1

GO classification annotation of 1732DMEGs. Table S7.2 KEGG pathways analysis of 1732 DMEGs.

Additional file 8. Table S8.1

Identification of differentially expressed transcription factors for DMEGs. Table S8.2 Result of TFs between two cultivars in treatment1. Table S8.3 Result of TFs between two cultivars in treatment2. Table S8.4 The correlation coefficient between the expression level and the methylation level of the gene encoding the transcription factor. Table S8.5 Result of kinases between two cultivars in treatment1. Table S8.6 Result of kinases between two cultivars in treatment2.

Additional file 9. Table S9.1

Data of RT-qPCR. Table S9.2 Data of RT-qPCR and RNA-seq.

Additional file 10. Table S10.1

Result of DMEGs involved in phosphate metabolism under treatment1. Table S10.2 Result of DMEGs involved in phosphate metabolism under treatment2. Table S10.3 Result of DMEGs involved in lipid signal under treatment1. Table S10.4 Result of DMEGs involved in lipid signal under treatment2. Table S10.5 Result of DMEGs involved in carbon fixation under treatment1. Table S10.6 Result of DMEGs involved in carbon fixation under treatment2. Table S10.7 Result of DMEGs involved in sugar metabolism under treatment1. Table S10.8 Result of DMEGs involved in sugar metabolism under treatment2. Table S10.9 Result of DMEGs encoded peroxidase and peroxisome enzymes under treatment1. Table S10.10 Result of DMEGs encoded peroxidase and peroxisome enzymes under treatment2. Table S10.11 Result of DMEGs involved in flavonoid biosynthesis under treatment1. Table S10.12 Result of DMEGs involved in flavonoid biosynthesis under treatment2. Table S10.13 Result of DMEGs involved in plant hormone signal transduction under treatment1. Table S10.14 Result of DMEGs involved in plant hormone signal transduction under treatment2.

Additional file 11. Figure S1.

Distribution of methylation levels in different genomic functional regions of two rapeseed cultivars. A, CG context; B, CHG context; C, CHH context.

Additional file 12. Figure S2.

Distribution of methylation levels in different transposable elements of two rapeseed cultivars. A, CG context; B, CHG context; C, CHH context.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 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 in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Zheng, G., Dong, X., Wei, J. et al. Integrated methylome and transcriptome analysis unravel the cold tolerance mechanism in winter rapeseed(Brassica napus L.). BMC Plant Biol 22, 414 (2022). https://doi.org/10.1186/s12870-022-03797-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-022-03797-1

Keywords

  • Winter rapeseed
  • Freezing stress
  • DNA methylation
  • Freezing tolerance
  • Transcriptome