Integration of transcriptome and proteome analysis reveals the mechanism of freezing tolerance in winter rapeseed

Winter rapeseed seedlings are susceptible to low temperature during overwintering in Northwest China, leading to reduced crops production. Freezing stress is one of the main environmental stresses in Northwest China from late autumn to early spring, an eventful period for overwinter survival rate of winter rapeseed. However, the molecular mechanism of freezing tolerance formation is still very backward in winter rapeseed. In this study, using a pair of freezing-sensitive and freezing-resistant cultivars NQF24 and NTS57, the exhaustive effects of freezing stress on freezing tolerance formation were evaluated by analyzing leaf at the levels of transcriptome, proteome, physiology and ultrastructure. There were 8497 and 7358 differentially expressed genes (DEGs) and 418 and 573 differentially abundant proteins (DAPs) identified in the leaf of NQF24 and NTS57 under freezing stress, respectively. Function enrichment analysis showed that most of the enriched DEGs and DAPs were associated with plant hormones signal transduction, fatty acid metabolism, ribosome, plant-pathogen interaction and secondary metabolites biosynthesis. Freezing tolerance is formed by enhanced signals transduction, increased the biosynthesis of protein and secondary metabolites, enhanced reactive oxygen species (ROS) scavenging, more osmolytes, lower lipid peroxidation, and stronger cell stability. These results can be taken as selection indicators in freezing tolerance breeding program in rapeseed.


Introduction
Freezing stress (< 0 °C) has profound and diverse effects on plant growth and development, it poses a serious threat to the geographical distribution of plant species (Ding et al. 2019a), it can also cause ice formation in the apoplast of plant tissues. The accumulation of intracellular ice physically disrupts the cell membrane , leading to severe cell dehydration. Freezing stress is common in the area of Northwest China which could not only impact on the quality of crops, but also reduce the production of crops, especially for winter crops grown in high latitudes and altitudes regions . Furthermore, freezing stress damaged the stability of proteins or protein complexes, and increased the accumulation of ROS and osmolytes during seed development Pu et al. 2019). Additionally, as we all know, intracellular organelles mainly include mitochondria, chloroplasts, nucleus, peroxisomes, endoplasmic reticulum, Golgi apparatus, endosomes, vacuoles, phagosomes, lysosomes, and secretory vesicles, which can not only transport and exchange materials with each other, but also work on specific cellular functions such as energy metabolism and transcription regulation . When plants are exposed to unfavorable environments, the stress-generated signals from all organelles are integrated to regulate stress responsive genes expression and other cellular activities to resume cellular homeostasis (Gong et al. 2020). Some previous studies in model plants have shown that the adaptation of the plant to freezing stress must change gene expression at the cellular level, which in turn alters the protein profiles Zhang et al. 2017). However, to date, the limited information is available about freezing-response genes or proteins in winter rapeseed.
Winter rapeseed (Brassica napus L.) is an important oilseed crop that is widely cultivated in Northwest China. However, its growth is easily persecuted to low temperature in Northwest China, which prevents winter rapeseed from overwintering and propagation, leading to severe reduction of crops yield . As a result, it is critical to acquire strong freezing resistance for winter rapeseed so as to overwinter safely. In recent years, the studies by Xu et al. (2018) and Zeng et al. (2018) suggested the impacts on the roots and leaves of two winter rapeseed cultivars after freezing treatment by isobaric tags for relative and absolute quantification. Pu et al. (2019) reported the effects on the freezing resistance formation in the leaves of two rapeseed cultivars after freezing treatment by high-throughput RNA sequencing (RNA-seq). However, all the above researches are based on traditional data-dependent acquisition (DDA) models, which might result in lower accuracy in the quantification of low abundance proteins. The data-independent acquisition (DIA) technique has an excellent detection rate, and it is considered as one of the most striking techniques by the researchers in Nature Methods (Egertson et al. 2013). Our recent study obtained a freezing-responsive molecular network after short-term successive freezing treatment by the sequencing-based transcriptome and DIA-based proteome (Wei et al. 2021).
In the present study, under long-term repeated freezing stress, on the one hand, we analyzed the physiological, biochemical and ultrastructural properties in the leaves of resistant and susceptible winter rapeseed cultivars; on the other hand, based on the results of DNA methylation analysis, the combined analysis of transcriptomics and proteomics was used to identify the freezing-responsive DEGs and DAPs in the leaves of two winter rapeseed cultivars; thereafter, the main metabolic pathways, cellular processes and differentially expressed transcription factors (TFs) associated with freezing stress were investigated; finally, we got a sophisticated molecular regulatory network of winter rapeseed in response to long-term freezing stress. It will be helpful to expedite the breeding and production of winter rapeseed cultivars with strong freezing tolerance.

