Temperature expression patterns of genes and their coexpression with LncRNAs revealed by RNA-Seq in non-heading Chinese cabbage

Background Non-heading Chinese cabbage (NHCC, Brassica rapa ssp. chinensis) is an important leaf vegetable grown worldwide. However, little is known about the molecular mechanisms underlying tolerance for extreme temperature in NHCC. The limited availability of NHCC genomic information has greatly hindered functional analysis and molecular breeding. Results Here, we conduct comprehensive analyses of cold and heat treatments in NHCC using RNA-seq. Approximately 790 million paired-end reads representing 136,189 unigenes with N50 length of 1705 bp were obtained. Totally, 14,329 differentially expressed genes (DEGs) were detected. Among which, 10 DEGs were detected in all treatments, including 7 up-regulated and 3 down-regulated. The enrichment analyses showed 25 and 33 genes were enriched under cold and heat treatments, respectively. Additionally, 10,001 LncRNAs were identified, and 9,687 belonged to novel LncRNAs. The expression of miRNAs were more than that of pri-miRNAs and LncRNAs. Furthermore, we constructed a coexpression network for LncRNAs and miRNAs. It showed 67 and 192 genes were regulated by LncRNAs under cold and heat treatments, respectively. We constructed the flowchart for identifying LncRNAs of NHCC using transcriptome. Except conducting the de novo transcriptome analyses, we also compared these unigenes with the Chinese cabbage proteins. We identified several most important genes, and discussed their regulatory networks and crosstalk in cold and heat stresses. Conclusions We presented the first comprehensive characterization for NHCC crops and constructed the flowchart for identifying LncRNAs using transcriptome. Therefore, this study represents a fully characterized NHCC transcriptome, and provides a valuable resource for genetic and genomic studies under abiotic stress. Electronic supplementary material The online version of this article (doi:10.1186/s12864-016-2625-2) contains supplementary material, which is available to authorized users.


Background
Nowadays, frequent occurrences of abnormal weather events have been observed all over the world, such as drought and extreme temperature. These stresses seriously impact plant growth and crops production [1,2]. Recently, several progresses have been made about the identification of stress-related genes, which potentially are able to increase the plant tolerance [3][4][5]. Understanding the molecular mechanism of the abiotic stresses response is important to improve tolerance using molecular techniques.
Until now, a large number of transcriptome sequencing projects have been conducted in many species. Genome-wide analyses have dramatically improved the efficiency of gene identification [16]. In Arabidopsis, about 30 % of the transcripts were related with abiotic stresses, and 2,409 genes played important roles in cold, salt, and drought stresses [24]. In chrysanthemum, 8,558 dehydration-responsive transcripts were detected using RNA-seq [25]. In wheat, about 2 % of the wheat genes were related with the cold stress [26]. In Populus and switchgrass, heat responsive genes were also identified by transcriptome sequencing [27,28]. In A. mongolicus, 9,309 up-regulated and 23,419 down-regulated genes were identified under cold stress [29].
Brassica rapa contains several subspecies, such as Chinese cabbage (B. rapa ssp. pekinensis), NHCC, and turnip (B. rapa ssp. rapa) [30,31]. The genome of Chinese cabbage had been sequenced, however, there is little information about the NHCC genome and gene dataset. Therefore, we conducted de novo assembly and gene annotation without prior genome information in this study. NHCC is one of the most important vegetables in China, and now is cultivated extensively worldwide. It is inevitable injured by low or heat stresses, which can directly lead to the production decrease and affect edible quality. The heat stress can affect the photosynthesis, and even induce the occurrence of several diseases, such as downy mildew, soft rot and virus diseases. The physiological change of temperature response mediated by several genes has been reported in model plants [32,33]. However, little is known about the temperature-regulated genes and the related pathways in NHCC.
In this study, we conducted the comprehensive characterization for NHCC using RNA-seq, and explored the effect of low and heat temperature on global change. We identified several most important genes in temperature response, and discussed their regulatory networks and crosstalk in cold and heat stresses. Using Illumina sequencing technology, we generated over 85 billion base of high quality sequence, and identified a larger number of differentially and specifically expressed transcripts. Furthermore, we also identified lots of LncRNAs, and constructed the coexpression network of LncRNAs and protein encoding genes using this transcriptome dataset.

