Tg and Sg C. japonica were selected based on latitude (Figure 1A), inferring differences in cold sensitivity. The local environmental temperature of Tg plants is lower than that of Sg plants, indicating that Tg might harbor cold-tolerant traits. Here, the two groups of plants were exposed to the same cold treatment of -5°C for 0, 1, 2, 3, 4, 5, 6, 7, 8, 12, and 24 h. Sg C. japonica suffered more damage than Tg C. japonica (Figure 1B and Supplementary Figure 1), based on malondialdehyde (MDA) and relative electrical conductivity measurements. Tg C. japonica had less osmotic and oxidative damage than Sg C. japonica (Figures 1C,D and Supplementary Table 1). Tg C. japonica had higher H₂O₂ catalase (CAT) activity than Sg C. japonica (Figure 1E). In contrast, a significant difference in the rate of change in electrical conductivity was detected between the two populations at 4 and 8 h (Figure 1F and Supplementary Table 1). The 4–24 h time interval revealed significant differences in MDA and CAT levels (Figures 1G,H). Thus, the transcriptome was sequenced using the 0 (CK), 4 (T1), 8 (T2), and 24 h (T3) cold treatments.

More than 168.23 Gb of clean data were collected from 24 samples with Q30 > 85.1% (Supplementary Table 2). The clean reads ranged from 20,404,743 to 31,731,765. Then, 49,394 unigenes were obtained using Trinity software. Furthermore, 5,855 proteins were quantified in the proteome. The Pearson’s correlation coefficients for the biological replicates of the different samples varied from 0.82 to 0.87 (Supplementary Table 2). These results indicate that the throughput and sequencing quality was sufficient to warrant further analysis.

We focused on the differences in transcriptome and proteome expression between the normal temperature treatments (control and unstressed). This approach further highlighted the significant differences in basic transcription. The Sg population was taken as the control. The screening condition was set to fold-change ≥ 2 and a false discovery rate (FDR) < 0.01 to identify differentially expressed genes (DEGs). In the Sg vs. Tg comparison, 816 genes were identified as DEGs; 384 genes were upregulated and 432 genes were downregulated. The screening conditions were set to fold change ≥ 1.5 and an FDR < 0.01 to identify differentially expressed proteins (DEPs). As results, 149 proteins were upregulated and 265 proteins were downregulated (Figure 2A).

The DEGs and DEPs were functionally enriched. The photosynthesis, ribosome-related genes, and proteins exhibited higher expression levels in Tg than in Sg C japonica (Figure 2B and Supplementary Table 3). Plant hormone signal transduction, plant-pathogen interactions, and biosynthesis of secondary metabolite pathways, such as flavonoid and anthocyanin biosynthesis-related genes and proteins, were expressed at higher levels in Sg than Tg C. japonica (Figure 2C and Supplementary Table 3). In addition, Gene Ontology (GO) enrichment analysis was performed, and terms such as photosynthesis-related “photosystem I and II,” “integral component of the membrane,” and “nucleic acid binding” were enriched in Tg (Figure 2D and Supplementary Table 4). Similarly, the regulation of DNA endoreduplication, serine/threonine kinase activity, and terpene synthase activity were enriched in Sg C. japonica (Figure 2E and Supplementary Table 4).

The transcriptomic data during the cold treatment showed that Tg and Sg C. japonica displayed similar transcriptional dynamics. Interestingly, in the T1 and T3 treatment groups, the number of DEGs in Tg C. japonica was lower than in the Sg. In the T2 treatment group, both Sg and Tg displayed the largest number of DEGs, as 7,114 and 7,742 genes were identified as DEGs, respectively (Figure 3A). Furthermore, the iTRAQ proteome of the T2 treatment group was identified. Relative to the control, 55 upregulated proteins and 27 downregulated proteins were detected in Tg C. japonica treated at −5°C for 8 h. Seventy-three upregulated proteins and 49 downregulated proteins were detected in Sg C. japonica during the same treatment. KEGG functional enrichment analysis showed that these DEGs and DEPs were enriched in the ribosome, photosynthesis-related, protein processing in the endoplasmic reticulum, plant hormone signal transduction, glycerophospholipid metabolism, and secondary metabolic pathways (Supplementary Table 5). In addition, these data suggest a broad restructuring of the proteome as a common response to cold stress.