Plant samples and freezing stress treatments
A pair of winter rapeseed cultivars freezing-resistant NTS57 (named NS) and freezing-sensitive NQF24 (named NF) used in the present study, were produced by Gansu Agricultural University (Wei et al. 2021). Seedlings of the two cultivars were grown in 5 L plastic pots for 5 weeks under constant 22 °C and 16 h /8 h (light/dark) cycle with 350 μmol m −2 s −1 flux density. Subsequently, potted plants were divided into two groups for freezing treatments. Briefly, the plants in freezing treatment group (Treatment, T) were transferred to an artificial climate room [− 4 °C/22 °C, 12 h/12 h (light/dark) cycle with 350 μmol·m −2 ·s −1 flux density] for 2 days. The plants in control group (Control, named C) at the equal growth stage were placed in a separate growth chamber at 22 ℃ for 2 days [12 h/12 h (light/dark) cycle with 350 μmol·m −2 ·s −1 flux density]. The third leaf after treatment was selected for all experiments with three biological replicates. Fresh samples were used for determination of net photosynthetic rate and transmission electron microscope (TEM) observations, while another set of samples were refrigerated in liquid nitrogen instantly and preserved at − 80 ℃ for total RNA and proteins extraction, and physiological and biochemical measurements.

Determination of net photosynthetic rate
Net photosynthetic rate was measured according to Ma et al. (2017) with some modifications. Photosynthesis was evaluated using a LI-COR 6400 portable gas analysis system (Lincoln, NE, USA) setting a built-in light source of 1000 μmol photons m −2 ·s −1 . Each sample was repeatedly determined six times during 10:00 am and 11:00 am to minimize the error.

TEM observations
Fresh rapeseed leaves from treatments and controls were used for TEM observations. Each sample was cut into small size of 2 mm 2 in the 2.5% glutaraldehyde fixative for 16 h at 4 °C, followed by a post-fixation with 1% OsO 4 for 2 h at room temperature. All fixed samples were dehydrated in a graded series of ethanol and embedded in Spurr resin. The ultrathin section 60-80 nm were finally stained with 2% uranium acetate saturated alcohol solution to avoid light for 10 min and 2.6% lead citrate to avoid CO 2 for 10 min, and then rinsed with ultra-pure water for 3 times. Observations were performed on an H-7650 electron microscope. High-resolution TEM images were produced as described in Toyooka et al. (2014).

Physiological and biochemical measurements
Rapeseed leaves (0.1 g) were used for measuring physiological and biochemical properties, including antioxidant enzymes [catalase (CAT), peroxidase (POD) and superoxide dismutase (SOD)] activities, lipid peroxidation level, and the contents of osmotic regulation substances (soluble sugar and protein). The antioxidant enzymes activities were measured using the kits from Beijing Solarbio Science & Technology following the manufacturer's protocols . The level of lipid peroxidation was evaluated by determining MDA content according to Wei et al. (2020). The samples were ground and extracted in 20% trichloroacetic acid buffer containing 0.67% thiobarbituric acid at 4 ℃, and then centrifuged at 12,000×g for 20 min. Absorbance of each sample supernatant was measured at 450, 532, and 600 nm. The contents of soluble sugar and soluble protein were measured as described by Song et al. (2015). For each sample, 0.5 g leaves were collected in a test tube and added 5 mL distilled water, the tube was incubated at 100 ℃ for 30 min, and then centrifuged at 10,000×g for 10 min, the supernatant was used for determining the content of sugar. The soluble protein was extracted from 0.2 g samples with 50 mM Tris-HCl (pH 8.8) containing 1 mM dithiothreitol, and the concentration was determined at 595 nm.

RNA extraction and RNA sequencing in rapeseed
The frozen leaf samples were ground in liquid nitrogen and total RNA were extracted using the TRIzol Reagent (Tiangen Biotech, China) according to the manufacturer's instructions. The sequencing library was constructed and sequenced in Gene Denovo Biotechnology Co. (Guangzhou, China) on an Illumina HiSeqTM 2500 platform. The transcriptomic raw data have been uploaded to the SRA database of NCBI (PRJNA685002).