Results and discussion
RNA sequencing and de novo assembly of NHCC transcriptome To obtain a global overview of NHCC transcriptome under different temperature treatments, we constructed and sequenced 15 RNA-Seq libraries, including cold treatments (4, 0 and -4°C), heat treatment (44°C), and normal condition (25°C). For each temperature, three samples as the biological replications were sequenced using Illumina HiSeq™ 2000. The base quality of reads was checked using FastQC (Additional file 1: Figure S1). We used relatively stringent criteria for quality control by removing the reads with adaptors and the low quality. Finally, 790,269,418 clean pair-end (PE) reads consisting of 71.12 billion nucleotides (nt) were obtained with an average GC content of 47.30 % (Table 1, Additional file  2: Table S1). After the first assembly, 1,596,012 contigs were obtained for all libraries, and the total length over 542.8 Mb (Table 1). The contigs were further joined into136,189 unigenes using paired-end information and gap filling process. The total length of all unigenes was 153.1 Mb, and the mean length of unigene was 1124 bp (Table 1, Additional file 2: Table S2). The PE sequencing not only increases the depth, but also improves de novo assembly efficiency. The N50 achieved 1705 bp, which was larger than most plants de novo assembled by RNA-Seq, such as radish (1095 bp), wax gourd (1132 bp), and celery (1088 bp) [34][35][36]. This phenomenon indicated that the high quality and accuracy of our assembled transcripts. Based on FRKM, we measured the correlation of three repeats for each temperature. The results showed that there was a good correlation among three repeats. The pearson's correlations of almost all comparisons were larger than 85 % (Fig. 1, Additional file 1: Figure S2).   Table  S3). The size distribution of BLAST-aligned coding sequence (89.31 %) and predicted proteins are analyzed (Additional file 1: Figure S3a,b). The remaining unigenes that did not match these databases were analyzed by three programs to predict coding regions. Finally, 2793, 2491, and 3119 coding sequences were predicted by ESTScan, CPC, and CNCI programs, respectively (Additional file 1: Figure S3c, Additional file 2: Table S4). The venn diagram showed that there were 684 coding sequences predicted by these three programs, so these genes were relatively reliable as coding genes (Additional file 1: Figure  S3d). A total of 105,217 coding transcripts were predicted in our study. Then we aligned these unigenes with the proteins of Chinese cabbage (E-value <10 -10 , identity >70 %). The results showed that 93,046 unigenes could align to the 3,2640 Chinese cabbage proteins (Fig. 2a). In addition, we found that over 70 % NHCC transcripts could match with more than 1 Chinese cabbage genes (Fig. 2b). This phenomenon might be caused by the genome duplication of B. rapa. The sequences without a homologous hit might represent novel genes in the genome, and some of them might be the specifically expressed in NHCC for temperature treatments. In addition, they also might be the non-coding, alternative transcription, lineage-specific or high allelic variant unigenes. For Nr annotations, 104,363 unigenes matched in this database (Additional file 2: Table S3). The result indicated that 89.20 % of the top hits showed strong homology with the E-value < 1E-15 (Additional file 1: Figure  S4a). The distribution pattern showed that 87.90 % of unigenes had a similarity higher than 60 % (Additional file 1: Figure S4b). The majority annotated unigenes were corresponded to the known plant genes, with 41.50 and 39.60 % matching with A. thaliana and A. lyrata, respectively (Additional file 1: Figure S4c). A total of 96,314 unigenes were assigned at least one GO term, and all GO terms were classified into three groups and further divided into 55 functional subgroups (Additional file 1: Figure S5). Overall, only 45,750 unigenes were assigned to COG classification (Additional file 2: Table  S3). Among 25 COG categories, the cluster for 'general functions prediction only' (36.84 %) represented the largest group, followed by 'Transcription' , and 'Replication, recombination and repair' (Additional file 1: Figure S6). To identify the biological pathways activated in NHCC, the assembled unigenes were annotated with KEGG. A total of 66,419 unigenes were significantly matched in this database, and were assigned to 128 KEGG pathways (Additional file 2: Table S3,5). The result showed that three largest pathway groups were metabolic pathways (ko01100, 21.12 %), biosynthesis of secondary metabolites (ko01110, 9.72 %), and plant-pathogen interactions (ko04626, 7.27 %). Following these three groups, the plant hormone signal transduction (ko04075) was about 6.60 % of all annotated genes. The level 1 of this pathway was 'Environmental Information Processing' , and the level 2 of it was 'Signal transduction' in the KEGG database.