We focused on a pairwise comparative analysis of the two populations under the same cold treatment to explore the differences in cold responsiveness. The pairwise comparative results suggested that 816, 694, 830, and 1,115 DEGs were identified at the CK, T1, T2, and T3 treatments, respectively (Figure 3B), and displayed larger differences between the populations than the control. These results indicate that marked interspecific variation and basal transcription were more conserved among the tested accessions than cold-related transcription. A total of 166 genes were shared by the two populations during cold stress, and 431 genes were specific to the T2 group (Supplementary Figure 2). The KEGG pathway enrichment analysis revealed numerous pathways related to carbon metabolism, biosynthesis of secondary metabolites (flavonoids and phenylpropanoids), and protein synthesis processing that were common and upregulated in Tg C. japonica, including protein processing in the endoplasmic reticulum, ribosome, and ribosome biogenesis, and photosynthesis and phenylpropanoid metabolism (Figure 3C and Supplementary Table 6).

An expression pattern analysis was performed across all treatments and populations using MaSigPro software (Figure 4). A total of 4,366 genes were divided into 15 clusters (R² ≥ 0.6, adjusted p < 0.05, Figure 4). Five of the 15 clusters included population-associated differences in the intensity of expression for 1,374 (31%) genes, independent of cold stress time. Cohort 1 included 2,638 genes with similar expression levels under cold stress, while cohort 2 included 1,515 (34%) genes with distinctly higher or lower expression in only a single population. Cluster 5 was the only cluster with higher Tg than Sg expression levels. The expression patterns of 213 (4%) genes were different or opposite (cohort 3), indicating that the expression differences could be transmitted across the domestication gradient.

To obtain the hub components associated with the differences in cold responsiveness, 3,119 genes with the largest contribution to 25% of the variance were screened from 12,475 DEGs for the weighted gene co-expression network analysis (WGCNA). The Power 8 network topology was selected to construct a hierarchical cluster tree. All genes were clustered into 12 modules (labeled with different colors, Figure 5A). All of the modules were highly correlated with the respective module. We associated each co-expression module with a cold treatment time in the two populations using Pearson’s correlation coefficient analysis (Figure 5B). The turquoise module was correlated with the control groups of both populations, but the black and red modules were correlated with the Sg-T2 and Sg-T3 groups, respectively (r = 0.59, p = 2e −3 and r = 0.65, p = 6e −4). The correlation between the pink module and the Tg-T2 group was the highest (r = 0.57, p = 3e −3), and the brown modules were correlated with the Tg-T1 group (r = 0.5, p = 3e −3). The transcriptional difference caused by the cold treatment was greater than that caused by the different populations. The difference among the populations under the cold treatment was greater than that under normal growth, indicating that C. japonica (Naidong) formed conservative and divergent evolution caused by environmental differences, and Tg C. japonica was more radical against cold stress.

