Characterization of cold stress responses in different rapeseed ecotypes based on metabolomics and transcriptomics analyses

The winter oilseed ecotype is more tolerant to low temperature than the spring ecotype. Transcriptome and metabolome analyses of leaf samples of five spring Brassica napus L. (B. napus) ecotype lines and five winter B. napus ecotype lines treated at 4 °C and 28 °C were performed. A total of 25,460 differentially expressed genes (DEGs) of the spring oilseed ecotype and 28,512 DEGs of the winter oilseed ecotype were identified after cold stress; there were 41 differentially expressed metabolites (DEMs) in the spring and 47 in the winter oilseed ecotypes. Moreover, more than 46.2% DEGs were commonly detected in both ecotypes, and the extent of the changes were much more pronounced in the winter than spring ecotype. By contrast, only six DEMs were detected in both the spring and winter oilseed ecotypes. Eighty-one DEMs mainly belonged to primary metabolites, including amino acids, organic acids and sugars. The large number of specific genes and metabolites emphasizes the complex regulatory mechanisms involved in the cold stress response in oilseed rape. Furthermore, these data suggest that lipid, ABA, secondary metabolism, signal transduction and transcription factors may play distinct roles in the spring and winter ecotypes in response to cold stress. Differences in gene expression and metabolite levels after cold stress treatment may have contributed to the cold tolerance of the different oilseed ecotypes.