Temperature-dependent gene expression patterns identified by RNA-Seq in NHCC
To view the gene expression, all genes were divided into three categories, including highly (FPKM >50), medium (5 < FPKM ≤50), and lowly (FPKM ≤5) expressed in each library. The results showed that most genes belonged to lowly expressed, followed by medium, and highly expressed (Additional file 1: Figure S7). The purple line shows the cumulative expressed gene number as the library number increased, and 134,980 genes were detected by all libraries.
To evaluate the temperature decrease course and temperature-dependent transcriptomic activities during cold-resistance process in NHCC, we performed a Fig. 2 The alignment analysis for NHCC RNA-Seq transcripts and Chinese cabbage (CC) proteins. a The summary of the aligned and un-aligned transcripts. b The distribution of the aligned transcripts of the NHCC temperature decrease course differential gene expression analysis by comparing any two adjacent cold treatments, using the higher temperature treatment as the denominator. There were 27 (3 3 ) possible patterns, including those that increased across all treatment boundaries, termed'up-up-up' (UUU); those that were similar across all boundaries, termed 'maintain-maintain-maintain' (MMM); and those that decreased across all boundaries, termed 'decrease-decrease-decrease' (DDD). The results showed that genes were non-randomly represented across all patterns, and the overall temperaturedependent patterns are analyzed (Fig. 3, Additional file 2: Table S6). Only few genes continuously decreased (DDD) or increased (UUU) in expression accompanying temperature decrease, and the number were 3 (0.002 %) and 11 (0.008 %), respectively. However, majority genes (79,852, 58.63 %) belonged to the MMM pattern, and the expression almost unchanged over the temperature decrease process. In addition to MMM, the MMU (15,675) and DMM (13,999) were also contained more genes than other expression patterns.

Differential expressed genes detection and compare them among each treatment
To identify the temperature respond genes, 14,329 DEGs were detected between each temperature-treated and control library (FC >2 and q-value >0.8) (Fig. 4). All DEGs were used for clustering analysis, and obtained a well cluster results (Fig. 5a). Three repeats for each temperature got together, and 0, 4 and 25°C formed one group, while -4 and 44°C formed another group. On the whole, all DEGs were divided into three groups, and defined as I, II, and III. Most DEGs in group I belonged to up-regulated genes under -4°C treatment, while similar in other four temperatures. In group II, the expression of most DEGs under 25 and -4°C was lower than that under other three temperatures. In group III, most DEGs had relatively low expression, except a few genes under 44°C.
To further survey the interaction of these treatments for DEGs, we constructed venn diagram using DEGs of each treatment. There were 5445, 7430, 8056, and 488 DEGs under 4, 0, -4, and 44°C, respectively. Among of them, 1075, 2865, 5254, and 286 belonged to each treatment-specific DEGs (Fig. 5b). Interestingly, we found 10 DEGs were detected in all treatments, including 7 up-regulated and 3 down-regulated genes (Additional file 1: Figure S8). The functional annotation showed that most of them belonged to the stresses related protein, such as LEA14 and KIN2 (Table 2). In addition, we also conducted qRT-PCR experiment to verify the accuracy of the RNA-Seq. The results showed that the expression trends of all genes were consistent with the RNA-Seq, and most genes were also significant differently expressed (p-value < 0.01) (Additional file 2: Table S7, Additional file 1: Figure S9). Among all DEGs,  Figure S10). Among these DEGs, the most treatment-specific up-regulated genes (4804) were detected under -4°C, and the most treatment-specific down-regulated genes (1377) were detected under 0°C (Additional file 1: Figure S8).