GO and KEGG analyses were performed for the gene modules mirroring the divergence in cold responsiveness. The red and black modules contained the specific genetic information for the response to cold stress in Sg C. japonica, and the brown and pink modules represented that for Tg C. japonica. The first 30 enriched GO terms (p < 0.05) showed that terms related to protein processes were enriched in the black and red modules, e.g., “response to unfolded protein,” “protein-folding chaperone,” and “protein storage vacuole.” In addition, many terms were related to the chloroplast, e.g., “chloroplast envelope” and “chloroplast.” Interestingly, the degree of enrichment of terms related to the response to a stimulus was higher in the brown and pink modules than in the red and black modules, e.g., “response to cold,” “response to jasmonic acid,” and “regulation of defense response.” In addition, “plasmodesma” and “DNA-binding transcription factor activity” were also enriched (Figure 5C and Supplementary Table 7). The KEGG enrichment analysis (p < 0.05) suggested that the ribosome, plant hormone signal transduction, photosynthesis-related, glycerophospholipid metabolism, and plant-pathogen interaction pathways were significantly enriched in Tg C. japonica, while protein processing in the endoplasmic reticulum and biosynthesis of secondary metabolites were significantly enriched in Sg C. japonica (Figure 5D and Supplementary Table 8), suggesting that the Tg and Sg proteomes’ response to cold stress is broad. Overall, these findings are consistent with the results obtained, demonstrating that our results are reliable.

Key candidate genes were screened from the transcriptome data based on standard gene names in the functional annotation of the pathways. Sixty-two and 64 DEGs were involved in the photosynthetic response to cold in Sg and Tg C. japonica, respectively (Supplementary Table 9). Thirty-two genes displayed different expression levels in the Sg vs. Tg comparison. Among these genes, six were specific to Tg, including CjFD3 (c153759.graph_c0), CjPSBR (c129109.graph_c0), CjFNR2 (c143865.graph_c0), CjATPC1 (c146541.graph_c0), CjPetA (c148557.graph_c2), and CjCF0 (c105617.graph_c0). Furthermore, 15 genes involved in carbon fixation during photosynthesis were identified as DEGs, e.g., glyceraldehyde-3-phosphate dehydrogenase (c117808.graph_c0) and ribulose bisphosphate carboxylase small chain (c122766.graph_c0) (Supplementary Table 10). To confirm the reliability of the transcriptome data, we selected five genes from these DEGs and validated them using quantitative real-time polymerase chain reaction (RT-qPCR). The expression trends in the RT-qPCR results were largely consistent with the transcriptome data (Supplementary Figure 3), suggesting that our data were reliable. These results indicate that these genes play an important role in photosynthesis regulated by local cold adaptation.

The differences in the transcriptomic and proteomic data were larger at the transcript level than at the protein level. Four photosynthesis-antenna proteins and five proteins involved in carbon fixation in photosynthetic organisms were identified as DEPs between Tg and Sg C. japonica under cold stress, e.g., glyceraldehyde-3-phosphate dehydrogenase B, transketolase, chloroplastic malate dehydrogenase, glyoxysomal, malate dehydrogenase, NADP-dependent malic enzyme-like, phosphoglycerate kinase, cytosolic, glyceraldehyde-3-phosphate dehydrogenase, photosystem I chlorophyll a/b-binding protein 3-1, and chlorophyll a-b binding protein CP29.1 (Supplementary Table 9). These results between Tg and Sg C. japonica at the translational level provide a stronger photosynthetic response to cold in Tg C. japonica than Sg C. japonica.

We focused on the divergence in ribosomes between Tg and Sg C. japonica during the cold treatment. Eighty-five genes displayed significantly different expression levels in response to cold stress between the two populations, including encoding large subunit ribosomal protein and small subunit ribosomal protein genes (Supplementary Table 11). Most of these genes were expressed at higher levels in Tg than Sg C. japonica; 56 DEPs involved in divergence in the ribosomal pathway between Tg and Sg C. japonica were screened in the proteome, e.g., 60S ribosomal protein L7a-1-like (c26982_f1p21_1112), 40S ribosomal protein S7 (c35197_f1p2_1028), 30S ribosomal protein S3 (chloroplast) (c16401_f1p4_4386), 50S ribosomal protein L27, chloroplastic (c74584.1_c1_orf5), and ribosomal protein S12 (chloroplast) (c75116.1_c0_orf7) (Supplementary Table 6). Consistent with the transcriptome results, only 30S ribosomal protein S5 and chloroplastic-like (c31055_f2p1_1207) were expressed at lower levels in Tg than Sg C. japonica (Supplementary Table 11). Taken together, the ribosome response to cold was stronger in Tg than Sg C. japonica, and protein processes played an important role.