Filtering of reads and gene expression analysis
Raw data were filtered to remove adapters, poly-N and low-quality reads to obtain clean reads. The high-quality paired-end reads were mapped to rapeseed reference genome (Brassica napus v2.0) using TopHat v2.0.3.12 (Kim et al. 2013). The levels of gene expression were quantified and normalized as fragments per kilobase per million mapped reads (FPKM), which can be diametrically used to identify DEGs in pair-wise comparisons (Trapnell et al. 2010). Genes with p-value < 0.001 and value of |log 2 FoldChange|≥ 2 were regarded as DEGs using the edgeR package (http:// www.rproje ct. org/).

Protein extraction and digestion
Total protein was extracted using the chilled acetone precipitation method according to Wei et al. (2020) with some changes. Briefly, samples were first ground into powder in liquid nitrogen, then dissolved in 2 mL lysis buffer containing 8 M Urea, 2% SDS, 1 × Protease Inhibitor Cocktail, followed by sonication on ice for 20 min and centrifugation at 11,000×g for 10 min at 4 ℃. The supernatant was pipetted 1 3 into a new centrifuge tube. Protein was precipitated with chilled acetone for 16 h at − 20 ℃, cleaned with acetone twice, re-dissolved in 8 M urea, and homogenized in ice for 3 min using an ultrasonic homogenizer. The homogenate was centrifuged at 11,000×g for 15 min at 4 ℃, then the supernatant was preserved and the protein concentration was gauged by the BCA Protein Assay Kit (Thermo Fisher, Waltham, MA). A total of 50 µg protein were suspended in 50 μL 50 mM ammonium bicarbonate, reduced by 10 mM dithiothreitol for 1 h at 55 °C and alkylated by 20 mM iodoacetamide for 1 h at 37 °C in the dark. Thereafter, each sample was precipitated with 500 ul chilled acetone for 16 h at − 20 ℃. The precipitates were washed three times with chilled acetone and resuspended in 50 mM ammonium bicarbonate. Finally, all samples were digested with sequencing-grade modified trypsin (Promega, Madison, WI) at a substrate/ enzyme ratio of 50:1 (w/w) at 37 °C for 16 h.

High pH reverse phase separation and LC-MS/MS analysis
Digested peptides from each sample were segregated on an Ultimate 3000 system (Thermo Fisher, USA) with a high pH reversed phase XBridge BEH C18 column (4.6 mm × 250 mm, 5 μm). High pH separation was performed using a linear gradient, starting from 5% B (20 mM ammonium formate in 80% acetonitrile, pH 10.0) to 45% B in 40 min. Peptides were eluted at a flow rate of 1 mL/min. For each sample, 10 fractions were convened and dried by a vacuum pump. All dried samples were saved at -80 ℃ until the next analysis.
The dried peptides were dissolved in 30 μL buffer D (0.1% formic acid and 0.1% formic acid in acetonitrile). LC-MS/MS analysis was operated on an Orbitrap Fusion Mass Spectrometer (Thermo Fisher, Germany) coordinated with a Nano ACQUITY UPLC system (Waters, Milford, MA). In short, 10 μL peptides were added to the trapping column (Acclaim PepMap C18, 100 μm × 2 cm) and then eluted on an analytical column (Acclaim PepMap C18, 75 μm × 25 cm) by a 120 min gradient of buffer D. The flow rate of column was retained at 500 nL/min, and the electrospray voltage was set to 2.1 kV.
For library generation, to acquire data in DDA mode, the Orbitrap Lumos was performed in a positive mode. MS spectra (350-1500 m/z) were collected at 1.2 × 10 5 resolution to attain an automatic gain control (AGC) target of 4.0 × 10 5 . The precursor ions were used for MS/MS detection with a collision energy of 32. Dynamic exclusion was set to 30 s to exclude all isotopes. MS/MS spectra were collected at 1.5 × 10 4 resolution to attain an AGC target of 5.0 × 10 4 . For data acquisition, we used DIA mode. MS spectra (350-1500 m/z) were collected at 1.2 × 10 5 resolution to attain an AGC target of 5.0 × 10 5 . MS/ MS spectra were collected at 3.0 × 10 4 resolution to attain an AGC target of 1.0 × 10 6 with a collision energy of 32. All 60 segments were selected for MS/MS. The acquisition window covered 1 m/z through 60 consecutive isolation windows.

Mass spectra data and protein quantification
Spectronaut Pulsar 11.0 software (Bruderer et al. 2015) was applied to analyze DDA MS/MS data, and the settings were as follows: digestion enzyme, trypsin; fixed modification, carbamidomethylation (C); variable modifications, oxidation (M); precursor mass tolerance, 20 ppm; fragment mass tolerance, 0.05 Da. The MS/MS data were searched from Brassica napus database (Brassica napus v2.0; Table S1), which were downloaded from the NCBI on May, 2018 with 123,465 sequences. False discovery rates of all peptides were set to 1%, which were evaluated by a reverse search sequence. Subsequently, raw data from each sample with the spectra library based on DIA MS/MS data were analyzed by Spectronaut Pulsar 11.0 with default parameters depending on the previously constructed DDA protein dataset. All proteins quantification normalization were performed by the local normalization with Pulsar software. The MS raw data have been submitted into a public iProX database (Ma et al. 2019a) with the accession number of IPX0002682000. Proteins with p-value < 0.001, peptides ≥ 2 and value of |log 2 FoldChange|≥ 1 were regarded as DAPs using the edgeR package.

Bioinformatics and statistics analysis
The correlation coefficient between three replicates was analyzed to assess the repeatability of the experimental results. Principal component analysis (PCA) was carried out to disclose the relationship between samples by R package (http:// www.r-proje ct. org/). Genes and proteins were annotated using the Gene Ontology (GO) database (http:// geneo ntolo gy. org/) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http:// www. kegg. jp/) to obtain their functions (Binns et al. 2009;Kanehisa et al. 2016). Enrichment analysis of DEGs and DAPs were conducted on the basis of GO database and KEGG database. All physiological data were performed using one-way analysis of variance and mean differences of three biological replicates were obtained by LSD test, p < 0.05 was considered for all experiments at a statistically significant level.

Physiological and biochemical responses to freezing stress
NF plants exhibited enhanced tissue damage than NS plants after freezing treatment (Fig. 1A). Compared to the controls, freezing stress provoked a significant reduction in the net photosynthetic rate of both rapeseed cultivars, but the net photosynthetic rate of rapeseed cultivar NF was significantly lower than that of cultivar NS under freezing treatment (Fig. 1B). In addition, TEM analysis was conducted in both cultivars after freezing treatment. Compared to the control, the chloroplasts in mesophyll cells were brutally collapsed in cultivar NF under freezing treatment, leading to thylakoid lamellae exposure, whereas, more plastoglobuli were observed in cultivar NS under freezing treatment (Fig. 1C). These results suggested that the freezing treatment resulted in a greater injury to the leaves of NF than that of NS.
The enzyme activities of CAT and SOD in the leaf of both cultivars were found to be significantly decreased after freezing treatment, but the CAT activity in cultivar NS was significantly higher than that of cultivar NF (Fig. 2AB). The POD activity and MDA content in the leaf of both cultivars were significantly increased after freezing treatment, but the content of MDA in cultivar NF was significantly higher than that of cultivar NS (Fig. 2CD), indicating that freezing stress resulted in enhanced lipid peroxidation in cultivar NF. Furthermore, no obvious changes in soluble sugar content were found in both cultivars, the content of protein was significantly increased in cultivar NS after freezing treatment (Fig. 2EF). These results indicated that enhanced ROS scavenging activities and increased osmolytes accumulation, and decreased lipid peroxidation might create stronger freezing resistance in cultivar NS under freezing stress.

Quality control analysis of transcriptome and proteome in winter rapeseed
In the current research, 80% transcripts were also detected at the protein level (Fig. 3A). The transcriptome in each sample was determined by RNA-Seq on the Illumina platform, and expression levels were estimated by FPKM, a total of 83,846 genes including 77,046 known genes and 6800 new genes were detected in all samples of developing leaves of both cultivars ( Fig. 3A; Table S1). A total of 29,798 unique peptides from 37,247 spectra, corresponding to 13,933 proteins, encoded by 11,147 genes in transcriptome, were identified in all samples of developing leaves of both cultivars ( Fig. 3A; Table S1). More down-regulated DEGs and more accumulated DAPs were discovered in both cultivars ( Fig. 3B; Table S3; Table S4). PCA results denoted that three biological replicates of NFC, NFT, NSC and NST had a good consistency (Fig. 3CD), and all replicates exhibited an extremely high overlap in proteome and transcriptome for each sample ( Figure S1). Moreover, the consistencies between three biological replicates of NFC, NFT, NSC and NST samples were analyzed by evaluating their correlation coefficients, all of which were higher than 0.9 ( Figures S2,   S3). These results ensured the reliability and accuracy of the proteome and transcriptome analysis.

Quantitative transcriptome and proteome analysis of leaves in both rapeseed cultivars under freezing treatment
To identify DEGs and DAPs that were induced by freezing treatment, the changes at transcription level and protein relative abundance in leaves of both cultivars were analyzed by calculating the relative level ratios of treatment to control. Compared to the corresponding control, a total of 18,117 (up, 7658; down, 10,459) genes and 16,978 (up, 6453; down, 10,525) genes were identified to be differentially expressed in leaves of cultivars NF and NS under freezing stress, respectively ( Fig. 4A; Table S2). Compared to the corresponding control, a total of 587 and 510 proteins were found to be accumulated in abundance in leaves of cultivars NF and NS under freezing stress, respectively, whereas 168 and 400 proteins were reduced in abundance in leaves of both cultivars under freezing stress ( Fig. 4B; Table S3). Both DEGs and DAPs numbers were significantly varied between both cultivars, Fig. 2 The physiochemical changes in leaves of winter rapeseed under freezing treatment. A CAT activity; B SOD activity; C POD activity; D MDA content; E Soluble sugar content; F Soluble protein content; one-way ANOVA followed by LSD 0.05 test was used for determination of significant differences suggesting that the dynamics of gene or protein expression regulation in both cultivars were very different.
A total of 9620 shared DEGs and 337 shared DAPs were owned in NF and NS, among the shared DEGs, 4411 were up-regulated and 5209 were down-regulated in NF, while 4443 were up-regulated and 5177 were downregulated in NS ( Fig. 4C; Table S4); among the shared DAPs, 288 were accumulated in abundance and 49 were reduced in NF, while 265 were accumulated in abundance and 72 were reduced in NS ( Fig. 4D; Table S4). Furthermore, 8497 DEGs and 418 DAPs were specifically identified only in NFC_NFT, while 7358 DEGs and 573 DAPs were specifically identified only in NSC_NST ( Fig. 4AB; Table S5), which are the molecular basis for the difference in freezing tolerance between both rapeseed cultivars. These unique DEGs and DAPs were considered as candidate genes and proteins associated with freezing resistance of winter rapeseed.

GO and KEGG analysis of DEGs and DAPs associated with freezing tolerance
All unique DEGs and DAPs in both cultivars under freezing stress were subjected to GO classification annotation ( Figure S4; Figure S5; Table S5) and KEGG enrichment analysis. Annotated DEGs and DAPs were divided into three major functional categories including biological process, cellular component and molecular function. In the cultivar NF, most DEGs and DAPs were enriched in the cellular process, metabolic process, single-organism process, response to stimulus and biological regulation in the biological process. The cell, cell part and organelle were the main cellular components enriched in most DEGs and DAPs. The binding and catalytic activity were the largest groups in the molecular function enriched by both DEGs and DAPs. Furthermore, top 20 KEGG enrichment pathways showed that ribosome (ko03010), linoleic acid metabolism (ko00591), biosynthesis of antibiotic (ko01130), metabolic pathways (ko01100), propanoate metabolism (ko00640), fatty acid biosynthesis (ko00061), alpha-linolenic acid metabolism (ko00592), fatty acid elongation (ko00062), and biosynthesis of secondary metabolites (ko01110) pathways were significantly (q < 0.05) enriched in DEGs ( Fig. 5A; Table S5), while endocytosis (ko04144), tryptophan metabolism (ko00380), fatty acid degradation (ko00071), riboflavin metabolism (ko00740), and glycine, serine and threonine metabolism (ko00260) pathways were significantly (q < 0.05) enriched in DAPs ( Fig. 5B; Table S5).

Identification of differentially expressed TFs in leaves of both rapeseed cultivars under freezing treatment
A total of 1518 and 1388 differentially expressed TFs were identified in transcriptome of NF and NS after freezing treatment, respectively (Table S6). In both rapeseed cultivars, most of the up-regulated TFs belonged to ethylene response factor (AP2/ERF), WRKY transcription factor (WRKY), NAC domain-containing protein (NAC), heat stress transcription factor (HSF), and protein TIFY (TIFY) families, while most of the down-regulated TFs belonged to bHLH (bHLH), zinc finger protein (C2C2), basic leucine zipper/ transcription factor TGA-like (bZIP), auxin response factor/growth-regulating factor (ARF/GFP), and transcription termination factor (mTERF) families (Fig. 7AB). However, only FHA, CSD and PBF families were identified and reduced in abundance in leaves of both cultivars under freezing treatment. Our results indicated that these TFs genes are important in response to freezing stress.

Discussion
It is acknowledged that freezing stress causes serious injuries to the plant, which is general in many winter crop growing areas around China (Zhu, 2016), especially for winter rapeseed grown in Northern China ). Some studies have suggested the effects of freezing stress on winter rapeseed resistance formation by the transcriptome, proteome and physiology Pu et al. 2019). However, the molecular mechanism and the formation of freezing tolerance in winter rapeseed are still poorly understood. The winter rapeseed cultivar NS can survive normally in winter in most northern areas from 38 degrees  north latitude to 42 degrees north latitude, and its tolerance at low temperature environments is better than other types of winter rapeseed, which contains abundant cold-resistant genes. Therefore, in this study, we focus on the levels of transcript, protein, physiology and biochemistry to elucidate the molecular basis of freezing tolerance in leaves of winter rapeseed freezing-tolerant NS and sensitive NF. It will be greatly helpful to plant molecular breeding under freezing stress.