The enrichment analyses revealed most DEGs related with cold and heat stresses
To understand the function of DEGs, we have conducted the GO enrichment analyses using all unigenes as background (Additional file 2: Table S8). Under 4°C, several cold related GO categories were significantly enriched, such as response to desiccation, response to cold, response to temperature stimulus, and cold acclimation. The photosynthesis, light harvesting, and translation categories were the top enrichment under 0°C. Of course, the cold related categories were also enriched at the top 10 categories. Under -4°C, the mainly enrichment categories were secondary metabolite biosynthetic process, S-glycoside metabolic process, glycosinolate, and glucosinolate metabolic process. However, we did not detect directly cold related categories as 4 and 0°C at the top 10 categories. This phenomenon indicated the regulatory mechanism might exit several differences between chilling (<7°C) and freezing (<0°C) temperatures, which was also consistent with the previous report [37]. Under 44°C, the mainly enrichment categories were response to high light intensity, heat acclimation, photosynthesis, and response to heat. In addition, the photosynthesis was also enriched under heat stress, which indicated that all temperature stresses could affect plant photosynthesis. We also analyzed genes belonged to the GO enrichment categories. Among all treatments, the most specific genes (2713) were found under -4°C, followed by 0°C (1115), 4°C (528), and 44°C (79). In addition, 90 genes were identified under all three cold treatments, and 5 genes were detected by all the cold and heat treatments (Additional file 1: Figure S11a).
In addition, we mapped DEGs to terms in KEGG database to identify significantly enriched pathways. Among the mapped pathways, 20, 11, 22, and 5 pathways were significantly enriched (Qvalue < 0.01) under 4, 0, -4, and 44°C treatments, respectively (Additional file 1: Figure  S12, Table S9). Notably, common enrichments were observed in photosynthesis pathway, metabolic pathway, and photosynthesis-antenna proteins pathway in all treatments. This results indicated that the cold and heat stresses affected the expression of genes involved in these pathways. Most enriched pathways were also detected by the previous reports, which partly reflected the accuracy of our results [24,27,29]. Interestingly, we found transcripts involved in protein processing in endoplasmic reticulum pathway were significantly enriched under 44°C treatment, while it did not enrich in cold treatment. This phenomenon indicated that this pathway might only play roles in heat resistance. We also   (44). Forty-five genes were identified under all three cold treatments, and 2 genes were detected by all the cold and heat treatments (Additional file 1: Figure S11b). Furthermore, we also surveyed the enrichment genes identified by combing the GO and KEGG databases. The results showed that 25 and 33 genes were enriched in the two databases under cold and heat treatments, respectively (Additional file 1: Figure S11c,d). These enriched genes will greatly enhance the potential utilization in cold and heat stresses of NHCC.

Identification of abiotic stresses related transcription factors from DEGs
Given that TFs have a major effect on the network of temperature-responsive genes, we also identified the temperature-inducible TFs. Overall, the number of Dehydrin, Chloroa_b-bind, p450, AP2, PSI_PsaH, and EF-hand was more than other TFs in three cold treatments (Additional file 1: Figure S13a,b,c,d,e,f ). However, many GST were identified under -4°C, while they were absent under 0 and 4°C. This phenomenon indicated that GST might play roles in the cold resistance below 0°C. Interestingly, HSP70 and HSP20 were identified under 0 and 4°C, indicating that there was a certain inherent association between cold and heat regulation. Under 44°C, HSP20 was significant enriched, which revealed that it played important roles in heat-resistance regulation (Additional file 1: Figure S13g,h). In all treatments, P450 and Chloroa_b-bind TFs were enriched, which indicated that cold and heat stresses had great impact on plant photosynthesis, and thus might affect crop yields. This revealed that they were related with plant photosynthesis, which was also consistent with previous reports [38,39]. Among all treatments, the most specific TFs were identified under -4°C (1781), followed by 0°C (775), and 4°C (313) in GO enrichment (Additional file 1: Figure S14a). In addition, 51 TFs were identified under all three cold treatments, and 2 genes were detected by all cold and heat treatments. Among all KEGG categories, the most specific TFs (819) were found under -4°C, followed by 0°C (368), and 4°C (128) (Additional file 1: Figure S14b). Twenty-nine genes were identified under all three cold treatments, and 1 genes were detected by all cold and heat treatments. Combing GO and KEGG enrichment analyses, 17 and 31 TFs were enriched only under cold and heat treatments, respectively (Additional file 1: Figure S14c,d, Table 3, Additional file 2: Table  S10).
AP2/ERF TFs mainly contained two subgroups, including CBF and DREB2. They interacted with DRE/ CRT cis-element and regulated ABA-independent gene expression. The CBF controlled many gene expression under several stresses, such as drought, salinity and freezing stresses. The DREB2 mainly affected gene expression under osmotic and heat stresses, while it slightly played role in cold stress [1]. Although the functions of these TFs were well explained, the interaction among them was rarely reported, especially in B. rapa. Therefore, we conducted correlation analysis for these TFs, and constructed the interaction network of them using expression values. A total of 38 unigenes were detected using BLAST alignment with ABFs, CBFs, and DREB2 of Arabidopsis (Additional file 2: Table S11). The pearson correlations coefficient (PCC) between two of these TFs were calculated using the expression value. Then the interaction network was constructed using part connections with the PCC larger than 80 %. Finally, this network contained 95 connections, including 82 positive and 13 negative connections (Fig. 6). This phenomenon revealed that most connections belonged to positive relationship among of these TFs. However, we noted that CL258.Contig16, a ABF TF, had negative connections with two DREB2A (CL13726.Contig1, CL13726.Contig2), CBF1 (CL1909.Contig9), and CBF3 (CL1909.Contig10).
To analyze the correlation of abiotic stresses and these TFs, we collected mainly cold and heat stresses related genes according to previous reports [1,12]. Then the candidate heat and cold related genes in NHCC were identified using BLAST alignment with the collected genes. We calculated PCC of these candidate genes and CBF or DREB2. The PCC values larger than 0.8 were selected to construct interaction network (Additional file 1: Figure S15). This network showed that there were more negative connections in DREB2 than that in CBF, which only contained four negative connections. Among all connections, we identified 228 transcripts, which had high PCC (>90 %) with DREB2, such as HSF, LEA, and MYB102. By CBFs, 96 transcripts were also identified, such as COR6.6, WD40, and ABF4. Moreover, 91 transcripts were detected by both DREB2 and CBF, including GRP7, P450, PP2C, and SRK2E (Fig. 7, Additional file 2: Table S12).

