Comparative Transcriptomics and Proteomics Analyses of Leaves Reveals a Freezing Stress-Responsive Molecular Network in Winter Rapeseed (Brassica rapa L.)

Winter rapeseed is susceptible to low temperature during winter in Northwest China, which could lead to a severe reduction of crop production. The freezing temperature could stress the whole plant, especially the leaf, and ultimately harm the survival rate of winter rapeseed. However, the molecular mechanism underlying freezing tolerance is still unclear in winter rapeseed. In this study, a comprehensive investigation of winter rapeseed freezing tolerance was conducted at the levels of transcript, protein, and physiology and biochemistry, using a pair of freezing-sensitive and freezing-resistant cultivars NQF24 and 17NTS57. There were 4,319 unique differentially expressed genes (DEGs) and 137 unique differentially abundant proteins (DAPs) between two cultivars identified in leaf under freezing stress. Function enrichment analysis showed that most of the enriched DEGs and DAPs were involved in plant hormone signal transduction, alpha-linolenic/linoleic acid metabolism, peroxisome, glutathione metabolism, fatty acid degradation, and secondary metabolite biosynthesis pathways. Based on our findings, it was speculated that freezing tolerance formation is caused by increased signal transduction, enhanced biosynthesis of protein, secondary metabolites, and plant hormones, elevated energy supply, greater reactive oxygen species scavenging, and lower lipid peroxidation as well as stronger cell stability in leaf under freezing stress. These results provide a comprehensive profile of leaf response under freezing stress, which have potential to be used as selection indicators of breeding programs to improve freezing tolerance in rapeseed.


INTRODUCTION
Cold stress is categorized into chilling stress (0-15°C) and freezing stress (<0°C) in crops, which imposes deleterious effects on plant survival, growth, and development, restricts the geographical distribution of plant species, and influences crop production (Chinnusamy et al., 2007;Ding et al., 2019). Freezing stress is one of the most important environmental stresses in Northwest China from late autumn to early spring, which could not only reduce the yield but also harm the quality of economically important crops, especially for winter crops grown in high-latitude and high-altitude regions . Moreover, freezing stress damages the stability of proteins or protein complexes and reduces the activities of reactive oxygen species (ROS)-scavenging enzymes during seed development .
Winter rapeseed (Brassica napus L.) is a fundamental and nutritious oilseed crop that is widely planted in Northwest China. However, the growth and the development of winter rapeseed are vulnerable to several factors, such as extremely low temperature, low level of precipitation, and heavy evaporation in the region, which affect the normal overwintering and propagation of rapeseed, leading to a reduction of yield and quality. Hence, there is a high demand to develop rapeseed cultivars with strong freezing tolerance to withstand the stress for overwinter safely.
So far, abundant studies have been focused on cold tolerance response in model plants, such as Arabidopsis and rice Zhang et al., 2017). However, the molecular mechanism underlying freezing response in winter rapeseed is far away lagging behind. Recently, Zeng et al. (2018) suggested the impacts on the roots of two rapeseed cultivars with different freezing tolerance under freezing stress by isobaric tags for relative and absolute quantification (iTRAQ). Xu et al. (2018) and Pu et al. (2019) expounded the effects on the freezing tolerance formation in the leaves of two rapeseed cultivars with different freezing tolerance under cold stress by iTRAQ and transcriptome sequencing, respectively. Nevertheless, most studies only investigated single omics response based on traditional datadependent acquisition (DDA) models. In the recent years, dataindependent acquisition (DIA) technology is becoming popular due to its extremely high detection rate, enabling the identification and quantification of almost all protein molecules in complex samples. In addition, the quantitative analysis of protein DIA is also rated as one of the most noteworthy technologies by the journal Nature Methods (Egertson et al., 2013).
In the present study, we comprehensively analyzed the physiological and structural characteristics of a pair of resistant and susceptible winter rapeseed cultivars under freezing stress, which include determination of the antioxidant enzyme activities, net photosynthetic rate, and microscopic observation. Meanwhile, the combination of RNA sequencing-based transcriptomics and DIA-based proteomics was employed to identify differentially expressed genes' (DEGs) and differentially abundant proteins' (DAPs) response to freezing tolerance between two winter rapeseed cultivars. Subsequently, the main metabolic pathways and cellular processes that participated in freezing stress were explored. The findings will expand the understanding of the molecular regulatory network by which winter rapeseed responds to cold stress.

Plant Samples and Freezing Stress Treatments
The two extension winter rapeseed cultivars in Northwest China, freezing-resistant cultivar 17NTS57 (NS, with more than 85% overwinter survival rate at −26°C, dark brown seed, 33.02% oleic acid, and 17.29% linoleic acid) and freezing-sensitive cultivar NQF24 (NF, with 0% overwinter survival rate below −10°C, dark brown seed, 36.98% oleic acid, and 16.73% linoleic acid), used in this study were provided by Gansu Agricultural University (Mi et al., 2021). The seedlings of the two cultivars were grown in plastic pots (5 L) filled with soil and normally maintained. At four-leaf stage, the potted plants were divided into two groups. The treatment group (treatment, T1) was transferred to a chamber [−4°C, 60% humidity, 12/12 h (day/ night) photoperiod, and 350 μmol m −2 ·s −1 irradiance] for 12 h, while the control group (control, T0) was maintained in normal condition [22°C, 60% humidity, 12/12 h (day/night) photoperiod, and 350 μmol m −2 ·s −1 irradiance]. The second leaf was sampled from the control and the treated plants, frozen in liquid nitrogen immediately, and stored at −80°C for RNA and protein extraction and determination of physiological and biochemical characteristics. The experiment was conducted in three independent biological replicates.

