Comparative Transcriptome Analyses Revealed Conserved and Novel Responses to Cold and Freezing Stress in Brassica napus L

Oil rapeseed (Brassica napus L.) is a typical winter biennial plant, with high cold tolerance during vegetative stage. In recent years, more and more early-maturing rapeseed varieties were planted across China. Unfortunately, the early-maturing rapeseed varieties with low cold tolerance have higher risk of freeze injury in cold winter and spring. Little is known about the molecular mechanisms for coping with different low-temperature stress conditions in rapeseed. In this study, we investigated 47,328 differentially expressed genes (DEGs) of two early-maturing rapeseed varieties with different cold tolerance treated with cold shock at chilling (4°) and freezing (−4°) temperatures, as well as chilling and freezing stress following cold acclimation or control conditions. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis indicated that two conserved (the primary metabolism and plant hormone signal transduction) and two novel (plant-pathogen interaction pathway and circadian rhythms pathway) signaling pathways were significantly enriched with differentially-expressed transcripts. Our results provided a foundation for understanding the low-temperature stress response mechanisms of rapeseed. We also propose new ideas and candidate genes for genetic improvement of rapeseed tolerance to cold stresses.

Abscisic acid (ABA) is a vital plant hormone that plays a key role in stress resistance during plant growth and development (Vishwakarma et al. 2017). It was reported that OST1/SnRK2E, a serine-threonine protein kinase in ABA core signaling pathway, acted upstream of CBFs to positively regulate freezing tolerance via phosphorylating ICE1 to prevent its 26S proteasome-mediated degradation by HOS1 (Ding et al. 2015). SnRK2E phosphorylated basic transcription factors 3 (BTF3) and BTF3-like factors, and facilitated their interactions with CBFs to promote CBF stability under cold stress .
Jasmonic acid (JA) is a lipid-derived plant hormone that plays an important role in the plant abiotic and biotic stresses (Wasternack and Song 2017;Wasternack and Strnad 2018). It was reported that JA signaling positively regulated the plant responses to freezing stress via the interaction between the repressors of jasmonate signaling proteins (TIFY/JAZ) and ICE1/2, thus suppressing the transcriptional activity of ICE proteins, thereby attenuating CBF genes expression (Hu et al. 2013).
Oil rapeseed is a winter biennial oil-seed crop, responds to vernalization and shows an excellent tolerance to cold stresses during vegetative stage. In recent years, more and more early-maturing rapeseed varieties were planted across China. Unfortunately, the earlymaturing rapeseed varieties with low cold tolerance have higher risk of freeze injury in cold winter and spring. Hence, it is vital to compare early-maturing rapeseed varieties tolerant to cold and evaluate molecular mechanisms that adapt to different low-temperature stress conditions.
In this study, RNA-seq technology was used to compare differentially expressed genes (DEGs) of two early-maturing rapeseed varieties with different cold tolerance, that were treated with cold shock (4°and 24°) temperatures, as well as chilling (4°) and freezing (24°) temperatures following cold acclimation and control conditions. Analyses identified conserved and novel signaling pathways and genes. Our results provided a foundation for understanding the low-temperature stress response mechanisms of rapeseed. We also propose new ideas and candidate genes for genetic improvement of rapeseed tolerance to cold stresses.
Total RNA extraction and sequencing The samples were sent to the sequencing cooperation of Novogene Co. Ltd (Bejing, China) for RNA isolation, examination and sequencing. HiSeq PE150 of Illumina was selected for the sample transcriptome sequencing.