Physiological and biochemical changes between both winter rapeseed cultivars after freezing treatment
Different abiotic stresses can extremely affect physiological and biochemical traits in different plants (Jamshidi Goharrizi et al. 2019, 2020a. Cold stress, including chilling (0-15 ℃) and freezing (< 0 ℃), is one of the main abiotic stresses that significantly affect the growth and agricultural productivity of plants . Chilling stress decreases membrane fluidity, destabilizes protein complexes, and impairs photosynthesis , whereas freezing stress causes more serious injuries to the plant. Freezing temperature triggers intracellular ice formation, which disrupts the cell membrane, causing a decrease in the water potential outside the cell, thus leading to severe cell dehydration, and may even create cell death Ding et al. 2019a). Over the past decades, it has been reported that ROS plays an important role in response to various abiotic stresses (Baxter et al. 2014). Under adverse environmental stresses, the ROS scavenging enzymes are necessary to maintain normal cellular redox homeostasis (Ray et al. 2012). For example, Jamshidi Goharrizi and colleagues proved that there is a negative correlation between ROS production and the enzymatic antioxidant defense system activity under different abiotic stresses. Also, they demonstrated that the level of enzymatic antioxidant defense system in more tolerant cultivars is more active than that of sensitive one (Jamshidi Goharrizi et al. 2020b. In this study, POD enzyme activity was significantly increased in both rapeseed cultivars under freezing stress, but the level of enzyme activity in cultivar NS was higher than that in cultivar NF (Fig. 2C). These results were broadly in line with Zeng et al. (2018) and Ma et al. (2019b). Some previous researches reported that cold-stressed arabidopsis roots showed the low electrolyte leakage and bad cell membrane integrity (Deng et al. 2015), the rapeseed leaves exhibited enhanced accumulation of MDA and soluble protein under freezing stress ). In the present study, freezing stress induced the production of MDA, soluble sugar and soluble protein in leaves of both cultivars, the level of lipid peroxidation in cultivar NS was significantly lower than that in cultivar NF (Fig. 2D), while the contents of soluble sugar and soluble protein in cultivar NS was higher than that in cultivar NF (Fig. 2EF). Taken together, high levels of ROS scavenging, osmolytes and lower level of lipid peroxidation in leaves of the winter rapeseed under freezing stress would furnish a redox homeostasis environment for plant growth and development, leading to the stronger freezing resistance in NS plants.