Evaluation of Physiological and Biochemical Characteristics
The measurement of net photosynthetic rate was carried out as per Chu et al. (2015), with minor modifications. After freezing stress, the third fresh leaves of both cultivars were used for photosynthesis evaluation, utilizing a LI-COR 6400 portable gas analysis system (Lincoln, NE, United States) setting a built-in light source of 1,000 μmol photons m −2 ·s −1 . Each selected sample was measured at least four times at 25°C. All determinations were operated at 9:00 am to 11:00 am to minimize the error.
A total of 0.2 g rapeseed leaves was used to measure the physiological and biochemical characteristics, including ROS-scavenging enzyme activities, lipid peroxidation level, free proline content, and relative electrolytic leakage. The ROS-scavenging enzyme activities were measured using kits from Beijing Solarbio Science and Technology, following the manufacturer's protocols . The samples were homogenized using 1 ml 0.05 mol/L phosphate-buffered saline extraction buffer (pH 7.8). The supernatant was centrifuged at 4°C for 15 min at 8,000 × g and used for catalase (CAT), superoxide dismutase (SOD), and peroxidase (POD) activity analysis using Thermo Multiscan FC. The level of lipid peroxidation was evaluated by determining the malondialdehyde (MDA) content according to Wei et al. (2020). The samples were fully ground and extracted in a buffer of 0.67% thiobarbituric acid in 20% trichloroacetic acid (TCA) at 4°C and then centrifuged at 12,000 × g for 20 min. The absorbance of each sample supernatant was measured at 450, 532, and 600 nm.
The proline measurement was executed as described by Zhang et al. (2018); 0.2 g rapeseed leaves from treatments and controls was ground in liquid nitrogen and extracted by the proline assay kit (Solarbio). After centrifugation at 12,000 × g for 20 min, the supernatants were used for proline analysis with a spectrometer. The relative electrolyte leakage (REL) was measured as described by Huang et al. (2017b); 0.1 g rapeseed leaves from treatments and controls was cut into 2 cm × 2-cm Frontiers in Plant Science | www.frontiersin.org segments and were transferred into a 50-ml centrifuge tube that was filled with 15 ml of deionized water. The centrifuge tube with fragments was shaken for 24 h at room temperature, and the primary conductivity (EL1) was measured with a conductivity meter. Subsequently, the leaves in the tube were heated at 121°C for 20 min to disrupt the tissues and release all the electrolytes into the solution completely. After the samples had cooled to room temperature, the second conductivity (EL2) was measured. The relative EL was calculated by using the equation: EL = (EL1/EL2) × 100%. All measurements were performed with at least three biological replications.

Transmission Electron Microscope Observations
For transmission electron microscope (TEM) observations, fresh rapeseed leaves were cut into 2 mm × 2-mm sections and fixed in 2.5% glutaraldehyde solution for 16 h at 4°C, followed by fixation with 1% KMnO 4 solution for 2 h. All the fixed samples were dehydrated in a graded series of ethanol and embedded in Spurr resin. Eighty-nanometerthick ultrathin sections were finally stained with uranyl acetate for 15 min and lead citrate for 10 min. Observations were operated on an H-7650 electron microscope. The construction of high-quality TEM images was manipulated as described by Toyooka et al. (2014).
Rapeseed Leaf RNA Extraction, RNA Sequencing, and Quantitative Real-Time PCR Total RNAs from four samples containing three biological replicates were ground in liquid nitrogen and extracted using TRIzol Reagent (Tiangen Biotech, China) according to the manufacturer's instructions. The library construction and sequencing were performed by Gene Denovo Biotechnology Co. (Guangzhou, China) on an Illumina HiSeqTM 2500 platform. After the Illumina sequencing, three replicates of raw sequences for each sample were filtered to generate clean reads for subsequent analysis. The expressions of 12 candidate freezingresponsive genes and proteins were analyzed by qRT-PCR. The qRT-PCR was performed according to the protocol of Wei et al. (2020). All the primers are listed in Supplementary Table S1. The relative quantification (2 −ΔΔCt ) of gene expression was evaluated using comparative cycle threshold method, and each sample was replicated three times.

Filtering of Reads and Gene Expression Analysis
Clean data were obtained by removing the reads containing adapters, reads containing poly-N, and low-quality reads from raw data. The high-quality paired-end reads from each sample were mapped to rapeseed reference genome by TopHat v2.0.3.12 as illustrated (Kim et al., 2013b). The gene expression levels were calculated and normalized as fragments per kilobase per million mapped reads, which can be directly used for identifying DEGs in pair-wise comparisons (Trapnell et al., 2010). Genes with a p-value <0.001 and a value of |log 2 foldchange| ≥2 by the edgeR package were assigned as DEGs. 1 The sequenced transcriptome raw data have been deposited to the Sequence Read Archive at NCBI with accession number PRJNA685002.