INTRODUCTION
As one of the most important environmental stresses, cold stress significantly affects plant growth, development and distribution, and it can be classified as chilling (<20 C) and freezing (<0 C) temperatures (Chinnusamy, Zhu & Zhu, 2007;Miura & Furumoto, 2013). Cell membrane stiffness, protein complex instability and photosynthesis weakness are detected under chilling stress, and much more serious injuries occur under freezing stress conditions. Many physiology and biochemistry processes can be influenced by cold stress, which is dependent on the plant species, stage of plant development, and length of the stress period. Many plant species, including Brassica napus L. (B. napus), have the ability to acclimate to the cold. However, this ability can be increased after exposure to chilling stress, and the capacity of acclimation varies greatly among species. Cold acclimation mechanisms have been fully explored and have expanded our understanding of Arabidopsis thaliana (A. thaliana), while information regarding the molecular basis of this process in B. napus is largely lacking despite their similar genetics.
In addition to the very sophisticated gene regulatory networks, plants have evolved various complex metabolite response patterns in the low temperature stress response (Kaplan et al., 2004). Increasing evidence has shown that gene expression and biochemical pathways are altered under cold stress conditions. Sugars, lipids and amino acids, for example, have been widely reported to be involved in the cold stress response, not only for protein synthesis process but as precursors of other key metabolites in the response to cold stress (Cook et al., 2004). As key functions of the metabolite response to cold stress, many more kinds of metabolites that are responsive to cold stress must be characterized in further studies.
Transcriptome sequencing technology has matured and can provide very efficient and high-throughput screening of low temperature response genes. In previous studies, genes involved in the cell wall and membrane stabilization, calcium signaling, hormones, TFs and mitogen-activated protein kinase (MAPK) pathways, and soluble sugar and protein biosynthesis and metabolic pathways were significantly altered after cold stress using RNA-Seq technologies (Wang et al., 2013a(Wang et al., , 2013bChen et al., 2014;Ren et al., 2014;Fu et al., 2016;Gong et al., 2018). Recent technological advances have made it possible to explore the molecular mechanisms underlying the response to low temperature stress using transcriptome and metabolomics techniques (Kaplan et al., 2007;Maruyama et al., 2014;Wu et al., 2016;Zhang et al., 2016). By integrating transcriptome and metabolome data, a broader and more comprehensive understanding of genes and metabolites co-regulation of low temperature stress in Arabidopsis has been performed (Kaplan et al., 2007). In rice, genes encoding proteins involved in sucrose-starch metabolism and the glyoxylate cycle were upregulated, and these changes were correlated with the accumulation of glucose, fructose and sucrose in rice after exposure to cold or dehydration, as revealed by joint metabolite, phytohormone, and gene transcript analyses (Maruyama et al., 2014). In A. thaliana, flavonoids are regarded as determinants of freezing tolerance and cold acclimation using 20 mutant lines based on qRT-PCR, LC-MS and GC-MS analysis (Schulz et al., 2016). Carboxylates, amino acids, carbohydrates and sugar alcohols were also made contributing to reduce freezing tolerance in the pgm mutant (Hoermiller et al., 2017). Joint RNA-Seq and metabolomics analyses have revealed a reactive oxygen species-dominated dynamic model response to chilling stress in two subspecies of Asian cultivated rice (Oryza sativa) with significant divergence in chilling tolerance (Zhang et al., 2016). Integrative analysis of transcripts and metabolites has revealed that many transcripts and metabolites are rearranged to cope with low temperature stress in Dendrobium officinale (Wu et al., 2016).
As the third largest oil crop in the world, rapeseed is the most important sources of vegetable oil in China. Rapeseed is the only overwintering oil crop in China. The research and breeding of cold-resistant varieties play an important role in the safety of edible oil in China. Until now, many valuable studies have been performed to examine the morphology, physiology and molecular biology of cold resistance (Burbulis et al., 2010;Chen et al., 2011), but few studies have addressed the cold resistance regulatory network mechanisms at the transcriptome and metabolome levels in rapeseed. ABA and IP3/Ca 2+ signal transduction pathways have been identified as key actors in response to low temperature in rapeseed (Xian et al., 2017). To our knowledge, the elucidation of candidate genes or the metabolite response to low temperature in rapeseed using metabolomics and transcriptomics analysis is still very lacking.
In this study, we aim to explore the molecular mechanisms in different rapeseed ecotypes at transcriptome and metabolome levels from a wide perspective and seek clues linking genes and metabolites in response to low temperature stress. Each rapeseed ecotype (spring and winter) containing five lines were used in our study. The spring rapeseed ecotype was more sensitive to low temperature than the winter rapeseed ecotype. Many cold response genes and metabolites were detected, and common ones had greater degrees of change in the winter than in the spring rapeseed ecotypes.

Plant materials and cold treatments
Five winter (N5, N7, N36, N49 and N60) and five spring ecotypes (M10, M11, N240, N244 and N249) of B. napus were used in this study. Seeds of each line were germinated and cultivated in a petri dish under normal conditions. After 10 days, the plants were moved to plastic pots in a greenhouse with a long-day photoperiod (28 C with 16 h light and 24 C with 8 h dark) until the plants grew to the fiveor six-leaf stage. Half plants of each line were grown under normal conditions and used as control plants. The remaining plants of each line were moved to 4 C for 12 h with the same photoperiod. Leaves of control and treatment plants were harvested at 0 and 12 h after cold stress. All leaves were immediately frozen in liquid nitrogen and stored at −80 C until use. Every sample consisted of leaves from more than five plants grown under the same conditions. Three biological replicates were performed, and 12 samples including control samples for winter (28W) and spring (28S) ecotypes, cold stress samples for winter (4W) and spring (4S) ecotypes were obtained and used for RNA-Seq and metabolic profiling.

RNA extraction and RNA-Seq
Total RNA of all samples was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. An Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA) and a spectrophotometer (NanoDrop, Wilmington, DE, USA) were used to detect the integrity, quality and quantity of total RNA. Then, high-quality RNA was sent to Genedenovo Biotechnology Co., Ltd. (Guangzhou, China) and sequenced on the Illumina sequencing platform. Briefly, mRNA was extracted using dynabeads oligo (dT) and fragmented by fragmentation buffer. Double-stranded cDNAs were synthesized using reverse-transcriptase and random hexamer primers, and the cDNA fragments were purified using a QIA quick PCR extraction kit. These purified fragments were washed with EB buffer for end reparation of poly (A) addition and then ligated to sequencing adapters. Following agarose gel electrophoresis and extraction of cDNA from gels, the cDNA fragments were purified and enriched by PCR to construct the final cDNA library. The cDNA library was sequenced on the Illumina sequencing platform (Illumina HiSeq TM 2500) using the paired-end technology.

