Skip to main content

Transcription-associated metabolomic profiling reveals the critical role of frost tolerance in wheat



Low temperature is a crucial stress factor of wheat (Triticum aestivum L.) and adversely impacts on plant growth and grain yield. Multi-million tons of grain production are lost annually because crops lack the resistance to survive in winter. Particularlly, winter wheat yields was severely damaged under extreme cold conditions. However, studies about the transcriptional and metabolic mechanisms underlying cold stresses in wheat are limited so far.


In this study, 14,466 differentially expressed genes (DEGs) were obtained between wild-type and cold-sensitive mutants, of which 5278 DEGs were acquired after cold treatment. 88 differential accumulated metabolites (DAMs) were detected, including P-coumaroyl putrescine of alkaloids, D-proline betaine of mino acids and derivativ, Chlorogenic acid of the Phenolic acids. The comprehensive analysis of metabolomics and transcriptome showed that the cold resistance of wheat was closely related to 13 metabolites and 14 key enzymes in the flavonol biosynthesis pathway. The 7 enhanced energy metabolites and 8 up-regulation key enzymes were also compactly involved in the sucrose and amino acid biosynthesis pathway. Moreover, quantitative real-time PCR (qRT-PCR) revealed that twelve key genes were differentially expressed under cold, indicating that candidate genes POD, Tacr7, UGTs, and GSTU6 which were related to cold resistance of wheat.


In this study, we obtained the differentially expressed genes and differential accumulated metabolites in wheat under cold stress. Using the DEGs and DAMs, we plotted regulatory pathway maps of the flavonol biosynthesis pathway, sucrose and amino acid biosynthesis pathway related to cold resistance of wheat. It was found that candidate genes POD, Tacr7, UGTs and GSTU6 are related to cold resistance of wheat. This study provided valuable molecular information and new genetic engineering clues for the further study on plant resistance to cold stress.

Peer Review reports


Wheat (Triticum aestivum L.) is one of the largest planted proportion crops globally, and is essential for human nutrition and animal feed. China is one of the countries where wheat is grown most widely in the world.. Wheat is often affected by different types of typically abiotic stresses such as cold injury and drought during growth and development, directly influencing wheat’s yield [1]. Climate-driven extreme temperatures affect the nutritional and reproductive growth of wheat, which in turn leads to a further decline in yield [2].. Cold is a climatic element that astrict crop growth and affects crop yield. Plants can respond to cold stress through series of physiological and biochemical changes [3]. Understanding the response of plants to cold stress is very important for crop molecular breeding [4,5,6,7]. Low temperature can cause irreparable damage to plants. Overstimulation of electron transport systems cause the enhanced concentration of reactive oxygen species (ROS), which affects membrane lipids, causes damage to proteins and nucleic acids, and ultimately results in cell death [8, 9]. When plants have appeared too low temperatures, the contents of metabolites, such as carbon [10], sugar [11], and protein [8], will be self-redistributed as the temperature decreases. Under cold stress, the generation of osmoprotective agents (trehalose, raffinose, proline, fructose, and betaine) can stabilize the protein and cell structure, which are widely found in plant species [12, 13]. Plant responses to cold stress can also protect cell metabolism by synthesizing late embryogenesis abundance (LEA) protein [14, 15], which stabilize cell structure by binding to DNA or RNA to form a compact hydrogen bond structure and reduce damage caused by abiotic stress [16].