Transcriptome assembly and annotation
Raw data (raw reads) of fastq format were first processed through in-house perl scripts. In this step, clean data/reads were obtained by removing reads containing adapter or, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30 and GC content for the clean data were calculated. All the downstream analyses were based on the clean data with high quality.
Reference genome and gene model annotation files were downloaded from genome website directly. Index of the reference genome was built using Hisat2 v2.0.4 and paired end reads were aligned to the reference genome using Hisat2 v2.0.4. Hisat2 was selected as the mapping tool because it can generate a database of splice junctions based on the gene model annotation file, thus providing better mapping result than other non-splice mapping tools. HTSeq v0.9.1 was used to count the reads mapped to each gene. FPKM of each gene was calculated based on the length of the gene and read counts mapped to the gene. FPKM considers the combined effects of sequencing depth and gene length for the read counts, and is currently the most commonly used method for estimating gene expression levels ).

Differential expression analyses of transcripts
Differential expression analysis of two conditions/groups (three biological replicates per condition) was performed using the DESeq R package (1.18.0). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution (K ij NB(m ij ,s ij 2 )). The resulting P-values were adjusted using the Multiple Hypothesis Testing for controlling false discovery rate. Genes with an adjusted P-value , 0.05 calculated by DESeq were assigned as differentially expressed.
GO and KEGG enrichment analyses of differentially expressed genes (DEGs) Gene ontology (GO) enrichment analysis of DEGs was implemented by the GOseq R package, in which gene length bias was corrected. GO terms with corrected P-value , 0.05 were considered significantly enriched in DEGs.
KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, at molecular level, especially large-scale datasets generated by Figure 1 The processes of different low-temperature treatments on two rapeseed (Brassica. napus L. cv. HX17 and HX58) seedlings.
genome sequencing and other high throughput technologies (http:// www.genome.jp/kegg/). KOBAS software was used to test the statistical significance of enrichment of DEGs in KEGG pathways.

qRT-PCR
To confirm the validity of the RNA-Sequencing data, we randomly selected nine DEGs for qRT-PCR analysis. 3 mg total RNAs/sample was used for cDNA biosynthesis with Superscript III reverse transcriptase (Invitrogen, San Diego, USA). The qRT-PCR was performed using the CFX96 Real-Time System (Bio-Rad, USA) with SYBR green (Bio-Rad, USA). The relative changes were calculated with 2-ΔCt and the BnaActin (BnaC05g34300D) from B. napus was used as the endogenous reference gene. The relative transcript level was determined and normalized using the reference level and averaged over the three technical replicates. Primers for the qRT-PCR (Table S10) were designed according to the cDNAs sequences using Primer Premier 5 (http://www.premierbiosoft.com/crm/jsp/com/pbi/crm/clientside/ ProductList.jsp).

Transcriptomic analyses
Two early-maturing semi-winter rapeseed varieties (HX17 and HX58) with different cold tolerance ( Figure S1) were selected for transcriptomic investigation by RNA-Seq following different low-temperature treatments.
Total RNA was extracted from the seedlings of HX17 and HX58 treated with different low-temperature conditions and sequenced using the Illumina HiSeq 2500 platform. A total of 296.67 Gb clean data (1,977,933,478 clean reads) were acquired after removal of adaptor sequences, ambiguous nucleotides and low-quality sequences (Table 1). Approximately 43.1-78.1 million clean reads were obtained from the libraries. The Q-score for 88.17-93.99% of reads was Q30, and 82.57-86.51% of the total reads were uniquely mapped to the reference genome (Table S1). The mapped sequences were assembled with Cufflinks, referenced against the annotated genome sequence of the tetraploid species B. napus provided by Index of /brassicanapus/data (http://www.genoscope.cns.fr/brassicanapus/ data/). RNA-Seq revealed that there was a total of 172,412 transcripts and 110,367 unigenes, and 10,578 genes on average for each sample, including 78,945 known genes and 26,233 novel genes on average (Table S2). Gene expression levels were estimated by fragments per kilo base of transcript per million fragments mapped (FPKM) (Langmead et al. 2009). Approximately 8.51-11.97% of total genes' FPKM was more than 15.0, and 59.26-65.50% of total genes' FPKM was lower than 1.0 (Table S3).
For the evaluation of DEGs' reliability and the filter of abnormal samples, Pearson correlation coefficient (PCC) analysis was conducted. The correlation analysis indicated that all three technical replicates replication samples showed more than 0.929 (R 2 ) of similarities (Table S4). Differentially expressed genes (DEGs) in rapeseed under different low-temperature treatments Differential expression analysis of treatments and control group was performed using the DESeq (Anders and Huber 2010). A threshold of fold change of $2 and p-adjusted q-value of , 0.005 was used for identifying DEGs.
A  (Figure 2 and Table S5). Statistically, DEGs in both HX17 and HX58 showed that FB had the most DEGs and CA had the least, in the order of FB . CB . FA. CA. HX58 presented with more DEGs than HX17 (Figure 2), consistent with the cluster analysis of DEGs (Figure 3 and Figure S2). In addition, the lower temperature correlated with more injuries, and HX58 was more susceptible to low-temperature stress ( Figure S1). DEGs under different stresses in both HX17 and HX58 were presented as CA vs. FA . CB vs. FB (Figure 2), confirming that the cold acclimation altered the responses to freezing stress in plant.
Pathway analyses of responses to different lowtemperature stresses in rapeseed Pathway analyses of DEGs against the Kyoto Encyclopedia of Genes and Genomes (KEGG) contributed to the understanding of gene functions (Kanehisa et al. 2008). As shown in Figure S3 and Table S6, there was no enriched pathway in both HX17 and HX58 treated with cold shock (CB and FB) compared to control condition (MB), indicating that cold shock stress caused serious injuries to the plant. Interestingly, there was no enriched pathway in HX58 acclimatized to cold (CA and FA), similar to cold shock stress. However, the plant hormone signal transduction pathway was enriched in HX17 acclimatized to cold (CA and FA), indicating that the plant hormone signal transduction was activated by cold acclimation in cold-tolerant rapeseed. Further analyses of the upregulated and downregulated DEGs enriched pathways are reported in Table S7. Pathways involved in primary metabolism and plant hormone signal transduction were enriched under different low-temperature stresses in both HX17 and HX58 (Fig. S3, Table  S6). Previous studies also reported changes in these two pathways in various crops and plants exposed to low-temperature stresses (Nievola et al. 2017;Qu et al. 2015;Zhang et al. 2017;Li et al. 2017;Bai et al. 2015;Zhang et al. 2014), suggested that low-temperature stresses were conservative to affect the normal metabolism of plants and the plant hormone signal transduction.
There were many reports that the plant hormone signal pathways played important roles in the plant responses to abiotic stresses (Colebrook et al. 2014;Kazan 2015;Gururani et al. 2015;Bielach et al. 2017;Verma et al. 2016;de Zelicourt et al. 2016;Khan et al. 2015;Sharma et al. 2017). The ABA and JA biosynthesis and signaling were induced by various abiotic stresses, and positively regulated the low-temperature tolerance in plants (Yoshida et al. 2014;Huang et al. 2017;Shinkawa et al. 2013;Hu et al. 2013;Du et al. 2013;Wang et al. 2016a). Hence, we analyzed the expression patterns of ABA and JA signaling genes. The results showed that among the 503 DEGs of the plant hormone signaling transduction (total of 1148 genes), 94 of them were the ABA signaling genes (total of 139 genes) and 43 of them were JA signaling genes (total of 61 genes)( Table S7). Most of the ABA signaling genes were changed under cold acclimation and cold shock treatments in both HX17 and HX58 but they presented different expression profiles (Figure 4). The ABA receptors pyrabactin resistance proteins/PYR-like proteins/regulatory components of ABA receptors (PYR/PYL/RCAR) PYL5 and PYL7 were both induced by cold and freezing treatments in HX17 and HX58, while the PYL1 and PYL9 were suppressed in all treatments. The PYR1 and PYL4 were only induced by cold treatment but not by freezing treatment. The PYL6 were induced by cold treatment and freezing treatment following cold acclimation. ABA co-receptor phosphatase 2C (PP2C) ABI1 and HAB1 were suppressed by all treatments, while HAB2 was upregulated. The Snf1related protein kinase 2 (SnRK2) SnRK2B and SnRK2D were induced by all low-temperature treatments, while SnRK2C was suppressed. Interestingly, SnRK2F was induced only by freezing treatment, similar to ABI5 expression, indicating that they played important roles in response to freezing in rapeseed. The Arabidopsis ABI5 was a bZIP transcription factor phosphorylated by SnRK2D and SnRK2I, and participated in ABA-regulated gene expression during seed development and subsequent vegetative stage by acting as the major mediator of ABA repression of growth (Skubacz et al. 2016;Nakashima et al. 2009;Brocard et al. 2002). It is possible that SnRK2F directly phosphorylated ABI5 in rapeseed. The SnRK2E/OST1 was induced only in HX17 but not in HX58, which may be a part reason for the better coldtolerance of HX17 than HX58.
It was reported that exogenous application of JA enhances plant freezing tolerance with or without cold acclimation, and overexpression of JAZ1/TIFY10A and JAZ4/TIFY6 suppressed the transcriptional function of ICE1, thereby attenuating the ICE-CBF-COR signaling pathway and freezing responses in Arabidopsis (Hu et al. 2013). As shown in Figure 3 Hierarchical cluster analysis of differentially expressed genes (DEGs) among all the 12 samples. Heatmap of DEGs among the 12 samples. Red indicates high expression, and blue indicates low expression. Color from red to blue represents descending log10 (FPKM + 1). Figure 5, the expressions of many JA signaling genes were altered in both HX17 and HX58 under low-temperature conditions. TIFY3 and TIFY6 were induced in both HX17 and HX58 under all lowtemperature conditions. Interestingly, TIFY7 was suppressed only under cold shock treatments, while the TIFY10A and TIFY11B were upregulated only by freezing treatments, similar to SnRK2F and ABI5 (Figure 4 and Figure 5). These changes suggest that they may be important regulators in the freezing response in rapeseed. Additionally, MYC2 was only induced by cold acclimation, but not by cold shock.
To test whether exogenous ABA and methyl-jasmonic acid (MeJA) could improve the rapeseed cold tolerance, 4-week-old rapeseed seedlings were sprayed with ABA (100 mM), MeJA (100 mM) and water before cold shock at -6°. The survival rate results demonstrated that exogenous ABA significantly improved the rapeseed freezing tolerance. On the contrary, exogenous MeJA reduced the rapeseed freezing tolerance ( Figure 6). The concentration of exogenous MeJA might be too high, affecting the rapeseed freezing tolerance. Other studies reported that treatment with exogenous JA induced leaf senescence and expression of senescence associated genes, while exogenous high concentration MeJA caused plant injures (Zhu et al. 2015). All the results indicated that a complex gene regulatory network composed of plant hormones signaling affected the low-temperature stresses in rapeseed, and appropriate concentrations of exogenous plant hormones could increase rapeseed's cold tolerance.
Identification of the transcriptional regulatory network responding to low-temperature stress in rapeseed Transcription factors (TF) play a central role in the gene regulatory networks that mediate various aspects of plant developmental processes and responses to environmental changes (Hoang et al. 2017;Riechmann and Ratcliffe 2000). The program iTAK was used to analyze DEGs and identified 7633 TFs. Of these, 3670 (48.1%) belong to 76 transcription factors and transcriptional regulator families. The 10 top families were MYB, AP2-EREBP, bHLH, Orphans, HB, bZIP, NAC, C2H2, WRKY and C3H (Table S8).
It was reported that CBFs were the core regulator in the coldresponse signaling in plants Guo et al. 2018;Ruelland et al. 2009;Wang et al. 2017). CBFs belong to the DREB1 proteins, including DRE1A/CBF3, DRE1B/CBF1, DRE1C/CBF2 and DRE1D/ CBF4 (Liu et al. 1998;Jaglo-Ottosen et al. 1998;Sakuma et al. 2002). We analyzed the ICE1, CBFs and some COR genes, and observed that COR15A and COR15B were significantly induced by all low-temperature treatments ( Figure 7A). Apart from BnaC03g71900D, all the CBF genes were induced by freezing stress more than chilling stress, regardless of cold acclimation ( Figure 7A), highlighting the importance of CBFs in rapeseed freezing tolerance. However, ICE1 was not induced by lowtemperature treatments, as opposed to AtICE1 (Chinnusamy et al. 2003).
The cold shock domain-containing proteins (CSP) act as RNA chaperones that destabilize mRNA secondary structures at low temperatures, and they are regulated by low-temperature stress (Sasaki and Imai 2011;Choi et al. 2015;Radkova et al. 2014;Chaikam and Karlson 2008). Overexpression of AtCSP2 resulted in decreased freezing tolerance despite cold acclimation (Sasaki et al. 2013). CSP functions are associated mainly with cold adaptation, but they are also involved in other biological processes under normal growth conditions (Sasaki and Imai 2011;Nakaminami et al. 2009;Park et al. 2009). All the 13 cold CSPs were regulated by the low-temperature treatments ( Figure 7B). As shown in Figure 5 Heatmap analysis of the genes in JA signal pathway. JAR1 (JASMONATE RESISTANT 1), Jasmonic acid-amido synthetase; COI1 (Coronatine-insensitive protein 1), Jasmonic acid receptor; TIFY, Jasmonate ZIM domain-containing protein; MYC2, bHLH transcription factor MYC2. The expression levels of each gene (FPKM values) in each sample, as indicated by different color rectangles. Red means high expression; Blue means low expression. Figure 6 Effect of ABA and MeJA on the freezing tolerance of ZS11 (Brassica. napus L. cv. Zhongshuang 11) seedlings. ZS11 seedlings were grown for 4-week under normal conditions at 20°C, and then sprayed with ABA (100 mM), MeJA (100 mM) or water (Mock) for 1 day. Then the seedlings were treated with -6°C for 2 days and recovered under normal conditions at 20°C for 7 days. Survival rates of seedlings were shown in brackets. Figure 7B, three CSP genes (BnaA07g02430D, BnaCnng42010D, and BnaAnng08360D) were suppressed by low-temperature treatments, and the others were induced by at least one treatment. BnaA01g01560D was induced by all low-temperature treatments. Four CSP (Novel01004, BnaC03g57650D, BnaA08g19010D, and BnaC07g11440D) were induced by cold shock (CB and FB) and chilling after cold acclimation (CA), but not freezing after acclimation (FA). BnaA09g43450D and Novel04090 were induced by cold shock (CB and FB) and chilling after cold acclimation (FA), but not freezing after acclimation (CA). All the results indicated that the CSP genes have a complex expression pattern in response to different low-temperature conditions.
Heat stress transcription factors (HSFs) play a crucial role in plant responses to high-temperature by regulating the expression of stressresponsive genes, such as heat shock proteins (Guo et al. 2016). Approximately half of HSFs were affected under low-temperature condition ( Table S8), suggesting that HSFs responded to high-temperature and low-temperature changes.
Plant-pathogen interaction pathways play important roles in low-temperature stresses in rapeseed As shown in Table S6-S7, the plant-pathogen interaction pathway was enriched in both HX17 and HX58 under low-temperature conditions, especially under the freezing stresses, consistent with expression of the plant-pathogen interaction pathway genes (Figure 8). Even though the plant-pathogen interaction pathway was consistently reported to be enriched under low-temperature treatments in various plants Qu et al. 2015;Yang and Huang 2018;Xu et al. 2016;Wang et al. 2016b;Tian et al. 2013;Du et al. 2015;Wang et al. 2015), the functions of this pathway in the plant cold responsiveness were neglected by researchers. The pathogenesis-related proteins (PR) were induced by cold acclimation (Kuwabara and Imai 2009;Hincha et al. 1997). Pretreatment of tomato fruit with MeJA or MeSA induced the synthesis of PR proteins, which led to increased chilling tolerance and resistance to pathogens (Ding et al. 2002). Overexpression of 3 PR genes (PR2, PR4, Glu) in Arabidopsis enabled the plant to tolerate freezing temperatures (Cabello et al. 2012). As shown in Figure 8b, PR1 and PR5 were induced by cold acclimation and cold shock. Overexpression of a rice PR5 gene in rapeseed enhanced the resistance to Sclerotinia sclerotiorum, a fungal pathogen that caused Sclerotinia stem rot in rapeseed (Aghazadeh et al. 2017). We propose that overexpression of PR genes in rapeseed might produce rapeseed germplasm with tolerance to cold and S. sclerotiorum.
Circadian rhythms play an important role in lowtemperature stresses in rapeseed As shown in Table S6-S7, the circadian rhythms pathway was enriched in the DEGs for both HX17 and HX58 under low-temperature conditions. It was also shown to respond to low-temperature conditions in several plants (Li et al. 2017;Sobkowiak et al. 2014;Liu et al. 2018;Abeynayake et al. 2015;Barah et al. 2013), and regulated CBFs and CORs gene expression (Edwards et al. 2006;Bieniawska et al. 2008). Almost all the genes associated with circadian rhythms showed differential expression (Figure 9). These included induction of the pseudoresponse regulators (APRR) APRR1/TOC1, APRR3 and APRR5 by chilling and freezing stresses in both HX17 and HX58, and suppression  of APRR7 and APRR9, consistent with the reports on Arabidopsis (Lee et al. 2005). A prr7-3 prr9-1 double mutant failed to entrain to temperature cycles (Salome and McClung 2005), while prr5/7/9 triple mutant plant was more tolerant to cold stresses through upregulation of CBF expression (Nakamichi et al. 2009). The Cholcone synthase (CHS), twin sister of FT (TSF), phytochrome A (PHYA) and transcription factor PIF3 were induced by chilling and freezing stresses. The circadian clock associated 1 (CCA1), early flowering 3 (ELF3), late elongated hypocotyl (LHY), phytochrome B (PHYB), suppressor of PHYA-105 1 (SPA1), repressor of UV-B photomorphogenesis 1 (RUP) and transcription factors HY5 and MYB75 were suppressed by chilling and freezing stresses. Cyclic dof factor 1 (CDF1) and cytochrome 1/2 (CRY1/2) were induced by cold acclimation, but were suppressed by cold shock. Gigantea (GI) and CHE/TCP21 were induced only by freezing stress.
Two phytochrome mutants phyB and phyD showed increased cold tolerance via the upregulation of CBF regulon (Franklin and Whitelam 2007). During the warm long-day condition growing season, the CBF pathway is actively suppressed by PHYB, PIF4, and PIF7 (Lee and Thomashow 2012). PHYA and PHYB function antagonistically to regulate cold tolerance via abscisic acid-dependent jasmonate signaling in tomato (Wang et al. 2016a). PIF3 acts as a negative regulator of plant cold acclimation by direct suppression of CBF expression .
HY5 transcription was regulated by low temperature through a CBFand ABA-independent pathway, and mutant hy5-1 provoked a significant reduction in the ability of Arabidopsis to cold acclimate (Catala et al. 2011). GI increased freezing tolerance via a CBF-independent pathway (Cao et al. 2005). CCA1/LHY bound to promoter of CBFs and promoted cold acclimation (Dong et al. 2011). COR27 and COR28 were direct targets of CCA1, meanwhile they bind to the chromatin of APRR1 and APRR5 and repress their expression (Li et al. 2016). Therefore, we could infer that the circadian rhythm pathways play a role in  Figure 11 qRT-PCR verification of 6 selected DEGs. Comparison of qRT-PCR (blue bar) with FPKM data (red line). The relative changes were calculated with 2-ΔCt and the rapeseed BnaActin gene was amplified as the reference gene. The relative qRT-PCR expression level (selected DEG/ BnaActin) is shown on the left y-axis. The relative transcript level was determined and normalized using the reference level and averaged over the three technical replicates. The FPKM from the RNA-Seq data are indicated on the right y-axis. freezing tolerance, through regulation of the circadian rhythm genes in rapeseed.

Vernalization in rapeseed
Although both cold acclimation and vernalization are responses to low temperature, the durations of cold exposure required to initiate these responses are quite distinct. The most upstream event of the vernalization pathway identified is the induction of vernalization insensitive 3 (VIN3) in response to low temperatures. Together with VIL1, VIN3 was required during vernalization for the modifications of FLC and FLM chromatin, it led to an epigenetically silenced state and acquisition of competence to flower (Sung and Amasino 2004).
We analyzed the expression pattern of all rapeseed VIN and FLC genes under different low-temperature conditions. As shown in Figure  10, only two VIN genes (BnaC07g42420D and BnaAnng10940D) were induced by cold acclimation and cold shock. BnaC01g08200D and BnaCnng18840D were suppressed by both cold acclimation and cold shock. BnaC03g68030D and BnaA08g13150D were induced only by cold acclimation but not by cold shock. The induction of a few VIN genes in response to low temperatures suggested that the vernalization pathway in early-maturing rapeseed was attenuated.
The expression patterns of nine rapeseed FLCs were different between HX17 and HX58 ( Figure 10). Two FLC genes were suppressed only in HX17 under cold-acclimation and cold-shock, while three FLC genes were suppressed only in HX58. BnaA10g22080D was suppressed in both HX17 and HX58 under cold shock stresses. All the results suggest that FLC genes vary extensively among the rapeseed populations.
Validation of gene expression patterns by qRT-PCR As shown in Figure 11, all the 4 candidate genes were induced by freezing stress strongly, with or without cold acclimation. Although the fold-changes in their expression detected by sequencing did not exactly match those detected by qRT-PCR, the detected expression patterns were mostly consistent for all the selected genes, confirming the reliability of the RNA-Seq results.
The differences in cold tolerance between HX58 and HX17 As shown in Figure 12, there are some genes involved in JA signaling, Ca 2+ -signaling, plant-pathogen interaction pathway, circadian rhythms pathway and flowering pathway with different expression patterns and levels between HX58 and HX17, and this may suggest that those genes' expression extent and amount lead the difference in cold tolerance between HX58 and HX17.

CONCLUSIONS
This study is the first report of the transcriptome data of two earlymaturing rapeseeds under different low-temperature stresses (cold acclimation and cold shock; chilling and freezing) using HiSeq PE150 of Illumina. A total of 47,328 DEGs were identified in rapeseed treated with different low temperatures. Further analysis of these DEGs showed that the low-temperature response was a complex process in rapeseed. Many KEGG pathways were enriched following low-temperature treatments, including the primary metabolism, plant hormone signal transduction, plant-pathogen interaction pathway, circadian rhythms and so on. These results provide a reference for our understanding of the rapeseed adaptation to different low-temperature environments and provide new insights into the molecular mechanisms of different low-temperature tolerance in rapeseed.

ACKNOWLEDGMENTS
This work was supported by Institutions of higher learning innovation ability enhancement (30552-520190100004) and National Natural Science Foundation of China (31801398). Figure 12 Heatmap analysis of the DEG genes between HX58 and HX17. The expression levels of each gene (FPKM values) in each sample, as indicated by different color rectangles. Red means high expression; Blue means low expression.