Identification and characterisation of NHCC LncRNA using RNA-seq
To identify potential LncRNA in NHCC, all sequences from NHCC transcriptome dataset were used. Based on previous reports [40][41][42], we designed the pipeline for LncRNA analyses. Finally, 10,001 LncRNAs were identified after a series filtering, including transcript length, coding potential, ORF size, and the exclusion of other Fig. 6 The interaction network for DREB2, CBF, and ABF TFs, which was constructed based on PCCs. The expression value of TFs at each treatment was used for calculating the PCCs. The blue lines represent the positive correlation, while the red lines represent the negative correlation ncRNA. Among these LncRNAs, 9,687 belonged to the novel LncRNAs. In addition, we also identified 50 pri-miRNA through comparing with Rfam and miRBase databases for comparative analyses.

LncRNA gene expression analyses and differently expressed LncRNAs identification
LncRNA is a group of endogenous RNAs that function as regulators of gene expression, which are involved in developmental and physiological processes [43]. They are longer than 200 bp, and several of them can also act as primary transcripts for the production of short RNAs [44]. We assessed the expression pattern under different temperatures using the expressed LncRNA (FPKM >0). A total of 2,236 LncRNAs were expressed in all the treatments. We observed that 73, 14, 107, 468, and 244 LncRNAs belonged to the temperature specificity for T25, T4, T0, TM4, and T44, respectively (Additional file 1: Figure S16). However, a hierarchical clustering of samples showed most LncRNA was low expression in each temperature (Additional file 1: Figure S17). In this study, 50 pri-miRNA, 10,001 LncRNAs, and 121,744 protein coding transcripts were identified. The average, maximum, and median expression values of these three type transcripts were calculated for comparative analyses. The results showed that the expression of protein coding transcripts were more than that of pri-miRNA and LncRNA (Fig. 8a,b,c, Additional file 2: Table S13). Most pri-miRNA and LncRNA had relatively low expression, which was consistent with the previous reports [45][46][47], indicating it was a common property of LncRNA. We also investigated the temperature specific expressed transcripts (SETs) among these three types transcripts. The results showed that 10.0 % pri-miRNA, 9.1 % LncRNA, and 2.6 % protein coding transcripts were detected as SETs. For pri-miRNA, 3 SETs were detected in 44°C, while no SET was found in 4 and 25°C (Fig. 8d). For LncRNA, the most SETs were identified in -4°C (468), followed by 44°C (244), and 0°C (107) (Fig. 8e). The similar trends were also found in the protein-coding transcripts, except the inverse of 0 and 25°C (Fig. 8f).
To further investigate the share or specific of these treatments for differently expressed LncRNAs (DELs), we conducted the venn diagram analyses. The results showed that 91, 418, 441, and 34 specifically expressed LncRNAs were identified under 4, 0, -4, and 44°C, Fig. 7 Model for DREB2 and CBF signal regulation in response to heat and cold stresses. The overlapping or specific abiotic stresses related genes were detected, which had high PCC (>90 %) with DREB2 or CBFs respectively (Additional file 1: Figure S18a). However, we not detected the share LncRNAs among all these treatments. Among of these DELs, most of them were downregulated, and there were the most down-regulated LncRNAs under 0°C (1103) than other three treatments (9~934) (Additional file 1: Figure S18b). Interestingly, the DEL number under three cold stresses (964~1327) was more than that of 44°C, which only contained 44 DELs. However, most LncRNAs (79.55 %) were upregulated under 44°C, while most LncRNAs belonged to the down-regulated under other three cold treatments (56.21 %~96.89 %).