The hormone network plays a key role in the plant’s response to stress. We carried out a targeted hormone quantification analysis specifically targeting gibberellic acid (GA), ABA, auxin, BRs, and zeatin (ZR) (Figures 6A–F). The contents of these hormones were higher in Sg than in Tg C. japonica. GA3 was downregulated in the two C. japonica populations under cold stress, whereas GA4 was upregulated initially and then downregulated in Tg under cold stress, but displayed the opposite trend (downregulated first and then upregulated) in Sg C. japonica (Figures 6A,B). The indole acetic acid (IAA) levels were similar. Surprisingly, ABA was downregulated first and then upregulated in leaves in the 24-h -5°C treatment in both populations. Interestingly, the ABA/IAA ratios decreased in Sg C. japonica under cold stress, and ABA/GAs ratios decreased in Tg C. japonica under cold stress (Figures 6G,H), while the ABA/GA ratio did not change in Sg C. japonica, suggesting that Tg C. japonica has higher adaptability to the cold (Supplementary Table 12).

Environmentally related differences in leaf morphology generally reflect the adaptive evolution of plants (Ordonez et al., 2009). The color and volume of leaves often mutate with a change caused by environmental stress (Givnish, 1978; Wright et al., 2005). No differences in leaf morphological traits were observed between the two populations of C. japonica under normal conditions. We cannot be sure whether this was affected by natural domestication because of the complexity of the climate conditions. Combining different types of data makes it possible to contrast the results and refine current hypotheses to provide a better understanding of the differences in cold tolerance.

As we speculated, the Tg plants suffered less osmotic and oxidative damage under the −5°C cold treatment than Sg plants. This finding is supported by the CAT activity results. The difference in latitude caused a large temperature difference, which promoted the difference in cold tolerance. The temperature is < 0°C for most of the winter in Qingdao, Shandong, and > 0°C for most of the winter in Zhoushan, Zhejiang (Supplementary Table 11). In our preliminary experiment, the same index was measured during the 4 and 0°C treatments, and no difference was observed between the two populations, indicating that plants evolve different cold-responsive mechanisms when populations grow under different cold stress intensities. Gene expression must change to modify physiological status and adapt to local environmental change; however, little is known regarding the molecular mechanism. Differences between populations of forsythia and Coffea mauritiana have been explored by transcriptome analyses (Garot et al., 2021; Li et al., 2021)Garot.

Variations in the transcriptome mirror genetic variations among species (von Wettberg et al., 2018). High genetic diversity has been reported in C. japonica (Li et al., 2013). Our results revealed significant differences in the transcript levels between the Tg and Sg populations under an unstressed (normal) condition. Genes related to photosynthesis contributed significantly to the changes in the transcriptome, with higher expression levels in the Tg than the Sg population, including photosystem II (PS II), photosystem I (PS I), ATP synthase, and chlorophyll a/b-binding proteins. The proteomic results were similar. These differences could reflect differences in light utilization efficiency or growth ability between Tg and Sg C. japonica. A study on differences in cold responsiveness in forsythia confirmed that the high latitude population has higher carboxykinase ATP activity (Li et al., 2021). The expression level of genes related to hormone signaling was lower in Tg than Sg C. japonica. Our hormone profile confirmed this result. The BR, ABA, ZR, and GA levels were lower in Tg C. japonica than those in Sg. Interestingly, the IAA level in Tg C. japonica was higher, and some genes related to the biosynthesis of auxins were upregulated. ABA and GA play opposite roles in regulating growth (Shi et al., 2021). Overall, Tg C. japonica displayed stronger viability than Sg C. japonica, which may reflect differences in cold tolerance affected by basic transcription.