Progress has been made in developing cold-response genes to increase tolerance to abiotic stresses in crops [17,18,19]. The response of plants to cold stress involves multiple levels of regulation, including molecular level, transcriptional regulation, metabolic changes and response pathways. Frost-resistant stress signaling and transcriptional regulation of wheat are complex, involving at least two major pathways (frost response and vernalization) and many other processes involving large numbers of genes (calcium-binding protein (CBP) [20, 21], Ethylene response element-binding protein (EREBP) [22, 23], cold-responsive/late embryogenesis abundant (COR/LEA) genes [24, 25], cold- and dehydration-responsive genes (CRT/DRE) [26], vernalization (VRN1, VRN2, VRN3) [27, 28], frost resistance loci (FR1, FR2) [29]. In the plant genome, about 7% of the coding sequences are associated with transcription factors (TFs) [30]. Wheat genomes contain a large number of transcription factors that play essential roles in biotic and abiotic stress biological processes, such as C-repeat-binding factor (CBF1, CBF2, and CBF3) [27, 31, 32], C-repeat (CRT) [26, 33]. Cold stresses cause cells to dehydrate, so various signaling pathways preventing cell dehydration have been studied, including MYB, WRKY, APETALA 2/ethylene response element-binding protein (AP2/EREBP), NAC, and alkaline leucine (Leu) zipper (BZIP) families [34, 35].

New approaches are needed to uncover the molecular mechanisms of wheat cold tolerance. Omics technology offers unique opportunities to study the regulatory mechanisms of complex networks related to cold hardiness [36]. With the speediness evolution of sequencing tools, transcriptome sequencing technology has been extensively applied to clarify gene function by high-throughput analysis of genomic information, and provided a potent way to systematically study the relationship between genotype and phenotype in transcriptional regulation [37]. For example, it has been widely used in research wheat stress resistance to reveal differentially expressed genes in different biological processes, such as cold [38,39,40], drought [41], heat [42], and salt tolerance [43]. However, low abundance and unknown transcripts are unable to distinguish [37]. Compared with conventional chemical analysis, metabolomics was promoted using liquid chromatography/mass spectrometry (LC/MS) by detecting lots of compounds [44]. Recently, the comprehensive analysis of metabolome and transcriptome has been widely used to study the regulatory network and correlation between genes and metabolites [35, 45]. Metabonomic, a concept first proposed by Nicholson et al., is a new omics following transcriptome, genomics and proteomic-related investigations [46].

Metabolomic can be used for qualitative and quantitative analysis of the small molecular metabolites in cells and organisms at a specific physiological period, and has been certified as a mighty instrumentality for studying gene modification and environmentally induced metabolite changes [47]. Currently, metabolomic analysis has been widely used to reveal metabolites profiles changes under low-temperature stress in many plants. Kaplan et al. [48] used non-target gas chromatography-mass spectrometry (GC-MS) to study the dynamic Arabidopsis thaliana changes under heat and cold stress, and discovered that cold shock had a more significant effect on metabolites. Cold acclimation analysis of two genotypes (cold sensitivity and cold tolerance) in perennial plants was studied by the nuclear magnetic resonance (NMR) method, and more than 40 metabolites were identified [49]. Metabolomics analysis of tillering node in winter wheat was conducted by gas chromatography time-of-flight mass spectrometry (GC-TOF/MS) during the overwintering period, and confirmed that monosaccharides, amino acids and lipids maintained carbon and nitrogen balance [8]. Zhao et al. [50] investigate changes in gene/metabolite activity in Jing 411 by integrated transcriptomics and metabolomics analyses, and indicated that critical pathways associated with ABA/JA signaling and proline biosynthesis played important roles in regulating cold tolerance in wheat.

Nevertheless, studies about the transcriptional and metabolic mechanisms underlying cold stresses in wheat are limited. This study investigated the influence of low-temperature stress on wheat seedlings by physiological, transcriptome and metabonomics analysis of Jimai 325 and its mutant MU-134. Widely Targeted Metabolome technology was applied to reveal the metabolic changes affected by different cold stresses. The findings suggest that candidate genes and metabolites involved in the essential biological pathways might regulate frost resistance in wheat, and a metabolic regulatory network was established in response to low-temperature stress. This study deepened the understanding of the molecular mechanism of wheat adaptive response to low-temperature stress and provided valuable information for wheat frost tolerance gene breeding.


Phenotypic variation under cold stress

To study the responses to cold stress, 28-days-old plants were exposed to cold stresses for progressive 7-day-period treatments. The cold stress significantly affected wheat’s normal growth and development and reduced the survival rate (Fig. 1A, B). To understand the effect of low temperature stress on wheat, several physiological characteristics between the mutant plant and Jimai325 were investigated. In our research, the activity of the peroxidase (POD), proline (Pro) content and soluble sugar contents significantly increased in Jimai325 (CK) and MU-134 when treated with cold stress compared with the untreated control plants (Fig. 1C, E, G). Catalase activity (CAT) and Superoxidase activity (SOD) were significantly increased in CK when treated with cold stress compared with the untreated control plants. Still, there was no significant difference in MU-134(Fig. 1D, H). Malondialdehyde content only significantly increased in MU-134 when treated with cold stress compared with the untreated control plants (Fig. 1F).

Fig. 1
figure 1

Cold Stress effects on plant development. A Shoot morphology. B Survival rate (%). C Proline content (ug/g). D, Catalase activity (μmoL*min*g)−1 FW. E, Soluble sugar contents (mg/g FW). F Malondialdehyde content (nmol/g FW). G Peroxidase activity (μmoL/(min*g) FW). H Superoxidase activity (U/(min*g) FW). The data are means ± SD(n ≥ 10). Different letters indicate significant differences between treatments and control plants from Student’s unpaired two-tailed t-test (P < 0.05). CK, Jimai 325; MU, MU-134

Global transcriptomic changes under cold stress in wheat

In order to identify the cold tolerance genes, RNA-seq was performed on MU and CK samples at normal and after low-temperature treatment, respectively. A total of 162.24 Gb nucleotides containing 1,081,481,144 clean reads were generated, with an average GC content of 55.66% and Q30 > 94% (Table S1). After denovo transcriptome assembly, a total of 66,152 unigenes were obtained to differentially expressed in 12 samples, using a fold change ≥1 and stringent FDR value < 0.05 were as thresholds (Table S2). High repeatability of sequencing data from 12 samples was proved by Pair-wise Pearson’s correlation coefficients (Fig. 2A). 23,070 DEGs, 9060 DEGs, 10,684 DEGs, and 23,338 DEGs were identified at CK1_vs_CK2, CK1_vs_MU1, CK2_vs_MU2 and MU1_vs_MU2, respectively (Fig. 2E). Up-regulated genes presented at CK1_vs_MU1(5556) and CK2_vs_MU2(6185) were significantly more than down-regulated genes (3504, 4499) (Fig. 2D). In contrast, the number of down-regulated genes at CK1_vs_CK2(11,572) and MU1_vs_MU2(12,611) were much more than those up-regulated (11,498, 10,727) ((Fig. 2C, Table S3). 895 specific genes were collectively found in DEGs of CK1_vs_CK2, CK1_vs_MU1, CK2_vs_MU2 and MU1_vs_MU2 (Fig. 2F).

Fig. 2
figure 2

Correlation and PCA analysis of all transcripts and identification of DEGs. A Correlation analysis of all samples. B PCA plot of transcriptome results. C Venn diagram of DEGs before and after cold stress in the same varieties; D, Venn diagram of DEGs between varieties before and after cold stress. Above and below the horizontal line are up-regulated and down-regulated genes, respectively; E Numbers of DEGs between different comparisons. F Venn diagram of DEGs among 4 comparison groups (FDR < 0.05 and FC ≥ 2). The numbers in parentheses showed percentages concerning the total genes. CK1, before treatment Jimai 325; CK2, cold treatment Jimai 325; MU1, before treatment MU-134; MU2, cold treatment MU-134

GO and KEGG analysis of the DEGs

To analyze the functions of the DEGs responding to cold between MU and CK, GO enrichment analysis was performed using the Goseq method in Blast2GO. Based on the GO annotations, 14,466 unigenes between CK1_vs_MU1, CK2_vs_MU2 were classified into 50 functional subgroups, including 38 in the biological process, 10 in molecular function and 2 in cellular component (Fig. 3A). The data showed that the major categories of biological processes are antibiotic degradation processes, hydrogen peroxide degradation processes, and flavonoid biosynthesis processes. Concerning the molecular function of the GO term, most DEGs were highly concerned with glucosidase activity, glutathione transferase activity, and beta-glucosidase activity. The Extrinsic component of plasma membrane and MCM complex were the significant categories annotated for the cellular component. The classification of GO categories was listed in Table S4.

Fig. 3
figure 3

GO (A) and KEGG (B) functional classifications of the annotated unigenes in wheat. The unigenes were distributed into three GO categories: biological process (a), cellular component (b), and molecular function (c). The unigenes were divided into five KEGG groups: cellular processes (a), environmental information processing (b), genetic information processing (c), metabolism (d) and organismal systems (e)

Further, 132 and 133 pathways of the DEGs linked with the two different treatment stages between MU and CK were exposed by KEGG pathway analysis, respectively. The summary of pathways confirmed in our research was listed in Table S5. The top 31 common KEGG pathways distinguished in the MU compared with CK were grouped into five categories, including genetic information processing, cellular processes, metabolism, environmental information processing, and organismal systems (Fig. 3B). The results revealed that the biosynthesis of secondary metabolites pathway was highly enriched, followed by the plant-metabolism interaction, plant hormone signal transduction, the MAPK signaling pathway - plant, Phenylpropanoid biosynthesis and Starch and sucrose metabolism. Therefore, the DEGs annotated as “biosynthesis of secondary metabolites pathway,” “plant-metabolism interaction” and “MAPK signaling pathway - plant” were further investigated in detail.

Transcription factors (TFs) are essential regulators in plants resistance to cold. A total of 183 significant differences TFs (FDR < 0.05 and FC ≥ 2 or FC < 0.05) were identified regulated wheat frost tolerance between MU and CK and divided into 49 categories, including 137 up TFs and 46 down TFs (Fig. 4A). ERF, B3, bHLH, C2H2, MYB, SET contain more transcription factors in different processing stages of MU and CK (Table S6). The expression levels of transcription factors related to plant frost resistance were analyzed, and which results showed the differential gene expression was significant in bHLH, MADS, MYB, and WRKY, although not always in the same direction. Concerning the general trends of TF gene expression among DEGs, bHLH were priority upregulated and WYB preferentially downregulated between MU and CK (Fig. 4B).

Fig. 4
figure 4

Transcription factor classification and expression profile. A Classification of transcription factors. B The expression profile of frost-resistance-related transcription factors. The numbers in the pie represent quantities. Expression was normalized; green represents low expression; red represents high expression

Untargeted metabolite profiles of wheat responding to cold stress

To comprehensively analyze the response of wheat to low-temperature stress, we simultaneously conducted nontargeted metabolome and transcriptome analysis on the same samples. Totally 3621 peaks were obtained in UPLC-MS/MS analysis and 650 metabolites were identified using the MS2 database (Table S7). In order to determine the different metabolites of CK and MU-134 under low-temperature stress, the metabolic profiles of 12 wheat samples were analyzed by principal component analysis (PCA). The results showed obvious separation between the two samples under control and cold treatment. The principal component analysis of metabolome data was stably consistent with the transcriptomic data PCA (Fig. 5A). The variable importance for projection (VIP) value > 1.5 and FC ≥ 2 or FC < 0.05 were used to define differentially expressed metabolites. 88 metabolites were differentially accumulated compared with MU before cold treatment, with 25 up-regulated and 63 down-regulated. After cold treatment, 132 differentially accumulated metabolites were detected compared with MU, with 56 up-regulated compared with MU and 76 down-regulated, respectively (Fig. 5B). Compared with MU, 54 differentially accumulated metabolites were commonly detected in CK before low-temperature stress and after. 16 specific metabolites were collectively found in differentially accumulated metabolites of CK1_vs_CK2, CK1_vs_MU1, CK2_vs_MU2 and MU1_vs_MU2 (Fig. 5C, Table S8).

Fig. 5
figure 5

Differentially accumulated metabolites of wheat in response to cold stress. A PCA plot of metabolomic results. B Numbers of differentially accumulated metabolites under different treatments. C Venn diagram of differentially accumulated metabolites (FDR < 0.05 and FC ≥ 2). D, K-means clustering analysis of the differentially accumulated metabolites into nine clusters according to their expression profile. The cluster names and the number of metabolites for each cluster are indicated. CK1, before treatment Jimai 325; CK2, cold treatment Jimai 325; MU1, before treatment MU-134; MU2, cold treatment MU-134

To identify metabolites with the accordance expression mode in clusters, we employed the K-means clustering algorithm, which groups the metabolites according to the similarity of the metabolome spectrum. A total of 9 clusters were confirmed (Fig. 5D), which can be divided into 5 categories: up-regulated in both CK1_vs_CK2 and MU1_vs_MU2 (class 2, 3, and 5), down-regulated in both CK1_vs_CK2 and MU1_vs_MU2 (class 6, 7, and 8), no change in CK1_vs_CK2 but upregulated in MU1_vs_MU2 (class 4), no change in CK1_vs_CK2 but down-regulated in MU1_vs_MU2 (class 9) and upregulated in CK1_vs_CK2 but no change in MU1_vs_MU2 (class 1). Among the three clusters with an overall decreasing trend of metabolite accumulation, cluster 6 metabolite was enhanced assembly in CK1 and CK2 but descended in MU1 and MU2 (Fig. 5D).

In order to identify the different metabolites between MU and CK in response to cold stress, all metabolite profiles were performed by hierarchical cluster analysis based on Pearson correlation. Within the control or cold treatment, the two samplings (MU and CK) could be detected a distinct separation. In addition, most of the differential metabolites were found to be flavonoids, termites and acids lipids (Fig. 6A). KEGG analysis of all differentially accumulated metabolites revealed the top 20 statistics of KEGG pathway enrichment. Metabolites participating in alpha-Linolenic acid metabolism, flavone and flavanol biosynthesis, Fructose and mannose metabolism, Plant hormone signal transduction and tryptophan metabolism were predominantly enriched in CK1_vs_MU1 (Fig. 6B). “flavone and flavanol biosynthesis” contained the greatest number of differentially accumulated metabolites, followed by “alpha-Linolenic acid metabolism.” Metabolites participating in flavone and flavanol biosynthesis, flavonoid biosynthesis, Benzoxazinoid biosynthesis, isoflavonoid biosynthesis and indole alkaloid biosynthesis were predominantly enriched in CK2_vs_MU2 (Fig. 6C). “flavone and flavanol biosynthesis” contained the greatest number of differentially accumulated metabolites, followed by “flavonoid biosynthesis. “.

Fig. 6
figure 6

Hierarchical cluster and KEGG enrichment analysis for all metabolome samples of frost resistance traits in wheat. A Hierarchical cluster analysis of all samples metabolite content. Blue represents low expression and red represents high expression. B KEGG pathway enrichment analysis based on the differentially accumulated metabolites in MU1_vs_CK1(B) and CK2_vs_MU2 (C). The top 20 statistics of KEGG pathway enrichment were shown. The dot color represents the size of the P-value. The smaller the q value, the closer the color is to red. The dot size represents the number of different genes contained in each pathway

We also analyzed the critical metabolites of CK1_vs_MU1 and CK2_vs_MU2, and found a total of 43 metabolites with significant differences, which could be divided into six categories, including 20 up-regulated metabolites and 23 down-regulated metabolites. Class alkaloids metabolites were the most abundant, followed by class flavonoids. p-Coumaroylputrescine has the highest difference multiple as high as 95.52(Fig. 7). Table 1 displayed a list of the 10 metabolites with the highest differential accumulation. The top upregulated metabolites were primarily alkaloids, amino acids and derivativ. Compared to CK2, p-Coumaroylputrescine showed the greatest fold change in MU2, with a log2 fold change (FC) (MU2/CK2) value of 6.58. Chlorogenic acid, the most downregulated metabolites, had the lowest log2 fold change (FC) (MU2/CK2) value of − 3.01 (Table 1).

Fig. 7
figure 7

The Heatmap of key metabolites in response to the cold treatment in MU and CK. Gray represents low fold change and red represents high fold change

Table 1 List of top 10 up-regulation and down-regulation metabolites in MU vs CK

Integrated analyses of transcriptomic and metabolomic changes involved in vital biological pathways

The metabolomics and transcriptome data were combined to further reveal the response of CK and MU to cold stress. We integrated analyzed all the data before and after cold stress treatment, focusing on analyzing related genes and metabolites between CK and MU after cold stress treatment. Comprehensive analysis of CK2 and MU2 showed that 12,511 genes were positively associated with 138 metabolites, including 9764 DEGs and 132 metabolites with different accumulations. The different accumulations metabolites and DEGs were enriched in the same KEGG pathways, including galactose metabolism (Fig. 8A). The DEGs were chiefly related to biotin metabolism, phenylalanine, tyrosine and tryptophan biosynthesis, glycolysis/gluconeogenesis, and starch and sucrose metabolism, suggesting that wheat may alleviate cold stress by a comprehensive series of energy metabolism and amino acid biosynthesis mechanisms. 290 metabolites and 27,308 genes in 6 clusters had consistent expression patterns by Kmeans integrated analysis of transcriptomics and metabolomics, indicating that gene expression was significantly correlated with metabolite accumulation (Fig. 8B).

Fig. 8
figure 8

Kegg functional classifications (A) and Kmeans (B) of the differentially expressed related genes and metabolites in response to the cold treatment in MU2 and CK2

Gene-metabolite interaction networks can help to understand functional relationships and identify novel regulatory elements. Spearman correlation method was used to calculate the correlation coefficient between DEGs and metabolites with significant differences related to the flavonol biosynthesis pathway, such as “phenylpropanoid biosynthesis,” “flavonoid biosynthesis” and “flavone and flavonol biosynthesis.” On the other hand, correlations were found for pathways related to sucrose metabolism and some amino acid biosynthesis (arginine and proline metabolism). 316 DEGs and 12 DAMs involved in flavonol biosynthesis pathway were executed Pearson correlation analysis. The result showed that 76 DEGs had strong positive correlation coefficient values (R2 > 0.9 and P-value < 0.05) with 11 metabolites, and 68 DEGs had strong negative correlation coefficient values (R2 < − 0.9 and P-value < 0.05) with 10 metabolites. For sucrose metabolism and some amino acid biosynthesis, 57 DEGs and 4 DAMs were carried out in Pearson correlation analysis. The result showed that 5 DEGs had strong positive correlation coefficient values (R2 > 0.9 and P-value < 0.05) with 3 metabolites, and 16 DEGs had strong negative correlation coefficient values (R2 < − 0.9 and P-value < 0.05) with 4 metabolites (Fig. 9). Further analysis revealed that the abundance of most metabolites was significantly positively correlated with the expression of four key genes, including coniferyl alcohol acyltransferase (CFAT), UDP-glycosyltransferase (UDPG), and tryptamine benzoyl transferase (TBT). Two genes negatively correlated with the abundance of most metabolites were flavonoid 3′ -monooxygenase (CYP75B4) and Eceriferum26(CER26).

Fig. 9
figure 9

Correlation network analysis of differential metabolites and differential expression genes related to frost tolerance in wheat. The red lines indicate positive correlation; the green line indicates negative correlation, the blue line indicates the correlation is not significant. Metabolites were marked red, structural genes were green. The size of the red circle represents the number of genes associated with the metabolite. The thickness of the ring frame of metabolite group was adjusted to indicate the differential multiple of metabolites. The line thickness between nodes represents the degree of correlation between two nodes

Flavonol biosynthesis pathway

The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathways were compared and analyzed between transcriptome and metabolome of wheat for frost resistance. Multiple pathways were enriched between CK and MU, including phenylpropanoid biosynthesis, flavonoid biosynthesis, flavone and flavonol biosynthesis, starch and sucrose metabolism, and arginine and proline metabolism.

The flavonol biosynthesis pathway includes three pathways: “phenylpropanoid biosynthesis,” “flavonoid biosynthesis”, and “flavone and flavonol biosynthesis” [51]. In order to further elucidate the specific role and clarify the regulatory mechanism of frost resistance at the seedling stage of wheat, we have combined analysed significantly different metabolites and significantly differentially expressed genes (DEGs) in the flavonol biosynthesis pathway. The results showed that 13 significantly difference metabolites and 14 key enzymes encoded by 83 markedly DEGs were mapped onto the phenylpropanoid/flavonoids metabolic pathway (ko00940, ko00941 and ko00944), as shown in Fig. 10. The expression of 59 genes was strongly up-regulated in 83 genes, including all genes that encoded CHI (chalcone isomerase), CHS (chalcone synthase), CYP73 (trans-cinnamate 4-monooxygenase), 4CL (4-coumarate--CoA ligase), PAL (phenylalanine ammonia-lyase) and CSE (caffeoyl shikimate esterase) (Table S9). In the phenylpropanoid pathway (ko00940), the DEGs PAL (TraesCS1B02G048100, TraesCS1B02G048200, TraesCS1B02G048300, TraesCS1B02G048400, TraesCS1B02G048500, TraesCS1B02G122800, TraesCS2A02G196400), 4CL (TraesCS2D02G150400), and CSE (TraesCS3D02G263200, TraesCS4D02G133100, TraesCS4A02G178300, TraesCS1B02G243700) were up-regulated, and the corresponding metabolites of cinnamic acid and caffeic acid ignificantly increased. The up-regulation or down-regulation of 22 DEGs encoding HCT (shikimate O-hydroxycinnamoyl transferase) and CYP98 (TraesCS7B02G241800) reduced the metabolites of p-Coumaroyl quinic acid and Chlorogenic acid metabolites. P-Coumaroyl coenzyme A (CoA) is converted to naringenin according to serial enzymes (CHS, CHI) in the “phenylpropanoid biosynthesis” pathway. Since then, afzelechin, prunin, luteolin, isovitexin, vitexin 2″-O-beta-L-rhamnoside, apigenin, vitexin, vitexin 2″-O-beta-D-glucoside, and trifolin metabolites were downregulated in the latter two pathways (ko00941, ko00944) through a series of enzymes, such as flavonoid 10′-monooxygenase (CYP75B1), flavone synthase II amino acid enzyme (CYP93B2_16), and flavonol-3-O-L-rhamnoside-7-O-glucosyltransferase (UGT73C6) (Fig. 10, Table S9).

Fig. 10
figure 10

Transcript and metabolic profiling of significantly differential metabolites and DEGs in flavonol biosynthesis pathway in wheat under cold treatment. A The map of integrative analysis in flavonol biosynthesis metabolites and their biosynthesis-related key enzymes. Red font: significantly up-regulated differential metabolites or DEGs encoding the corresponding enzyme; green font: significantly down-regulated differential metabolites or DEGs encoding the corresponding enzyme; Blue font: significantly up/down-regulated DEGs encoding the corresponding enzyme; Black font: the enzyme encoding genes and differential metabolites detected but with no significant difference. Black arrow lines: ko00940; Blue arrow lines: ko00941; Red arrow lines: ko00944. DEGs and DAMs were investigated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) to map to the possible KEGG pathway maps for the biological interpretation of systemic functions ( B Heatmap of differential metabolites accumulation flavonol biosynthesis pathway. The color bar presented at the top right represents the level of metabolites, where red indicates the metabolites with a higher level, and blue indicates metabolites with a lower level. C Heatmap of DEGs enriched flavonol biosynthesis pathway. Gene expression was scaled using normalized FPKM for the mean value of three biological replicates. A color bar was presented at the top right, and the colors from blue to red indicate low to high expression. 4CL, 4-coumarate--CoA ligase; ANS, anthocyanidin synthase; CHI, chalcone isomerase; CHS, chalcone synthase; CYP73, trans-cinnamate 4-monooxygenase; CYP75B1, flavonoid 3′-monooxygenase; CYP93B2_16, flavone synthase II; CYP98, 5-O-(4-coumaroyl)-D-quinate 3′-monooxygenase; FG2, Kaempferol-3-O-rutinoside; FLS, flavonol synthase; HCT, shikimate O-hydroxycinnamoyltransferase; PAL, phenylalanine ammonia-lyase; UGT73C6, flavonol-3-O-L-rhamnoside-7-O-glucosyltransferase; CSE, caffeoyl shikimate esterase

Sucrose and amino acid biosynthesis pathway

To better analyze the cooperation between differentially expressed genes (DEGs) and metabolites, we also chose some DEGs and metabolic associated with sugar and amino acid pathways to draw the sketch map directing the difference between CK and MU after cold treatment (Fig. 11). The results showed that 12 significantly different metabolites and 19 key enzymes encoded by 43 markedly DEGs were mapped onto the sucrose and amino acid biosynthesis pathway (ko00051, ko00052, ko00270, ko00330 and ko00380), as shown in Fig. 11. In the 43 DEGs, the expression of 21 genes was strongly up-regulated, including all genes that encoded alpha-galactosidase (GLA), sucrose-phosphate synthase (SPS), sucrose-6-phosphatase (SPP), glucose-6-phosphate isomerase (GPI), L-asparaginase (ansA), glutamate 5-kinase/glutamate-5-semialdehyde dehydrogenase (proA), 4-aminobutyrate-pyruvate transaminase (POP2), arginine decarboxylase (ADC). The expression of 20 genes was strongly down-regulated, including all genes that encoded inositol 3-alpha-galactosyltransferase (GOLS), endoglucanase (cel), xylulokinase (xylB), D-ribulokinase (rbtK), pyruvate kinase (pyk) and cysteine synthase (cysE) (Table S10).

Fig. 11
figure 11

Transcript and metabolic profiling of significantly differential metabolites and DEGs through sugar and amino acids metabolic pathways in wheat under cold treatment. A The map of integrative analysis in sugar and amino acids metabolites and their biosynthesis-related key enzymes. Red font: significantly up-regulated differential metabolites or DEGs encoding the corresponding enzyme; green font: significantly down-regulated differential metabolites or DEGs encoding the corresponding enzyme; Blue font: significantly up/down-regulated DEGs encoding the corresponding enzyme; Black font: the enzyme encoding genes and differential metabolites detected but with no significant difference. DEGs and DAMs were investigated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) to map to the possible KEGG pathway maps for the biological interpretation of systemic functions ( B Heatmap of differential metabolites accumulation in sugar and amino acids biosynthesis pathways. The color bar presented at the top right represents the level of metabolites, where red indicates the metabolites with a higher level, and blue indicates metabolites with a lower level. C Heatmap of DEGs enriched in sugar and amino acids biosynthesis pathways. Gene expression was scaled using normalized FPKM for the mean value of three biological replicates. A color bar was presented at the top right, and the colors from blue to red indicate low to high expression. SPP, sucrose-6-phosphatase; ADC, arginine decarboxylase; rbtK, D-ribulokinase; ansA, L-asparaginase; GOLS, inositol 3-alpha-galactosyltransferase; ahcY, adenosylhomocysteinase; xylB, xylulokinase; HK, hexokinase; pyk, pyruvate kinase; proA, glutamate 5-kinase; SPS, sucrose-phosphate synthase; cysE, serine O-acetyltransferase; GLUL, glutamine synthetase; GLA, alpha-galactosidase; GPI, glucose-6-phosphate isomerase; aroK, shikimate kinase; cel, endoglucanase; RFS, raffinose synthase; POP2, 4-aminobutyrate---pyruvate transaminase; cysK, cysteine synthase

In the galactose metabolism pathway (ko00052), the DEGs RFS (raffinose synthase, TraesCS7A02G236500, TraesCS5D02G271300) were up-regulated, and the level of corresponding metabolites of raffinose markedly increased. Meanwhile, the DEGs of SPS (TraesCS3B02G461800, TraesCS3D02G013000), SPP (TraesCS1A02G088600, TraesCS1B02G107600), GLA (TraesCS5A02G173200, TraesCS1B02G182000) and GPI (TraesCS5A02G247200) were significantly up-regulated resulted in a increase in D-glucose and D-glucose-6p metabolites. Whereas, the ribitol metabolite was observed significantly reduce in the CK vs. MU comparison due to rbtk (TraesCS1B02G077300, TraesCS1B02G142000 and TraesCS1B02G146300) down-regulated. In the cysteine and methionine metabolism pathway (ko00270), the down-regulated of DEGs cysE (TraesCS4A02G042700) and ahcY (adenosylhomocysteinase, TraesCS2B02G517500, TraesCS2B02G538600, TraesCS2B02G538700, TraesCS2B02G538800) caused the homocysteine significantly decline. Strikingly, the homocysteine level decreased more than 7685-fold in the CK2 vs. MU2 comparison, whereas only a 3-fold decrease was observed in the CK1 vs. MU1 comparison. In the arginine and proline metabolism pathway (ko00330), the up-regulation of DEGs encoding proA (TraesCS3B02G395900, TraesCS3A02G363700 and TraesCS3D02G357200) and POP2 (TraesCS7A02G334200) resulted in a increase in proline metabolites. Since then, p-Coumaroylputrescine, Feruloylagmatine and p-Coumaroylagmatine metabolites were up-regulated through a series of enzymes, such as ansA (TraesCS2B02G016000) and ADC (TraesCS1B02G018200). Significantly, the level of p-Coumaroylagmatine decreased more than 30-fold in the CK2 vs. MU2 comparison (Fig. 11, Table S10).

Validation of differentially expressed proteins at the transcription level

To verify the reliability of expression levels in RNA-seq data, quantitative real-time PCR (qRT-PCR) was applied to validate the twelve differentially expressed frost-resistance related genes (Fig. 12A). Correlation analysis further exhibited a significant correlation (correlation coefficient of 0.94) between RNA-Seq and qRT-PCR (Fig. 12B). The qRT-PCR results corroborated that the expression trends of twelve DEGs were highly consistent with the transcriptome data and supported the credibility of RNA-seq data. UGT73C6 (TraesCS3B02G144800), UGTs (TraesCS3A02G124800) and GSTU6 (TraesCS4A02G455000) showed downregulation trend after receiving cold treatment, and the expression level of MU was significantly lower than CK. However, HCT (TraesCSU02G166000), CYP75B1 (TraesCS7B02G310900) and ADC (TraesCS1B02G018200) were verified up-regulated after cold treatment, and the expression level of MU was significantly higher than that of CK. After cold treatment, POD2 (TraesCS1B02G355200), Tacr7 (TraesCS2A02G460800, TraesCS2B02G482600) and GSTU6 (TraesCS7A02G034500) showed a sharp up-regulation trend of CK, but the expression levels of MU were very low before and after cold treatment. The relative expression of POD72 (TraesCS3D02G185000, TraesCS3A02G180200) in MU1 was found significantly higher than that in CK1 and CK2. After cold treatment, the expression of these two DEGs was up-regulated in CK and decreased extremely to close to zero in MU. We believe that candidate genes POD, Tacr7, UGTs and GSTU6 were related to frost resistance of wheat and could be utilized as transgenic or molecular marker-assisted selection to promote frost resistance of wheat. We will proceed genetic transformation of candidate genes in wheat to verify further.

Fig. 12
figure 12

The relative expression levels of twelve selected DEGs were compared by RNA-seq and qRT-PCR. A Expression patterns of the 12 genes involved in the Cold stress of wheat. B correlation analysis based on RNA-Seq and qRT-PCR data. The line chart shows the gene expression level from the transcriptome (FPKM); The qRT-PCR expression levels were calculated as a ratio relative to the level of expression of CK1, which was set as 1. Bars indicate means ± standard deviations (SDs) of at least three independent biological replicates


Large numbers of DEGs were found in leaves between the wild-type and mutants before and after cold treatment. GO enrichment analysis shows that these DEGs were related to large-scale biological processes and molecular functions, such as hydrogen peroxide catabolic, extrinsic component of plasma membrane, and flavonoid biosynthetic process, which were expected, but also related to other pathways due to interaction between elements. Compared with MU1, the number of DEGs in CK1 was 9060 (5556 up-regulated and 3504 down-regulated). In contrast, cold treatment had the greatest effect on wheat transcriptome, with 10,684 DEGs changed (6185 up-regulated and 4499 down-regulated; Fig. 3E). Many genes were found in catabolic, extrinsic component of plasma membrane and flavonoid biosynthetic process before and after cold treatment. These results indicated that the frost resistance of wild-type and mutant wheat in this experiment was extremely related to plasma membrane, oxidative decomposition and flavonoids, which was consistent with previous research results [52,53,54,55], but also involved in sugar and energy metabolism [56]. In addition, these DEGs were mapped to KEGG pathways. The four largest KEGG categories were “Biosynthesis of Secondary metabolites,” “MAPK signaling Pathway-Plant,” “plant hormone signal” and “Phenylpropanoid biosynthesis” (Fig. 3B). These results were also consistent with previous studies on biological processes [57, 58]. However, there were also many genes enriched in the “plant-interaction” pathway, because low temperature significantly changed the membrane permeability of plants leading to a decrease in plant resistance to pathogens [59].

Transcription factors (TFs) also play an important role in low-temperature stress response of plants, such as DREB [60], WARY [61], NAC [62], ICE [63] and MAD [64]. Several transcription factors induced by low temperature stress were also identified, such as TraesCS2D02G011700 (WRKY), TraesCS3B02G416400 (MYB) and TraesCS7B02G152800 (bHLH). These transcription factors induced by low-temperature stress may play a critical role in the frost resistance mechanism of wheat. Several key transcription factors directly regulate frost resistance of wheat in different ways. However, the number of transcription factors regulating plant frost tolerance in specific plants is minimal, and the relationship between transcription factors has not been thoroughly studied. In this study, transcriptome sequencing data provided the basis for distinguishing regulatory networks and novel regulatory factors involved in frost resistance of wheat (Fig. 4). As the second-largest transcription factor in plants, bHLH transcription factor plays an essential role in plant growth and development, physiological metabolism and stress response. It mainly regulates the expression of related genes through the interaction between specific amino acids and target genes. NtbHLH123 transcription factor in tobacco enhances tolerance to cold stress through regulatory pathways and reactive oxygen homeostasis [65]. MYB and bHLH families are believed to regulate plant secondary metabolism and stress responses according to combinatorial interactions of MYB and bHLH proteins or multi-protein complexes formed of different subgroups of bHLH proteins [66, 67]. In this study, TraesCS7B02G152800 (bHLH) and TraesCS3B02G416400 (MYB) were highly induced by cold stress. This discovery is a living example of cooperation between different transcription factors to regulate plant stress tolerance. The expression of TraesCS7B02G152800 was significantly up-regulated under cold stress, and the expression of TraesCS3B02G416400 was significantly down-regulated under cold stress, which indicated that the two genes were related to the difference of frost tolerance between two wheat varieties.

Metabolomics can define and quantify the small molecular weight metabolites in biological cells, and has been widely used to study the accumulation patterns of metabolites. In addition, the changes of genes and metabolites in different tissues were revealed by using this technique. Metabolites are intermediate and final products that regulate plant growth and development. The total number of metabolites in plants is as high as 200,000, making them an ideal target for studying biosynthesis regulation [68]. In recent years, the genetic mechanism of plant metabolite changes has been studied [69, 70]. We detected 650 metabolites that could be annotated in MS2 database using the extensive targeted metabolomics method, and the number of detected metabolites was higher than that of previous studies on frost resistance of wheat [71]. The four largest KEGG categories associated with the expression of differential metabolites were “flavone and flavonol biosynthesis,” “flavonoid biosynthesis,” “isoflavonate biosynthesis” and “plant hormone signal transduction.” These results are consistent with previous studies on plant frost resistance biological processes [57]. Ten metabolites were detected with significant difference in frost resistance at seedling stage, of which 8 were increased and 2 were decreased. They fall into three main categories:alkaloids, mino acids and derivativ, and phenolic acids. Notably, the levels of two alkaloids, p-coumaroyl putrescine and p-coumaroyl cadaverine, were significantly increased by 95 and 59 fold, respectively. Yu et al. [72] found that p-coumaroyl glucoside is one of the primary differential accumulative metabolites of maize frost tolerance. Therefore, we believe that the accumulation of p-coumaroyl may be a metabolite formed during plant resistance to low-temperature stress.

Comparative analysis of metabolomics and transcriptome can facilitate interpreting biological processes and mechanisms involved in regulatory responses. In the process of low-temperature stress, plants will mobilize a significant number of factors for defense, especially metabolites. These reactions first require ample energy supply. Therefore, when low-temperature stress occurs, the energy metabolism of plants increases significantly [73]. Most energy of plants comes from respiration, and the substrate of respiration is glucose. The regulation of enzymes involved in sucrose synthesis and decomposition is more complex at low temperatures. Sucrose phosphate synthase (SPS) and sucrose-6-phosphatase (SPP) catalyze sucrose synthesis and their activity reflects the ability of sucrose biosynthesis pathway. The research by Guy et al. [74] showed low-temperature stress can improve the activity of SPS in spinach leaves. In this study, both SPS and SPP were up-regulated, which may be caused by cold stress in order to compensate for the effects of low temperature. In our study, both wild-type and mutant glucose content was higher after cold treatment than before. Interestingly, the glucose content of the wild-type was 2.4 times higher than that of the frost-sensitive mutant, and the corresponding phosphate isomerase (GPI) gene was also up-regulated 4.1 times. One of the reasons may be that plants regulate osmotic potential and enhance cell water retention through carbohydrate accumulation. Second, other protective substances and energy are generated through glucose metabolism. Thus, lines with poor frost resistance may exhibit lower sugar accumulation. Raffinose family (RFOs), as the second group of sugars in plants, plays an essential role in plant growth and development and abiotic stress resistance. Studies show that sugar beet produce an oligosaccharide raffinose to protect tissue from freezing damage in low-temperature environment [75]. Overexpression of AnGolS1 in tomato can increase the contents of inositol galactoside, sucrose and glucose, and improve frost resistance of tomato [76]. In this study, up-regulated expression of raffinose synthase (RFS) resulted in increased accumulation of raffinose metabolites.

Amino acids are essential for protein synthesis, and previous studies have shown that amino acid metabolism plays an integral role in enhancing plant abiotic stress tolerance. Proline is an osmotic regulator involved in plant response to abiotic stress. It can also induce the expression of cold stress response and enhance the cold tolerance of wheat [71]. In our study, under cold treatment, the expression levels of ProA and POP2 genes were up-regulated, leading to significantly higher proline content of CK and MU than before treatment, which was also consistent with the results of physiological indicators determined by us. AnsA and ADC significantly up-regulated p-Coumaroyl putrescine, Feruloyl agmatine, and p-Coumaroyl agmatine metabolite accumulation. The increase of proline, P-Coumaroyl putrescine, Feruloyl agmatine, and P-Coumaroyl agmatine reduced cold stress damage on wheat. These results strongly suggest that the accumulation of osmoprotectants (soluble carbohydrates, proline and amino acids) contributes to the improvement of cold resistance of plants [77].

Environmental stress can induce biosynthesis of secondary metabolites, such as phenylpropanoid, flavonoids, and flavone and flavonol biosynthesis [78]. The phenylpropanoid biosynthetic pathway is one of the most critical metabolite pathways in plants, producing a large number of secondary metabolites, such as lignin and flavonoids [79]. Lignification of plant tissues is a stress perception and signal transduction mechanism that can prevent freezing injury and cell collapse [78]. Recent studies have shown that the expression of structural genes, such as Chalcone synthase (CHS) and chalcone isomerase (CHI),were induced by cold stress in phenylpropanoid pathway. Therefore, the accumulation of flavonoids and lignin can adapt plants to low-temperature environment [80]. This study found that some frost-response genes involved in the phenylpropanoid metabolism biosynthesis were expressed higher after low-temperature treatment than before. For example, PAL is a rate-limiting enzyme that catalyzes the deamination of phenylalanine to cinnamic acid and is the entry point for synthesizing all phenylalanine. Compared with before the treatment, the transcription level of 7 PALs was higher after treatment, which allowed more carbon flux to enter the phenylpropanoid pathway, leading to increased accumulation and better frost tolerance. As a branch product of phenylpropanoid metabolism, flavonoids are essential and affect basic metabolism and stress resistance in plants. When plants suffer adversity, the excessive accumulation of oxygen free radicals caused by metabolic disorders in the body is removed by a large number of flavonoids synthesized by stress, and the antioxidant activity of plants is enhanced [81]. It has been confirmed that flavonoid accumulation in Arabidopsis and Spinach was positively correlated with plant stress resistance [82, 83]. The expression levels of genes related to flavonoid synthesis were significantly up-regulated in the samples treated with low temperature, confirming that flavonoid synthesis enhanced the low-temperature tolerance of wheat leaves and reduced the damage caused by low temperature. However, flavonoid metabolite content was no change in this study, and only the decrease of metabolites of Naringenin, Afzelechin and Prunin in the flavonoid pathway was detected. In the present study, we also found that Vitexin 2″-O-beta-L-rhamnoside, Vitexin,Vitexin 2″-O-beta-D-glucoside after cold treatment metabolite accumulation were decreased. Kirakosyan [84] found that the metabolites content of Vitexin 2 “-O-rhamnoside in hawthorn increased 1.68 times after cold treatment, which may be caused by different crops and materials.

The frost resistance related genes found in this study provide reference for regulating frost resistance of wheat and have potential agricultural application value. However, the functions of most of the DEGs or DAMs found in this study remain largely unknown, and we will continue to study the functions of these genes and metabolites.


In summary, the differences between frost-resistant material CK and frost-sensitive mutant MU were analyzed by physiological, metabolomics and transcriptomics. The underlying mechanisms between gene expression and metabolite biosynthesis were discussed and suggested that the biosynthesis of flavonol, sucrose and amino acid may play an essential role in wheat resistance to cold. This research provides valuable molecular information for frost resistance mechanisms in wheat, and will significantly promote the selection of frost-resistant varieties by biomarker-assisted selection. The critical regulatory genes of flavonol, sucrose and amino acid biosynthesis and their roles in frost resistance need to further explore and verify.

Materials and methods

Plant growth conditions and stress treatments

The MU-134 mutant utilized in this study was isolated from EMS (ethyl methanesulfonate) mutated population (3219 lines) originated by Dr. Liangjie Lv from a high yield, drought-resistant cultivar Jimai 325 cultivated by the Institute of Cereal and Oil Crops, Hebei Academy of Agriculture and Forestry Sciences (HAAFS). All plants were planted within the experimental plot at HAAFS in Shijiazhuang, China. Jimai 325 (CK) and MU-134 (M7 generation) individuals were planted in parallel and collected.

Seeds, Jimai 325 and mutant MU-134, were sterilized in 10% (v/v) H2O2 for 25 min and rinsed with distilled water, then germinated on paper towels soaked with distilled water for 2 d. When the roots were about 1.5 cm long, uniform seedlings were transplanted into a growth chamber with a temperature of 22 °C/18 °C, a 14/10 hour light/dark cycle (light for 14 hours from 6 a.m. to 8 p.m.) and 65% humidity. 28-days-old seedlings were transferred to a pre-cooled chamber with a gradually decreasing temperature (1 °C per hour, from 22° to 6°). After cold acclimation, plants were placed at 6° during the day / -5° during the night for 7 days. Throughout the experiment, Water was normally irrigation with the same soil moisture as the water capacity of field. 60 individuals from 6 pots were selected from each material in each period. 12 individuals were used for sequence and physiological analysis. Finally, the survival rate (measured as the percentage of total seedlings that survived in relation to the total number tested) of 16 plants per replicate were measured to evaluate cold tolerance with three replicates. All experiments were performed simultaneously in the light cycle and repeated at least three times.

Assays of physiological parameters

The total soluble sugar content was calculated according to the standard curve and measured by anthrone colorimetry method with glucose as the standard [85] 0.2 g seedling tissue was ground with 1 mL distilled water, heated for 10 min at 100 °C, cooled, centrifuged at 8000 g for 10 min, and diluted to 10 mL volume with distilled water. The reaction mixture of 0.3 mL includes 0.04 mL extract, 0.04 mL distilled water, 0.02 mL mixed reagent (1 g anthrone and 50 mL ethyl acetate) and 0.2 mL 98% (w / w) H2SO4. The mixture was heated at 95 °C for 10 minutes, and the absorbance was determined at 620 nm using a Scanning spectrophotometer (Agilent8453e, CA, USA). Malondialdehyde (MDA) content was determined by the thiobarbituric acid reactive substances (TBARS) method using the difference of absorbance at 532 nm and 600 nm [86]. According to the method of WST-8 [87], WST-8 and superoxide anion (O2) reaction catalyzed by xanthine oxidase (XO) can produce water-soluble formazan dyes. The activity of superoxide dismutase (SOD) is negatively correlated with the amount of black paint, so the activity of SOD can be calculated by colorimetric analysis of WST-8 products. Peroxidase (POD) enzyme activity can be catalyzed for H2O2 oxidation of specific substrates by the peroxidase activity detection kit (BC0090, Beijing SolarBio Science & Technology Co., Ltd., Beijing, China), with characteristic light absorption detected at 470 nm [88]. Catalase (CAT) was determined by Li et al. [89] method using catalase kit (S0051, Beyotime Biotechnology, Shanghai, China). According to the method of proline (Pro) content detection kit (BC0295, Beijing SolarBio Science & Technology Co., Ltd., Beijing, China), Pro was extracted with sulfosalicylic acid (SA). Approximately 200 mg seedling tissue was homogenized in 2 mL sulphosalicylic acid (3%), followed by centrifugation at 13000 g for 15 min at 4 °C. The extract (0.5 mL) was transferred to a new microcentrifuge tube and mixed with acid ninhydrin and acetic acid. The reaction mixture was boiled in a water bath at 100 °C for 30 min, cooled at 4 °C for 30 min and thoroughly mixed with 1 mL toluene. Finally, the 0.2 mL absorbance was measured at 520 nm using a spectrophotometer (PowerWave XS2) [90].

RNA-seq and bioinformatics analysis

As three biological replicas of cold stress treatments and control, twelve fully expanded leaves of 12 individual plants were analyzed by RNA-Seq. Trizol reagents (Invitrogen, CA, USA) was used to extract total RNA. According to the manufacturer’s instructions (Illumina), 100 ng RNA is used to construct the RNA-seq library and synthesize the first-strand cDNA using SuperScript II reverse transcriptase (Invitrogen, CA, USA). After the second-strand cDNA was synthesized and linked, the enriched and purified cDNA fragment was sequenced on the RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA).

The sequenced raw readings are processed into clean readings by filtering low-quality readings. The clean reads were localized to the reference genome of wheat ( using HISAT2 [91]. The R package, edgeR v3.16, was employed to identify the differentially expressed genes (DEGs) between cold stress and control samples. DEGs of cold-treated plants relative to control samples were detected with a P < 0.05 and |log2 ratio ≥ 2|.

as the threshold. Genes were annotated by Nr, Nt, Pfam, KOG/COG, Swiss-Prot, KO, and GO databases. The central and highly connected metabolites of specific modules were identified using Cytoscape v3.9.0 visualizing [92]. Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of DEGs were generated using R, based on the hypergeometric distribution. Plant transcription factors were predicted using ITAK software ( [93]. Venny diagrams are drawn using Venny 2.1 on the website ( [94]. Raw and processed RNA-Seq data files were deposited in SRA ( under the following accession numbers PRJNA800586(

Metabolomic profiling

Twelve fully unfolded leaves from 12 individual plants of three biological repeats were analyzed by metabonomics profiling under stress treatment and control. Each sample (25 mg) was feed into an EP tube containing an internal standard (L-2- Chlorophenylalanine, 2 μg/mL) of 1000 μ L extract solution (acetonitrile: methanol: water = 2:2:1). Samples were incubated (− 40 °C, 1 h) and centrifuged (10,000 rpm, 15 min) after homogenization and sonication cycle. The supernatant (500 μL) was dried in a vacuum concentrator and ultrasonic degradation 200 μL of 50% acetonitrile. 75 μL centrifuged supernatant was transferred to a fresh glass vial for LC/MS analysis. The quality control (QC) sample was made by mixing the supernatants of all samples in equal quantities. Metabolomic analyses were performed using LC-MS as described by Law et al. [95]. Samples were analyzed on an Agilent 7890B gas chromatography system coupled to an Agilent 5977A MSD system (Agilent Technologies Inc., CA, USA). MS raw data files were peak deconvolution, alignment and integration processed by R package XCMS (version 3.2). The metabolome resulting data were identified in MS2 database and analyzed in MetaboAnalyst 4.0 online [96]. The differential metabolites were selected based on the combination of a statistically significant threshold of variable influence on projection (VIP) values obtained from the OPLS-DA model and p values from the normalized peak areas from different groups, where metabolites with VIP > 1.0 and p < 0.05 were considered as differential metabolites. Hierarchical clustering (Euclidean distance) was performed with MeV4.9 to explore the pattern of metabolite abundances. The ratio of stress-treated plants to control plants was used to calculate the fold-change of each metabolite. As the significance threshold, the variable influence projection (VIP) value of the OPLS-DA model and the p-value of two-tailed Student’s t-test on the normalized peak area were used to select the differential metabolites of different groups, in which the metabolites with VIP > 1.0and p < 0.05were regarded as differential metabolites. To explore the distribution of metabolite abundance, MeV4.9 was carried out to analyze hierarchical clustering (Euclidean distance).

Verification of candidate genes by quantitative real-time PCR (qRT-PCR)

Quantitative real-time PCR (qRT-PCR) was utilized to validate the twelve differentially expressed unigenes transcript expression levels using the same RNA samples used for sequencing. The Specific primers were designed using Primer Premier 6.0 and are listed in Supplementary Table S11. The housekeeping gene was β-actin. According to the user manual, the first cDNA strand fragments were synthesized from total RNA using the PrimeScript™ RT Master Mix kit (Takara, Japan). The qRT-PCR was performed on an ABI7500 real-time fluorescent quantitative PCR instrument (Applied Biosystems, Foster City, USA). Three replicates were performed for each sample. The relative expression levels of the candidate genes were measured using the comparative threshold cycle (2-ΔΔCT) method [41].

Statistical analysis

All data analyses and graphics drawings were accomplished in R 4.1.2. T The statistical significance was calculated by Student’s t-test. The gene expression values of FPKM were scaled to Z-score to draw the heatmap of transcriptome and metabolome. Graphics and Heatmaps were plotted with the ggplot and the heatmap packages of R language respectively. Venn diagrams were drawn with the Interactivenn website ( Each experiment was repeated three times.

Availability of data and materials

The datasets used and analysed during the current study are available in the NCBI Bioproject repository, [PRJNA800586]. The experimental materials, Jimai 325 and MU-134, were bred and mutated by the Institute of Cereal and Oil Crops, Hebei Academy of Agriculture and Forestry Sciences (HAAFS). Seeds were obtained from the Institute of Cereal and Oil Crops of the Hebei Academy of Agriculture and Forestry Sciences (HAAFS). All databases in this study are available to the public.


  1. Liu H, Yang Y, Liu D, Wang X, Zhang L. Transcription factor TabHLH49 positively regulates dehydrin WZY2 gene expression and enhances drought stress tolerance in wheat. BMC Plant Biol. 2020;20(1):259.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  2. Hassan M, Xiang C, Farooq M, Muhammad N, Yan Z, Hui X, et al. Cold stress in wheat: plant acclimation responses and management strategies. Front Plant Sci. 2021;12:676884.

    PubMed  PubMed Central  Article  Google Scholar 

  3. Knight H, Knight MR. Abiotic stress signalling pathways: specificity and cross-talk. Trends Plant Sci. 2001;6(6):262–7.

    CAS  PubMed  Article  Google Scholar 

  4. Jia Y, Ding Y, Shi Y, Zhang X, Gong Z, Yang S. The cbfs triple mutants reveal the essential functions of CBFs in cold acclimation and allow the definition of CBF regulons in Arabidopsis. New Phytol. 2016;212(2):345–53.

    CAS  PubMed  Article  Google Scholar 

  5. Wang X, Wang H, Liu S, Ferjani A, Li J, Yan J, et al. Genetic variation in ZmVPP1 contributes to drought tolerance in maize seedlings. Nat Genet. 2016;48(10):1233–41.

    CAS  PubMed  Article  Google Scholar 

  6. Jiang B, Shi Y, Zhang X, Xin X, Yang S. PIF3 is a negative regulator of the CBF pathway and freezing tolerance in Arabidopsis. Proc Natl Acad Sci U S A. 2017;114(32):E6695.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. Shi Y, Ding Y, Yang S. Molecular regulation of CBF signaling in cold acclimation. Trends Plant Sci. 2018;23(7):623–37.

    CAS  PubMed  Article  Google Scholar 

  8. Bao Y, Yang N, Meng J, Wang D, Cang J. Research on the adaptability of winter wheat Dongnongdongmai 1 (Triticum aestivum L.) to overwintering in alpine regions. Plant Biol. 2021:23:445–55.

  9. Szalai G, Pap M, Janda T. Light-induced frost tolerance differs in winter and spring wheat plants. Plant Physiol. 2009;166(16):1826–31.

    CAS  Article  Google Scholar 

  10. Maruyama K, Takeda M, Kidokoro S, Yamada K, Sakuma Y, Urano K, et al. Metabolic pathways involved in cold acclimation identified by integrated analysis of metabolites and transcripts regulated by DREB1A and DREB2A. Plant Physiol. 2009;150(4):1972–80.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. Dong S, Diane M, Beckles C. Dynamic changes in the starch-sugar interconversion within plant source and sink tissues promote a better abiotic stress response. J Plant Physiol. 2019;234:80-93.

  12. Krasensky J, Jonak C. Drought, salt, and temperature stress-induced metabolic rearrangements and regulatory networks. J Exp Bot. 2012;63(4):1593–608.

    CAS  PubMed  Article  Google Scholar 

  13. Janska AA, Marsík P, Zelenková S, Ovesná J. Cold stress and acclimation-what is important for metabolic adjustment. Plant Biol. 2010;12(3):395–405.

    CAS  PubMed  Article  Google Scholar 

  14. Battaglia M, Olvera-Carrillo Y, Garciarrubio A, Campos F, Covarrubias AA. The enigmatic LEA proteins and other hydrophilins. Plant Physiol. 2008;148(1):6–24.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  15. Tolleter D, Jaquinod M, Mangavel C, Passirani C, Saulnier P, Manon S, et al. Structure and function of a mitochondrial late embryogenesis abundant protein are revealed by desiccation. Plant Cell. 2007;19(5):1580–9.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  16. Koubaa S, Brini F. Functional analysis of a wheat group 3 late embryogenesis abundant protein (TdLEA3) in Arabidopsis thaliana under abiotic and biotic stresses. Plant Physiol Biochem. 2020;156:396–406.

    CAS  PubMed  Article  Google Scholar 

  17. Tran LS, Nishiyama R, Yamaguchi-Shinozaki K, Shinozaki K. Potential utilization of NAC transcription factors to enhance abiotic stress tolerance in plants by biotechnological approach. GM Crops. 2010;1(1):32–9.

    PubMed  Article  Google Scholar 

  18. Wani S, Singh N, Devi T, Haribushan A, Jeberson S. Engineering abiotic stress tolerance in plants: extricating regulatory gene complex; 2013.

    Google Scholar 

  19. Aghaei M, Jaafar M-A, Farhad E, et al. Identification of CBF14 and NAC2 genes in aegilops tauschii associated with resistance to freezing stress. Appl Biochem Biotechnol. 2015;176:1059–70.

    PubMed  Article  CAS  Google Scholar 

  20. Orvar BL, Sangwan V, Omann F, Dhindsa RS. Early steps in cold sensing by plant cells: the role of actin cytoskeleton and membrane fluidity. Plant J. 2010;23(6):785–94.

    Article  Google Scholar 

  21. Sangwan V, Foulds I, Singh J, Dhindsa RS. Cold-activation of Brassica napus BN115 promoter is mediated by structural changes in membranes and cytoskeleton, and requires Ca2+ influx. Plant J. 2010;27(1):1–12.

    Article  Google Scholar 

  22. Balti I, Benny J, Perrone A, Caruso T, Abdallah D, Salhi-Hannachi A, et al. Identification of conserved genes linked to responses to abiotic stresses in leaves among different plant species. Funct Plant Biol. 2020;48(1):54–71.

    PubMed  Article  CAS  Google Scholar 

  23. Le M, Pagter M, Hincha D. Global changes in gene expression, assayed by microarray hybridization and quantitative RT-PCR, during acclimation of three Arabidopsis thaliana accessions to sub-zero temperatures after cold acclimation. Plant Mol Biol. 2015;87(1–2):1–15.

    CAS  PubMed  Article  Google Scholar 

  24. Yokota H, Iehisa JC, Shimosaka E, Takumi S. Line differences in Cor/Lea and fructan biosynthesis-related gene transcript accumulation are related to distinct freezing tolerance levels in synthetic wheat hexaploids. J Plant Physiol. 2015;176:78–88.

    CAS  PubMed  Article  Google Scholar 

  25. Kobayashi F, Takumi S, Handa H. Identification of quantitative trait loci for ABA responsiveness at the seedling stage associated with ABA-regulated gene expression in common wheat. Theor Appl Genet. 2010;121(4):629–41.

    CAS  PubMed  Article  Google Scholar 

  26. Xu Z, Ni Z, Liu L, Nie L, Li L, Chen M, et al. Characterization of the TaAIDFa gene encoding a CRT/DRE-binding factor responsive to drought, high-salt, and cold stress in wheat. Mol Gen Genomics. 2008;280(6):497–508.

    CAS  Article  Google Scholar 

  27. Babben S, Schliephake E, Janitza P, Berner T, Keilwagen J, Koch M, et al. Association genetics studies on frost tolerance in wheat (Triticum aestivum L.) reveal new highly conserved amino acid substitutions in CBF-A3, CBF-A15, VRN3 and PPD1 genes. BMC Genomics. 2018;19(1):409.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  28. Zhu J, Pearce S, Burke A, See DR, Skinner DZ, Dubcovsky J, et al. Copy number and haplotype variation at the VRN-A1 and central FR-A2 loci are associated with frost tolerance in hexaploid wheat. Theor Appl Genet. 2014;127(5):1183–97.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  29. Motomura Y, Kobayashi F, Iehisa JC, Takumi S. A major quantitative trait locus for cold-responsive gene expression is linked to frost-resistance gene Fr-A2 in common wheat. Breed Sci. 2013;63(1):58–67.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  30. Udvardi MK, Kakar K, Wandrey M, Montanari O, Murray J, Andriankaja A, et al. Legume transcription factors: global regulators of plant development and response to the environment. Plant Physiol. 2007;144(2):538–49.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  31. Jung W, Seo Y. Identification of novel C-repeat binding factor (CBF) genes in rye (Secale cereale L.) and expression studies. Gene. 2019;684:82–94.

    CAS  PubMed  Article  Google Scholar 

  32. Sieber AN, Longin CF, Leiser WL, Wurschum T. Copy number variation of CBF-A14 at the Fr-A2 locus determines frost tolerance in winter durum wheat. Theor Appl Genet. 2016;129(6):1087–97.

    CAS  PubMed  Article  Google Scholar 

  33. Takumi S, Koike A, Nakata M, Kume S, Ohno R, Nakamura C. Cold-specific and light-stimulated expression of a wheat (Triticum aestivum L.) Cor gene Wcor15 encoding a chloroplast-targeted protein. J Exp Bot. 2003;54(391):2265–74.

    CAS  PubMed  Article  Google Scholar 

  34. Qin Y, Tian Y, Liu X. A wheat salinity-induced WRKY transcription factor TaWRKY93 confers multiple abiotic stress tolerance in Arabidopsis thaliana. Biochem Biophys Res Commun. 2015;464(2):428–33.

    CAS  PubMed  Article  Google Scholar 

  35. Guo Q, Li X, Niu L, Jameson PE, Zhou W. Transcription-associated metabolomic adjustments in maize occur during combined drought and cold stress. Plant Physiol. 2021;186(1):677–95.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  36. Alotaibi F, Alharbi S, Alotaibi M, Al Mosallam M, Motawei M, Alrajhi A. Wheat omics: classical breeding to new breeding technologies. Saudi J Biol Sci. 2021;28(2):1433–44.

    CAS  PubMed  Article  Google Scholar 

  37. Yan J, Qian L, Zhu W, Qiu J, Lu Q, Wang X, et al. Integrated analysis of the transcriptome and metabolome of purple and green leaves of Tetrastigma hemsleyanum reveals gene expression patterns involved in anthocyanin biosynthesis. PLoS One. 2020;15(3):e0230154.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  38. Zheng X, Shi M, Wang J, Yang N, Wang K, Xi J, et al. Isoform sequencing provides insight into freezing response of common wheat (Triticum aestivum L.). Front Genet. 2020;11:462.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  39. Kalapos B, Dobrev P, Nagy T, Vitamvas P, Gyorgyey J, Kocsy G, et al. Transcript and hormone analyses reveal the involvement of ABA-signalling, hormone crosstalk and genotype-specific biological processes in cold-shock response in wheat. Plant Sci. 2016;253:86–97.

    CAS  PubMed  Article  Google Scholar 

  40. Díaz M, Soresi DS, Basualdo J, Cuppari SJ, Carrera A. Transcriptomic response of durum wheat to cold stress at reproductive stage. Mol Biol Rep. 2019;46:2427–45.

    PubMed  Article  CAS  Google Scholar 

  41. Lv L, Zhang W, Sun L, Zhao A, Zhang Y, Wang L, et al. Gene co-expression network analysis to identify critical modules and candidate genes of drought-resistance in wheat. PLoS One. 2020;15(8):e0236186.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  42. Peng LN, Xu YQ, Wang X, Feng X, Zhao QQ, Feng SS, et al. Overexpression of paralogues of the wheat expansin gene TaEXPA8 improves low-temperature tolerance in Arabidopsis. Plant Biol (Stuttg). 2019;21(6):1119–31.

    CAS  PubMed  Article  Google Scholar 

  43. Liu P, Xu ZS, Pan-Pan L, Hu D, Chen M, Li LC, et al. A wheat PI4K gene whose product possesses threonine autophophorylation activity confers tolerance to drought and salt in Arabidopsis. J Exp Bot. 2013;64(10):2915–27.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  44. Buscher JM, Czernik D, Ewald JC, Sauer U, Zamboni N. Cross-platform comparison of methods for quantitative metabolomics of primary metabolism. Anal Chem. 2009;81(6):2135–43.

    CAS  PubMed  Article  Google Scholar 

  45. Hao Z, Liu S, Hu L, Shi J, Chen J. Transcriptome analysis and metabolic profiling reveal the key role of carotenoids in the petal coloration of Liriodendron tulipifera. Hortic Res. 2020;7:70.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  46. Nicholson JK, Lindon JC, Holmes E. 'Metabonomics': understanding the metabolic responses of living systems to pathophysiological stimuli via multivariate statistical analysis of biological NMR spectroscopic data. Xenobiotica. 1999;29(11):1181–9.

    CAS  PubMed  Article  Google Scholar 

  47. Xiaojian W, Kangfeng C, Guoping Z, Fanrong Z. Metabolite profiling of barley grains subjected to water stress: to explain the genotypic difference in drought-induced impacts on malting quality. Front Plant Sci. 2017;8:1547.

    Article  Google Scholar 

  48. Kaplan F, Kopka J, Sung DY, Zhao W, Popp M, Porat R, et al. Transcript and metabolite profiling during cold acclimation of Arabidopsis reveals an intricate relationship of cold-regulated gene expression with modifications in metabolite content. Plant J Cell Mol Biol. 2010;50(6):967–81.

    Article  CAS  Google Scholar 

  49. Hyacinthe LG, Fontaine JX, Roland M, Jérme P, François M, Françoise G, et al. NMR‐based Metabolomics to Study the Cold‐acclimation Strategy of Two Miscanthus Genotypes[J]. Phytochemical Analysis Pca. 2017;28:58-67.

  50. Zhao Y, Ke X, Li J, Li S, Zhang S, Yang X. Integrated transcriptomics and metabolomics analyses provide insights into cold stress response in wheat. Crop J. 2019;7(6):857–66.

    Article  Google Scholar 

  51. Kan L, Liao Q, Chen Z, Wang S, Ma Y, Su Z, et al. Dynamic transcriptomic and Metabolomic analyses of Madhuca pasquieri (Dubard) H. J. Lam during the post-germination stages. Front Plant Sci. 2021;12:731203.

    PubMed  PubMed Central  Article  Google Scholar 

  52. Takahashi D, Li B, Nakayama T, Kawamura Y, Uemura M. Plant plasma membrane proteomics for improving cold tolerance. Front Plant Sci. 2013;4:90.

    PubMed  PubMed Central  Google Scholar 

  53. Guo X, Liu D, Chong K. Cold signaling in plants: insights into mechanisms and regulation. J Integr Plant Biol. 2018;60(9):745–56.

    PubMed  Article  Google Scholar 

  54. Chan Z, Wang Y, Cao M, Gong Y, Mu Z, Wang H, et al. RDM4 modulates cold stress resistance in Arabidopsis partially through the CBF-mediated pathway. New Phytol. 2016;209(4):1527–39.

    CAS  PubMed  Article  Google Scholar 

  55. Schulz E, Tohge T, Zuther E, Fernie AR, Hincha DK. Flavonoids are determinants of freezing tolerance and cold acclimation in Arabidopsis thaliana. Sci Rep. 2016;6:34027.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  56. Tian Y, Peng K, Bao Y, Zhang D, Meng J, Wang D, et al. Glucose-6-phosphate dehydrogenase and 6-phosphogluconate dehydrogenase genes of winter wheat enhance the cold tolerance of transgenic Arabidopsis. Plant Physiol Biochem. 2021;161:86–97.

    CAS  PubMed  Article  Google Scholar 

  57. Ahmed N, Park J, Jung H, Hur Y, Nou I. Anthocyanin biosynthesis for cold and freezing stress tolerance and desirable color in Brassica rapa. Funct Integr Genomics. 2015;15(4):383–94.

    CAS  PubMed  Article  Google Scholar 

  58. Hwarari D, Guan Y, Ahmad B, Movahedi A, Min T, Hao Z, et al. ICE-CBF-COR signaling Cascade and its regulation in plants responding to cold stress. Int J Mol Sci. 2022;23(3):1549.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  59. Pła ek A, Rapacz M, Hura K. Relationship between quantum efficiency of PSII and cold-induced plant resistance to fungal pathogens. Acta Physiol Plant. 2004;26(2):141–8.

    Article  Google Scholar 

  60. Zong JM, Li XW, Zhou YH, Wang FW, Wang N, Dong YY, et al. The AaDREB1 transcription factor from the cold-tolerant plant Adonis amurensis enhances abiotic stress tolerance in transgenic plant. Int J Mol Sci. 2016;17(4):611.

    PubMed Central  Article  CAS  Google Scholar 

  61. Zhang Y, Yu H, Yang X, Li Q, Ling J, Wang H, et al. CsWRKY46, a WRKY transcription factor from cucumber, confers cold resistance in transgenic-plant by regulating a set of cold-stress responsive genes in an ABA-dependent manner. Plant Physiol Biochem. 2016;108:478–87.

    CAS  PubMed  Article  Google Scholar 

  62. Zhuo X, Zheng T, Zhang Z, Zhang Y, Jiang L, Ahmad S, et al. Genome-wide analysis of the NAC transcription factor gene family reveals differential expression patterns and cold-stress responses in the Woody Plant Prunus mume. Genes (Basel). 2018;9(10):494.

    Article  CAS  Google Scholar 

  63. Zuo ZF, Kang HG, Park MY, Jeong H, Sun HJ, Song PS, et al. Zoysia japonica MYC type transcription factor ZjICE1 regulates cold tolerance in transgenic Arabidopsis. Plant Sci. 2019;289:110254.

    CAS  PubMed  Article  Google Scholar 

  64. Saha G, Park JI, Jung HJ, Ahmed NU, Kayum MA, Chung MY, et al. Genome-wide identification and characterization of MADS-box family genes related to organ development and stress resistance in Brassica rapa. BMC Genomics. 2015;16:178.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  65. Yang X, Luo Y, Bai H, Li X, Tang S, Liao X, et al. DgMYB2 improves cold resistance in chrysanthemum by directly targeting DgGPX1. Hortic Res. 2022;9:uhab028.

    PubMed Central  Article  Google Scholar 

  66. Zhao Q, Xiang X, Liu D, Yang A, Wang Y. Tobacco transcription factor NtbHLH123 confers tolerance to cold stress by regulating the NtCBF pathway and reactive oxygen species homeostasis. Front Plant Sci. 2018;9:381.

    PubMed  PubMed Central  Article  Google Scholar 

  67. Chezem WR, Clay NK. Regulation of plant secondary metabolism and associated specialized cell development by MYBs and bHLHs. Phytochemistry. 2016;131:26–43.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  68. Guo J, Wu Y, Wang G, Wang T, Cao F. Integrated analysis of the transcriptome and metabolome in young and mature leaves of Ginkgo biloba L. Ind Crop Prod. 2020;143:111906.

    CAS  Article  Google Scholar 

  69. He L, Liu X, Wu Z, Teng N. Transcriptome and metabolome analyses provide insights into the Stomium degeneration mechanism in lily. Int J Mol Sci. 2021;22(22):12124.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  70. Ji X, Ren J, Zhang Y, Lang S, Wang D, Song X. Integrated analysis of the metabolome and transcriptome on anthocyanin biosynthesis in four developmental stages of Cerasus humilis Peel coloration. Int J Mol Sci. 2021;22(21):11880.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  71. Bao Y, Yang N, Meng J, Wang D, Fu L, Wang J, et al. Adaptability of winter wheat Dongnongdongmai 1 (Triticum aestivum L.) to overwintering in alpine regions. Plant Biol (Stuttg). 2021;23(3):445–55.

    CAS  Article  Google Scholar 

  72. Yu T, Zhang J, Cao J, Li X, Li S, Liu C, et al. Metabolic insight into cold stress response in two contrasting maize lines. Life (Basel). 2022;12(2):282.

    CAS  Google Scholar 

  73. Yang L, Wang Z, Zhang A, Bhawal R, Li C, Zhang S, et al. Reduction of the canonical function of a glycolytic enzyme enolase triggers immune responses that further affect metabolism and growth in Arabidopsis. Plant Cell. 2021;34:1745–67.

    Article  Google Scholar 

  74. Guy CL, Huber JL, Huber SC. Sucrose phosphate synthase and sucrose accumulation at low temperature. Plant Physiol. 1992;100(1):502–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  75. Keller I, Mudsam C, Rodrigues CM, Kischka D, Zierer W, Sonnewald U, et al. Cold-triggered induction of ROS- and Raffinose metabolism in freezing-sensitive taproot tissue of sugar beet. Front Plant Sci. 2021;12:715767.

    PubMed  PubMed Central  Article  Google Scholar 

  76. Liu Y, Zhang L, Meng S, Liu Y, Zhao X, Pang C, et al. Expression of galactinol synthase from Ammopiptanthus nanus in tomato improves tolerance to cold stress. J Exp Bot. 2020;71(1):435–49.

    CAS  PubMed  Article  Google Scholar 

  77. Wang A, Li J, Al-Huqail AA, Al-Harbi MS, Ali EF, Wang J, et al. Mechanisms of chitosan nanoparticles in the regulation of cold stress resistance in Banana plants. Nanomaterials (Basel). 2021;11(10):2670.

    CAS  Article  Google Scholar 

  78. Vfmab C, Qlb D, Zza C, Fc A, Xs A, Zwc E, et al. Proteomic variation in Vitis amurensis and V. vinifera buds during cold acclimation. Scientia Horticulturae, 263.

  79. Dong N, Lin H. Contribution of phenylpropanoid metabolism to plant development and plant-environment interactions. J Integr Plant Biol. 2021;63(1):180–209.

    CAS  PubMed  Article  Google Scholar 

  80. Wang X, Liu Y, Han Z, Chen Y, Huai D, Kang Y, et al. Integrated Transcriptomics and metabolomics analysis reveal key metabolism pathways contributing to cold tolerance in Peanut. Front Plant Sci. 2021;12:752474.

    PubMed  PubMed Central  Article  Google Scholar 

  81. Fürtauer L, Weiszmann J, Weckwerth W, Nägele T. Dynamics of plant metabolism during cold acclimation. Int J Mol Sci. 2019;20(21):5411.

    PubMed Central  Article  CAS  Google Scholar 

  82. Schulz E, Tohge T, Zuther E, Fernie A, Hincha D. Flavonoids are determinants of freezing tolerance and cold acclimation in Arabidopsis thaliana. Sci Rep. 2016;6:34027.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  83. Watanabe M, Ayugase J. Effect of low temperature on flavonoids, oxygen radical absorbance capacity values and major components of winter sweet spinach (Spinacia oleracea L.). J Sci Food Agric. 2015;95(10):2095–104.

    CAS  PubMed  Article  Google Scholar 

  84. Kirakosyan A, Seymour E, Kaufman P, Warber S, Bolling S, Chang S. Antioxidant capacity of polyphenolic extracts from leaves of Crataegus laevigata and Crataegus monogyna (hawthorn) subjected to drought and cold stress. J Agric Food Chem. 2003;51(14):3973–6.

    CAS  PubMed  Article  Google Scholar 

  85. Liu L, Cang J, Yu J, Wang X, Huang R, Wang J, et al. Effects of exogenous abscisic acid on carbohydrate metabolism and the expression levels of correlative key enzymes in winter wheat under low temperature. Biosci Biotechnol Biochem. 2013;77(3):516–25.

    CAS  PubMed  Article  Google Scholar 

  86. Zhang ZY, Liu HH, Sun C, Ma QB, Bu HY. A C2H2 zinc-finger protein OsZFP213 interacts with OsMAPK3 to enhance salt tolerance in rice. J Plant Physiol. 2018;229:100–10.

    CAS  PubMed  Article  Google Scholar 

  87. Liu J, Tang X, Wang Y, Zang Y, Zhou B. A Zostera marina manganese superoxide dismutase gene involved in the responses to temperature stress. Gene. 2016;575(2 Pt 3):718–24.

    CAS  PubMed  Article  Google Scholar 

  88. Wang Y, Liang C, Meng Z, Li Y, Abid MA, Askari M, et al. Leveraging atriplex hortensis choline monooxygenase to improve chilling tolerance in cotton. Environ Exp Bot. 2019;162:364–73.

    CAS  Article  Google Scholar 

  89. Li XC, Xing YZ, Jiang X, Qiao J, Tan HL, Tian Y, et al. Identification and characterization of the catalase gene pycat from the red alga pyropia yezoensis (bangiales, rhodophyta). J Phycol. 2012;48(3):664–9.

    CAS  PubMed  Article  Google Scholar 

  90. Vieira S, Silva TM, Glória M. Influence of processing on the levels of amines and proline and on the physico-chemical characteristics of concentrated orange juice. Food Chem. 2010;119(1):7–11.

    CAS  Article  Google Scholar 

  91. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12(4):357–60.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  92. Bader GD, Hogue CW. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 2003;4(1, article 2):2.

    PubMed  PubMed Central  Article  Google Scholar 

  93. Zheng Y, Jiao C, Sun H, Rosli HG, Pombo MA, Zhang P, et al. iTAK: A program for genome-wide prediction and classification of plant transcription factors, transcriptional regulators, and protein kinases. Mol Plant. 2016;9(12):1667–70.

    CAS  PubMed  Article  Google Scholar 

  94. Oliveros JC. Venny. An interactive tool for comparing lists with Venn's diagrams. 2007-2015,

  95. Law SR, Chrobok D, Juvany M, Delhomme N, Lindén P, Brouwer B, et al. Darkened leaves use different metabolic strategies for senescence and survival. Plant Physiol. 2018;177:00062.02018.

    Article  CAS  Google Scholar 

  96. Jasmine C, Othman S, Carin L, Iurie C, Shuzhao L, et al. MetaboAnalyst 4.0: towards more transparent and integrative metabolomics analysis. Nucleic Acids Res. 2018;W1:W486–94.

    Google Scholar 

Download references


We appreciated that the Hebei Academy of Agriculture and Forestry Sciences for permission to conduct wheat research.


This research was supported by Hebei Scientific and Technological Innovation Team of Modern Wheat Seed Industry (21326318D), HAAFS Science and Technology Innovation Special Project (2022KJCXZX-LYS-2), and Key R & D project of Hebei Province (20326313D). Funders did not participate in the design of the study, analysis of the results and writing of the manuscript, but provided financial support for the manuscript.

Author information

Authors and Affiliations



Conceptualization, LL; Data curation, CD and AZ; Formal analysis, LL and AZ; Investigation, LL, CD and YL; Methodology, LL and YZ; Project administration, HL and XC; Supervision, HL and XC; Validation, LL, YL and YZ; Writing – review & editing, LL. All authors reviewed the manuscript.

Corresponding authors

Correspondence to Hui Li or Xiyong Chen.

Ethics declarations

Ethics approval and consent to participate

The collection and experimental research of wheat material involved in this study were approved by Hebei Academy of Agriculture and Forestry Sciences (HAAFS) in China and comply with the guidelines and legislation in HAAFS.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing 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.

Sequencing output statistics of 12 samples. Table S2. Detailed information of identified differentially expressed genes. Table S3. DEGs identified in all comparison groups. Table S4. GO enrichment analysis of the DEGs. Table S5. KEGG enrichment analysis of the DEGs. Table S6. Significant differences in transcription factors between CK and MU. Table S7. Differentially expressed metabolites identified in all comparison groups. Table S8. Venn diagram of differentially expressed metabolites were found in all comparison groups. Table S9. Significantly differential metabolites and DEGs annotated in phenylpropanoid/flavonoid biosynthesis pathways. Table S10. Significantly differential metabolites and DEGs annotated in sugar and amino acids metabolic pathways. Table S11. Primers of genes used in RT-qPCR.

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 The Creative Commons Public Domain Dedication waiver ( 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

Lv, L., Dong, C., Liu, Y. et al. Transcription-associated metabolomic profiling reveals the critical role of frost tolerance in wheat. BMC Plant Biol 22, 333 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Wheat (Triticum aestivum)
  • Frost-resistant
  • Transcriptome
  • Metabolome
  • Flavonoid biosynthesis
  • Sucrose and amino acid biosynthesis