Construct the coexpression network between LncRNAs and protein-coding genes
We constructed a coexpression network for LncRNAs and protein-coding genes according to the previously proposed method [40]. For each treatments and for each pair of genes (LncRNA or protein-coding), we computed the PCCs of expression patterns using the expression values. We found that approximately 57.59 % were positive connections, and 42.40 % were negative connections (Additional file 1: Figure S19a). Among all connections, the PCCs of 29.27 % were between -0.4 to -0.2, and followed by 26.29 % connections were between 0.8 to 1. Furthermore, we identified 65,568,352 connections between protein-coding and protein-coding genes, among which, 42.11 % were positive connections, and 57.89 % were negative connections (Additional file 1: Figure  S19b). However, the opposite result was found between LncRNA and LncRNA connections, and the values were 50.12 and 49.88 % for positive and negative connections, respectively.
To be more accurate and intuitive showed the relationship between the LncRNAs and protein-coding genes, we selected the connections with the high correlation (|PCC| > 0.95) to construct the interaction network. Overall, the whole network constituted by these connections was divided into 8 clusters, including 3 large networks and 5 relative small networks (Additional file 1: Figure S20a). In the cluster 1, about 67 % connections with the PCC > 0.95, and~33 % with the PCC = 1. The PCCs of all connections in the cluster 2 were larger than 0.95, but less than 1. Only three connections with PCC less than -0.95 were located in the cluster1 and cluster3. Most connections (362,213) in the networks belonged to positive, and only 636 connections were negative correlation (Additional file 1: Figure S20b). This phenomenon was also found in LncRNA vs LncRNA and protein-coding vs protein-coding genes. Among all of these three types, 99.72 % connections were positive, and 0.28 % were negative connections.
To further analyze the correlation between LncRNA and protein-coding genes under temperature stresses, we annotated the function of the target genes. A total of 67 target genes were regulated by LncRNAs under all three cold treatments, and comparing them with Arabidopsis (Additional file 2: Table S14). The annotation showed that most of them belonged to the cold respond proteins, such as CBF1, COR6.6, and LEA14. Similar, we Fig. 8 The comparative analyses of LncRNA, pri-miRNA, and protein-coding genes expression pattern. a, b, c The average, maximum, and median expression analysis of LncRNA, pri-miRNA, and protein coding transcripts in NHCC using boxplot. d, e, f The temperature specific expressed transcripts of LncRNA, pri-miRNA, and protein coding transcripts in NHCC identified 192 target genes of LncRNAs under 44°C treatment, and most of them belonged to the heat respond genes, such as HSP, LTP, and CBF4 (Additional file 2: Table S15). Furthermore, we also found several target genes of LncRNAs under cold and heat stresses, such as KIN2. This phenomenon indicated that they might play important roles in the corsstalk between cold and heat stresses responses.