Despite differences in cold resistance and gene expression levels, both C. japonica populations displayed their largest transcriptional responses after 8 h of cold treatment. Combining functional enrichment of the transcriptome and proteome suggested that a range of conserved pathways responded to the cold. Cold stress has negative impacts on photosynthesis (Batista-Santos et al., 2011). We identified key responders that may account for the changes in photosystem-related transcription under cold stress, such as oxygen-evolving enhancer protein (OEE)-like and chlorophyll-a apoprotein P700, which were upregulated. OEE1 specifically stimulates chloroplast fructose bisphosphatase, promotes carbon fixation (Charles and Halliwell, 1981; Heide et al., 2004), and protects the light-harvesting chlorophyll-binding protein against proteolysis depending on the chlorophyll-a apoprotein P700 (Eichacker et al., 1996). Furthermore, ferredoxin was upregulated in our results. Overexpression of ferredoxin-like protein in Oryza sativa promotes abiotic stress tolerance by reducing the accumulation of reactive oxygen species (Huang et al., 2020). These results further support the conclusion that photosynthesis is a core C. japonica response to cold. In addition, most genes related to the carbon fixation pathway were upregulated, while PEPC was downregulated in Tg and Sg plants. Our previous results suggest that these genes were upregulated and played an important role in C. japonica at 4°C (Fan et al., 2021a), indicating that there is divergence in the regulatory mechanisms among the responses to different cold stressors. Overexpression of the PEPC gene increases the net photosynthetic rate and reduces O₂-associated inhibition of photosynthesis under warm conditions (Zhang et al., 2014). Some studies have shown that upregulating PEPC may be important during cold acclimation (not freezing temperatures) (Sosińska and Kacperska-Palacz, 1979; Vu et al., 1995).

Circadian rhythms merit attention, including the input, central oscillator, and output pathways (Srivastava et al., 2019). The circadian system perceives temperature cues and regulates the interaction between temperature and the internal system of the plant (Srivastava et al., 2019). Most of the key members of the RESPONSE REGULATOR were downregulated in Tg and Sg C. japonica, e.g., PRR7, LHY, and TOC1 (Supplementary Table 7). Some transcriptome analyses have shown that the circadian clock regulates the maximum gene response to abiotic stress in Arabidopsis (Covington et al., 2008; Hubbard et al., 2009), rice (Duan et al., 2014), and barley (Habte et al., 2014). The circadian clock regulates the expression of the cold-responsive gene C-repeat binding factor 1 (CBF1) to provide resistance to cold (Maibam et al., 2013). PHYTOCHROME-INTERACTING FACTOR 3 (PIF3) acts as a negative regulator of freezing tolerance by binding to the CBF gene promoter (Jiang et al., 2017). In addition, overexpression of TOC1 (key member of the RESPONSE REGULATOR) leads to stress, as it delays stomatal closure (Legnaioli et al., 2009).