Protein Extraction and Digestion
Protein extraction was conducted using the method of cold TCA/acetone precipitation according to Wei et al. (2020), with minor modifications. Plant leaves (0.2 g) from four samples containing three biological replicates were ground to powder in liquid nitrogen and then dissolved in 2 ml lysis buffer containing 7 M urea, 2% SDS, 1 × protease inhibitor cocktail (Roche Ltd. Basel, Switzerland), followed by sonication on ice for 20 min and centrifugation at 12,000 rpm for 10 min at 4°C. The supernatant was transferred into a fresh tube. For each sample, the proteins were precipitated with ice-cold acetone at −20°C overnight. The precipitations were cleaned with acetone three times and re-dissolved in 8 M urea and homogenized for 3 min in ice using an ultrasonic homogenizer. The homogenate was centrifuged at 12,000 rpm for 15 min at 4°C. The supernatant was collected, and the protein concentration was measured by the BCA Protein Assay Kit (Thermo Fisher Scientific, Waltham, MA). Ten, 50 μg proteins extracted from leaves was suspended in 50 μl of 50 mM ammonium bicarbonate, reduced by adding 10 ul of 100 mM dithiothreitol for 1 h at 55°C and alkylated by adding 5 ul of 20 mM iodoacetamide for 1 h at 37°C in the dark. Subsequently, protein samples were precipitated using 300 ul pre-chilled acetone at −20°C overnight. The precipitates were washed twice with cold acetone and resuspended in 50 mM ammonium bicarbonate. The samples were digested with sequencing-grade modified trypsin (Promega, Madison, WI) at an enzyme/substrate ratio of 1:50 (w/w) for 16 h at 37°C.

High-pH Reverse-Phase Separation and Liquid Chromatography With Tandem Mass Spectrometry Analysis
The digested peptides were separated on an Ultimate 3,000 system (Thermo Fisher Scientific, United States) with a high-pH reverse-phase XBridge BEH C18 column (4.6 mm × 250 mm, 5 μm) from Waters. The peptides were eluted at a flow rate of 1 ml/min. Buffer A consisted of 20 mM ammonium formate in water (pH 10.0), and buffer B consisted of 20 mM ammonium formate in 80% v/v acetonitrile (pH 10.0). High pH separation was performed using a linear gradient, starting from 5% B to 45% B in 40 min. For each sample, 10 fractions were collected and dried in vacuum centrifugation. All the samples were stored at −80°C until further analysis.
The peptides were dissolved in 30 μl buffer which consisted of 0.1% formic acid in water and 0.1% formic acid in acetonitrile (buffer D). Liquid chromatography with tandem mass spectrometry analysis was carried out on an Orbitrap Fusion Lumos mass spectrometer (Thermo Fisher Scientific, Germany) coupled with a Nano ACQUITY UPLC system (Waters Corporation, Milford, MA). Briefly, 10 μl of peptides was loaded onto a 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 column flow rate was maintained at 500 nl/min, and electrospray voltage of 2.1 kV was used.
For library generation, the Orbitrap Fusion Lumos was operated in positive mode to acquire data in DDA mode. Mass spectrometry (MS) spectra (350-1,500 m/z) were collected at 1.2 × 10 5 resolution to reach an automatic gain control (AGC) target of 4.0 × 10 5 . The maximum injection time was set to 50 ms. The precursor ions were used for MS/MS detection using a normalized 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 reach an AGC target of 5.0 × 10 4 with a maximum injection time of 35 ms. For data acquisition, DIA was used. MS spectra (350-1,500 m/z) were collected at 1.2 × 10 5 resolution to reach an AGC target of 5.0 × 10 5 . The maximum injection time was set to 50 ms. MS/MS spectra were collected at 3.0 × 10 4 resolution to reach an AGC target of 1.0 × 10 6 with a collision energy of 32. A total of 60 segments were selected for MS/MS using normalized collision energy of 32. The acquisition window covered 1 m/z through 60 consecutive isolation windows.

Mass Spectra Data and Protein Quantification
Spectronaut Pulsar 11.0 (Biognosys AG) software (Bruderer et al., 2015) was used to analyze and process the DDA MS/MS data with the following settings: trypsin as the digestion enzyme, carbamidomethylation (C) as the fixed modification, oxidation (M) as the variable modifications, 20 ppm as precursor mass tolerance, and 0.05 Da as fragment mass tolerance. The MS/ MS data were searched from the database of B. napus, which was downloaded from the NCBI (Brassica napus v2.0, 123,465 entries; Supplementary Table S1). False discovery rate of peptides was dynamically set at 1%, which was calculated by a reverse search sequence. Subsequently, the raw data of every sample with the spectra library based on DIA MS/MS data were analyzed by Spectronaut Pulsar 11.0 with the default parameters depending on the previously constructed DDA protein dataset, which included protein sequences (19,394), precursors (49,118), and peptides (39,809). All protein quantification normalizations were performed by local normalization with the Pulsar software. The mass spectrometry data have been submitted into a public iProX database  with accession number IPX0002682000. Proteins with a value of p < 0.001, and peptides with a value of p ≥ 2 and a value of |log 2 foldchange| ≥1 by the edgeR package were filtered as DAPs (see Footnote 1).

Bioinformatics and Statistical Analysis
The correlation coefficient between three replicates was calculated to evaluate the repeatability of the experimental results between samples. Principal component analysis (PCA) was performed to reveal the structure/relationship of samples by R package gmodels (see Footnote 1). The proteins and genes were annotated against Gene Ontology (GO) database and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases (Binns et al., 2009;Kanehisa et al., 2016). 2,3 Enrichment analysis was performed based on GO and KEGG databases. Physiological data were subjected to one-way analyses of variance, and the mean differences of three replicates were compared by least significant difference test.