Conclusion
In this study, approximately 790 million paired-end reads representing 136,189 unigenes with a total length of 153.1 Mb were obtained. Only few genes were DDD or UUU expression patterns, and majority genes belonged to MMM pattern during the temperature decrease process. 14,329 DEGs were detected between at least one treatment and control library. Among which, 10 DEGs were identified in all treatments, including 7 up-regulated and 3 down-regulated genes. The enrichment analyses demonstrated that most temperature related categories were discovered under cold and heat treatments. Among the enrichment categories, 25 and 33 genes were identified in both CO and KEGG databases under cold and heat treatments, respectively.
Totally, 10,001 LncRNAs were identified from NHCC transcriptome dataset, and 9,687 belonged to novel LncRNAs. The analyses indicated the expression of protein coding transcripts were higher than that of pri-miRNA and LncRNA. We constructed a coexpression network for LncRNAs and protein-coding genes. A total of 67 and 192 target genes were regulated by LncRNAs under three cold and heat treatments, respectively. Furthermore, we also identified several shared target genes of LncRNAs under cold and heat treatments, which indicated that they might play important roles in the corsstalk between cold and heat stresses.
In conclusion, we conduct comprehensive analyses for cold and heat stresses in NHCC using RNA-seq, and identified numerous differentially and specifically expressed transcripts. Many important genes and TFs response to treatment stress were detected, and their crosstalk between cold and heat stress responses was discovered. In addition, we also identified large number of LncRNAs, and constructed the coexpression network of LncRNAs and protein encoding genes. This study provides a platform for elucidating physiologic responses to low and high temperature in B. rapa.

Plant materials, growth conditions, and treatments
The NHCC advanced inbred line, 'Suzhouqing' , was used in this study. The surface-sterilized seeds were grown in pots containing a soil: vermiculite mixture (3:1) in a controlled-environment growth chamber programmed for 16/8 h at 25/18°C for day/night. Seedlings at the five-leaf stage were transferred to growth chambers set at -4, 0, 4°C as cold treatments, 25°C as control, and 44°C as heat treatment under 4 h for RNA-Seq and qRT-PCR. Three samples of each treatment were generated from different batches of plants for three biological replicates. All leaf samples collected from control and treated plants were washed with distilled water, immediately frozen in liquid nitrogen, and stored at −80°C for RNA extraction.

RNA extraction for transcriptome sequencing and RT-PCR validation
The RNA was isolated from leaves using RNA kit (Tiangen, Beijing, China) according to manufacturer's instructions. RNA samples were treated with RNase free DNase I to avoid DNA contamination. The RNA was reverse transcribed into cDNA using Prime Script RT reagent Kit (TaKaRa, Kyoto, Japan). The cDNA libraries were constructed using an mRNA-seq assay with a fragment length range of 200 bp (±25 bp). Finally, the library was sequenced for paired-end reads of 90 bp using Illumina HiSeq™ 2000 platform, which was performed by the Beijing Genomics Institute (BGI) (http://www.genomics.cn/ index). For qRT-PCR, the actin gene (AF111812) was used as an internal control to normalize the expression level of the target gene. Primer 5.0 designed the specific primers according to gene sequences. The qRT-PCR assays were performed with three biological and technical replicates. Each reaction was performed in 20-μL reaction mixtures containing a diluted cDNA sample as template, SYBR Premix Ex Taq (2×) (TaKaRa, Kyoto, Japan) and gene-specific primers. qRT-PCR was performed according to our previous report [48]. The comparative Ct value method was adopted to analyze the relative gene expression. RNA expression levels relative to actin gene were calculated as 2 -ΔΔCT according to a previous analysis [48,49].