Some genes and pathways in our transcriptome and proteome results suggest the specific response to cold stress in Tg and Sg C. japonica. ABA and GA are crucial hormone antagonists (Finkelstein, 2013; Shu et al., 2018). The ABA–PYR/PTL/RCAR complex is formed by ABA interacting with PYR/PYL (Rodríguez-Gacio Mdel et al., 2009). PYR/PYL was upregulated in Tg and Sg C. japonica. Interestingly, PP2C was downregulated in Tg, indicating the dynamics of ABA content. GID1 interacts with GA to promote the accumulation of the GA-GID1-Della protein, which stimulates plant growth (Sun et al., 2010). GID1, GID2, and DELLA were upregulated in Tg C. japonica, suggesting that GA synthesis increased during the cold treatment. Two DELLA proteins were downregulated in Sg C. japonica and displayed a complex expression trend. Our hormone profiling confirmed that GA content was upregulated in Tg C. japonica, but downregulated in Sg C. japonica. In addition, PIF3, a negative responder to cold stress (Jiang et al., 2017), was downregulated. A previous study suggested that the balance between ABA and GA regulates the development and response to stimuli in plants (Wang et al., 2006). A high level of ABA with high ratios of ABA to GA3 result in low-temperature dormancy in C. sinensis; however, a high IAA level with high ratios of IAA to ZR cause summer dormancy (Pan et al., 2000). The ABA/GA ratio remains low during the formation of tender shoots in C. sinensis during winter (Dai et al., 2021). Our hormone data revealed a similar trend, but there was divergence, as GA4 accumulated in Tg C. japonica, which caused the ABA/GA ratio under cold stress to be lower in Tg C. japonica than in Sg C. japonica, indicating that ABA was suppressed by accumulating GA in Tg C. japonica to maintain growth under cold stress. Moreover, the ABA and IAA ratio also supported this conclusion, e.g., the balance in hormones led to the divergence in cold tolerance between Tg and Sg C. japonica, but the underlying mechanism is unclear and worthy of further exploration.

The transcriptional and translational responses to cold stress were most evident in Tg C. japonica, and were strongly biased toward ribosome and protein processes. Seven ribosomal protein genes and 55 ribosomal proteins were specifically induced under cold stress in Tg C. japonica. Overrepresented ribosomal proteins (e.g., 40S ribosomal protein S7 and large subunit ribosomal protein L22e) are implicated in protein translation (Barkan et al., 2007) and folding (Kaur et al., 2015). In Arabidopsis, the p70 ribosomal S6 kinase homologous gene may function in the adaptation of cells to cold conditions (Mizoguchi et al., 1995). Cold temperatures affect protein translation by ubiquitinating ribosomal proteins (Gong et al., 2020). Overall, targeted changes in selected ribosomal proteins may be an integral part of the difference in freezing tolerance. In addition, plastid ribosome proteins are important components of chloroplast biogenesis. The encoding gene (TCD11) mutant of ribosomal small subunit protein S6 prevents the assembly of ribosomes in chloroplasts under low temperatures (Wang et al., 2017). Thus, these genes may regulate photosynthesis-related pathways. F-type H + -transporting ATPase was specifically induced in Tg C. japonica. The transgenic yeast overexpressing ThVHAc1 exhibits increased tolerance to these abiotic stressors (Wang et al., 2020). Therefore, photosynthesis is a potential pathway regulating cold tolerance.

Although Tg C. japonica was more cold-tolerant, Sg C. japonica expressed many pathways in response to the cold. We observed transcripts and proteins related to primary and secondary metabolism that differed across the populations, with higher expression of phenylpropanoid and tryptophan biosynthetic genes and proteins in Sg C. japonica than in Tg C. japonica. Tryptophan promotes membrane association for the glycosyltransferase involved in the abiotic stress response (Ge et al., 2011). In our results, the higher expression level of glucosyltransferase (c140084.graph_c0) in Sg C. japonica supports this conclusion. Tryptophan protects plants from the effects of salt (Kahveci et al., 2021). In addition, a phenylpropanoid-like pathway may play a role in PsCA-enhanced chilling tolerance in cucumber (Wang B. et al., 2021), and GhANN1 improves salinity tolerance by regulating the phenylpropanoid pathway in cotton (Zhang et al., 2021).

The two populations of C. japonica were obtained from Qingdao (36°19′N, 120°67′E) and Zhoushan, China (30°05′N, 122°21′E) in October 2019. We selected C. japonica plants with more than six true leaves at similar developmental stages. The seedlings were transferred to the greenhouse at the Research Institute of Subtropical Forestry, Chinese Academy of Forestry (30°06′N, 119°96′E) for 3 months of normal culture (plastic pots: top/bottom diameter 18/9 cm, height 10 cm; temperature: 25 ± 2°C; relative humidity: 60%; soil: peat soil 80% + river sand 20%, pH: 6.0 ± 0.5; natural light). Seedlings of similar growth status from the two regions were divided into two groups. Each group was treated in a −5°C incubator for 24 h. A piece of leaf was collected every hour, physiological analyses were performed, and the sample was immediately snap-frozen in liquid N₂, and stored at −80°C until RNA or protein extraction. Each time point contained three biological replications.