Physiological and Biochemical Responses to Freezing Stress
After the freezing treatment, NF plants appeared more severely wilting than NS plants ( Figure 1A). Freezing stress caused a significant decrease in the net photosynthetic rates of both rapeseed cultivars in comparison with the normal condition, whereas under freezing stress, the net photosynthetic rate of rapeseed cultivar NS was significantly (p < 0.05) higher than that of cultivar NF. Interestingly, if we do the relative values on net photosynthesis, we discovered that the net photosynthetic rate decreased for approximately 50% in both cases ( Figure 1B). Through TEM analysis, the chloroplasts in mesophyll cells were found to be disintegrated harshly in rapeseed cultivar NF after freezing stress, whereas more plastoglobuli were found in cells of rapeseed cultivar NS under freezing stress ( Figure 1C). These results indicated that the freezing stress exerted greater damage to the leaf of rapeseed cultivar NF than to that of cultivar NS.
The freezing stress caused a significant increase in the enzyme activities of CAT and POD in the leaf of both rapeseed cultivars. However, the enzyme activities of CAT, SOD, and POD in the leaf of cultivar NS maintained a higher level than those of cultivar NF under freezing stress (Figures 2A-C). Furthermore, the level of plasma membrane lipid peroxidation in the leaf of both cultivars was evaluated by measuring MDA and REL, respectively. Both MDA content and REL level were significantly (p < 0.05) increased in the leaf of both cultivars under freezing stress, but the levels in cultivar NF were significantly (p < 0.05) higher than those in cultivar NS (Figures 2D-F), indicating that the freezing treatment resulted in more severe plasma membrane damage in the leaf of cultivar NF than in that of cultivar NS. Moreover, free proline was significantly (p < 0.05) increased only in cultivar NS under freezing stress ( Figure 2E). These results suggested that a high level of ROS scavenging capacity and a low level of lipid peroxidation might contribute to enhance freezing tolerance in winter rapeseed cultivar NS under freezing stress.

Quality Control Analysis of Rapeseed Transcriptome and Proteome
Overall, transcripts were detected for 80.0% of the proteins ( Figure 3A). The transcriptome analysis detected a total of 81,833 genes, including 75,034 known genes and 6,799 new genes in both rapeseed cultivars ( Figure 3A; Supplementary Table S1).
For the proteome analysis, a total of 30,471 unique peptides were identified from 38,159 spectra, which corresponded to 13,921 proteins and were encoded by 11,141 genes in transcriptome ( Figure 3A; Supplementary Tables S1,S2). The PCA indicated that three biological replicates of NFT0, NFT1, NST0, and NST1 had good conformity (Figures 3B,C), and all biological replicates showed a high overlap in transcriptome and proteome, respectively (Supplementary Figure S1). In addition, a repeatability analysis between three biological replicates of NFT0, NFT1, NST0, and NST1 samples showed that their correlation coefficient was greater than 0.9 (Supplementary  Figures S2, S3). These results indicated a high level of reliability of the transcriptome and proteome analyses.

Quantitative Transcriptomic and Proteomic Analysis of Rapeseed Leaves Under Freezing Stress
To identify DEGs and DAPs that were induced by freezing treatment, the changes of transcription level and protein relative abundance in leaves of both cultivars were analyzed through comparison of T1 to T0. Compared to the corresponding controls, a total of 9,669 (up, 4,669; down, 5,000) genes and 10,544 (up, 4,049; down, 6,495) genes were identified to be differentially expressed in the leaves of rapeseed cultivars NS and NF under freezing stress, respectively ( Figure 4A; Supplementary Table S3). In comparison with controls, a total of 392 and 166 proteins were found to be accumulated in abundance in the leaves of cultivars NS and NF under freezing stress, respectively, whereas 161 and 190 proteins were reduced in abundance in cultivars NS and NF under freezing stress, respectively ( Figure 4A; Supplementary Table S4). It is worth noting that the number of DEGs and DAPs was significantly changed in both cultivars, implying a greatly diverse dynamics of gene or protein expression regulation in both cultivars. In addition, more accumulated proteins were found in rapeseed cultivar NS than in NF, and a small number of DAPs were shared between both cultivars ( Figure 4B).

Comparative Analysis Between Protein Abundance and Gene Expression Levels
The relationship between the levels of transcription and protein would be evaluated. DAPs and DEGs were compared between both rapeseed cultivars at different freezing treatment time points (T0 and T1), respectively. A total of 4,625 shared DEGs and 40 shared DAPs were jointly owned in NFT0_NST0 and NFT1_NST1. Among the shared DEGs, 2,608 were up-regulated and 2,017 were down-regulated in NFT0_NST0, while 2,706 were up-regulated and 1,919 were down-regulated in NFT1_NST1 ( Figure 4C; Supplementary Table S5). Among the shared DAPs, 25 were accumulated in abundance and 15 were reduced in both NFT0_NST0 and NFT1_NST1, and the accumulation trend of all proteins is exactly consistent (Figure 4D; Supplementary Table S5). Furthermore, 4,319 DEGs and 137 DAPs were identified only in NFT1_NST1 but not in NFT0_ NST0 (Figures 4C,D; Supplementary Table S6), which were the molecular basis for the difference in freezing tolerance between both cultivars, and they were regarded as unique genes/proteins. In addition, 21 DAPs corresponding to DEGs were found in the common portion in transcriptome and proteome changes between both cultivars during freezing stress ( Figure 4E; Supplementary Table S7), and the concordance changes between the gene and its encoded protein were analyzed using Pearson's correlation test. The results indicated a significant positive correlation between the fold changes of DEGs and the corresponding DAPs, with a Pearson of 0.415 in NFT1_NST1 ( Figure 4F). These results suggested that DEGs and their encoded DAPs play a collaborative role in the freezing tolerance of winter rapeseed.

