Integration of Metabolome and Transcriptome Reveals the Major Metabolic Pathways and Potential Biomarkers in Response to Freeze-Stress Regulation in Apple (Malus domestica)

Abstract

Freezing stress is the main factor affecting the normal growth and distribution of plants. The safe overwintering of a perennial deciduous plant is a crucial link to ensuring its survival and yield. However, little is known about the molecular mechanism of its gene regulation metabolites as related to its freeze-tolerance. In order to enhance our comprehension of freeze-tolerance metabolites and gene expression in dormant apple trees, we examined the metabolic and transcriptomic differences between ‘Ralls’ and ‘Fuji’, two apple varieties with varying degrees of resistance to freezing. The results of the freezing treatment showed that ‘Ralls’ had stronger freeze-tolerance than ‘Fuji’. We identified 302, 334, and 267 up-regulated differentially accumulated metabolites (DAMs) and 408, 387, and 497 down-regulated DAMs between ‘Ralls’ and ‘Fuji’ under −10, −15, and −20 °C treatment, respectively. A total of 359 shared metabolites were obtained in the upward trend modules, of which 62 metabolites were associated with 89 pathways. The number of up-regulated genes accounted for 50.2%, 45.6%, and 43.2% of the total number of differentially expressed genes (DEGs), respectively, at −10, −15, and −20 °C. Through combined transcriptome and metabolome analysis, we identified 12 pathways that included 16 DAMs and 65 DEGs. Meanwhile, we found that 20 DEGs were identified in the phenylpropanoid biosynthesis pathway and its related pathways, involving the metabolism of p-Coumaroyl-CoA, 7, 4′-Dihydroxyflavone, and scolymoside. These discoveries advance our comprehension of the molecular mechanism underlying apple freeze-tolerance and provide genetic material for breeding apple cultivars with enhanced freeze-tolerance.

1. Introduction

Freezing stress refers to the adverse effects on plants caused by exposure to freezing temperatures, which can damage plant tissues and reduce plant growth and productivity [1,2]. Freezing stress can have significant effects on fruit trees, particularly those that are not adapted to cold temperatures [3,4]. Exposure to freezing temperatures can damage fruit tree tissues, including the bark, buds, flowers, and fruits, resulting in reduced growth, yield, and even the death of the tree [4,5]. Specifically, freezing stress can cause damage to the cell structure and membranes of fruit tree tissues, leading to cell dehydration, rupture, and death [6]. This can result in various symptoms, such as the browning and wilting of leaves, flower and fruit drop, the cracking and deformation of fruits, and the dieback of branches [7,8,9]. The severity of the effects of freezing stress on fruit trees depends on several factors, including the duration and intensity of exposure, the stage of growth, and the species and cultivar of the tree [10,11,12]. Some fruit tree species, such as apples, pears, and peaches, are more tolerant to freezing stress than others, like citrus and avocado trees. Plants have various mechanisms to adapt to these stress conditions, including changes in gene expression, the accumulation of osmoprotectants, and modification of membrane fluidity, among others [13,14,15]. However, repeated or prolonged exposure to freezing temperatures can still cause significant damage to plants, particularly in crops and other economically important species.

Metabolome and transcriptome are two important levels of biological information that can provide insights into the molecular mechanisms underlying freezing stress in plants. The field of metabolomics is concerned with the investigation of small molecules, also known as metabolites, which participate in a range of biochemical processes within plants [16,17]. Freezing stress can cause changes in the metabolic profile of plants, such as the accumulation of compatible solutes, sugars, and proteins, which can help to protect cells from freezing damage [10,14,18]. Metabolomics can be used to identify these metabolites and understand their functions in regard to freeze-tolerance [16,17]. Transcriptomics, on the other hand, is the study of the full set of transcripts, or RNA molecules, in a cell or tissue. These transcripts represent the genes that are being actively expressed at a given time and can provide insights into the molecular pathways and regulatory networks involved in the response to freezing stress [19,20,21]. Transcriptomics can be used to identify genes that are up-regulated or down-regulated in response to freezing stress and to understand their functions in freeze-tolerance [22,23]. Combining metabolomics and transcriptomics can provide a more comprehensive view of the molecular mechanisms underlying freezing stress in plants [24,25]. This approach, known as metabolome-transcriptome analysis, can help to identify key metabolic and genetic pathways involved in the response to freezing stress and to develop strategies for improving freeze-tolerance in crops and other economically important plants.

The molecular mechanisms underlying the response to freezing stress involve a complex network of signaling pathways and gene expression changes. During cold acclimation, plants undergo a series of physiological and biochemical changes that increase their tolerance to freezing stress [22,26,27]. This process involves the upregulation of genes encoding for proteins involved in various protective mechanisms, such as the accumulation of osmoprotectants (e.g., sugars and amino acids), the synthesis of anti-freeze proteins, and the activation of antioxidant enzymes [28,29]. The expression of these genes is controlled by a number of transcription factors and regulatory proteins, which are themselves regulated by multiple signaling pathways. These pathways include the calcium-mediated signaling pathway, the mitogen-activated protein kinases (MAPK) signaling pathway, the abscisic acid (ABA) signaling pathway, and the INDUCER OF CBF EXPRESSION 1 (ICE)-C-REPEAT BINDING FACTOR (CBF) signaling pathway, among others [30,31,32,33].

Apple is a fruit tree species with a large production area and abundant yield worldwide, with an excellent potential presence on the international market. Apple trees yield less and have a limited cultivation range under low temperature stress; therefore, much research is focused on understanding the response of apple trees to such stress. At present, studies on the molecular mechanism of apple tree response to freezing stress mainly focus on the molecular model of transcription factors, but the analysis of the molecular mechanism of apple tree frost resistance based on multi-omics analysis is rarely carried out. In this study, we first compared the freeze-tolerance of annual branches from different varieties of apple trees and determined their limit temperatures for freezing resistance. Transcriptome and metabolome analysis identified potential critical pathways and genes involved in apple response to freezing stress. The study provides insights into the genetic factors underlying tolerance to freezing in apples and can aid in the development of new crop varieties that are more resilient and productive, which is crucial for global fruit security and sustainability.

2. Materials and Methods

2.1. Plant Materials

Two apple cultivars (‘Ralls’ and ‘Fuji’) were selected for this study. ‘Ralls’ apples were grown in the apple orchard area in Yong’an village (40.44 N, 115.51 E), Zhangjiakou City, and ‘Fuji’ apples were grown in the apple cultivation area of Yuanfang village (38.43 N, 113.89 E), Shijiazhuang City, China. The plant materials used in the experiment were all from 5-year-old apple trees. Annual branches were collected as experimental material on 15 December 2021. All plant materials were stored at a constant temperature of −10 °C for 7 d.

2.2. Measurement of Physiological Indexes and Recovery Culture of Plant Materials

Plant materials with constant growth potential were treated at low temperatures. The materials treated at −10 °C were used as the control group. Plant materials in the experimental group were treated at −15 °C, −20 °C, −25 °C, and −30 °C for 12 h, respectively. The physiological parameters of plant samples were evaluated, and each sample was subjected to three biological replicates. The content of soluble sugar was determined using the anthrone method, and the content of soluble protein was determined via the colorimetric method with Coomassie Brilliant Blue G−250 as the colorimetric liquid [26]. We utilized a conductivity meter to determine the relative electrolyte leakage (REL) of plants treated under freezing conditions, as demonstrated previously [34]. The rate of temperature change was 4 °C per hour, whether the temperature was decreased or increased. The freeze-damaged plant materials were cultured in MS medium for 7 d in an artificial climate room (16 h light/8 h dark cycles, 24 °C, and 70% relative humidity), and then the leaf bud germination rate was measured.

2.3. Metabolomic and Bioinformatic Analysis

We selected the freeze-treated annual apple branches for metabolome detection, based on the physiological index data of the branches under low temperature treatment. ‘Ralls’ materials in the experimental group were treated at −10 °C, −15 °C, −20 °C, −25 °C, and −30 °C for 12 h, respectively. ‘Fuji’ materials in the experimental group were treated at −10 °C, −15 °C, and −20 °C for 12 h, respectively. In this study, we conducted six biological replicates to increase the reliability of the data.