The leaves were used to measure the oxidation index (relative electric conductivity) and osmotic index (MDA and CAT activity) to evaluate the degree of divergence in cold tolerance between Tg and Sg C. japonica. The CAT activity and MDA content were measured with the Malondialdehyde Content Detection kit (SolarbioBC0020, Beijing, China) and the Catalase Activity Detection kit (Solarbio^®BC0200), following the manufacturer’s instructions. The cold-treated leaves for the relative electric conductivity assay were placed in 10 mL of deionized water and recorded as S0. The leaves were shaken at 20°C for 1 h, and S1 was determined. After boiling in water for 30 min, the samples were shaken at 20°C to cool, and S2 was determined. (S1–S0)/(S2–S0) was used to measure electrolyte leakage. The average change in the MDA contents, CAT activity, and relative electrical conductivity during each 2 h interval was calculated.

All 24 samples (0, 4, 8, and 24 h cold treatment times of the two populations, with three biological replications each) of total RNA were extracted. Leaves were ground in liquid N₂. Flour (100 mg) was suspended in 2.0 mL of a buffer containing RNAiso Plus (Takara^®9108Q, Dalian, China), centrifuged at 12,000 rpm for 8 min at 4°C, and the supernatant was collected. One mL of 25:24:1 (v/v/v) phenol: chloroform: isoamyl alcohol mixture was added and centrifuged as described above to collect the supernatant. The nucleic acids were precipitated by adding 0.5 volumes of isopropyl alcohol, and ethanol was used to wash the nucleic acids. The colloidal solid was dried under ambient conditions for 15 min and dissolved in RNase-free water.

The TrueSeq RNA Sample Prep kit (Illumina Technologies, San Diego, CA, United States) was used to construct the library and to sequence the RNA of the 24 samples. mRNA was enriched using the oligo dT-bead method. After fragmentation, second-strand cDNA sample synthesis, end-repair, and ligation of the NEBNext Adaptor with a hairpin loop structure were completed. The 150–400 bp insert fragment was selected and amplified by polymerase chain reaction (PCR). The libraries were sequenced based on the Illumina HiSeq™ 2,000 platform. The raw data were stored in the NCBI database (Accession number: PRJNA769967). Adaptor sequences and low-quality reads (padj < 10) in the raw data were removed to obtain clean reads, and Trinity software (v. 2.4.0) (Grabherr et al., 2011) was used to assemble the transcriptome based on the left.fq and right.fq files. The Gene Ontology database (Ashburner et al., 2000); the Kyoto Encyclopedia of Genes and Genomes database (Kanehisa et al., 2004); the KOG (Koonin et al., 2004), the COG [(Tatusov et al., 2000); the Pfam (Finn et al., 2013) and the NR (Deng et al., 2006) databases were used to annotate the genes. Cufflinks program (Trapnell et al., 2010) was used to calculate the gene expression levels based on the FPKM values (Fragments per Kilobase of transcripts per Million mapped reads). The DEGs were identified via the DESeq (v1.10.1)]. P-values were corrected using the Benjamini and Hochberg approach for controlling the false discovery rate (FDR) (Love et al., 2014).