Differentially Expressed Transcription Factors Under Freezing Stress
A total of 883 and 962 differentially expressed transcription factors (TFs) were identified as corresponding to the leaf transcriptome of cultivars NS and NF under freezing stress, respectively (Supplementary Table S9). In both cultivars, most of the up-regulated TFs belonged to ethylene response factor (AP2/ERF), NAC domain-containing protein (NAC), WRKY transcription factor (WRKY), heat stress transcription factor (HSF), TAZ domain-containing protein (TAZ), protein TIFY (TIFY), and zinc finger protein (C2C2) families, while most of the down-regulated TFs belonged to the transcription factor bHLH (bHLH), zinc-finger homeodomain protein (zf-HD), squamosa promoter-binding-like protein (SBP), transcription repressor OFP (OFP), and trihelix transcription factor (Trihelix) families (Figures 6A,B). Contrary to the findings of transcriptome, we did not find many differentially expressed TFs in proteome. There were only two and six TFs identified in the leaf proteome of cultivars NS and NF under freezing stress, respectively, and seven out of eight TFs were reduced in abundance. These results indicated the complicated dynamic changes between gene expression regulation and protein expression regulation. In addition, a TF analysis of unique DEGs in NST1_NFT1 was executed. Generally, in line with the above-mentioned studies, it was shown that most TFs classified into the AP2/ERF, NAC, WRKY, HSF, and TIFY families were up-regulated in cultivar NS under freezing stress compared to cultivar NF ( Figure 6C; Supplementary Table S9). Our results suggested that these TFs in tolerant winter rapeseed cultivar NS are more susceptible to activation under freezing stress, thereby regulating the expression of downstream genes in response to freezing stress.

RNA-Seq and DIA Quantitation Validation by qRT-PCR
The qRT-PCR analysis was used to validate the reliability of our transcriptome and proteome data. Based on high fold change, 12 freezing-responsive DEGs or DAPs were selected as targets. Among them, 11 out of 12 were found to be consistent at both the mRNA and RNA-Seq or protein levels under freezing stress in rapeseed cultivar NS compared to cultivar NF. Besides that, one gene showed an inconsistency between the mRNA and RNA-Seq levels (Figure 7). These results indicated that the qRT-PCR expression patterns were in good agreement with the RNA-Seq and DIA quantitation, and the RNA-Seq and DIA quantitation results were reliable in the present study.

DISCUSSION
It is generally believed that chilling and freezing stresses lead to a considerable decline in productivity of winter rapeseed, which is common in many winter rapeseed production areas around China (Zhu, 2016). Leaves, as the primary organ for photosynthetic accumulation products during cold acclimation, are crucial for overwinter survival rate and the reconstruction of above-ground organ after the revival stage once it perceives cold signaling, followed rapidly with response to cold stress, which would increase the cold tolerance and maintain its function. Some previous studies have also investigated the effects of freezing stress on winter rapeseed leaf tolerance formation by proteomics Zeng et al., 2018), transcriptome, and physiology (Pu et al., 2019). However, the current unraveling of the intricate mechanism of the winter rapeseed in response to freezing stress remains elusive. The NS winter rapeseed can survive normally in winter in most northwestern areas from 38° to 42° north latitude, and its tolerance at low temperature environments is better than those of other types of winter rapeseed, which contain abundant cold-resistant genes. Therefore, in the present study, the thorough effects of cold stress on cold tolerance formation were investigated through the analysis of cold-tolerant NS and sensitive NF rapeseed leaf at the levels of transcript, protein, and physiology and biochemistry. It will be of great importance to elucidate the exhaustive mechanism of winter rapeseed for plant molecular breeding under cold stress.

Physiological, Ultrastructural, and Photosynthetic Changes Between Both Winter Rapeseed Cultivars Under Freezing Stress
Cold stress can trigger a series of comprehensive physiological and biochemical events, leading to the synthesis of many substances or protective proteins in plants (Kaplan et al., 2007). Freezing stress provokes intracellular ice formation, which results in membrane lipid peroxidation and cell structural damage . For example, Gao et al. (2019) suggested that the level of lipid peroxidation in jojoba was significantly increased after cold treatment. Similarly, low electrolyte leakage and bad cell membrane integrity were exhibited in cold-stressed Arabidopsis roots (Deng et al., 2015). In this study, freezing stress gives rise to an increase in MDA content and REL level in the leaves of both cultivars, but the level of lipid peroxidation in cultivar NF was significantly higher than that in cultivar NS (Figures 2D,F). Coincidentally, the TEM micrographs exhibited a stronger cell ultrastructure stability in cultivar NS plants, compared with NF, after the freezing treatment ( Figure 1C). Additionally, consistent with proteomic results, some beta tubulins said to contribute to cytoskeleton modifications were found to be accumulated in abundance (Supplementary Table S6), which implied that freezing stress induced rearrangements in the tubulin cytoskeleton, leading to the formation of a more stable cell structure (Livanos et al., 2014). In the past two decades, it has become increasingly evident that ROS plays a pivotal signal-regulating role in response to diverse 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). In the present study, ROS-scavenging enzyme activities were significantly increased in both cultivars under freezing stress, but the levels of enzyme activities in cultivar NS were higher than those in cultivar NF, and the enzyme activities of SOD were significantly different between NS and NF (Figure 2). Overall, these findings are in accordance with the findings reported by Zeng et al. (2018). Taken together, the enhanced ROS scavenging and the lower level of lipid peroxidation accompanied with more stable cell morphology in the leaves of winter rapeseed under freezing treatment would maintain a normal redox environment for leaf growth and development, thus leading to the stronger freezing tolerance in NS plants.