Data filtering and de novo assembly
Raw reads generated by Illumina Hiseq™ 2000 were initially processed to get clean reads through the following three steps. i) Remove reads with adaptors contamination; ii) Discard reads with ambiguous sequences "N" larger than 5 %; iii) Remove low quality reads, which contained more than 20 % Q <20 bases [50]. In addition, we used FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) to check and visualize the quality of RNA-seq reads (Additional file 1: Figure S1). After filtering, all clean reads were assembled using a de novo assembly software Trinity [51]. Firstly, clean reads with a certain length of overlap were combined to generated contigs. Then, the paired-end reads were realigned to contigs to obtain unigene, which could identify different contigs in the same transcript and ensure the interval among these contigs. The contigs in one transcript were assembled by Trinity and gained the sequence not being extended on either end, which defined as unigene [34]. Then, the TGICL program was used to delete redundant unigene and further assembled all unigenes to form a single set of non-redundant unigenes [52].
Gene expression quantification and differential expression analyses RNA-Seq reads were aligned to the assembled transcripts using TopHat pipeline with the built-in Bowtie mapping program [53]. The expression of all unigenes was estimated by calculating read density as 'fragments per kilobase of exon per million mapped reads' (FPKM) [54]. The DEGs (FC >2, q-value >0.8) between normal and stress-treated conditions were identified using NOISeq (http://www.bioconductor.org/) [55]. GO enrichment analyses were performed using Blast2GO [56]. The temperature-dependent gene expression patterns were analyzed according to the previous report [57]. The comparisons were made between two adjacent temperatures-that is, 25°C vs 4°C, 4°C vs 0°C, and 0°C vs -4°C. A gene with FC > =2 was grouped into 'up' pattern, a gene with FC < =0.5 was grouped into 'decrease' , and the remaining genes were grouped into 'maintain'. Therefore, a gene was grouped to 1 out of 27 patterns, ranging from up-up-up (UUU), maintain-maintain-maintain (MMM), to decrease-decrease-decrease (DDD).

Functional annotation and classification of the transcripts
All assembled transcripts were annotated with the publicly available protein databases, including Nr (http:// www.ncbi.nlm.nih.gov), GO (http://www.geneontology.org), COG (http://www.ncbi.nlm.nih.gov/COG), Swiss-Prot protein (http://www.expasy.ch/sprot), and KEGG (http:// www.genome.jp/kegg) databases using BLAST (E-value <10 -5 ). Then, the best alignments were used to decide sequence direction and to predict coding regions of the unigenes. ESTScan software was used to decide sequence direction and coding regions when a unigene unaligned to none of the above databases [58]. WEGO software was used to conduct GO classification for understanding the distribution of gene function [59]. The unigenes were also aligned to COG database to predict and classify possible functions. In addition, KEGG was used to annotate the pathway of the unigenes.

LncRNA detection
To de novo detect LncRNAs using RNA-seq, we developed a flowchart according to previous reports with slightly modification (Additional file 1: Figure S21) [40,41,60,61]. We applied several filters to ensure reliability of LncRNAs. Firstly, all unigenes (136,189) were annotated using BLAST (E-value <10 -5 ) alignment with NR,NT,Swiss-Prot,KEGG,COG, and GO databases. There were 14,445 unigenes un-annotated by any protein databases above mentioned. Among which, 14,300 unigenes selected for transcripts greater than or equal to 200 bp. Secondly, Coding-Non-Coding Index (CNCI, http://www.bioinfo.org/software/cnci) was applied on all candidate unigenes in order to distinguish proteincoding and non-coding sequences [62]. The unigenes with score <0 were defined as non-coding. In addition, Coding Potential Calculator (CPC, http://cpc.cbi.pku.e du.cn) was also used for identifying all candidate transcript models in order to assess their coding potential by a second independent method [63]. In order to extract potential non-coding transcripts with a high reliability from our dataset, all transcripts with a score (CPC < -1) were retained as potential non-coding. By combining these two methods, 10,930 unigenes were identified as potential non-coding RNAs. Thirdly, we discard transcripts with an ORF greater than 100 amino acids by ORF-Predictor (E-value <10 -5 ), and 813 unigenes were discard.
Finally, 10,117 candidate transcripts were identified and compared against several non-coding RNA databases, including Rfam, miRBase, NONCODE with designated threshold value (E-value <10 -5 , identity >90 %) by BLAST [64][65][66]. Candidate transcript models with known protein motifs were discarded. We obtained 50 pri-miRNA sequences by comparing with the miRBase and Rfam databases. To identify LncRNA, we filtered other non-coding RNAs through comparing with Rfam and NONCODE databases. Finally, 10,001 unigenes were identified as LncRNA based on a series of analyses above mentioned. Furthermore, among which, 9,687 belonged to the novel novel LncRNAs by comparing with NONCODE databases.

Statistical analysis
The differential expression levels of genes under the different temperatures were clustered using Cluster program (http://bonsai.hgc.jp/~mdehoon/software/cluster), and visualized using Tree View software (http://jtree view.sourceforge.net/) [67]. Using in-house Perl script, PCCs were calculated for correlation studies, including the three repeats correlation, LncRNA and target mRNA correlation, and key candidate mRNA-mRNA correlation. The coexpression interaction networks were constructed using Cytoscape (http://www.cytoscape.org/) according to PCC [68]. The numbers of specific and common genes were plotted using Venn diagram in R package [69].

Ethics and consent to participate
This article does not contain any human or animals data performed by any of the other committee.