Proteins were extracted according to Rao et al. (2021). A BCA kit (Beyotime, Shanghai, China) was used to measure the protein concentrations, following the manufacturer’s instructions. Briefly, 0.5 M TEAB was used to dissolve the peptides, and a TMT kit (Thermo Fisher Scientific, Waltham, MA, United States) was used to label the peptides. The peptides were graded through a high pH high-performance liquid chromatography system (Agilent 300 Extend C18; Agilent Technologies Inc., Santa Clara, CA, United States). The peptides were dissolved in mobile phase A (0.1% formic acid and 2% acetonitrile) and separated on the EASY-nLC 1,000 (ThermoFisher). Mobile phase B was 0.1% formic acid and 90% acetonitrile. The liquid chromatography-electrospray ionization-tandem mass spectrometry analysis was performed according to Li J. et al. (2017). The MS/MS data were searched against our transcriptome database using MaxQuant (v1.5.2.8) (Cox and Mann, 2008). The following parameters were used for the search: Trypsin/P was the cleavage enzyme, allowing up to two missing cleavages. The mass tolerance of the precursor ions for the first search and main search was 20 and 5 ppm, respectively. The mass tolerance for the fragment ions was 0.02 Da. Carbamidomethyl on Cys was specified as a fixed modification, and methionine oxidation, protein N-terminal acetylation, and deamidation (NQ) were specified as variable modifications. The FDR was adjusted to 1% and the minimum score for the peptides was > 40. For the DEPs analysis, A t-test with p < 0.05, fold change > 1.5, or < 0.667 was considered significantly upregulated or downregulated, respectively.

The GO annotated proteome was derived from the UniProt-GOA database,¹ and InterProScan² was used to supplement the functional annotation. The proteins were classified by GO annotation based on biological processes, cellular components, and molecular functions. KAAS was used to annotate the KEGG protein database description, and the pathways were mapped in KEGG.

The functional enrichment analyses of the transcripts and proteins were performed using KOBAs (Mao et al., 2005). The R (v4.0.5) package WGCNA was used to perform the coexpression analysis using genes that contributed 25% or more of the variance (Zhang and Horvath, 2005), parameters: weighted network = unsigned, power = 8, the module genes, and Pearson’s correlation values were used to identify the correlation between the modules and the cold treatment groups. Cytoscape (v3.5.1) was used to visualize the gene networks. The expression trend analysis was performed using the R package MaSigPro (Nueda et al., 2014).

The candidate genes related to the photosynthesis, ribosome, and plant hormone signaling pathways were screened. RT-qPCR of the 10 genes was carried out using the same RNA used for transcriptomic sequencing, with the same operational approaches as described before (Fan et al., 2021b). The NADPH gene was used as the reference gene (Pan et al., 2020). The relative expression level was determined by the 2^–ΔΔCt method (Livak and Schmittgen, 2001). The expression level of Tg-CK was set to 1, and that of other samples was determined relative to Tg-CK. The primer sequences are shown in Supplementary Table 14.

Leaves were ground in liquid nitrogen, and 10 ml of methanol (80%, v/v) was added to extract the hormones. The mixture was centrifuged at 3,600 rpm for 11 min at 4°C, the supernatant was collected, and passed through a Chromosep C18 column (C18 Sep-Park Cartridge, Waters Corp., Milford, MA, United States) and dried. An enzyme-linked immunosorbent assay (ELISA) was used to quantify the hormones according to Chen et al. (2021). Briefly, the mouse antigens and antibodies and the horseradish peroxidase substrate were produced at China Agricultural University (Beijing, China). Coating buffer (1.5 g/L Na₂CO₃, 2.93 g/L NaHCO₃, and 0.02 g/L NaN₃) and the antigens were incubated at 37°C for 4 h for ZR, GAs, BR, and ABA and 24 h for the IAA assay at 4°C. The plates were washed with PBS. An ELISA reader (model EL310, Bio-TEK, Winooski, VT, United States) was used to check color. Each group contained three biological replicates and three technical repetitions.

The statistical analysis was performed via SPSS (v21.0) software (IBM Corp., Chicago, IL, United States). The Student’s t-test was used to calculate differences in the mean values. A p-value < 0.05 was considered significant. The data are presented as mean ± standard deviation of three independent experiments.