We referred to previous research methods to extract metabolites from annual apple branches [35]. The liquid chromatography-mass spectrometry system used for metabolomics analysis consists of a Waters Acquity I-Class PLUS ultra-high-performance liquid chromatography coupled with a Waters Xevo G2-XS QTOF high-resolution mass spectrometer. The chromatographic column used was the Acquity UPLC HSS T3 column (1.8 um 2.1 × 100 mm) purchased from Waters [36]. The raw data collected using MassLynx V4.2 were processed through Progenesis QI software for peak picking and alignment. The identification was based on Progenesis QI software with an online METLIN database, public databases, and a custom-built library by BGI. Theoretical fragmentation recognition was performed with a mass deviation of 100 ppm for precursor ions and 50 ppm for fragment ions [37,38]. Metabolites were detected by searching against the METLIN database and a self-constructed sequencing platform database. Additionally, theoretical fragments were recognized based on parent ion mass deviation <100 ppm and fragment ion mass deviation <50 ppm [36]. We normalized the original peak area information with the total peak area and carried out subsequent analysis. The quality control samples and the repeatability of the samples within the group were evaluated using Spearman correlation analysis and principal component analysis (PCA). The classification and pathway information on compounds in the study were searched using Lipidmaps, the Kyoto Encyclopedia of Genes and Genomes (KEGG), and Human Metabolome Database (HMDB) databases [39,40]. The significance (p-value) of the difference for each compound was determined by conducting t-tests based on grouping information to calculate and compare fold change differences. Differential metabolites were identified by integrating the difference multiples, p-values, and VIP values from the OPLS-DA model using screening criteria of fold change (FC) > 2, p-value < 0.05, and VIP ≥ 1. The significance of the differential metabolites in KEGG pathway enrichment was determined by conducting a hypergeometric distribution test. We used the omicShare online tool (https://www.omicshare.com/tools/Home/Soft/trend) (accessed on 10 October 2022) for trend analysis, with parameter settings of p-value < 0.05 and fold change >2.

2.4. Transcriptomics and Bioinformatics Analysis

The samples used for transcriptome sequencing were consistent with the metabolome sequencing samples. The TaKaRa MiniBEST Plant RNA Extraction Kit was utilized to extract the total RNA from the samples, and their concentration and purity were evaluated through a NanoDrop 2000 spectrophotometer. Biomarker Technologies Co., Ltd. (Beijing, China) finished constructing the cDNA library and sequencing the transcriptome. Filtered data obtained from the Illumina high-throughput sequencing platform excluded reads with joints, N ratio exceeding 10%, and mass values Q ≤ 10 accounting for over 50% of the entire read. The quality score was determined to estimate base-calling error probability using the filtered data [41]. We used the high-quality genome of the apple anther-derived homozygous line HFTH1 as a reference genome [42]. We utilized the HISAT2 software to promptly and precisely match clean reads against the reference genome, thus acquiring the positioning data of reads on the reference genome. StringTie software was used to align the aforementioned reads, and the resulting transcriptome was employed for further analysis [43]. The main steps of gene expression calculation are as follows: (1) count the number of fragments (or reads) that map to each gene in the transcriptome; (2) normalize this count by the length of the gene and the total number of fragments in the sample; (3) multiply the normalized count by one million to obtain a per-million scaling factor; (4) express the final result as FPKM units, which represent the estimated number of fragments per kilobase of gene length per million total fragments in the sample [44]. The DESeq2 program was utilized to conduct differential expression analysis among sample groups [45]. The selection criteria included a corrected false discovery rate (FDR) that was less than 0.05 and an absolute fold change (FC) greater than 2 (i.e., |log₂(FC)| > 1), while also accounting for the probability (p-value) of hypothesis testing. Genes with a p-value lower than 0.05 were deemed to be differently expressed [46].

2.5. Detection of Relative Gene Expression

Relative expression levels were measured in the annual branches of ‘Ralls’ and ‘Fuji’ treated at −10 °C, −15 °C, and −20 °C. Total RNA was isolated with the RNAprep Pure Polysaccharide Polyphenol Plant Total RNA Extraction Kit from TianGen Biotech (Beijing) Co., Ltd. (Beijng, China), and the cDNA library was constructed using the PrimeScript™ II 1st Strand cDNA Synthesis Kit from Takara Biomedical Technology (Beijing) Co., Ltd., Beijing, Chinas. qPCR primers were designed using Integrated DNA Technologies (IDT) based on the CDS sequence of the target gene. The primer and its sequence information can be found in Supplementary Table S1. The TB Green intercalating fluorescent dye method was utilized for real-time PCR. The reaction system and conditions followed the instructions for the TaKaRa TB Green™ Premix Ex Taq™ II (Tli RNaseH Plus) Kit. The apple Actin gene was employed as the endogenous control, and the experiment was conducted with three biological replicates and three technical replicates. The 2^−ΔΔCt method was utilized for determining the relative expression levels of genes [47].

3. Results

3.1. Physiological Index Evaluation of Apple under Freezing Stress

We employed a randomized factorial design to assess the impact of freezing damage on both ‘Ralls’ and ‘Fuji’ apples. Overall, the soluble sugar and protein content of ‘Ralls’ was higher than that of ‘Fuji’ (Figure 1A,B). The soluble sugar content in ‘Ralls’ and ‘Fuji’ showed a trend of decreasing after increasing with decreasing temperature. The soluble sugar content in ‘Fuji’ and ‘Ralls’ reached the highest at −15 °C and −25 °C, respectively (Figure 1A). The variation of soluble protein content was not as obvious as that of soluble sugar content with decreasing temperature. The soluble protein content of ‘Fuji’ and ‘Ralls’ reached the highest at −25 °C, which were 2.87 mg/g and 2.81 mg/g, respectively (Figure 1B). When the temperature was lower than −15 °C, the relative electrolyte leakage of apples increased rapidly, especially ‘Fuji’ (Figure 1C and Figure S1). At the same time, we found that the germination rate of ‘Fuji’ in the recovery culture process after −20 °C treatment was zero. ‘Ralls’ leaf buds could not germinate in the recovery culture process after low temperature treatment below −30 °C (Figure 1D and Figure S2). In summary, the ultimate low temperatures that ‘Fuji’ and ‘Ralls’ could tolerate were about −20 °C and −30 °C, respectively. ‘Ralls’ were more resistant to freezing than ‘Fuji’.

3.2. Metabolome Profiling of Forty-Eight Apple Samples

Two apple cultivars with differential freeze tolerance (‘Ralls’ and ‘Fuji’) were selected for this study. The freeze-tolerant (‘Ralls’) material was treated at −10, −15, −20, −25, and −30 °C, and the freezing sensitive (‘Fuji’) material was treated at −10, −15, and −20 °C, respectively. In this study, we detected 1483 metabolites across forty-eight samples (Table S2). Spearman’s rank correlation indicated the reliability of biological duplication of the sample, and further confirmed the reliability of the obtained metabolites (Figure S3). We performed PCA on all samples in two and three dimensions. The contributory ratio of the three principal components was 32.63%, 13.93%, and 9.01%, respectively, and the cumulative contributory ratio was 55.57% (Figures S4 and S5). These samples were mainly divided into two groups, ‘Ralls’ and ‘Fuji’, of which the ‘Fuji’ samples were more distributed on three dimensional levels (Figures S4 and S5).

3.3. Identification and Functional Enrichment of Differentially Accumulated Metabolites

To achieve the precise identification of differentially accumulated metabolites (DAMs), we initially conducted PCA analysis and OPLS-DA on comparison groups. PCA and OPLS-DA results showed that the samples in the pairwise comparison group were clearly divided into two clusters (Figure S6). Here, we identified 302, 334, and 267 up-regulated DAMs and 408, 387, and 497 down-regulated DAMs between ‘Ralls’ and ‘Fuji’ under −10, −15, and −20 °C treatment, respectively (Figure 2A,B, Figures S7 and S8). We found that the number of down-regulated DAMs was more than the number of up-regulated DAMs, and the total number of DAMs decreased with the decrease in temperature (Figure 2A). The DAMs showed similar patterns of metabolite changes within the groups, further proving the reliability of the data (Figure 2C, Figures S9 and S10). DAMs among the three pairwise comparison groups were significantly associated with 112, 124, and 103 pathways, respectively, of which 83 pathways were shared (Table S3 and Figure S11). We found that some pathways had only one associated metabolite, while others had multiple associated metabolites (Table S3 and Figure 2D). For example, seven metabolites were identified in the alph-Linolenic acid metabolism pathway between ‘Ralls’ and ‘Fuji’ under −10 °C treatment, including trans-2-Enoyl-OPC6-CoA, 12-OPDA, norlinolenic acid, trans-2-Enoyl-OPC8-CoA, 3-Oxo-OPC6-CoA, colnelenic acid, and (+)-7-Isojasmonic acid CoA (Figure 2D).

3.4. Trend of Metabolite Content under Low Temperature Treatment

We conducted a trend analysis of metabolite response to different low temperature treatments in order to investigate the change pattern of metabolites. Five significant trend change modules were identified in ‘Ralls’ under low temperature treatment of −10, −15, −20, −25, and −30 °C (Figure 3A). Meanwhile, two significant trend change modules were identified in ‘Fuji’ under low temperature treatment of −10, −15, and −20 °C (Figure 3B). Module 19 of ‘Ralls’ and module 7 of ‘Fuji’ showed an upward trend, including 435 and 1049 metabolites, respectively (Figure 3C). Here, we obtained 359 shared metabolites in the upward trend modules, of which 62 metabolites were associated with 89 pathways (Figure 3C and Table S4). The pathways significantly associated with these metabolites involved cutin, suberine and wax biosynthesis, flavone and flavonol biosynthesis, fatty acid metabolism, biosynthesis of unsaturated fatty acids, biosynthesis of terpenoids and steroids, and biosynthesis of plant hormones (Table S4). Some metabolites were identified in multiple KEGG pathways, for example, arachidonate was identified in 23 pathways and chorismate was identified in 17 pathways (Figure 3D). These metabolites were associated with an average of 3.4 pathways, which also meant that these metabolites were involved in multiple biological processes. The change pattern analysis results of ten metabolites, including 3-Dehydroshikimate, chorismate, and arachidonate, were consistent with the trend analysis results, showing that the accumulation of metabolites showed an increasing trend with the decrease in temperature (Figure 3E).

3.5. Transcriptome Profiling and DEGs Analysis

We obtained 171.1 Gb of clean bases in 24 samples for gene expression analysis. The Q30 of each sample data reached more than 94%, with an average of 95.3%. The average CG content of all samples was 46.33%, and 91.2% of reads were aligned to the apple reference genome (Table S5). Taken together, we obtained high-quality transcriptome data, which could be used in subsequent studies.

Here, 1517, 1522, and 1640 DEGs were identified under different freezing stress conditions between ‘Ralls’ and ‘Fuji’ (Figure 4A). The number of up-regulated genes accounted for 50.2%, 45.6%, and 43.2% of the total number of genes, respectively, under −10, −15, and −20 °C. As the temperature decreased, the up-regulated gene count was gradually reduced while the count of down-regulated genes increased. (Figure 4A). We identified 240 shared up-regulated genes under different freezing stress (Figure 4B). These shared genes were divided into eight clusters based on gene expression patterns (Figure S12). Overall, these genes showed similar expression patterns in both ‘Ralls’ and ‘Fuji’. The one hundred genes contained in Cluster 8 had low expression level, and one gene contained in Cluster 1 had the highest expression level. The difference of gene expression in Cluster 3 (10 genes), Cluster 4 (34 genes), Cluster 5 (20 genes), and Cluster 7 (47 genes) was most obvious between ‘Ralls’ and ‘Fuji’ (Figure 4C,D and Figure S10).

To confirm the dependability of RNA-seq, RT-qPCR was employed to measure the relative expression levels of genes in Cluster 3. The findings indicated that the relative expression levels of all genes were comparable to those obtained from RNA-seq (Figure 4D and Figure 5). The expression of these genes was notably increased under freezing stress in both cultivars, albeit with minor variations in the degree of change. The expressions of the CYP92C6 (HF19026), C1H46 (HF06976), and NAC071 (HF37215) genes in ‘Ralls’ were significantly higher than those in ‘Fuji’ at 15 °C. The expression of BETV6 (HF21117) gene in apple was increased by 8.4- and 7.3-fold, and the expression of BETV6 gene in ‘Fuji’ was increased by 3- and 2.5-fold, respectively, at −15 °C and −20 °C (Figure 5). According to the gene expression level, the expression of these genes was related to the cultivars and treatment temperature in apple. Based on RNA-seq and relative expression levels of genes, we inferred that these genes belong to the broad spectrum of freezing stress response genes in apple.

3.6. GO and KEGG Analysis of DEGs

To understand the potential functions of up-regulated genes in response to freezing stress, we performed GO and KEGG enrichment analyses. These DEGs were mainly annotated into 17 cell components, 11 molecular functions, and 17 biological process terms (Table S6). In the biological process, more than 50 genes were enriched in five terms, including metabolic process (GO:0006591), cellular process (GO:0048522), multi-organism process (GO:0051704), biological regulation (GO:0065007), and response to stimulus (GO:0048583) terms (Table S6). Binding term (GO:0005488), which belongs to molecular function, had the most genes enriched, including 333 genes in the −10 °C treatment group, 256 genes in the −15 °C treatment group, and 234 genes in the −20 °C treatment group. Based on the KEGG database, we enriched 80 pathways based on these shared up-regulated genes, including the plant–pathogen interaction (ko04626), MAPK signaling pathway—plant (ko04016), inositol phosphate metabolism (ko00562), the fatty acid metabolism (ko01212), and the plant hormone signal transduction (ko04075) (Table S7).

3.7. Transcriptomic and Metabolomic Analysis of Apple under Freezing Stress

We identified 12 pathways, that included 16 DAMs and 65 DEGs, through combined transcriptome and metabolome analysis (Table 1). L-Serine (pos_557) was shared in the biosynthesis of amino acids pathway, cyanoamino acid metabolism pathway, glycine, serine, and threonine metabolism pathway, respectively. Meanwhile, 3, 3, and 12 DEGs were found in the biosynthesis of amino acids pathway, cyanoamino acid metabolism pathway, glycine, serine and threonine metabolism pathway, respectively. 1,3-Diaminopropane (pos_1927) was shared in the arginine and proline metabolism pathway and glycine, serine, and threonine metabolism pathway, respectively. Based on the protein sequences encoded by these genes, 55 orthologous proteins were identified, among which 29 proteins had interactions in the STRING database (Table S8 and Figure 6). These interacting proteins were associated with 12 metabolic pathways. The PSY protein belonging to the carotenoid biosynthesis pathway interacted with CAS1, CAMS1, and SMT2 proteins in the steroid biosynthesis pathway. We found that proteins related to the phenylpropanoid biosynthesis pathway interacted with multiple other pathway-related proteins. Meanwhile, the phenylpropanoid biosynthesis pathway was also an upstream essential link of the flavone and flavonol biosynthesis pathway and flavonoid biosynthesis pathway. Scolymoside and 7,4′-Dihydroxyflavone metabolites were identified as DAMs in flavone and flavonol biosynthesis pathway and flavonoid biosynthesis pathway, respectively, involving 2 and 7 DEGs. The integration of DAMs in multiple pathways showed that the changes in DAMs were regulated by multiple genes in response to freezing stress, which further regulated the freeze-tolerance of plants.

3.8. Analysis of Phenylpropanoid Biosynthesis Pathway in Apples under Low Temperature Treatment

The combined metabolome and transcriptome analysis revealed the significant involvement of the phenylpropanoid biosynthesis pathway in apple’s adaptation to freezing stress. Here, p-Coumaroyl-CoA was identified as the DAM and a connecting metabolite for both phenylpropanoid and flavonoid biosynthesis pathways. The identification of 7,4’-Dihydroxyflavone metabolites via comparative metabolome analysis confirmed their involvement in the flavonoid biosynthesis pathway, which subsequently leads to downstream flavone and flavonol biosynthesis (Figure 7A). Meanwhile, 11, 7, and 2 DEGs were identified in phenylpropanoid biosynthesis, flavonoid biosynthesis and flavone and flavonol biosynthesis pathways (Table 1). Under normal growth conditions, we observed significantly higher expression levels of these genes in ‘Ralls’ compared to ‘Fuji’. The expression level of the F25L23.50 gene (HF40459), belonging to the phenylpropanoid biosynthesis pathway, was not detected in ‘Fuji’, but showed a trend of increasing first and then decreasing with the decrease in treatment temperature in ‘Ralls’. The differentially expressed multiples of T8E3.5 (HF00746), UGT73B4 (HF29367) and F7K24.190 (HF16814) genes in phenylpropanoid biosynthesis pathway were negatively correlated with temperature. C4H (HF39985), PCBER1 (HF42225), F22H5.18 (HF07700), and PMAT1 (HF33401) genes also showed a similar trend of differentially expressed multiples (Figure 7B,C and Figure S13). In addition, the expression of other genes showed different differential expression multiples under different freezing conditions. The differentially expressed multiples of these genes at different temperatures were likely to be the key factors for the difference in freeze tolerance between ‘Fuji’ and ‘Ralls’.

4. Discussion

In the field of fruit trees, according to the statistics of apple varieties, there are at least 7500 kinds of apple varieties worldwide, of which the cultivar ‘Ralls’ and cultivar ‘Fuji’ have excellent fruit taste and are very popular among people. Apple tree varieties differ in frost resistance, as shown by several studies [20,48,49,50]. In this study, it was found that the ultimate low temperatures that ‘Fuji’ and ‘Ralls’ could tolerate were about −20 °C and −30 °C, respectively. Our results have important guiding significance for the cultivation of ‘Fuji’ and ‘Ralls’. At the same time, the difference in freeze tolerance between ‘Fuji’ and ‘Ralls’ makes for good test material to study the freeze tolerance of apples.

Metabolomics are used to study the metabolic changes that occur in plants or organisms when they are exposed to freezing stress. Cold/freezing-stress-induced metabolomic responses have been extensively investigated in Arabidopsis, and more recently in crops, grasses, and medicinal plants [51,52,53]. Cold/freezing stress can cause an increase in various metabolites, such as proline, TCA-cycle intermediates, GABA, citrulline, ascorbate, putrescine, polyamines lipids, and sugars (glucose, fructose, inositol, etc.) [18,54,55,56,57]. In the identification of metabolites of two cultivars under different low temperature stress, we found that the amount of DAMs was related to apple cultivar and temperature. Similar changes are also found in the differential metabolites of Arabidopsis, wheat, and maize [52,56,57]. The DAMs show significant differences in plant response to cold/freezing stress at the compound level. Therefore, identifying and analyzing DAMs under cold/freezing stress can help understand the biochemical pathways involved in the response to this stress and potentially identify biomarkers for breeding or engineering cold/freezing tolerant plants.

Transcriptome has been widely used in the study of plant freeze tolerance. Transcriptome analysis can provide information on DEGs, which can be useful for identifying the key regulatory pathways involved in the response to freezing stress. The identification of genes involved in plant freeze tolerance may also lead to the development of markers that can be used to select genotypes with better freezing. In recent years, transcriptome analysis has often been combined with other omics to analyze plant freeze tolerance [19,26,58,59]. Combined transcriptome and metabolome analysis is also a strategy to truly reveal plant freeze tolerance. Metabolome and transcriptome studies indicate that the collaboration of flavonoids, abscisic acid, and nitric oxide enhances the freeze tolerance of Liriope spicata [24]. Metabolome and transcriptome are combined to discover indigo biosynthesis in Phaius flavus flowers when subjected to freezing treatment [60]. By analyzing the transcriptome and metabolome together, we can identify the specific genes that regulate metabolites. The ADC1 gene has a crucial function in enabling potato to endure freezing temperatures following cold acclimation, as demonstrated by the integration of transcriptome and metabolome data that showed its involvement in the putrescine pathway [23]. In response to cold/freezing stress, certain genes are activated for the production of enzymes that induce biosynthesis of substances such as callose, fermentation products, phospholipids, starch, sugars, flavonoids, protein amino acids, GABA, and terpenoids. Meanwhile, other genes are repressed to inhibit gluconeogenesis, amino acids, photorespiration, ethylene, folic acid, betaine, fatty acid, sulfate assimilation, chlorophyll biosynthesis, and brassinosteroids [22].

Freezing stress responses are regulated by a complex interplay of genes and environmental factors, and involve a variety of physiological pathways [61,62,63,64]. In this study, we identified 12 pathways that included 16 DAMs and 65 DEGs through combined transcriptome and metabolome analysis. Phenylpropanoid biosynthesis is a metabolic pathway that produces compounds such as flavonoids, lignins, and other phenolic compounds in plants [25,65,66,67]. In this study, 20 genes involved in the ABA pathway and related pathways showed significantly up-regulated expression in ‘Ralls’ apples. Through low-temperature stress treatment of apple trees, we found that 11, 7, and 2 DEGs were identified in the phenylpropanoid biosynthesis, flavonoid biosynthesis, and flavone and flavonol biosynthesis pathways, respectively, involving the metabolism of p-Coumaroyl-CoA, 7,4′-Dihydroxyflavone, and scolymoside. Freezing stress can induce the synthesis of these compounds through the activation of specific genes involved in the pathway, providing protection against damage from cold/freezing temperatures.

5. Conclusions

In conclusion, we used freeze-sensitive and freeze-tolerant apple cultivars as experimental materials to comprehensively analyze the molecular mechanisms underlying freeze tolerance in apples by combining transcriptome and metabolome analysis. The results of the freezing treatment showed that ‘Ralls’ had stronger freeze tolerance than ‘Fuji’. We integrated metabolome and transcriptome analyses to identify 12 pathways involved in apple response to freezing stress, which encompassed 16 DAMs and 65 DEGs. Moreover, it was found that the phenylpropanoid biosynthesis pathway and its related pathways are the key pathways for the difference in freeze tolerance between ‘Ralls’ and ‘Fuji’. These discoveries advance our comprehension of the molecular mechanism underlying apple freeze tolerance and furnish genetic materials for breeding apple cultivars with enhanced freeze tolerance.

Author Contributions

Conceptualization, Y.Y. and P.L.; methodology, Y.W.; software, P.L.; validation, W.L., J.L. and P.L.; formal analysis, Y.Y.; investigation, Y.W.; resources, W.L.; data curation, Y.Y.; writing—original draft preparation, Y.Y. and P.L.; writing—review and editing, Y.Y. and P.L.; visualization, P.L.; supervision, Y.Y.; project administration, Y.Y.; funding acquisition, Y.Y. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The data information is publicly accessible at the NCBI under BioProject PRJNA977804 and the figshare—Digital Science under https://doi.org/10.6084/m9.figshare.23264861 (accessed on 2 Jun 2023).

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Physiological and biochemical indices of apple under freezing stress and recovery culture after freezing treatment. (A). Changes soluble sugar content of apple under freezing treatment at 4 (CK), −10, −15, −20, −25, −30, −35, and −40 °C. (B). Changes soluble protein content of apple under freezing treatment at 4 (CK), −10, −15, −20, −25, −30, −35, and −40 °C. (C). The REL of apple under freezing treatment at −15, −20, and −25 °C. (D). The germination percentage of apple under freezing treatment at −15, −20, and −25 °C. fs, ‘Fuji’; gg, ‘Ralls’. The single asterisk (*) and double asterisk (**) represent the significance level determined by t-test as p < 0.05 and p < 0.01, respectively.

Figure 2. Identification of differential metabolites between ‘Ralls’ and ‘Fuji’. (A). The number of differential metabolites between ‘Ralls’ and ‘Fuji’ under different low temperatures. The letters gg represent ‘Ralls’, and the letters fs represent ‘Fuji’. The numbers 1, 2, and 3 represent the low temperature treatments of −10, −15, and −20 °C, respectively. (B). The volcano plot of differential metabolites between ‘Ralls’ and ‘Fuji’ under −10 °C treatment. (C). Clustering heat map of differential metabolites between ‘Ralls’ and ‘Fuji’ under −10 °C treatment. Quantitative values of metabolites after Z-score standardization. (D). Partial network diagram of KEGG enrichment for differential metabolites between ‘Ralls’ and ‘Fuji’ under −10 °C treatment. The figure displays pathways represented by light yellow nodes, and the small connected nodes are specific metabolites identified within these pathways.

Figure 3. Trends of metabolite content in ‘Ralls’ and ‘Fuji’ under low temperature treatment. (A). The trend of metabolite content in ‘Ralls’ under low temperature treatment of −10, −15, −20, −25, and −30 °C. (B). The trend of metabolite content in ‘Fuji’ under low temperature treatment of −10, −15, and −20 °C. (C). Venn diagram of the cluster with continuous increase in metabolites. Green and red represent ‘Ralls’ and ‘Fuji’, respectively. (D). Number of pathways significantly associated with metabolites. (E). Partial shared metabolite change patterns in ‘Ralls’ and ‘Fuji’ under low temperature treatment. The abundance values are processed with log10.

Figure 4. Identification and analysis of DEGs under freezing stress. (A). The number of DEGs between ‘Ralls’ and ‘Fuji’ under different low temperatures. (B). Venn diagram of up-regulated genes. (C). The number of genes in shared gene cluster analysis. The letter C represents cluster group. (D). The level of gene expression in the Cluster 3 (C3). The letters gg represent ‘Ralls’, and the letters fs represent ‘Fuji’. The numbers 1, 2, 3, 4, and 5 represent the low temperature treatments of −10, −15, −20, −25, and −30 °C, respectively.

Figure 5. The relative expression level of genes under artificial low temperature treatment. The letter gg represents ‘Ralls’, and the letter fs represents ‘Fuji’. The single asterisk (*) and double asterisk (**) represent the significance level determined by t-test as p < 0.05 and p < 0.01, respectively.

Figure 6. Integrating metabolome and transcriptome to identify the interaction network of key gene expression proteins. Proteins with the same color and pattern markings were annotated to the same pathways. The pathway names were indicated by the full name of the adjacent or graphic note. The linear connected proteins represented protein–protein interactions.

Figure 7. Integrating phenylpropanoid biosynthesis, flavonoid biosynthesis and flavone and flavonol biosynthesis pathways in apples under low temperature treatment. (A). Schematic diagram of the association between those pathways. The DAMs identified in this study are shown in red. The pathway name is shown in purple. (B). Differential expression multiples of phenylpropanoid biosynthesis pathway-associated genes. (C). Differential expression multiples of flavonoid biosynthesis pathway-associated genes. The letters gg represent ‘Ralls’, and the letters fs represent ‘Fuji’. The numbers 1, 2, and 3 represent the low temperature treatments of −10, −15, and −20 °C, respectively.

Table 1

DAMs and DEGs identified in combination with metabolome and transcriptome.

Pathway	Metabolites ID 1	Gene ID 2
Propanoate metabolism	pos_3699	HF37733; HF39626
Phenylalanine metabolism	pos_3975	HF32403; NewGene_4237; NewGene_9458
Fatty acid metabolism	pos_2511; pos_4795	HF08643; HF37733; NewGene_7563
Arginine and proline metabolism	pos_1927; pos_4010	HF04894; HF27116; HF41275
Carotenoid biosynthesis	pos_4978	HF05780; HF05786; HF12739; HF23511; HF34033; HF06845
Flavone and flavonol biosynthesis	pos_4009	HF33398; HF33401
Flavonoid biosynthesis	pos_3214	HF00746; HF06720; HF07700; HF15427; HF21107; HF39985; HF42225
Phenylpropanoid biosynthesis	pos_2056	HF00746; HF06720; HF09713; HF16285; HF16814; HF29367; HF31060; HF39985; HF40459; NewGene_1559; NewGene_9073
Steroid biosynthesis	pos_5705	HF02832; HF16819; HF19558; HF20971; HF28860; HF28861; HF35942
Histidine metabolism	pos_2179	HF03481; HF29862
Cyanoamino acid metabolism	pos_557	HF32157; HF40459; HF44479
Glycine, serine, and threonine metabolism	pos_1927; pos_557	HF32157; HF32319; HF44479
Tryptophan metabolism	pos_865	HF09172; HF10553; HF10561; HF19130; HF22784; NewGene_6625; NewGene_9458
Biosynthesis of amino acids	pos_557	HF02826; HF02828; HF03481; HF04894; HF12391; HF27116; HF29862; HF32157; HF36209; HF38329; HF44479; NewGene_963
Purine metabolism	pos_2642; pos_3608	HF01558; HF04086; NewGene_495
Tyrosine metabolism	pos_3561	HF09884; HF32403; NewGene_4237; NewGene_9458
Cutin, suberine, and wax biosynthesis	pos_4739; pos_4973	HF05413; HF07533; HF22738; HF27176

¹ The metabolites corresponding to metabolites ID are shown in Table S2. ² The genes were consistent with the reference genome, and NewGene was the newly detected gene obtained by transcriptome sequencing.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/metabo13080891/s1. Figure S1: The relative electrolyte leakage of apple under freezing treatment at 4 (CK), −10, −15, −20, −25, −30, −35, and −40 °C. fs, ‘Fuji’; gg, ‘Ralls’; Figure S2: The germination percentage of apple under freezing treatment at 4 (CK), −10, −15, −20, −25, −30, −35, and −40 °C. fs, ‘Fuji’; gg, ‘Ralls’; Figure S3: The spearman rank correlation analysis of sample metabolites. fs, ‘Fuji’; gg, ‘Ralls’; Figure S4: Principal component analysis of two dimensions. fs, ‘Fuji’; gg, ‘Ralls’; Figure S5: Principal component analysis of three dimensions. fs, ‘Fuji’; gg, ‘Ralls’; Figure S6: Principal component analysis and orthogonal projections to latent structures—discriminant analysis of pairwise comparison groups. A, B, and C. Principal component analysis on pairwise comparison groups. D, E, and F. Orthogonal projections to latent structures—discriminant analysis of pairwise comparison groups. fs, ‘Fuji’; gg, ‘Ralls’; Figure S7: The volcano plot of DAMs between ‘Ralls’ and ‘Fuji’ under −15 °C treatment; Figure S8: The volcano plot of DAMs between ‘Ralls’ and ‘Fuji’ under −20 °C treatment; Figure S9: Clustering heat map of DAMs between ‘Ralls’ and ‘Fuji’ under −15 °C treatment. Quantitative values of metabolites after Z-score standardization. fs, ‘Fuji’; gg, ‘Ralls’; Figure S10: Clustering heat map of DAMs between ‘Ralls’ and ‘Fuji’ under −20 °C treatment. Quantitative values of metabolites after Z-score standardization. fs, ‘Fuji’; gg, ‘Ralls’; Figure S11: KEGG annotation of DAMs under different low temperatures. fs, ‘Fuji’; gg, ‘Ralls’; Figure S12: Expression patterns of shared up-regulated genes. The letter gg represents ‘Ralls’, and the letter fs represents ‘Fuji’. The number 1, 2, and 3 represent the low temperature treatment of −10, −15, and −20 °C, respectively; Figure S13: Differential expression multiples of flavone and flavonol biosynthesis pathway-associated genes. The letter gg represents ‘Ralls’, and the letter fs represents ‘Fuji’. The number 1, 2, and 3 represent the low temperature treatment of −10, −15, and −20 °C, respectively; Table S1: Primer list and their sequence information; Table S2: Information and quantification of all metabolites in forty-eight samples; Table S3. KEGG annotation of DAMs under different low temperatures; Table S4: KEGG annotations of metabolites with the same trend; Table S5: RNA-seq data statistics and quality control; Table S6: GO enrichment of up-regulated genes between ‘Ralls’ and ‘Fuji’ under different low temperatures; Table S7: KEGG enrichment of shared up-regulated genes; Table S8. Identification and annotation of orthologous sequences of DEG encoded proteins.

References

1. Takahashi, D.; Uemura, M.; Kawamura, Y. Freezing Tolerance of Plant Cells: From the Aspect of Plasma Membrane and Microdomain. Survival Strategies in Extreme Cold and Desiccation: Adaptation Mechanisms and Their Applications; Iwaya-Inoue, M.; Sakurai, M.; Uemura, M. Springer: Singapore, 2018; pp. 61-79. [DOI: https://dx.doi.org/10.1007/978-981-13-1244-1_4]

2. Li, J.; Zhang, Z.; Chong, K.; Xu, Y. Chilling tolerance in rice: Past and present. J. Plant Physiol.; 2022; 268, 153576. [DOI: https://dx.doi.org/10.1016/j.jplph.2021.153576]

3. Villouta, C.; Workmaster, B.A.; Bolivar-Medina, J.; Sinclair, S.; Atucha, A. Freezing stress survival mechanisms in Vaccinium macrocarpon Ait. terminal buds. Tree Physiol.; 2020; 40, pp. 841-855. [DOI: https://dx.doi.org/10.1093/treephys/tpaa028]

4. Charrier, G.; Ngao, J.; Saudreau, M.; Ameglio, T. Effects of environmental factors and management practices on microclimate, winter physiology, and frost resistance in trees. Front. Plant Sci.; 2015; 6, 259. [DOI: https://dx.doi.org/10.3389/fpls.2015.00259]

5. Neuner, G. Frost resistance in alpine woody plants. Front. Plant Sci.; 2014; 5, 654. [DOI: https://dx.doi.org/10.3389/fpls.2014.00654]

6. Webb, M.S.; Uemura, M.; Steponkus, P.L. A Comparison of Freezing Injury in Oat and Rye: Two Cereals at the Extremes of Freezing Tolerance. Plant Physiol.; 1994; 104, pp. 467-478. [DOI: https://dx.doi.org/10.1104/pp.104.2.467] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12232096]

7. Maurya, J.P.; Bhalerao, R.P. Photoperiod- and temperature-mediated control of growth cessation and dormancy in trees: A molecular perspective. Ann. Bot.; 2017; 120, pp. 351-360. [DOI: https://dx.doi.org/10.1093/aob/mcx061] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28605491]

8. Bratzel, F.; Turck, F. Molecular memories in the regulation of seasonal flowering: From competence to cessation. Genome Biol.; 2015; 16, 192. [DOI: https://dx.doi.org/10.1186/s13059-015-0770-6] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26374394]

9. Ambroise, V.; Legay, S.; Guerriero, G.; Hausman, J.-F.; Cuypers, A.; Sergeant, K. The Roots of Plant Frost Hardiness and Tolerance. Plant Cell Physiol.; 2020; 61, pp. 3-20. [DOI: https://dx.doi.org/10.1093/pcp/pcz196]

10. Hincha, D.K.; Zuther, E. Plant Cold Acclimation and Freezing Tolerance. Plant Cold Acclimation: Methods and Protocols; Hincha, D.K.; Zuther, E. Springer: New York, NY, USA, 2014; pp. 1-6. [DOI: https://dx.doi.org/10.1007/978-1-4939-0844-8_1]

11. Morales, M.; Munné-Bosch, S. Oxidative Stress: A Master Regulator of Plant Trade-Offs?. Trends Plant Sci.; 2016; 21, pp. 996-999. [DOI: https://dx.doi.org/10.1016/j.tplants.2016.09.002]

12. Tylewicz, S.; Petterle, A.; Marttila, S.; Miskolczi, P.; Azeez, A.; Singh, R.K.; Immanen, J.; Mähler, N.; Hvidsten, T.R.; Eklund, D.M. et al. Photoperiodic control of seasonal growth is mediated by ABA acting on cell-cell communication. Science; 2018; 360, pp. 212-215. [DOI: https://dx.doi.org/10.1126/science.aan8576]

13. Friedrich, T.; Faivre, L.; Bäurle, I.; Schubert, D. Chromatin-based mechanisms of temperature memory in plants. Plant Cell Environ.; 2019; 42, pp. 762-770. [DOI: https://dx.doi.org/10.1111/pce.13373] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29920687]

14. Kaplan, F.; Guy, C.L. β-Amylase Induction and the Protective Role of Maltose during Temperature Shock. Plant Physiol.; 2004; 135, pp. 1674-1684. [DOI: https://dx.doi.org/10.1104/pp.104.040808]

15. Rohde, P.; Hincha, D.K.; Heyer, A.G. Heterosis in the freezing tolerance of crosses between two Arabidopsis thaliana accessions (Columbia-0 and C24) that show differences in non-acclimated and acclimated freezing tolerance. Plant J.; 2004; 38, pp. 790-799. [DOI: https://dx.doi.org/10.1111/j.1365-313X.2004.02080.x]

16. Jang, C.; Chen, L.; Rabinowitz, J.D. Metabolomics and Isotope Tracing. Cell; 2018; 173, pp. 822-837. [DOI: https://dx.doi.org/10.1016/j.cell.2018.03.055]

17. Johnson, C.H.; Ivanisevic, J.; Siuzdak, G. Metabolomics: Beyond biomarkers and towards mechanisms. Nat. Rev. Mol. Cell Biol.; 2016; 17, pp. 451-459. [DOI: https://dx.doi.org/10.1038/nrm.2016.25] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26979502]

18. Cook, D.; Fowler, S.; Fiehn, O.; Thomashow, M.F. A prominent role for the CBF cold response pathway in configuring the low-temperature metabolome of Arabidopsis. Proc. Natl. Acad. Sci. USA; 2004; 101, pp. 15243-15248. [DOI: https://dx.doi.org/10.1073/pnas.0406069101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15383661]

19. 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. [DOI: https://dx.doi.org/10.3389/fpls.2021.664311]

20. Zhou, P.; Li, X.; Liu, X.; Wen, X.; Zhang, Y.; Zhang, D. Transcriptome profiling of Malus sieversii under freezing stress after being cold-acclimated. BMC Genom.; 2021; 22, 681. [DOI: https://dx.doi.org/10.1186/s12864-021-07998-0]

21. Han, Z.; Xu, X.; Zhang, S.; Zhao, Q.; Li, H.; Cui, Y.; Li, X.; Wang, L.; Chen, S.; Zhao, X. Transcriptomics Profiling of Acer pseudosieboldianum Molecular Mechanism against Freezing Stress. Int. J. Mol. Sci.; 2022; 23, 14676. [DOI: https://dx.doi.org/10.3390/ijms232314676] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36499002]

22. Kaplan, F.; Kopka, J.; Sung, D.Y.; Zhao, W.; Popp, M.; Porat, R.; Guy, C.L. 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.; 2007; 50, pp. 967-981. [DOI: https://dx.doi.org/10.1111/j.1365-313X.2007.03100.x]

23. Kou, S.; Chen, L.; Tu, W.; Scossa, F.; Wang, Y.; Liu, J.; Fernie, A.R.; Song, B.; Xie, C. The arginine decarboxylase gene ADC1, associated to the putrescine pathway, plays an important role in potato cold-acclimated freezing tolerance as revealed by transcriptome and metabolome analyses. Plant J.; 2018; 96, pp. 1283-1298. [DOI: https://dx.doi.org/10.1111/tpj.14126]

24. Peng, Z.; Wang, Y.; Zuo, W.T.; Gao, Y.R.; Li, R.Z.; Yu, C.X.; Liu, Z.Y.; Zheng, Y.; Shen, Y.Y.; Duan, L.S. Integration of Metabolome and Transcriptome Studies Reveals Flavonoids, Abscisic Acid, and Nitric Oxide Comodulating the Freezing Tolerance in Liriope spicata. Front. Plant Sci.; 2021; 12, 764625. [DOI: https://dx.doi.org/10.3389/fpls.2021.764625] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35154173]

25. Sun, S.; Fang, J.; Lin, M.; Hu, C.; Qi, X.; Chen, J.; Zhong, Y.; Muhammad, A.; Li, Z.; Li, Y. Comparative Metabolomic and Transcriptomic Studies Reveal Key Metabolism Pathways Contributing to Freezing Tolerance under Cold Stress in Kiwifruit. Front. Plant Sci.; 2021; 12, 628969. [DOI: https://dx.doi.org/10.3389/fpls.2021.628969]

26. Li, P.; Zheng, T.; Li, L.; Liu, W.; Qiu, L.; Ahmad, S.; Wang, J.; Cheng, T.; Zhang, Q. Integration of chromatin accessibility and gene expression reveals new regulators of cold hardening to enhance freezing tolerance in Prunus mume. J. Exp. Bot.; 2023; 74, pp. 2173-2187. [DOI: https://dx.doi.org/10.1093/jxb/erad027]

27. Schieber, M.; Chandel, N.S. ROS Function in Redox Signaling and Oxidative Stress. Curr. Biol.; 2014; 24, pp. R453-R462. [DOI: https://dx.doi.org/10.1016/j.cub.2014.03.034]

28. Takagi, T.; Nakamura, M.; Hayashi, H.; Inatsugi, R.; Yano, R.; Nishida, I. The Leaf-Order-Dependent Enhancement of Freezing Tolerance in Cold-Acclimated Arabidopsis Rosettes is not Correlated with the Transcript Levels of the Cold-Inducible Transcription Factors of CBF/DREB1. Plant Cell Physiol.; 2003; 44, pp. 922-931. [DOI: https://dx.doi.org/10.1093/pcp/pcg117] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/14519774]

29. Dai, F.; Huang, Y.; Zhou, M.; Zhang, G. The influence of cold acclimation on antioxidative enzymes and antioxidants in sensitive and tolerant barley cultivars. Biol. Plant.; 2009; 53, pp. 257-262. [DOI: https://dx.doi.org/10.1007/s10535-009-0048-5]

30. Guo, X.; Liu, D.; Chong, K. Cold signaling in plants: Insights into mechanisms and regulation. J. Integr. Plant Biol.; 2018; 60, pp. 745-756. [DOI: https://dx.doi.org/10.1111/jipb.12706] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30094919]

31. Shi, Y.; Ding, Y.; Yang, S. Molecular Regulation of CBF Signaling in Cold Acclimation. Trends Plant Sci.; 2018; 23, pp. 623-637. [DOI: https://dx.doi.org/10.1016/j.tplants.2018.04.002] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29735429]

32. Takezawa, D.; Komatsu, K.; Sakata, Y. ABA in bryophytes: How a universal growth regulator in life became a plant hormone?. J. Plant Res.; 2011; 124, pp. 437-453. [DOI: https://dx.doi.org/10.1007/s10265-011-0410-5]

33. Zhu, J.-K. Abiotic Stress Signaling and Responses in Plants. Cell; 2016; 167, pp. 313-324. [DOI: https://dx.doi.org/10.1016/j.cell.2016.08.029] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27716505]

34. Li, H.; Ding, Y.; Shi, Y.; Zhang, X.; Zhang, S.; Gong, Z.; Yang, S. MPK3- and MPK6-Mediated ICE1 Phosphorylation Negatively Regulates ICE1 Stability and Freezing Tolerance in Arabidopsis. Dev. Cell; 2017; 43, pp. 630-642. [DOI: https://dx.doi.org/10.1016/j.devcel.2017.09.025] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29056553]

35. Dunn, W.B.; Broadhurst, D.; Begley, P.; Zelena, E.; Francis-McIntyre, S.; Anderson, N.; Brown, M.; Knowles, J.D.; Halsall, A.; Haselden, J.N. et al. Procedures for large-scale metabolic profiling of serum and plasma using gas chromatography and liquid chromatography coupled to mass spectrometry. Nat. Protoc.; 2011; 6, pp. 1060-1083. [DOI: https://dx.doi.org/10.1038/nprot.2011.335]

36. Wang, J.; Zhang, T.; Shen, X.; Liu, J.; Zhao, D.; Sun, Y.; Wang, L.; Liu, Y.; Gong, X.; Liu, Y. et al. Serum metabolomics for early diagnosis of esophageal squamous cell carcinoma by UHPLC-QTOF/MS. Metabolomics; 2016; 12, 116. [DOI: https://dx.doi.org/10.1007/s11306-016-1050-5]

37. Shoko, T.; Maharaj, V.J.; Naidoo, D.; Tselanyane, M.; Nthambeleni, R.; Khorombi, E.; Apostolides, Z. Anti-aging potential of extracts from Sclerocarya birrea (A. Rich.) Hochst and its chemical profiling by UPLC-Q-TOF-MS. BMC Complement. Altern. Med.; 2018; 18, 54. [DOI: https://dx.doi.org/10.1186/s12906-018-2112-1]

38. Qi, D.; Brownridge, P.; Xia, D.; Mackay, K.; Gonzalez-Galarza, F.F.; Kenyani, J.; Harman, V.; Beynon, R.J.; Jones, A.R. A Software Toolkit and Interface for Performing Stable Isotope Labeling and Top3 Quantification Using Progenesis LC-MS. OMICS A J. Integr. Biol.; 2012; 16, pp. 489-495. [DOI: https://dx.doi.org/10.1089/omi.2012.0042] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22888986]

39. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res.; 2000; 28, pp. 27-30. [DOI: https://dx.doi.org/10.1093/nar/28.1.27]

40. Wishart, D.S.; Feunang, Y.D.; Marcu, A.; Guo, A.C.; Liang, K.; Vázquez-Fresno, R.; Sajed, T.; Johnson, D.; Li, C.; Karu, N. et al. HMDB 4.0: The human metabolome database for 2018. Nucleic Acids Res.; 2018; 46, pp. D608-D617. [DOI: https://dx.doi.org/10.1093/nar/gkx1089]

41. Ewing, B.; Hillier, L.; Wendl, M.C.; Green, P. Base-Calling of Automated Sequencer Traces UsingPhred. I. Accuracy Assessment. Genome Res.; 1998; 8, pp. 175-185. [DOI: https://dx.doi.org/10.1101/gr.8.3.175]

42. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods; 2015; 12, pp. 357-360. [DOI: https://dx.doi.org/10.1038/nmeth.3317]

43. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.-C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol.; 2015; 33, pp. 290-295. [DOI: https://dx.doi.org/10.1038/nbt.3122]

44. Florea, L.; Song, L.; Salzberg, S.L. Thousands of exon skipping events differentiate among splicing patterns in sixteen human tissues. F1000Res; 2013; 2, 188. [DOI: https://dx.doi.org/10.12688/f1000research.2-188.v1]

45. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol.; 2014; 15, 550. [DOI: https://dx.doi.org/10.1186/s13059-014-0550-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25516281]

46. Robinson, M.D.; McCarthy, D.J.; Smyth, G.K. edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics; 2010; 26, pp. 139-140. [DOI: https://dx.doi.org/10.1093/bioinformatics/btp616]

47. Livak, K.; Schmittgen, T. Analysis of Relative Gene Expression Data using Real-Time Quantitative PCR. Methods; 2002; 25, pp. 402-408. [DOI: https://dx.doi.org/10.1006/meth.2001.1262] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11846609]

48. Liang, Y.; Wang, S.; Zhao, C.; Ma, X.; Zhao, Y.; Shao, J.; Li, Y.; Li, H.; Song, H.; Ma, H. et al. Transcriptional regulation of bark freezing tolerance in apple (Malus domestica Borkh.). Hortic. Res.; 2020; 7, 205. [DOI: https://dx.doi.org/10.1038/s41438-020-00432-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33328456]

49. Pramsohler, M.; Hacker, J.; Neuner, G. Freezing pattern and frost killing temperature of apple (Malus domestica) wood under controlled conditions and in nature. Tree Physiol.; 2012; 32, pp. 819-828. [DOI: https://dx.doi.org/10.1093/treephys/tps046] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22628198]

50. Chen, X.; Li, S.; Zhang, D.; Han, M.; Jin, X.; Zhao, C.; Wang, S.; Xing, L.; Ma, J.; Ji, J. et al. Sequencing of a Wild Apple (Malus baccata) Genome Unravels the Differences Between Cultivated and Wild Apple Species Regarding Disease Resistance and Cold Tolerance. G3 Genes Genomes Genet.; 2019; 9, pp. 2051-2060. [DOI: https://dx.doi.org/10.1534/g3.119.400245] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31126974]

51. Xu, Y.; Fu, X. Reprogramming of Plant Central Metabolism in Response to Abiotic Stresses: A Metabolomics View. Int. J. Mol. Sci.; 2022; 23, 5716. [DOI: https://dx.doi.org/10.3390/ijms23105716]

52. Sun, C.X.; Gao, X.X.; Li, M.Q.; Fu, J.Q.; Zhang, Y.L. Plastic responses in the metabolome and functional traits of maize plants to temperature variations. Plant Biol.; 2016; 18, pp. 249-261. [DOI: https://dx.doi.org/10.1111/plb.12378]

53. Ghassemi, S.; Delangiz, N.; Asgari Lajayer, B.; Saghafi, D.; Maggi, F. Review and future prospects on the mechanisms related to cold stress resistance and tolerance in medicinal plants. Acta Ecol. Sin.; 2021; 41, pp. 120-129. [DOI: https://dx.doi.org/10.1016/j.chnaes.2020.09.006]

54. Kaplan, F.; Kopka, J.; Haskell, D.W.; Zhao, W.; Schiller, K.C.; Gatzke, N.; Sung, D.Y.; Guy, C.L. Exploring the Temperature-Stress Metabolome of Arabidopsis. Plant Physiol.; 2004; 136, pp. 4159-4168. [DOI: https://dx.doi.org/10.1104/pp.104.052142]

55. Wienkoop, S.; Morgenthal, K.; Wolschin, F.; Scholz, M.; Selbig, J.; Weckwerth, W. Integration of Metabolomic and Proteomic Phenotypes: Analysis of Data Covariance Dissects Starch and RFO Metabolism from Low and High Temperature Compensation Response in Arabidopsis thaliana. Mol. Cell. Proteom.; 2008; 7, pp. 1725-1736. [DOI: https://dx.doi.org/10.1074/mcp.M700273-MCP200] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18445580]

56. Mazzucotelli, E.; Tartari, A.; Cattivelli, L.; Forlani, G. Metabolism of γ-aminobutyric acid during cold acclimation and freezing and its relationship to frost tolerance in barley and wheat. J. Exp. Bot.; 2006; 57, pp. 3755-3766. [DOI: https://dx.doi.org/10.1093/jxb/erl141]

57. Korn, M.; Gärtner, T.; Erban, A.; Kopka, J.; Selbig, J.; Hincha, D.K. Predicting Arabidopsis Freezing Tolerance and Heterosis in Freezing Tolerance from Metabolite Composition. Mol. Plant; 2010; 3, pp. 224-235. [DOI: https://dx.doi.org/10.1093/mp/ssp105]

58. Wei, J.; Shen, Y.; Dong, X.; Zhu, Y.; Cui, J.; Li, H.; Zheng, G.; Tian, H.; Wang, Y.; Liu, Z. DNA methylation affects freezing tolerance in winter rapeseed by mediating the expression of genes related to JA and CK pathways. Front. Genet.; 2022; 13, 968494. [DOI: https://dx.doi.org/10.3389/fgene.2022.968494] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36061187]

59. Janmohammadi, M.; Zolla, L.; Rinalducci, S. Low temperature tolerance in plants: Changes at the protein level. Phytochemistry; 2015; 117, pp. 76-89. [DOI: https://dx.doi.org/10.1016/j.phytochem.2015.06.003]

60. Zhang, Y.M.; Su, Y.; Dai, Z.W.; Lu, M.; Sun, W.; Yang, W.; Wu, S.S.; Wan, Z.T.; Wan, H.H.; Zhai, J. Integration of the metabolome and transcriptome reveals indigo biosynthesis in Phaius flavus flowers under freezing treatment. PeerJ; 2022; 10, e13106. [DOI: https://dx.doi.org/10.7717/peerj.13106]

61. Zheng, G.; Dong, X.; Wei, J.; Liu, Z.; Aslam, A.; Cui, J.; Li, H.; Wang, Y.; Tian, H.; Cao, X. Integrated methylome and transcriptome analysis unravel the cold tolerance mechanism in winter rapeseed(Brassica napus L.). BMC Plant Biol.; 2022; 22, 414. [DOI: https://dx.doi.org/10.1186/s12870-022-03797-1]

62. Liu, Q.; Kasuga, M.; Sakuma, Y.; Abe, H.; Miura, S.; Yamaguchi-Shinozaki, K.; Shinozaki, K. Two Transcription Factors, DREB1 and DREB2, with an EREBP/AP2 DNA Binding Domain Separate Two Cellular Signal Transduction Pathways in Drought- and Low-Temperature-Responsive Gene Expression, Respectively, in Arabidopsis. Plant Cell; 1998; 10, pp. 1391-1406. [DOI: https://dx.doi.org/10.1105/tpc.10.8.1391]

63. Raza, A.; Charagh, S.; Garcia-Caparros, P.; Rahman, M.A.; Ogwugwa, V.H.; Saeed, F.; Jin, W. Melatonin-mediated temperature stress tolerance in plants. GM Crops Food; 2022; 13, pp. 196-217. [DOI: https://dx.doi.org/10.1080/21645698.2022.2106111] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35983948]

64. Liu, Q.; Ding, Y.; Shi, Y.; Ma, L.; Wang, Y.; Song, C.; Wilkins, K.A.; Davies, J.M.; Knight, H.; Knight, M.R. et al. The calcium transporter ANNEXIN1 mediates cold-induced calcium signaling and freezing tolerance in plants. EMBO J.; 2021; 40, e104559. [DOI: https://dx.doi.org/10.15252/embj.2020104559] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33372703]

65. Xu, K.; Zhao, Y.; Gu, J.; Zhou, M.; Gao, L.; Sun, R.-X.; Wang, W.-W.; Zhang, S.-H.; Yang, X.-J. Proteomic analysis reveals the molecular mechanism underlying the cold acclimation and freezing tolerance of wheat (Triticum aestivum L.). Plant Sci.; 2022; 318, 111242. [DOI: https://dx.doi.org/10.1016/j.plantsci.2022.111242]

66. Wang, X.; Liu, Y.; Han, Z.; Chen, Y.; Huai, D.; Kang, Y.; Wang, Z.; Yan, L.; Jiang, H.; Lei, Y. et al. Integrated Transcriptomics and Metabolomics Analysis Reveal Key Metabolism Pathways Contributing to Cold Tolerance in Peanut. Front. Plant Sci.; 2021; 12, 752474. [DOI: https://dx.doi.org/10.3389/fpls.2021.752474]

67. Schulz, E.; Tohge, T.; Winkler, J.B.; Albert, A.; Schäffner, A.R.; Fernie, A.R.; Zuther, E.; Hincha, D.K. Natural Variation among Arabidopsis Accessions in the Regulation of Flavonoid Metabolism and Stress Gene Expression by Combined UV Radiation and Cold. Plant Cell Physiol.; 2021; 62, pp. 502-514. [DOI: https://dx.doi.org/10.1093/pcp/pcab013] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33544865]