Comparison of the differences in freezing tolerance between both winter rapeseed cultivars
Numerous TFs have been reported in plants that monitor the expression of C-repeat binding transcription factor (CBF) by combining corresponding cis-elements in their promoters under cold stress. After exposure to freezing stress, CBF proteins rapidly recognize the promoter regions of downstream cold-regulated genes to activate their expression, thereby enhancing freezing tolerance (Jia et al. 2016;Zhao et al. 2016). For instance, bHLH transcription factor ICE1, calmodulin-binding transcription activator 3 (CAMTA3), brassinazole-resistant 1 (BZR1), heat shock transcription factor C1 (HSFC1), and MYB88 /MYB124 transcription factors positively regulate the CBF expression at the transcriptional level that contribute to cold stress Park et al. 2015;Xie et al. 2018). In contrast, MYB15, phytochrome interacting factor 3 (PIF3), PIF4, PIF7 and ethylene sensitivity 3 (EIN3) negatively modulate the transcriptional repression of CBF (Shi et al. 2012;Kim et al. 2017;Jiang et al. 2017). In the present study, most of AP2/ERF, WRKY, NAC, HSF, and TIFY transcription factors were up-regulated in both cultivars, while transcription factors bHLH, bZIP, ARF/GFP, OFP, and mTERF were down-regulated ( Fig. 7AB; Table S6). In addition, we found that AP2/ERF and HSF genes have differences in DNA methylation levels (Table 1), indicating that these transcription factor genes may defend against freezing stress through epigenetic modifications. Interestingly, these differentially expressed TFs encoded by genes were not found at the protein level, implying that those transcriptional moderators might regulate the gene expression just at the transcription level in contribution to the freezing resistance of winter rapeseed under freezing stress.
Furthermore, some protein kinases including calciumdependent protein kinases (CDPKs), CBL-interacting serine/ threonine-protein kinases (CIPKs), mitogen-activated protein kinases (MAPKs), receptor-like kinases (RLKs), serine/ threonine-protein kinases (SPK), and SNF1-related protein kinases (SnRKs) have been characterized to be crucial regulators of cold-stress responses in plant Zhao et al. 2017;Ding et al. 2019b). In addition, a putative cold sensor chilling-tolerance divergence 1 (COLD1), was reported to mediate a cold-sensing calcium channel in rice, leading to the activation of cold-regulated genes Shi et al. 2018). A similar result was obtained in the present study, a total of 7 CDPKs, 4 CIPKs, 13 MAPKs, 81 RLKs, 1 SnRK, and 42 SPKs were up-regulated in the cultivar NF under freezing stress, whereas 3 CDPKs, 3 CIPKs, 95 RLKs, 2 SnRKs, and 37 SPKs were down-regulated; a total of 9 CDPKs, 8 CIPKs, 9 MAPKs, 22 RLKs, 1 SnRK and 30 SPKs were up-regulated in the cultivar NS under freezing stress, whereas 9 CDPKs, 6 CIPKs, 11 MAPKs, 124 RLKs, 8 SnRKs and 49 SPKs were down-regulated (Table S5). Combined with the differences in DNA methylation levels of CDPK and CIPK genes (Table 1), a popular explanation is that freezing treatment might induce DNA methylation modification of these genes, thereby trigger the activation of other freezing-responsive genes, which needs further study. Ribosomal proteins (RPs) are well-known for their role in maintaining the stability of the ribosomal complex, which performs mRNA directed protein synthesis (Moin et al. 2016). A drastic lowering of the ambient temperature entails extensive reprogramming of the gene expression pattern (He et al. 2021). In this study, compared to the control, all DAPs enriched in the ribosome pathway were accumulated in abundance in NF under freezing stress ( Fig. 5B; Table S5), implying that the protein synthesis was enhanced, which might be beneficial to cope with freezing stress in NF.
Photosynthesis is the main driving force for the growth and biomass production in plants, it can supply the energy and carbon for the biosynthesis of necessary organic compounds (Nowicka et al. 2018). In this study, the photosynthesis pathway was significantly enriched by DEGs and DAPs in NS, and all ferredoxins (FDX) and chlorophyll a/b-binding proteins (CBPs) enriched in this pathway were reduced in transcription and abundance ( Fig. 6; Table S5), which can recruit photons and transfer the absorbed energy towards the photosynthetic reaction centre (Palm et al. 2018), suggesting that the photosynthesis activity was impaired under freezing stress. Consistent with transcriptome and proteome results, the measured net photosynthetic rates were decreased in NS and NF under freezing stress, but it was significantly higher in the former than that in latter (Fig. 1B), which might attribute to more stable chloroplasts ultrastructure in the leaves of NS.
Plant hormones, including abscisic acid (ABA), auxin, brassinosteroid (BR), cytokinin, ethylene, gibberellin, jasmonic acid (JA) and salicylic acid, combine endogenous substances with environmental signals to regulate plant growth, development and defense (Huang et al. 2017). In many instances, plants respond to environmental stresses by producing amount of ABA, ethylene and JA (Huang et al. 2016;Hu et al. 2017). Similar results were obtained in this study, a total of 115 and 104 unique DEGs in NF and NS were respectively involved in plant hormone signal transduction pathway (Table S5), among them, 10 and 7 DEGs were up-regulated in NF and NS, respectively, and participated in ABA signaling, encoding abscisic acid-insensitive protein (ABI), SPK and protein phosphatase 2C (PP2C), whereas 6 and 10 DEGs encoding ABI, SPK and abscisic acid receptor (ABR) were down-regulated in NF and NS, respectively. Five out of 48 DEGs were up-regulated in NF, and 15 out of 39 DEGs were up-regulated in NS. They are involved in auxin signaling, encoding auxin-induced protein (AIP), auxin-responsive protein (ARP), auxin transporter protein (ATP), and GH3 auxin-responsive promoter (GH3). Therefore, we surmised that there may be more than one ABA and auxin metabolic pathway related to freezing tolerance, and one of them was enhanced. Furthermore, 14 upregulated DEGs and 7 down-regulated DEGs were found in NF and NS, respectively, and gone in for ethylene signaling, encoding ethylene insensitive protein (EIN), EIN3-binding F-box protein (EBF), ethylene-responsive transcription factor (ERF) and ethylene receptor (ER). Correspondingly, we found that two key enzymes of ethylene biosynthesis, 1-aminocyclopropane-1-carboxylate synthases (ACS) and 1-aminocyclopropane-1-carboxylate oxidases (ACO), were up-regulated in NF and down-regulated in NS. This may be due to feedback regulation in NS leading to reduced ethylene metabolism. Five DEGs were up-regulated in NF and engaged in JA signaling, encoding TIFY proteins and transcription factor MYC2, which was a key regulator for the activation of JA under cold stress (Kazan, 2015). These are consistent with what has been found in the previous research, which indicated that high levels of plant hormones could improve the freezing resistance (Pu et al. 2019;Wei et al. 2021). In addition, one DEG was down-regulated in NF, encoding protein BZR, which is a plant-specific positive regulator of BR signaling . BRs can be perceived by the brassinosteroid insensitive1 (BRI1; She et al. 2011), and BRI1 interacted with BRI1-associated receptor kinase and positively regulated BZR1 . Interestingly, in this study, 2 BRI1 genes were found and they were down-regulated in NF, however, we found only one down-regulated BZR2 gene in NS. Therefore, we speculate that BZR1 could combine with BRI1 to negatively regulate BR signaling in NF under freezing stress. Apart from that, numerous studies have reported that BRs could improve photosynthesis and ROS scavenging enzymes activities (Xia et al. 2018;Yan et al. 2020). Overall, these findings are in accordance with the results of photosynthesis and antioxidant enzymes (Figs. 1B, 2). Simultaneously, other DEGs identified in both cultivars in our study mainly represent some TFs, which are essential for plants under abiotic stresses (Li et al. 2019;Frank et al. 2020). Taken together, these results suggested that plant hormones play profound roles in coping with freezing stress of rapeseed.
Plant lipoxygenases (LOXs) are a kind of multifunctional fatty acid dioxygenases that catalyze peroxidation of linolenic acid and linoleic acid. The synthesis of jasmonic acid begins with the peroxidation of α-linolenic acid by lipoxygenase. It was reported that biotic and abiotic stresses could cause the liberation of α-linolenic acid, which is an important precursor for the biosynthesis of jasmonic acid (Yan et al. 2013;Zhao et al. 2014), which interacted with other hormones signaling to regulate leaf tolerance to cold stress . In this study, all unique LOXs participated in the alpha-linolenic acid metabolism were found to be down-regulated in NF, whereas they were up-regulated in NS ( Fig. 5A; Table S5). An interesting finding was that the alpha-linolenic acid metabolism was enriched by DAPs in NF, and all DAPs were accumulated in abundance ( Fig. 5B; Table S5). These results suggested that there are differences in temporal and spatial pattern between the expression of gene and protein. Beyond those, combined with JA signaling resulting from DEGs, implying that increased jasmonic acid signal might regulate the expression of downstream freezing-responsive genes.
Plant secondary metabolites often have no fundamental role in the maintenance of plant life processes, but they are crucial for the plant to interrelate to its environment for adaptation and defense (Ramakrishna et al. 2011). In the present study, there were 7 and 6 unique up-regulated DEGs in NF and NS, respectively, encoding caffeic acid O-methyltransferase (CCMT), shikimate O-hydroxycinnamoyltransferase (SHCT) and S-adenosylmethionine-dependent methyltransferase (SAMTs), implicated in the secondary metabolism ( Fig. 5B; Table S5), which were proclaimed in many plant species (Ranjan et al. 2020). It is in agreement with those reported by Zeng et al. (2018) and Pu et al. (2019). It is interesting to note that some DNA methylation differences in the CCMT gene were found under freezing inducement (Table 1), which is particularly important when investigating the epigenetic freezing-reponsive mechanism of winter rapeseed. A further novel finding was that the plant-pathogen interaction pathway was significantly enriched by DEGs in NS, and most of DEGs enriched in this pathway were up-regulated ( Fig. 6A; Table S5), implying that its defense system was strengthened after freezing treatment. These results suggested that accumulated secondary metabolites and enhanced plant-pathogen interaction together contribute to the formation of freezing tolerance of rapeseed cultivar NS under freezing stress.
It is acknowledged that heat shock proteins (HSPs) regarded as the target genes of HSFs, are beneficial to defense against cold stress in grasses (Wang et al. 2016). Similarly, Lin et al. (2019) reported that some HSPs and chaperone proteins (CPs) were induced by cold stress. Scarpeci et al. (2008) indicated that ROS induce enhanced abiotic stress tolerance in plants through direct interaction with HSFs. In the present study, some up-regulated HSFs, CPs and HSPs were concomitant with a high level of POD activity ( Fig. 2C; Table S5), showing that the HSFs might interact with ROS signals to regulate the expression of CPs and HSPs under freezing stress. Furthermore, several upregulated DEGs encoding protein disulfide isomerase (PDI) associated with protein synthesis were found in NS, whereas they were down-regulated in NF (Table S5). In addition, we discovered 10 unique down-regulated DEGs in NF, which encode mitochondria-localized ATP synthase subunit (ASS) and cytochrome c oxidase (CCO), and 2 unique up-regulated CCOs in NS (Table S5). As we all know, ASS is the structural basis of proton transportation and energy generation in mitochondria (Klusch et al. 2017). CCO is the last respiratory complex of the electron transmit chain in mitochondria and is responsible for transferring electrons to the final acceptor oxygen in the respiratory metabolism (Dahan et al. 2014). The down-regulated ASS and CCO indicated the decrease of energy requirement in the developing leaf of NF under freezing stress, which might be one of the major reasons why it was more freezing-sensitive than the cultivar NS.

