Regulation of phenylpropanoid metabolism during moderate freezing and post-freezing recovery in Dendrobium ofﬁcinale

ABSTRACT Freezing injury is a major environmental limiting factor affecting plant yield and geographical distribution. Little is known about the regulation of phenylpropanoid metabolism during freezing (FT) and post-freezing recovery (FR). Dendrobium ofﬁcinale is a valuable medicinal plant and hypersensitive to low temperature, especially freezing in winter. Here, we used metabolomics and transcriptomics analyses to reveal the regulation of of phenylpropanoid metabolism during FT and FR in D. ofﬁcinale. Most anthocyanin and flavonol signiﬁcantly increased during FR, including quercetin 3-O-glucoside (Q3G) and cyanidin 3-O-sophoroside (C3S) increased 104.6-fold and 66.8-fold in FR vs. CK. Transcriptomics analyses revealed the expression levels of flavonoid 3′- hydroxylase (F3`H) and ﬂavonol synthase (FLS), which were related to an increase in quercetin derivatives. Six phytohormone-related GO terms were enriched during FT and FR. Four genes were significantly affected by FT and FR, including MYB97, MYB39, fructosyltransferase, and zinc finger transcript factor. Four GA response-related genes, 2-oxoglutarate-dependent dioxygenases, were specially induced during freezing stress. Moreover, jasmonic acid (JA) and salicylic acid (SA) were significantly increased during FT, and only SA decreased during FR. This study provided understanding of phenylpropanoid metabolism under freezing stress and will be useful for enhancing freezing tolerance via molecular breeding in D. officinale.


Introduction
Cold stress is a major environmental stress that limits yield and quality of crops, ornamental plants and herbs. Cold stress, which can be classified as chilling (0-5°C) and freezing (< 0°C) stresses, affects water and nutrient uptake, membrane rigidification, destabilizes protein activity, and impairs chloroplast structure (Li et al. 2008;Miura and Furumoto 2013;Hu et al. 2019). Freezing stress causes more serious injuries to tissue cells and leads to the formation of ice crystals in the extracellular space (Bredow and Walker 2017). Extracellular ice has a very low water potential and can reduce cell volume, leading to membrane damage (Nievola et al. 2017). Accumulation of reactive oxygen species (ROS) often occurs during the early phase of freezing injury. If not scavenged, ROS can cause drastic lipid changes, leading to perturbations in membrane structure and functioning (Min et al. 2014).
Phenylpropanoid compounds are a large class of secondary metabolites that play an important role in plant responses to biotic and abiotic stresses. Phenylpropanoids are derived from chalcones, which are converted to other flavonoid classes, including anthocyanins, flavones, and flavonols (Dixon and Paiva 1995). The accumulation of anthocyanins is taken as a general indicator of plant stress exposure (Schulz et al. 2016). Abiotic stresses strongly increases flavonoid content and the expression of flavonoid biosynthesis genes, which are directly associated with oxidative tolerance in various plant species (Nakabayashi et al. 2014;Ahmad et al. 2019;Zhan et al. 2019). Antioxidative activity assays in vitro also showed that glycosidic flavonoids had strong antioxidative activity (Nakabayashi et al. 2014). However, the protective mechanism of flavonoids during freezing was not only the scavenging of ROS. Reduced flavonoid content indeed impaired leaf freezing tolerance, but freezing treatment of plants was performed in the dark to exclude ROS formation in the chloroplast (Steponkus 1984). Thus, further research is needed to reveal the protective mechanisms that link flavonoid biosynthesis to plant stress tolerance.
D. officinale has been widely cultivated and is becoming an important economic crop for health care in China, but the yield of D. officinale is greatly affected by low temperatures under nature conditions (Wu et al. 2016). The study is necessary to determine the optimum temperature range for D. officinale growth and the temperature threshold for cold acclimation induction. The active medicinal components of D. officinale primarily include polysaccharides, alkaloids, terpenes and flavonoids (Jin et al. 2016). Our previous studies indicated that the primary flavonoids of D. officinale were delphinidin and quercetin derivatives (Zhan et al. 2020). Wu et al reported changes of metabolites and genes under cold stress (0°C) in D. officinale, but the response of secondary metabolism was not clear (Wu et al. 2016). Recently, expression profiling of the SET DOMAIN GROUP family and glycolipid synthase family were identified in abiotic stress Zhan et al. 2021). However, studies on abiotic stress in D. officinale are rarely reported, especially during freezing stress. In this study, metabolomic and transcriptome analyses are integrated to investigate the regulation of phenylpropanoid metabolism during moderate freezing (FT) and post-freezing recovery (FR). These results will facilitate the identification of coldresistant candidate genes in D. officinale and provide new insights into the molecular mechanisms related to freezing tolerance in plants.

Physiological differences during FT and FR in D. officinale
Freezing damage seriously affected the cultivation of D. officinale in winter. To investigate the overall effect of freezing damage on D. officinale, we monitored the responses to cold stress under a gradual temperature decrease ( Fig S1). Firstly, we compared the electrolyte leakage at different temperatures. The electrolyte leakage displayed an increase from 13.2% to 80.5%, and showed a 4.6-fold increase at −6°C, and subsequently kept stabilized ( Fig S2a). The electrolyte leakage significantly increased at 4°C compared with that at 8°C during three days ( Fig S2b). Moreover, the electrolyte leakage significantly decreased during FR compared with that during FT (Figure 1a). FT and FR significantly increased anthocyanin, ROS and soluble sugar contents compared with CK ( Figure 1b-d). ROS, as a cold tolerance indicator, showed a twofold increase at 3 h of FT compared with 0 h and subsequently remained stable. Considering the above results, we decided to take −6°C for 3 h as FT and 8°C for 24 h as FR for transcriptomic and metabolomic analyses.
Overview of the metabolic profiles in response to CA, FT and FR To explore the metabolic response to the cold stimulus in D. officinale, we examined plant metabolic profiles under CK, FT, and FR using nontargeted metabolomics as described previously (Zhan et al. 2020). After quality validation, 9,022 potential metabolites with annotations were identified in CK, FT, and FR (Table S1). Principal component analysis (PCA) was used to analyze the rate of contribution of the first two primary components, which were 44.71% and 17.48%, respectively, the three treatments were obviously separated and formed a cluster (Figure 2a). Hierarchical cluster analysis showed that the CK, FT, and FR metabolite profiles of the six independent biological replicates clustered together ( Figure 2b). Furthermore, a total of 1,952 metabolites were mapped to the KEGG Pathway. The largest number of metabolites belonged to the categories 'biosynthesis of other secondary metabolites' (251 metabolites), followed by 'metabolism of terpenoids and polyketides' (151 metabolites), 'amino acid metabolism' (192 metabolites) and 'carbohydrate metabolism' (230 metabolites) (Figure 2c, Table S2).
KEGG pathway enrichment analysis of differentially accumulated metabolites (DAMs) KEGG pathway enrichment analysis was further performed on the DAMs during FT and FR (Figure 3a). Among these pathways, five pathways were significantly enriched (p < 0.05) in the two comparisons, including 'porphyrin and chlorophyll metabolism,' 'flavonoid biosynthesis,' 'flavone and flavonol biosynthesis,' 'phenylpropanoid metabolism,' and 'anthocyanin biosynthesis. ' We performed an HCA heat map to investigate distinct clusters of phenylpropanoid profiles during FT and FR ( Figure 3b and Table S3). Nine quercetin and cyanidin derivatives were significantly increased in FR vs. CK, including quercetin 3-O-glucoside (104.6-fold), cyanidin 3-O-sophoroside (66.8-fold), leucodelphinidin (4.3-fold), and dihydroquercetin (21.6-fold). Several lignins and phenylpropanoids were significantly decreased in FT vs. CK comparison, such as sinapoyl_aldehyde (4.7-fold), naringenin (5.3-fold), 4-coumarate (9.0fold). Substantial changes at FT and FR were also seen for compounds connected to chlorophyl degradation ( Figure  3c), including pheophorbide_a approximately increased 7.0-fold and pheophytin_a approximately decreased 7.0fold in the two comparisons. In addition, citrate was increased 2.4-fold in FT vs. CK and isoprene was increased 5.6-fold in FR vs. CK. TCA intermediates as injury-related metabolites were found to accumulate during short-term (3 h) freezing stress (Min et al. 2020). Low temperature also affected photosystem II activity and flavonols could alleviate it (Mattila et al. 2020). These results suggested that plants tried to relieve freezing damage by synthesizing a broad collection of metabolites.

Differentially expressed gene (DEG) analysis in CK, CA, FT and FR
To explore the gene expression patterns at CK, FT, and FR, RNA-seq analysis was performed with three independent replicates of each treatment. Using the paired-end HiSeq4000 sequencing platform, 60.25 Gb of sequence data were yielded. The GC percentages ranged from 45% to 46%. The percentage of valid reads ranged from 92.94% to 95.47%. The probability of incorrect base calling was used to evaluate the sequencing quality by the value of Q30, and a high proportion of Q30 (97%) for each library showed a high quality of RNA-seq (Table S4). To obtain an overview of the transcriptomic variation, PCA was performed, and the explained values of PC1 and PC2 were 42.96 and 15.72%, respectively ( Fig S3).
Pairwise comparisons of CK, FT, and FR, respectively identified 8,561 DEGs in FT vs. CK, 8,645 DEGs in FR vs. CK, and 1,812 DEGs in FR vs. FT, with a significant difference in their expression levels with the threshold of FDR (false discovery rate) ≤ 0.05 ( Figure 4a). The numbers of DEGs identified in each comparison are shown in a Venn diagram (Figure 4b). A total of 6,002 DEGs were identified in FT vs. CK and FR vs. CK. Based on hierarchical clustering, the expression levels of the DEGs were divided into two clusters ( Figure 4c). Furthermore, 66 GO terms and 63 GO terms were significantly enriched (adjusted p-value < 0.05) in FT vs. CK and FR vs. CK, respectively (Table S5). The majority of the DEGs were significantly enriched in secondary metabolic processes, photosynthesis processes, cold response processes, and phytohormone response processes. For instance, the top of 'biological process' GO subcategories (adjusted pvalue < 0.001) were related to photosynthesis in the two comparisons ( Figure 4d). However, carbohydrate metabolism was only significantly enriched in FT vs. CK, and 'ion transmembrane transport' was only significantly enriched in FR vs. CK.

Combined transcriptome and metabolome analysis revealed the phenylpropanoid metabolic pathways in CA, FT and FR
To further identify the functions of the DEGs and their related biological processes, KEGG pathway enrichment analysis was performed. 20 pathways were significantly enriched (p < 0.05) in FT vs. CK and FR vs. CK (Figure 5a and b). 16 pathways were significantly enriched in the two comparisons, such as 'glycosyltransferases,' 'cytochrome P450,' 'transporters,' and 'flavonoid biosynthesis' (Fig S4). We further identified the numbers of phenylpropanoid pathway-related DEGs and 38 DEGs were divided into two clusters (Fig S5 and Table S6). However, the phenylpropanoid pathway is a sophisticated metabolic pathway, and two important branches were analyzed in our present study: the lignin and flavonoid metabolic pathways. The coexpression network analysis of DEGs and DAMs of phenylpropanoid pathways was conducted (pearson correlation coefficient > 0.8 or < − 0.8, p value < 0.05; Table  S7). The visualized network of phenylpropanoid biosynthesis showed that 40 pairs had a positive correlation, and 20 pairs were negatively correlated. The visualized network of flavonoid biosynthesis showed that 8 pairs had a positive correlation, and 3 pairs were negatively correlated. (Figure 6a).
To provide a systematic view on the variation of phenylpropanoid biosynthesis pathways, the metabolome and transcriptome data were integrated according to correlation (Figure 6b). Phenylpropanoid biosynthesis was initiated from phenylalanine, cinnamic acid and 4-coumarate, which were all decreased during FT. In the lignin biosynthesis pathway, caffeate, coniferyl alcohol, and sinapoyl aldehyde decreased, while ferulate increased during FT and FR. Most identified genes related to lignin biosynthesis showed a reduction during freezing stress. Among them, three peroxidase (POD) genes significantly decreased in three comparisons, including POD (LOC110112233) decreased 700.2fold and 500.1-fold in FT and FR compared to that in CK, respectively. Lignan is also an important active compound in herbs, which shares a common upstream pathway with lignin and branches after the synthesis of coniferyl alcohol, so the medicinal value of D. officinale requires further exploration. In the flavonoid biosynthesis pathway, the metabolites showed a significantly improved production during FT and FR. In particular, quercetin 3-O-glucoside (Q3G) and cyanidin 3-O-sophoroside (C3S) increased 104.6-fold and 66.8-fold in FR vs. CK, respectively. The increased levels of dihydroquercetin and quercetin were associated with significant increases in mRNA, including F3`H (LOC110115941) and FLS (LOC110097028) increased 13.7-fold and 34.5-fold in FR vs. CK, respectively. The results showed that the cyanidin and quercetin derivatives principally accumulated during FR.
To test the precision of the RNA-Seq data, RT-PCR was performed on 8 DEGs selected from phenylpropanoid biosynthesis, including HCT, CCOMT, FLS, CYP84A1, POD, CAD, CHI, and F3`H. All 8 selected genes showed similar expression patterns in RNA-seq data (Fig S6), which indicated that the RNA-seq data were reliable.

Phytohormone signaling associated with freezing stress
Six phytohormone-related GO terms were selected from Table S5 (Figure 7a). 'Response to jasmonic acid,' 'response to salicylic acid,' and 'response to gibberellin' were significantly enriched (p < 0.05) in FT vs. CK and FR vs. CK. 'Response to cytokinin' and 'response to auxin' were significantly enriched in FR vs. CK, while 'response to brassinosteroid' was significantly enriched in FT vs. CK. We further analyzed the expression of enriched DEGs in the two comparisons and the top 10 of DEGs were shown in the volcano plot (Figure 7b-c and Table S8). MYB39, FT and ZC3H were significantly down-regulated in the two comparisons, while MYB97 was significantly up-regulated in the two comparisons. However, several DEGs were associated with specific cold treatment, such as EGL1 belonging to the bHLH family that was up-regulated in FT vs. CK, and UROD and PRS1 were up-regulated in FR vs. CK. It suggested these DEGs might play a crucial role in the response mechanism of cold stress.
Jasmonic acid (JA) interacts salicylic acid (SA) to regulate leaf senescence and tolerance to cold stress (Hu et al. 2017). JA and its precursor (linoleate) levels were significantly upregulated during FT and FR, while SA levels were significantly up-regulated during FT (Figure 7d). SA accumulation could indicate the degree of injury during freezing stress (Min et al. 2020). However, SA levels decreased during FR and remained at the same level in CK. These results suggested that JA and SA played different roles during freezing stress.

Discussion
Freezing injury is one of the most universal abiotic stresses affecting plant development during winter. The accumulation of intercellular ice could physically damage the cell membrane, causing a decrease in the water potential in the extracellular space and leading to severe cell dehydration (Steponkus 1984). The plants regulate their osmotic potential through the accumulation of solutes in the interior of their cells during cold stress, such as soluble sugars (Zuther . Soluble sugars, including glucose, fructose, and sucrose, significantly accumulated under cold treatment (at 0°C during 20 h) in D. officinale (Wu et al. 2016). Moreover, soluble sugars are generally considered to be an antioxidant to defend against anoxic injury (Couée et al. 2006). ROS accumulation is directly associated with sugar accumulation to adapt to the ill-effects of environmental stress (Roitsch 1999). In our study, soluble sugars and ROS levels were significantly induced during FT and FR (Figure 1). Glucose and sucrose levels were significantly increased during cold stress, but sucrose levels were decreased during recovery (Fig S7). Reductions in sucrose levels were also found during recovery after freezing treatment in cold acclimated leaves (Vyse et al. 2020). The phenomenon might be caused by a higher demand for carbohydrates for respiration during recovery.
Previous work had implicated the role of supercoolingpromoting substances in the defence against freezing stress, including antifreeze protein, phenylpropanoid, flavonoid, terpenoid, tannin, phenol (Fujikawa et al. 2018). In the present study, most flavonoids and anthocyanins were induced during FR. Among them, Q3G and C3S accumulated 104.6fold and 66.8-fold in FR compared to CK, respectively (Figure 3b). Q3G as a novel anti-ice nucleation substance had been identified from deep supercooling tree xylem cell (Kasuga et al. 2008). However, not much work has been done to report cyanidin glycosides as supercooling-promoting substances. Due to the antioxidant activity of flavonoids, the role of flavonoids during freezing stress might be in the scavenging of ROS (Wang et al. 2006;Schulz et al. 2016). Q3G and cyanidin 3-O-glucoside (C3G) exhibited higher ROS-scavenging activity than kaempferol 3-Oglucoside in vitro (Nakabayashi et al. 2014). We found total anthocyanin content and ROS levels significantly increased during FT and FR (Figure 1d and e). Metabolome data also showed that Q3G and cyanidin glycoside levels significantly increased during cold stress. These results suggested that flavonoid biosynthesis was associated with the ROS accumulation. However, we found that the accumulation of Q3G and cyanidin glycosides in FT vs. CK were lower than in FR vs. CK (Figure 6b). Previous work indicated that freezing stress of plants was performed in the dark to exclude ROS formation in the chloroplast, and a lower light intensity also reduced ROS accumulation during cold acclimation (Schulz et al. 2016). It suggested flavonoid accumulation was not only related to the scavenging of ROS, but also to enhanced freezing tolerance in plants. The role of flavonoids during freezing stress deserves to be studied further.
The biosynthesis of flavonoids via the phenylpropanoid pathway is controlled by a few key enzymes. In the present study, no obvious increases in early phenylpropanoid biosynthetic genes were found in FT vs. CK and FR vs. CK. Integrated analyses showed that increased levels of dihydroquercetin and quercetin were associated with significant increases in mRNA of F3`H and FLS during cold stress (Figure 6b). F3`H and FLS played important roles in the response to cold stress. For instance, the expression profiling of F3`H from Antarctic moss was influenced by diverse abiotic stresses including cold stress (Liu et al. 2014). A tomato F3H-like gene stimulated flavonoid biosynthesis in response to chilling stress, and FLS levels were higher in the overexpression plants than in the WT plants (Meng et al. 2015). Furthermore, lots of DEGs were classified into three kinds of GO terms, including photosynthesis, carbohydrate and cell wall metabolism (Figure 4d). In Arabidopsis, freezethaw cycles caused photosynthesis and cell wall lesions that led to tissue death (Vyse et al. 2020). The flow of individual metabolites with mRNA, including caffeate, coniferyl alcohol, sinapoyl alcohol, POD, HCT, F5H, and CCOMT, suggested that lignin precursors were significantly downregulated during cold treatment (Figure 6b). Lignin is a component of the plant cell wall, and it is a branched polymer of phenylpropanoid compounds. Lignification and changes in lignin content are complex during cold stress. Lignin content increased at low temperatures and participated in cell wall modification by strengthening it (Le Gall et al. 2015). However, some studies indicated that no change or decrease in the levels of lignin and its precursors were observed in plants during low temperatures (Moura et al. 2010). Cold acclimation (8°C) could alleviate chilling injury (0°C) in peach fruit by enhancing activity of cell wall hydrolases and reducing transcript levels of POD . Low temperature (4°C) also suppressed the activities of lignin biosynthesis enzymes, including PAL, CAD, and POD (Lwin et al. 2020). The ternary complex of bHLH-MYB-AP2 repressed the transcription of 4CL to reduce lignification of loquat fruit at 0°C (Xu et al. 2019). These findings suggested that the regulation of phenylpropanoid metabolism was effective in enhancing freezing tolerance and cold acclimation in plants. shows the numbers of the DEGs in the two comparisons. c. Hierarchical cluster analysis of 6,002 DEGs. d. The Go enrichment analysis of DEGs in the two comparisons. The x-axis represents the fold enrichment of GO, while the y-axis represents -log 10 adjusted p-value. The dot sizes represent the number of differentially enriched genes. The 'biological process' GO subcategory is labeled with a cut-off (adjusted p-value < 0.001).
Plant hormones play critical roles in plant adaptation to biotic and abiotic stresses. GO enrichment analysis showed that plenty of DEGs were significantly enriched in the JA and SA (Figure 7a). JA plays a positive role in modulating tolerance to cold stress. Blocking JA biosynthesis and signaling pathways resulted in hypersensitivity to freezing stress (Hu et al. 2017). Exogenous methyl jasmonate (MeJA) treatment had been shown to increase endogenous JA levels and cold tolerance in plants (Zhao et al. 2013;Kazan 2015). Our results also showed that JA and its precursor (linoleate) levels significantly increased during freezing stress (Figure 7f). The SA pathway is better known as a response to plant defense. Short-term cold stress activated disease resistance through activation of the SA pathway (Wu et al. 2019). SA accumulation in response to freezing stress showed a biphasic response in spinach leaves, increasing first and then decreasing, and again increasing (Min et al. 2020). SA levels increase in Arabidopsis that were exposed to low temperatures for one week. However, previous research found that SA biosynthesis did not contribute to freezing tolerance, but had a major role in the induction of 'defense response' genes (Kim et al. 2013). Our data revealed SA levels significantly increased in FT and then decreased in FR (Figure 7f). These suggested that SA played an important role during short-term freezing stress. In addition, gibberellin (GA) pathways also were enriched in the two comparisons ( Figure  7a). GA regulates various developmental processes in plants, including involvement in plant responses to cold. GAdeficient mutant (ga1) showed increased freezing tolerance, while DELLA knockout mutant (gai-t6 rga-24) was hypersensitive to freezing stress (Achard et al. 2008;Richter et al. 2010). Previous study showed that DELLA proteins fine-tuned JA signaling pathway through interaction with JA ZIM-domain 1 (JAZ1) protein (Hou et al. 2010). JA delayed GA-mediated DELLA protein degradation, and the della mutant was less sensitive to JA for growth inhibition (Yang et al. 2012). However, the mechanism of JA and GA crosstalk remains unknown under freezing stress. We further selected the DEGs from the 'response to gibberellin,' and found that most DEGs were down-regulated during cold stress (Fig S8). Out of these, seven 2-oxoglutarate-dependent dioxygenases (2OGD) genes were identified and five 2OGDs were specially induced during freezing stress. Based on functional annotations in the NCBI database, four 2OGDs were potentially involved in the biosynthesis of GA (Fig S9). Members of the 2OGD superfamily were involved in various metabolic processes, such as phytohormone biosynthesis (Kawai et al. 2014). The tomato 2OGD gene (SlF3HL) had been reported to positively regulate chilling stress tolerance by modulating JA biosynthesis (Hu et al. 2019). However, AtDMR6, which had its highest sequence homology with SlF3HL, was involved in the hydrolysis of SA (Zeilmaker et al. 2015). Thus, the identification of 2OGD and phytohormone crosstalk contributes to dissecting the molecular mechanisms involved in tolerance to freezing stress.

Plant materials and freezing treatment
Two-year-old cultivated D. officinale were collected from a growth chamber at Zhejiang University, Hangzhou, China. The growth conditions were set at 25 ± 2°C with a light/ dark cycle of 12/12 h, fluorescent lighting of 80 μmol photons m −2 s −1 , and 65%-75% relative humidity as a control treatment (Zhan et al. 2020). D. officinale is an epiphytic plant. To ensure air permeability, the compositions of soil for D. officinale were mainly pine barks, ceramic pieces, perlites, and small amounts of nutritious soil. Freezing treatment followed the protocol described by Li et al. (2008), but without lighting. Briefly, plants were placed at 0°C and then subjected to a gradual drop from 0 to −6°C within 3 h. The ice chips were placed onto the soil to induce crystallization and prevent supercooling. The temperature was held at −6°C for 0, 3, 9, 16, and 24 h. For post-freezing recovery, the temperature was gradually raised to 8°C within 12 h without lighting and then held at 8°C for 12 h under a light of 30 μmol photons m −2 s −1 before sampling. All samples were immediately frozen in liquid nitrogen and stored at −80°C.
Measurement of soluble sugar, ROS, anthocyanin content, and electrolyte leakage Soluble sugar content was measured using an anthrone-sulfuric acid method with glucose as the standard (Zhan et al. 2017). Briefly, pre-weighed leaves were ground in liquid nitrogen and then incubated in water at 95°C for 30 min. After centrifugation at 12,000 rpm, 0.8 ml of supernatant was incubated with 0.2 ml anthrone-sulfuric acid (0.3 g/ml) at room temperature for 1 min, and the absorbance was recorded at 620 nm. Soluble sugars were separated by thinlayer chromatography [developing solvent composed of chloroform: methanol (60: 40, by volume) to measure sucrose and glucose. The sugars were visualized by α-naphthol and measured by the anthrone-sulfuric acid method with glucose as the standard. We used NBT (nitro blue tetrazolium, Sangon, Shanghai, China) to detect ROS production according to a method described previously (Zhan et al. 2019). Anthocyanin content was determined according to the methods described previously (Zhan et al. 2020). Electrolyte leakage measurements were described previously (Wu et al. 2016;Vyse et al. 2020). In brief, leaf tissues were incubated in 10 mL of deionized water and shaken at 25°C for 2 h. The conductivity was measured at first using a conductively meter (DDSJ308A, Shanghai, China). Afterward, the samples were boiled for 30 min and then left to cool to room temperature and carried out a second measurement. The ratio between the two values was calculated and control values were subtracted to obtain the final electrolyte leakage values.

Nontargeted metabolomic analyses
The leaf samples were harvested from different cold treatments for metabolome analysis (six biological replicates). The tissue samples were ground into a fine powder in liquid nitrogen, and then extracted with precooled methanol. Metabolite profiling was conducted using UPLC-ESI-MS/ MS. For nontargeted metabolomic analyses, metabolite extraction, UPLC-ESI-MS/MS analysis and bioinformatic analysis were performed according to our previous work (Zhan et al. 2020). In brief, identified metabolites were annotated using the KEGG Compound database (http://www. kegg.jp/kegg/compound/), annotated metabolites were then mapped to the KEGG Pathway database (http://www.kegg. jp/kegg/pathway.html). For DAMs selection, significantly regulated metabolites between groups were determined by log 2 fold-change (FC) > 1 or log 2 FC < −1 and with statistical significance (p value < 0.05). The DAMs were then enriched in the KEGG pathway using a MBROLE website (http://csbg. cnb.csic.es/mbrole2/).

Transcriptomic analyses
Total RNA samples of 10 μg from leaves (three biological replicates) were prepared. The methods of RNA quantity, cDNA library preparation, and transcriptomic analysis were the same as our previously published work (Zhan et al. 2020). Expression levels for each unigene were calculated as the FPKM (fragments per kilo bases of exons for per million mapped reads) using the software package StringTie (Pertea et al. 2015). The DEGs were screened with criteria: FDR ≤ 0.05, log 2 fold-change (FC) > 1 or log 2 FC < −1 and with statistical significance (p value < 0.05). The DEGs were also subjected to Gene Ontology (GO) enrichment analysis using the GOseq software and KEGG pathway enrichment analysis by KOBAS 2.0 (Young et al. 2010;Xie et al. 2011). The RNA-seq data have been submitted to the BIG Data Center of the Chinese Academy of Sciences (http://bigd.big.ac.cn) with accession number CRA003229 and CRA005177.

Quantitative real-time PCR analysis
Total RNA was isolated from leaf samples using TransZol reagent (TransGen Biotech, Beijing). First-strand cDNA was reverse transcribed using the TIANscript RT Kit according to the manufacturer's instructions (TransGen Biotech, Beijing). Quantitative real-time PCR analyses were performed using a SYBRGreen qPCR kit (TransGenBiotech) with a MyiQ system (Bio-Rad) as described previously (Zhan et al. 2017). The primers for amplification are listed in Table S9.