Sequence processing and analysis
After removal of adapter and low quality sequences, clean reads were obtained and mapped to the B. napus reference genome (http://www.genoscope.cns.fr/brassicanapus/). They were then assembled using TopHat 2.0.0 and DEGseq. FPKM (fragments per kilobase of exon per million mapped fragments) was used to calculate gene expression levels, and differentially expressed genes (DEGs) were screened with Cuffdiff using the following two criteria: (i) false discovery rate P-value correction of <0.05 and (ii) |log2 (fold change)| > 1 (Trapnell et al., 2010). Based on blastp analysis, homologs of TFs and hormone genes in B. napus were screened from our previous study (Jian et al., 2018), and possible TF and plant hormone genes among these DEGs were also obtained. In addition, all sequencing data were uploaded to NCBI under accession number PRJNA530291.

qRT-PCR analysis
Here, 1 mg of total RNA was used for cDNA synthesis using a Transcriptor First-Strand cDNA Synthesis Kit (Roche, Basel, Switzerland). qRT-PCR was performed on a Bio-Rad CFX96 Real-time System with SYBR Ò Green PCR Supermix (CA, USA). Each reaction of 20 mL contained 10 mL of SYBR Green premix (2X) (Roche, Basel, Switzerland), 0.5 mL of each primer (10 mM), 7.0 mL of H 2 O and 2 mL of cDNA. Three biological replications with three technical replicates were used for each reaction. The following program for qPCR was used: 98 C for 30 s, followed by 40 cycles of 98 C for 10 s and 60 C for 30 s. BnActin7 and Bn26S gene were used as a control. All primers used in this study are shown in Table S1.

Metabolite extraction
The metabolite extraction protocol is based on previous studies with a few modifications (Wang et al., 2018;Ying et al., 2020). The cryopreserved biomaterial samples were removed and vacuum freeze-dried. The dried sample was ground in MM 400 (Retsch, Haan, Germany) at 30 Hz for 1.5 min, 100 mg of the powder was weighed, and 1.0 ml 70% methanol containing 0.1 mg/l lidocaine was used as an internal standard for overnight extraction at 4 C. After extraction, the sample was centrifuged at 10,000×g for 10 min, the supernatant was absorbed, and the sample was filtered with a microporous membrane (0.22-mm pore size), and stored in the injection bottle for LC-MS analysis. The quality control (QC) sample was prepared by mixing the sample extracts and used to analyze the reproducibility of the samples under the same treatment method.
Mass spectrometry information was measured on API 4500 QTRAP LC/MS/MS system both in positive and negative mode, the main parameters of the linear ion trap and triple quadruple rod were as follows: electrospray ion source (electrospray ionization, ESI) temperature 550 C, mass spectrometry voltage 5.5 kV, curtain gas (curtain gas, CUR) 25 psi. The collision-induced ionization (collision-activated dissociation, CD) parameter was set to high. In the qualitative analysis of the metabolite, the isotopic signals, repetitive signals containing K + ions, Na + ions and NH4 + ions were removed. Metabolites were analyzed with reference to internal database, MassBank (http://www.massbank.jp/), KNAPSAcK (http://kanaya.naist.jp/KNApSAcK/), HMDB (http://www.hmdb.ca/) (Wishart et al., 2013), MoTo DB (https://omictools.com/moto-db-tool) and METLIN (http://metlin.scripps.edu/index.php) (Zhu et al., 2013), as well as other available mass spectrometry public databases . In the triple quad (QQQ), each ion pair was scanned and detected according to the optimized declustering voltage (declustering potential, DP) and collision energy (collision energy, CE). The m/z range was set between 50 and 1,000. After mass spectrometry analysis, the original data in wiff format could be opened and browsed using Analyst 1.6.1, which could be used for qualitative and relative quantitative analysis. Cinematographic peak areas were determined and corrected to ensure the accuracy of the metabolite relative quantitative analysis.

Multivariate statistical analysis
Principal component analysis (PCA) carried out by SPSS for Windows (Version 18.0; SPSS Inc., Chicago, IL, USA), partial least squares discriminant analysis (PLS-DA) and orthogonal partial least-squares discriminant analysis (OPLS-DA) were performed by R software package (http://bioconductor.org/packages/release/bioc/html/ropls.html) with data from 12 samples to observe differences in metabolic composition between two B. napus ecotypes after cold stress. The loading plot helped to identify the metabolites that contributed the most to the changes in metabolite patterns between the comparison groups. The variables far from the origin in the transverse coordinate direction contributed more to the distinction between the two groups of samples. In addition, the variables in a certain position of the load graph often provided an important contribution to the samples in the same position on the score graph. We combined the multivariate statistical analysis of the VIP value of OPLS-DA and the univariate statistical analysis of the T test P value to screen significantly different metabolites among different comparison groups. The thresholds used to determine a significant difference were VIP ≥ 1 and T-test P < 0.05.

RESULTS
Transcriptome analysis for winter and spring ecotypes of B. napus in response to cold stress To detect transcriptome differences in response to cold treatment between different oilseed ecotypes, we performed an RNA-Seq analysis of each five winter and spring ecotypes of B. napus using Illumina-based 2 × 150-bp paired-end read sequencing. In total, 598,775,190 clean reads with 587,373,972 high-quality reads were obtained from 12 libraries (Table 1). In our study, all high-quality reads were mapped to the assembled B. napus genome of "Darmor-bzh" (Chalhoub et al., 2014) after mapping to Ribosome. Finally, more than 80.83-82.23% of the clean reads in 12 libraries could be mapped on the reference genome (Table S2).

Identification of DEGs in response to cold stress
In this study, 25,461 DEGs with 5,853 upregulated and 19,607 downregulated genes between the control and cold-treated spring ecotypes were detected. For winter ecotypes, there were 5,912 upregulated and 22,600 downregulated DEGs ( Fig. 2A; Table S4). Interestingly, among the upregulated genes, 3,343 were commonly detected in two ecotypes, and 2,569 and 2,510 genes were specific to winter and spring ecotypes, respectively ( Fig. 2B; Table S4). For downregulated genes, 13,703 genes were common in both winter and spring ecotypes, and 8,897 and 5,904 genes were specific to winter and spring ecotypes ( Fig. 2C; Table S4). In addition, most of the DEGs (67.4% and 71.2% upregulated genes, and 78.2% and 76.5% downregulated genes in winter and spring ecotypes, respectively) had a |fold change| ≤ 8 and <4% DEGs had a |fold change| ≥ 1,024 ( Fig. 2D; Table S4).
As critical regulators, TFs make great contributions in response to various (a) biotic stresses in plants. In this study, 2,346 TFs genes were identified as DEGs in two ecotypes, with 1,855 and 1,585 genes in the winter and spring ecotypes, respectively. Among them, 1,094 genes including 300 upregulated and 794 downregulated genes, were commonly detected in both two ecotypes, of which 761 included 204 upregulated and 557 downregulated genes and 491 included 180 upregulated and 311 downregulated genes in the winter and spring ecotypes, respectively (Fig. 3A). These 2,346 TF genes were unevenly distributed in 53 TF families, of which the top four families, namely bHLH, ERF, MYB and WRKY, made up the majority in all groups. Interestingly, most of the bHLH family members were downregulated in both two ecotypes after cold stress (Table S5). Plant hormones are not only involved in plant development, but they also have crucial functions in abiotic stresses. In total, 434 genes, including 159 upregulated and 275 downregulated, involved in hormone metabolism or signal transduction process were identified. In addition, 409 genes, including 127 upregulated and 282 downregulated, and 393 genes, including 110 upregulated and 283 downregulated, were screened in winter and spring ecotypes after cold stress, respectively (Fig. 3B). Of them, 268 genes, including 78 up-and 190 downregulated, were commonly detected in both ecotypes after cold stress (Fig. 3B). These hormone genes, ABA, auxin and ET genes made up the majority (Table S6).
In Arabidopsis, 320 genes were identified as key factors involved in low or high temperature stresses. In this study, 568 (234 upregulated and 334 downregulated) genes in total and 284 (138 upregulated and 146 downregulated) genes in common were detected in two ecotypes of B. napus after cold stress. In the winter ecotype, 436 genes, including 190 upregulated and 246 downregulated, and 416 genes, including 182 upregulated and 234 downregulated, were screened in winter and spring ecotypes, respectively (Fig. 3C). Among them, key genes encoded calcineurin B-like protein 1 (CBL1), cold regulated 15b   (Table S7). When plants perceive a low external temperature, they transmit the signal from the cell membrane to the cell, which causes a change in the calcium ion concentration on the cell membrane and triggers a phosphorylation cascade of kinases. The calcium concentration increased to promote the upregulation of ICE1 or ICE1-like, CAMTA1-3 and ZAT12, which directly or indirectly promoted the upregulation of CBFs, followed by adaption to low temperature stress through RAP2.1, 2.6 and ZAT10 (Fig. 4A). In this process, the two ecotypes of rapeseed displayed similar expression patterns, but the degree of change was slightly different (Fig. 4B). Several of them were confirmed by qRT-PCR analysis (Fig. 5). High correlation coefficients were detected between RNA-Seq and qRT-PCR technologies, suggesting that the RNA-Seq data were reliable.

Functional analysis of DEGs involved in the response to cold stress
To classify the functions of the DEGs in the two ecotypes after cold stress, GO was used according to biological process, cellular component and molecular function (Fig. S2). Among all DEG comparisons, "metabolic process" and "cell process", "cell part" and "cell" and "catalytic activity" and "binding" were dominant in "biological process", "cellular component" and "molecular function", respectively. In addition, DEGs identified in the two ecotypes after cold stress were used for KEGG pathway analysis to explore their potential functions. In total, 27 and 44 significantly enriched pathways were detected among the up-and downregulated genes, respectively ( Fig. 6; Table S8). The upregulated genes that were commonly detected in the two ecotypes were mainly involved in primary metabolism, including amino acid metabolism (arginine, proline, phenylalanine, cysteine and methionine metabolism), carbohydrate metabolism (amino sugar and nucleotide sugar, inositol phosphate, starch and sucrose, glycolysis/ gluconeogenesis, ascorbate and aldarate and butanoate metabolism) lipid metabolism (glycerophospholipid and glycerolipid metabolism) and secondary metabolism (flavonoid biosynthesis, phenylpropanoid biosynthesis, flavone and flavonol biosynthesis and carotenoid biosynthesis), suggesting that this reprograming of primary and secondary metabolism might represent a response to cold stress ( Fig. 6; Table S8). In addition, we observed signal transduction (phosphatidylinositol signaling system), transport and catabolism (endocytosis) and environmental adaptation (circadian rhythm-plant and plant-pathogen interaction). Moreover, plant hormone signal transduction was only detected in the winter B. napus ecotype and glutathione metabolism and stilbenoid, diarylheptanoid and gingerol biosynthesis were only detected in the spring ecotype ( Fig. 6; Table S8), suggesting that hormone signaling and secondary metabolism were involved in the cold stress response in different B. napus ecotypes. Conversely, among the downregulated genes commonly detected in the two ecotypes were pathways involved in amino acid metabolism (arginine biosynthesis, valine, leucine and isoleucine biosynthesis, lysine biosynthesis and valine, leucine and isoleucine degradation, phenylalanine, tyrosine and tryptophan biosynthesis and alanine, aspartate and glutamate metabolism), carbohydrate metabolism (propanoate metabolism), lipid metabolism (steroid biosynthesis, fatty acid elongation, fatty acid biosynthesis and alpha-linolenic acid metabolism), nucleotide metabolism (purine metabolism) and glycan biosynthesis and metabolism ( Fig. 6; Table S8). In addition, energy metabolism genes (photosynthesis, sulfur metabolism and photosynthesis-antenna proteins) also seemed to be downregulated. Plant-pathogen interaction pathways were only detected in the spring B. napus ecotype after cold stress ( Fig. 6; Table S8). In the winter B. napus ecotype, eight pathways were also detected: glucosinolate biosynthesis, RNA degradation, microbial metabolism in diverse environments, glycerophospholipid metabolism, folate biosynthesis, riboflavin metabolism, brassinosteroid biosynthesis and regulation of autophagy ( Fig. 6; Table S8), suggesting that different pathways were involved in different ecotypes after cold stress.
To further explore the molecular mechanisms underlying the response to cold stress in the two B. napus ecotypes, the software MapMan was used to classify the DEGs involved in the metabolic pathways (Fig. 7). As shown in Fig. 7, genes involved in photosynthesis and circadian rhythm were mostly suppressed, while cell wall, lipids, and biosynthesis of secondary metabolites were mostly activated. All DEGs involved in different metabolic pathways are summarized in Tables S9 and S10.

Combined analysis of gene expression and metabolite data
To screen key genes and metabolites involved in the response to cold stress, joint data from gene expression and metabolism using pairwise correlations and O2PLS were assessed. A pairwise correlation analysis of all DEGs and DAMs was performed, and the top 50 DEGs and top 50 DAMs are displayed in Fig. 9; Table S13. Two latent variables were identified using O2PLS analysis, with 98.0-98.3% of the total variation in the transcript dataset and 97.4-98.3% of the total variation in the metabolite dataset (Table S14).
To visualize the importance of the elements in the two omics, a loading plot of different omics was drawn for all transcriptional and metabolic data (Tables S15 and S16). As shown in Figs. 9A-9C and 9B-9D, the loading coefficient thresholds were higher for the metabolomics dataset than for the transcriptomic dataset, and the top 10 genes and metabolites are represented as a red circular plot (Figs. 9A-9C and 9B-9D). In the winter B. napus ecotype, four metabolites (D-(+)-sucrose, nicotinic acid adenine dinucleotide, syringic acid O-feruloyl-O-hexoside and dihydromyricetin) were increased, and the rest of them (MAG (18:3) isomer2, methylmalonic acid, adenine, L-threonine, N-acetyl-L-tyrosine and kynurenic acid) were decreased (Table 2). Interestingly, all ten genes were downregulated, two of which were involved in protein posttranslational modifications, three in the cell, one in amino acid metabolism, one in RNA regulation, one in glycolysis and two encoded unknown proteins (Table 3). In the spring B. napus ecotype, four metabolites (Coniferin, 3-O-p-coumaroyl quinic acid O-hexoside, LysoPC 18:1 (2n isomer) and L-phenylalanine) were decreased and six metabolites (MAG (18:3) isomer1,  (Table 2). Among the top 10 genes, nine were upregulated and one was downregulated with unknown functions. Nine upregulated genes were found to encode MYB111, RAV1, CARBONIC ANHYDRASE 1, FASCICLIN-LIKE ARABINOGALACTAN PROTEIN 8, NICOTIANAMINE SYNTHASE 3, UDP glucosyl and glucoronyl transferases, anthocyanins, anthocyanin 5-aromatic acyltransferase, eukaryotic aspartyl protease family protein and one unknown function protein (Table 3).

DISCUSSION
As one major adverse environmental factor, low temperature decreases crop production worldwide and has become a great concern for crop breeders (Khodakovskaya et al., 2006;Chinnusamy, Zhu & Zhu, 2007). Plants have evolved fine regulatory mechanisms in response to low temperature stress, including changes in gene expression and metabolites (Maruyama et al., 2014). In this study, alterations of genes and metabolites after cold stress in both winter and spring B. napus ecotypes were detected using RNA-Seq and LC-MS analysis, respectively.

ABA signal transduction components have critical roles in response to cold stress
Increasing evidence has revealed the types of mechanisms involved in cold stress response.
As an important stress phytohormone, ABA participates in regulating diverse abiotic stresses. COR genes are precisely regulated by both ABA-dependent and -independent pathways. Hormone homeostasis may be altered under cold stress compared with normal conditions (Lee, Henderson & Zhu, 2005). ABA levels are increased after cold stress in potato (Chen, Li & Brenner, 1983). Plant freezing tolerance can be enhanced by exogenous application of ABA at a normal temperature, while the functions of ABA in the cold regulation process remain elusive.

Signal transduction in response to cold stress
The signal transduction pathway makes great contributions in response to various environmental stresses, including cold stress (Janska et al., 2010;D'Angeli & Altamura, 2016;Tanou et al., 2017). As an important messenger, Ca 2+ plays a crucial role in stress signaling (Ye et al., 2013). The Ca 2+ level increases rapidly under low temperature stress and activates various signaling pathways in response to cold stress and cold-induced gene regulatory processes (Monroy & Dhindsa, 1995;Dodd, Kudla & Sanders, 2010;Reddy et al., 2011;Shi, Ding & Yang, 2015). Calcium signals are sensed by three major family proteins in higher plants: calmodulins (CaMs), calcium-dependent protein kinases (CDPK) and calcineurin B-like (CBL) proteins (Sheen, 1996;Luan, 2009;Boudsocq & Sheen, 2013;Miura & Furumoto, 2013). The MAPKs cascade and receptor-like protein kinases (RLKs) also contribute in the perception of signals in the external environment and have important roles in environmental stress responses (Morris & Walker, 2003;Lehti-Shiu et al., 2009). The calmodulin-like OsMSR2 gene is significantly induced by low temperature, and its over-expression in Arabidopsis enhances the salt and drought tolerance ability of the plants (Xu et al., 2011). Five calmodulin genes, two CDPK genes and one CBL gene in Camellia sinensis were identified to take part in signal transduction after cold stress using RNA-Seq analysis (Wang et al., 2013b). To determine whether these gene families were involved in cold response in B. napus, we explored these genes in the list of DEGs. The results showed that 135 genes, including 44 MAPK cascade genes, 35 CDPKs, 38 CaMs, 16 CIPKs and 2 Ca 2+ -ATPase genes, were differentially expressed, with 68 upregulated and 67 downregulated in winter B. napus ecotype after cold stress. In the spring B. napus ecotype, 117 genes, including 38 MAPK cascade genes, 44 CDPKs, 19 CaMs, 15 CIPKs and 1 Ca 2+ -ATPase genes, were differentially expressed, with 34 upregulated and 83 downregulated after cold stress (Table S17). Among them, 54 genes (including 25 MAPK cascade genes, 26 CDPKs, 14 CaMs, 13 CIPKs and 1 Ca 2+ -ATPase gene) were commonly detected, with 27 upregulated and 52 downregulated in both winter and spring B. napus ecotypes. Interestingly, two genes, encoding CBL1 (BnaA01g08510D) and ACA2 (BnaA08g15730D), showed opposite expression patterns between the two B. napus ecotypes after cold stress. These two genes were downregulated in the winter B. napus ecotype but upregulated in the spring B. napus ecotype after cold stress (Table S17). Our results provide further evidence supporting the key functions of Ca 2+ binding proteins and the MAPK pathway in response to cold stress in B. napus.
Several COR genes are not directly regulated by the CBF/DREB family (Fowler & Thomashow, 2002;Kreps et al., 2002). WRKY and NAC TFs also play key roles in response to cold stress (Hu et al., 2008;Yokotani et al., 2013;Yang et al., 2015;Zhuo et al., 2018). An involvement of the WRKY gene family in the cold stress response has been suggested in Arabidopsis (Lee, Henderson & Zhu, 2005), rice (Ramamoorthy et al., 2008), soybean (Zhou et al., 2008), Vitis vinifera , cucumber (Zhang et al., 2016) and cassava (Wei et al., 2016). In our study, 150 (128 in the winter and 115 in the spring ecotype) genes were differentially regulated after cold stress, with 47 up-and 103 downregulated. Among them, 26 up-and 46 downregulated were detected in both winter and spring B. napus ecotypes after cold stress. Key members, such as BnWRKY8s, 23s, 28s, 31s, 42s, 48s and 123s, were upregulated, while BnWRKY7s, 51s, 70s and 74s were downregulated after cold stress in winter and/or spring B. napus ecotypes. Key members of NAC genes involved in plant cold responsiveness have been reported, such as TaNAC2, OsNAC5 and SsNAC23 (Hu et al., 2008;Song et al., 2011). Thirty-two and twenty-eight NAC genes were upregulated after cold stress in the winter and spring B. napus ecotypes, respectively, and nine of them were common to the two ecotypes. Many other kinds of TFs, such as bHLH, C2H2, MYB, MYC, MADS and b-ZIP families, were also detected as DEGs after cold stress (Table S5), suggesting that these TFs are also involved in the low temperature regulatory process (Chinnusamy et al., 2003;Yang, Dai & Zhang, 2012;Peng et al., 2013). These up-and downregulated TFs suggest very complex mechanisms of low temperature regulation.

CONCLUSIONS
In this study, a total of 25,460 DEGs and 41 differentially expressed metabolites (DEMs) in spring oilseed ecotype and 28,512 DEGs and 47 DEMs in winter oilseed ecotype were identified by transcriptome and metabolome analyses after cold stress. In addition, more specific genes and metabolites were detected in winter ecotypes. The combined data indicate that lipid, ABA, secondary metabolism, signal transduction and TFs could be involved in the complex regulation in the spring and winter ecotypes in response to cold stress. This study provides new evidence and insights into how different ecotypes of B. napus respond to cold stress from metabolomics and transcriptomics levels. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Grant Disclosures
The following grant information was disclosed by the authors: Talent

Competing Interests
Kun Lu is an Academic Editor for PeerJ.

Author Contributions
Hongju Jian conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Ling Xie performed the experiments, analyzed the data, prepared figures and/or tables, and approved the final draft. Yanhua Wang performed the experiments, prepared figures and/or tables, and approved the final draft. Yanru Cao performed the experiments, prepared figures and/or tables, and approved the final draft. Mengyuan Wan performed the experiments, analyzed the data, prepared figures and/or tables, and approved the final draft. Dianqiu Lv conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft. Jiana Li conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft. Kun Lu conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft. Xinfu Xu conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft. Liezhao Liu conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft.

Data Availability
The following information was supplied regarding data availability: Data is available at the NCBI: PRJNA530291.

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.8704#supplemental-information.