Comparison of Effects of Freezing Stress on Freezing Tolerance Formation Between Both Winter Rapeseed Cultivars
A multitude of transcription factors have been reported in numerous higher plants, which monitor the expression of C-repeat binding transcription factor (CBF) by combining the corresponding cis-elements in their promoters under cold stress. After exposure to freezing stress, the CBF proteins rapidly recognize the promoter regions of downstream coldregulated genes to activate their expression that enhances freezing tolerance (Jia et al., 2016;Zhao et al., 2016). For instance, bHLH transcription factor ICE1, calmodulin-binding transcription activator 3, brassinazole-resistant1 (BZR1), heat shock transcription factor C1 (HSFC1), and MYB88 /MYB124 transcription factors positively regulate the CBF expression at the transcriptional level, which contributes to cold stress (Kim et al., 2013a;Park et al., 2015;Xie et al., 2018). In contrast, MYB15, phytochrome interacting factor3 (PIF3), PIF4, PIF7, and ethylene sensitivity3 negatively modulate the transcriptional repression of CBF (Shi et al., 2012;Jiang et al., 2017;Kim et al., 2017). In the present study, most of the AP2/ERF, NAC, WRKY, HSF, TAZ, TIFY, and C2C2 transcription factors were up-regulated in both rapeseed cultivars, while the transcription factors bHLH, zf-HD, SBP, OFP, and Trihelix were down-regulated (Figures 6A,B). Interestingly, these differentially expressed TFs encoded by genes were not found at the protein level, implying that those TFs might regulate the gene expression just at the transcription level, contributing to the freezing tolerance of winter rapeseed under freezing stress. Furthermore, over the past few years, several protein kinases including mitogen-activated protein kinases (MAPKs), Ca 2+ -dependent protein kinases (CDPKs), calcium/calmodulin-regulated receptor-like kinases (CRLKs), calcineurin-B-like interacting protein kinases (CIPKs), and receptor-like kinases (RLKs) have also been characterized to be crucial regulators of cold stress responses in plant Zhao et al., 2017;Ding et al., 2019). In addition, a putative cold sensor chillingtolerance divergence1 was reported to mediate a cold-sensing calcium channel in rice, leading to the activation of coldregulated genes Shi et al., 2018). A similar pattern of results was obtained in the present study: a total of nine MAPKs, five CDPKs, six CIPKs, and some RLKs/ CRLKs were up-regulated in cultivar NS under freezing stress compared to cultivar NF (Supplementary Table S6). These results suggested that the plasma membrane-located sensors might perceive the signal by the freezing treatment and initiate the Ca 2+ signaling pathway in rapeseed. Following Ca 2+ influx to activate the plasma membrane, some protein kinases interact with TFs and trigger downstream freezingresponsive gene expressions. Phytohormones, including abscisic acid (ABA), jasmonic acid (JA), ethylene (ETH), gibberellin, auxin (IAA), cytokinin, salicylic acid (SA), and brassinosteroid (BR), combine endogenous substances with environmental signals to regulate plant growth and development as well as defense (Huang et al., 2017a). In many instances, plants respond to environmental stresses by producing amounts of ABA, BR, ETH, JA, and SA (Huang et al., 2016;Hu et al., 2017;Wu et al., 2019). Shibasaki et al. (2009) found that cold stress inhibited the expression and transport of auxin-responsive related genes in root. Similar results were obtained in this study; a total of 87 unique DEGs in NST1_NFT1 were enriched in plant hormone signal transduction pathway (Figure 5A; Supplementary Table S6). Among them, 10 DEGs were up-regulated and involved in ABA signaling, encoding abscisic acid-insensitive proteins, serine/threonine-protein kinase SnRK2 (SnRK), and protein phosphatase 2C (PP2C), whereas three DEGs encoding abscisic acid receptors were down-regulated; four DEGs were up-regulated and engaged in ETH signaling, encoding ethylene-responsive transcription factor2 and ethylene response sensor2, while six DEGs encoding ethylene-responsive transcription factors 1/15 and ethylene-insensitive protein were down-regulated. Interestingly, we found that two key enzymes of ethylene biosynthesis, 1-aminocyclopropane-1carboxylate synthases and 1-aminocyclopropane-1-carboxylate oxidases, were up-regulated. Therefore, we speculate that there may be at least two different ABA and ethylene metabolic pathways associated with freezing tolerance, and one of them was enhanced. Furthermore, 19 DEGs were up-regulated and participated in JA signaling, encoding jasmonic acid synthetase and TIFY proteins; 28 DEGs (19 up-regulation and nine down-regulation) were found and contributed to auxin signaling, encoding auxin-responsive proteins, auxin transporter proteins (ATP), and GH3 auxinresponsive promoter. These are consistent with what has been found in previous research, which suggested that high auxin expression levels could improve freezing tolerance (Pu et al., 2019). Moreover, one DEG was up-regulated, encoding protein BZR1, which is a plant-specific positive regulator of BR signaling (Yin et al., 2005;. BRs are perceived by brassinosteroid insensitive1 (BRI1), which cooperates with BRI1-associated receptor kinase1 and positively modulates BZR1 (She et al., 2011;Fang et al., 2019). Fortunately, a leucine-rich repeat receptor-like kinase was found in this study, encoded by BRI1 gene, which was up-regulated. Therefore, we speculated that BZR1 and BRI1 work together to actively mediate BR signaling in response to freezing stress. Apart from that, a large number of studies have shown that BRs can enhance photosynthesis and the activities of antioxidant enzymes (Zhou et al., 2014;Xia et al., 2018). A similar conclusion was reached by the determination of photosynthesis and ROS-scavenging enzymes (Figures 1B, 2). Simultaneously, other DEGs (12 up-regulation and three down-regulation) representing 13 transcription regulators, one gibberellin receptor, and one xyloglucan endotransglucosylase/hydrolase, which plays a crucial role in plant cell wall remodeling (Witasari et al., 2019), were also identified in our study. Taken together, these results indicated that phytohormones play critical roles in response to the freezing stress of rapeseed.
Plant lipoxygenases are a kind of fatty acid dioxygenase with diverse functions and catalyze the peroxidation of polyunsaturated fatty acids such as linolenic and linoleic acids. It was reported that biotic and abiotic stresses cause the release of α-linolenic acid, which fascinated particular attention because it represents an essential precursor of jasmonic acid biosynthesis (Yan et al., 2013;Zhao et al., 2014). The synthesis of jasmonic acid begins with the peroxidation of α-linolenic acid by lipoxygenase, the product of which is converted to 12,13-epoxyoctadecatrienoic acid (12,13-EOT) by the allene oxide synthase (AOS). Subsequently, the 12,13-EOT is processed into 12-oxophytodienoic acid (OPDA) by the allene oxide cyclase (AOC), and OPDA is reduced by the OPDA reductase. In this study, all specific DEGs in NST1_NFT1 involved in the alpha-linolenic acid metabolism were found to be up-regulated. Among them, there are nine lipoxygenases, four AOC, two AOS, and seven 12-oxophytodienoate reductases ( Figure 5A). Another promising finding was that the alphalinolenic acid metabolism was also significantly enriched by DAPs, and all DAPs were accumulated in abundance ( Figure 5B; Supplementary Table S6). Beyond those and combined with JA signaling resulting from DEGs, these results implied that enhanced jasmonic acid synthesis might regulate the expression of downstream freezing-responsive genes under freezing stress.
Peroxisomes can generate some retrograde signals such as ROS and many metabolites that are important for stress responses (Ng et al., 2014). Glutathione is well known for its ROS scavenging function in biotic stress . In this study, six unique DAPs in NST1_NFT1 belonging to the peroxisome pathway were significantly enriched after the freezing treatment ( Figure 5B; Supplementary Table S6). Under freezing stress, four catalases with accumulated abundance in cultivar NS compared to cultivar NF were consistent with the determination of CAT enzyme activity ( Figure 2A); nevertheless, two peroxidases with reduced abundance in cultivar NS compared to cultivar NF were contrary to the determination of POD enzyme activity ( Figure 2C; Supplementary Table S6). This divergence might be attributed to posttranslational modifications (Park et al., 2018). A further novel finding is that the glutathione metabolism pathway was significantly enriched by DAPs, and most of the glutathione S-transferases (GST) enriched in this pathway were accumulated in abundance. In addition, one glutathione peroxidase (GPX) gene was also found to be up-regulated. Our results demonstrated that ROS partook in response to rapeseed freezing tolerance.
Plant secondary metabolites often have no fundamental role in the maintenance of plant life processes, but they are important for the plant to interrelate to its environment for adaptation and defense (Ramakrishna and Ravishankar, 2011). Phenylpropanoid compounds congregate a large family of secondary metabolites that originate from phenylalanine and are beneficial to defense against various stresses in plants (Liu et al., 2015). In the present study, there were five unique abundant accumulated DAPs in NST1_NFT1 encoded by caffeic acid O-methyltransferase, shikimate O-hydroxycinnamoyltransferase, and 3-dehydroquinate dehydratase/ shikimate dehydrogenase and implicated in phenylpropanoid biosynthesis ( Figure 5B; Supplementary Table S6). It is in agreement with those reported by Zeng et al. (2018) and Pu et al. (2019). In addition, we also found some unique up-regulated DEGs encoding S-adenosylmethionine-dependent methyltransferase, which are involved in the secondary metabolism of many plant species (Ranjan et al., 2020). These results suggested that the accumulated secondary metabolites contribute to the freezingstressed rapeseed cultivar NS under freezing stress compared to the cultivar NF.
It has been reported that heat shock proteins (HSPs), regarded as the target genes of HSFs, are correlated with cold stress tolerance in grasses (Wang et al., 2016). Similarly, Wang et al. indicated that a large amount of HSPs and chaperone proteins (CPs) was induced by cold stress (Wang et al., 2004). Scarpeci et al. (2008) suggested that ROS, through a direct interaction with HSFs, mediate cellular signals during abiotic stress. In the present study, some up-regulated HSFs, HSPs, and CPs are accompanied with high levels of antioxidation enzyme activities (Figure 2), thus implying that the HSFs might interact with ROS signals to regulate the expression of HSPs and CPs during cold stress. Except the HSPs and CPs, some other unique DEGs encoding cell division control protein, endoplasmin, and protein disulfide isomerase involved in protein synthesis were discovered in NST1_NFT1, and they were up-regulated. In addition, we found some unique up-regulated DEGs in NST1_NFT1, which encode some mitochondria-localized ATP synthase subunit (ASS), cytochrome c oxidase (CCO), and persulfide dioxygenase (PDO). It is known that ASS is the structural basis of proton transportation and energy generation in the mitochondria (Klusch et al., 2017). CCO is the last respiratory complex of the electron transport chain in mitochondria and is responsible for transferring electrons to the final acceptor oxygen in respiratory metabolism (Dahan et al., 2014). The up-regulated ASS and CCO insinuated the increase of energy requirement during leaf development under freezing stress.

A Possible Freezing Stress-Responsive Molecular Network in Winter Rapeseed
In the present study, based on the transcriptomics and proteomics, a freezing stress-responsive molecular network was generated from the most of significantly enriched freezing-responsive genes and proteins in developing rapeseed leaves (Figure 8). This network covers five foremost functional components, including balancing between ROS production and scavenging, increased signal transduction and transcriptional regulation, enhanced phytohormone (ETH, ABA, IAA, JA, and BR) biosynthesis, accelerated biosynthesis of proteins and secondary metabolites, abundant energy supply, and impaired photosynthesis. Such a molecular network allows us to further understand the principle of freezing tolerance in freezing-treated FIGURE 8 | A presumed model for freezing tolerance mechanism in developing rapeseed leaves under freezing stress. The freezing-responsive genes and proteins up-regulated are marked by "↑," and those down-regulated are marked by "↓." The enhanced metabolic pathways and cellular processes are marked with the red "↑," while the reduced ones are marked with the blue "↓." Dashed arrows indicate possible regulation. rapeseed leaves, and some of them can be used as selection index in freezing tolerance breeding program in rapeseed.
When rapeseed plants are exposed to freezing stress, some plasma membrane-located sensors mobilize the signals such as Ca 2+ , ROS, and plant hormones by signal transduction pathways; following these signal influxes, their related protein kinases are activated to interact with TFs and trigger the expression of downstream freezing-tolerant genes. Under persistent freezing stress, the balance between ROS production and scavenging in the cells of rapeseed plants is disrupted. To eliminate the oxidative damage resulting from increased ROS production in mesophyll cells, the rapeseed plant leaves would strive to maintain cellular redox homeostasis by elevating the biosynthesis of antioxidation enzymes (GST, CAT, SOD, and GPX; Figure 2); under freezing stress, the rapeseed plant leaves increase the expression of key enzymes (ASS, CCO, and PDO) in the energy metabolism pathway in response to freezing stress; under freezing stress, the rapeseed plant leaves also heighten the Ca 2+ signal pathway by up-regulation of MAPK, CDPK, SnRK2, CIPK, RLK, and CRLK. Subsequently, these up-regulated protein kinase genes incorporate with TFs to cause changes in metabolism pathways (enhanced synthesis of proteins and secondary metabolites and increased plant hormone synthesis), suggesting that the production of plant hormones and secondary metabolites may be positively contributing to freezing tolerance. However, these findings need to be demonstrated by further studies.
Furthermore, freezing-tolerant cultivar NS possesses higher ROS scavenging capacity, more stable cell morphology, more ETH, JA, IAA, and ABA production, and stronger photosynthesis and energy supply than freezing-sensitive cultivar NF under freezing stress, which may be the major reasons why cultivar NS is more freezing tolerant than cultivar NF.

CONCLUSION
Under freezing stress, the rapeseed cultivar more tolerant to freezing stress would produce stronger signals (Ca 2+ -mediated, ROS, ABA, ETH, JA, IAA, and BR), metabolic pathways (photosynthesis, protein biosynthesis, and secondary metabolite synthesis), and ROS scavenging capacity, more stable cell morphology, and a lower level of lipid peroxidation in leaves.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: The sequenced transcriptome raw data have been deposited to the SRA at NCBI with the accession number of PRJNA685002. The mass spectrometry data have been submitted into a public iProX database  with the accession number of IPX0002682000.

AUTHOR CONTRIBUTIONS
The work presented here was carried out in collaboration among all authors. JW executed the omics data analysis and wrote the manuscript. GZ cultivated all plant samples. XY polished and revised this paper. SL performed the analysis of physiological and biochemical characteristics. XD and XC participated in the determination of physiological and biochemical characteristics. HL and XF were engaged in the measurement of photosynthesis. JJ and WM carried out the protein data analysis. ZL designed the experiments. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We are grateful to Guangzhou Genedenovo Biotechnology Co., Ltd., for assisting in the sequencing and bioinformatics analysis.