A possible freezing-responsive molecular network in developing winter rapeseed leaves
In this study, a freezing-responsive molecular network was produced in developing rapeseed leaves based on the results of transcriptome and proteome (Fig. 8). This network includes three important functional components: first, rapeseed plant leaves enhance Ca 2+ signal pathway by upregulation of CDPK, CIPK, MAPK, RLK, SnRK and SPK under freezing stress. Subsequently, these up-regulated protein kinases genes corporate with TFs to trigger the expression of downstream freezing-tolerant genes. Last, those activating genes cause some changes in metabolic pathways (enhanced plant hormones biosynthesis, ROS balancing and plant-pathogen interaction, accelerated biosynthesis of proteins and secondary metabolites, decreased photosynthesis), showing that the production of plant hormones, antioxidation enzymes and secondary metabolites might positively contribute to freezing tolerance.

Conclusions
After freezing treatment, the more tolerant rapeseed cultivar in response to freezing stress would fabricate stronger signals (Ca 2+ , ROS, ABA, Auxin, BR, Ethylene and JA), metabolic pathways (protein biosynthesis, plant-pathogen interaction and secondary metabolites synthesis) and ROS scavenging capacity, more osmolytes (soluble sugar and protein) and lower level of lipid peroxidation in leaves, as well as more stable cell morphology.