Comparative Transcriptome Analysis Reveals Heat-Responsive Genes in Chinese Cabbage (Brassica rapa ssp. chinensis)

Chinese cabbage (Brassica rapa ssp. chinensis) is an economically and agriculturally significant vegetable crop and is extensively cultivated throughout the world. Heat stress disturbs cellular homeostasis and causes visible growth inhibition of shoots and roots, severe retardation in growth and development, and even death. However, there are few studies on the transcriptome profiling of heat stress in non-heading Chinese cabbage. In this study, we investigated the transcript profiles of non-heading Chinese cabbage from heat-sensitive and heat-tolerant varieties “GHA” and “XK,” respectively, in response to high temperature using RNA sequencing (RNA seq). Approximately 625 genes were differentially expressed between the two varieties. The responsive genes can be divided into three phases along with the time of heat treatment: response to stimulus, programmed cell death and ribosome biogenesis. Differentially expressed genes (DEGs) were identified in the two varieties, including transcription factors (TFs), kinases/phosphatases, genes related to photosynthesis and effectors of homeostasis. Many TFs were involved in the heat stress response of Chinese cabbage, including NAC069 TF which was up-regulated at all the heat treatment stages. And their expression levels were also validated by quantitative real-time-PCR (qRT-PCR). These candidate genes will provide genetic resources for further improving the heat-tolerant characteristics in non-heading Chinese cabbage.


INTRODUCTION
The species Brassica rapa includes various vegetable crops, such as turnip, field mustard, and Chinese cabbage. Non-heading Chinese cabbage (B. rapa ssp. chinensis) is an economically and agriculturally significant vegetable crop that is cultivated extensively worldwide. Non-heading Chinese cabbage originated from China and has a long cultivation history (Song et al., 2014a). The adaptable growth temperature for Chinese cabbage ranges from 18 to 22 • C and its production is usually impaired by heat stress in many regions . Heat stress, triggered by high environmental temperature, can affect plant performance, leading to severe retardation in vegetative growth, yield depression and even death (Caers et al., 1985;Barnabás et al., 2008;Song et al., 2014b). One of the phenotypes is leaf etiolation and bleaching with clearcut inhibition of photosynthetic activity (Wang L. et al., 2011). Based on global climate model analysis, the predictions suggest that global warming and extreme heat events will threaten food safety by reducing crop production in the future (Battisti and Naylor, 2009;Rosenzweig et al., 2014). Therefore, discovery of the genes related to heat tolerance and investigating the molecular mechanism play important roles in genetic improvement of crops.
Signal transduction components, transcription factors (TFs) and proteins associated with the metabolism of stress-generated reactive oxygen species (ROS) are mainly responsive to the high temperature (Grover et al., 2013). Identification of heat responsive genes from suitable genotypes can give some insights into the heat-tolerance mechanism. Transcript profiling of two Chinese cabbage (B. rapa ssp. pekinensis) inbred lines showed that many genes are affected by high temperatures including heat shock proteins (HSPs), genes associated with membrane leakage and enzymes involved in ROS homeostasis . Previous studies have also shown that genes involved in oxidative stress, protection of proteins, programed cell death, biotic stress responses and metabolism were differentially expressed under high temperatures .
In recent years, a great deal of attention has been paid to the elucidation of the mechanisms of heat-tolerance for breeding heat-resistant cultivars of Chinese cabbage and other important crops. In Chinese cabbage, heat-responsive miRNA and nat-siRNAs were identified, and some of these small RNA were upregulated under heat stress (Yu et al., , 2013. Another non-coding small RNAs, chloroplast small RNAs (csRNAs), were also reported to be highly sensitive to heat stress and results showed that high temperature suppresses the production of some csRNAs (Wang L. et al., 2011). HSPs have also been assumed to play a central role in the heat stress response and in acquired thermotolerance in plants (Kotak et al., 2007). Using microarray, the transcript profiles of two Chinese cabbage inbred lines were studied and some heat responsive genes, such as heat shock proteins and MYB41, were identified . Lately, a total of 9687 novel lncRNAs were identified and 192 genes were regulated by these lncRNAs in Non-heading Chinese cabbage NHCC under heat treatment. However, the responsive genes to heat stress in non-head Chinese cabbage remain not fully understood.
Recently, the whole genome of B. rapa was sequenced and annotated, providing a solid foundation to study the expression of heat stress responsive genes (Wang X. W. et al., 2011). Highthroughput sequencing methods such as Illumina SOLEXA, ABI SOLiD, and Roche 454 have observably increased the efficiency and reduced the cost of sequencing, making the study of transcriptomes easier and more feasible. RNA sequencing (RNA-Seq), a high-throughput sequencing method, is developed to analyze the transcriptome either with or without genome information. It's an efficient tool to promise simultaneous estimation of abundance and new transcript discovery (Cloonan et al., 2008). By RNA-Seq, researchers can obtain almost all of the expressed genes, especially genes with very low abundance. Therefore, genes with abundant expression differences and interesting pathways can be analyzed exhaustively. Compared with microarray, RNA-seq has very low, if any, background signal, and has also been shown to be highly accurate for quantifying expression levels (Wang et al., 2009). RNA-Seq has been applied for transcriptome analysis in many plant species under abiotic stress, including carnation (Wan et al., 2015), banana (Yang et al., 2015) and Brassica juncea (Bhardwaj et al., 2015).
In the present study, genome-wide analysis of gene expression in leaves of Chinese cabbage under heat stress was performed using RNA-seq based on the Illumina HiSeq2000 platform. The aim of the study is to identify heat responsive genes and provide new insight into the heat tolerance in non-heading Chinese cabbage.

Plant Materials
A heat sensitive variety "GHA" (S) and a tolerant variety "XK" (T) of non-heading Chinese cabbage (B. rapa ssp. chinensis) were used in this study. Sterile seeds were sown in pots and germinated in a growth chamber. For heat treatments, 3-week-old seedlings were grown at 37 • C (high temperature) for 1, 6, 12, and 18 h. For control samples (CK), the leaf discs were incubated at 25 • C. Then, leaves of the two varieties under CK and heat treatments were collected, frozen immediately in liquid nitrogen, and stored at −80 • C for use.

RNA Isolation, cDNA Library Construction and Illumina Sequencing
Total RNA was isolated from the leaves using Trizol Reagent (Invitrogen, USA) according to the manufacturer's instructions. The extracted RNA was treated with DNase I (Promega, USA) to remove the contaminated DNA. The RNA quality and purity were varified by Nanodrop 2000 and electrophoresis on 1.0% agarose gels. Using poly-T oligo-attached magnetic beads, mRNAs were purified from the total RNA. Then, the mRNAs were fragmented and cDNA was synthesized using random hexamer, DNA polymerase I and RNase H. The double-stranded cDNAs were purified and ligated to adaptors for Illumina paired-end sequencing. The quality and quantity of the library was verified using an Agilent 2100 Bioanalyzer and ABI real time RT-PCR, respectively. The cDNA libraries were sequenced using the Illumina HiSeq2000 platform by the Beijing Genomics Institute (BGI).

Sequence Data Analysis and Annotation
Raw reads in fastq format were firstly filtered and the adaptor sequences and low quality reads were removed. Stringent filtering criterion was carried out to minimize the effects of sequencing errors during the assembly. Briefly, bases with Phred quality score lower than 20 and reads length short than 50 bp were discarded. Reads of 70% bases in one read having high quality scores (≥20) FIGURE 1 | Volcano plots of the transcriptome between "GHA" and "XK" at different heat treatment stages. In the volcano plot, statistical significance (log 10 of p-value; Y-axis) has been plotted against log 2 -fold change (X-axis). Frontiers in Plant Science | www.frontiersin.org were filtered. The genome sequence of B. rapa retrieved from the Brassica Database (BRAD) (http://brassicadb.org/brad/) was used as the reference database (Wang X. W. et al., 2011). All the clean reads were then mapped to the reference genome using TopHat v 1.4.0 (Trapnell et al., 2010). The transcripts abundance was normalized by the FPKM (fragments per kilobase of exon per million fragments mapped) using Cufflinks (Trapnell et al., 2012).

Identification of Differentially Expressed Genes (DEGs)
The annotation of differentially expressed genes (DEGs) was performed by analyzing syntenic relationships between B. rapa and Arabidopsis thaliana. The syntenic orthologs of DEGs in A. thaliana were retrieved from BRAD (http://brassicadb.org/brad/). B. rapa genes sharing the same A. thaliana ortholog were defined as a paralogous set. FDR (false discovery rate) value less than 0.01 and |log 2 (fold change)|≥2 were used to recognize the significance of the gene expression difference. After normalization, hierarchical clustering and k-means clustering analysis of the expression patterns were performed by Mutiexperimental Viewer v4.7 (Saeed et al., 2003).

GO and KEGG Enrichment Analysis
To identify putative biological functions and pathways for the DEGs, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) database were searched for annotation. The results of GO annotation were submitted to WEGO for GO classification (Ashburner et al., 2000). Enrichment analysis of GO and KEGG was used by the AgriGO and KOBAS2.0 packages, respectively (Du et al., 2010;Xie et al., 2011). GO functional enrichment and KEGG pathway enrichment analysis were also tested at a significance cutoff of 0.05 false discovery rate (FDR).

Validation of RNA-seq by Real Time RT-PCR
For the real time quantitative RT-PCR, 1 µg of total RNA was used to synthesize the cDNA using the RevertAid First Strand cDNA Synthesis Kit (Thermo Scientific). Quantitative PCR was performed using the FastStart Universal SYBR Green Master (Roche) according to the manufacturer's instruction on the StepOne Plus Real time PCR Platform (Applied Biosystems). The qRT-PCRs were carried out with the following protocol: 95 • C for 10 min, followed by 40 cycles of 95 • C for 15 s, and at 60 • C for 60 s. The BrActin was used as the internal control, which has been shown to be one of the reference genes in Chinese cabbage. After the amplification, the melting curve was determined for specific product. Three independent biological replicates for each sample and three technical replicates for each biological replicate were analyzed. Significant differences of gene expression level between "GHA" and "XK" were evaluated using a student's t -test. All of the primers for the selected genes are listed in Table S1.

Availability of Supporting Data
The data sets supporting the results of this article are included in the Supplementary Online Materials. The sequencing raw data of this article have been deposited in a SRA database at the NCBI under the accession number: SRP064703.

Overview of the RNA Sequencing
In total, 10 libraries of Chinese cabbage leaves were created, from plants exposed to normal (CK) and high temperature (37 • C) for 1, 6, 12, and 18 h. Totally, 36.8G raw reads were generated ( Table 1). After filtering the reads with low quality, the proportion of clean reads were all above 92% for each library. Of the total clean reads from the ten samples, 40.5-43.28% were perfect match, 26.08-28.15% had no more than five base mismatches, 67.24-70.89% were mapped to unique genome locations and 29.11-32.76% were unmapped reads ( Table 1). All clean reads were aligned to the reference B. rapa genome (v1.5) and a total of 43537 predicted B. rapa genes were annotated.
Differential Expression Profiling of Heat Treatment between "GHA" and "XK" To investigate the gene expression patterns between heat sensitive variety "GHA" and heat tolerant variety "XK", FPKM values were calculated for each sample to normalize the expression under different conditions. Volcano plots were constructed to see how many transcripts were significantly regulated during heat treatment for different time periods (Figure 1). The center of the volcano is a line at which fold change is zero, while both sides of the line are indicating down-regulation (negative vales) and up-regulation (positive values), respectively. And the significantly DEGs are represented by red dots with |log 2 (fold change)|≥2 and FDR value less than 0.001. The results showed that many transcripts were up-regulated in "GHA" at CK and at 1 h after heat treatment (Figures 1A,B). However, more transcripts were up-regulated in "XK" for heat treatment at 18 h ( Figure 1E). A total of 625 DEGs between "GHA" and "XK" were detected by RNA-seq in the present study. Heatmap of these DEGs at all the stages revealed that many DEGs were up-regualted in "XK" (Figure 2A). Clustering analysis showed that large abundant of DEGs were sepcific expressed in each stage (Figures 2B-F Frontiers in Plant Science | www.frontiersin.org Figure S1). Some of the DEGs were induced in "XK" under room temperature (CK) (Cluster 5) and 1 h heat treatment (Cluster 2) (Figures 2B,C). However, other DEGs of clusters 9, 18, and 15 were highly expressed in "GHA" at 1, 6, and 12 h after heat treatment, respectively (Figures 2C-E).
Comparative transcriptome analysis between "GHA" and "XK" showed that many genes were up or down regulated. Then, Venn diagrams illustrated the number of genes uniquely expressed in each sample or genes that were shared between one or more other samples (Figure 3). A total of 625 genes were differentially expressed in the five samples ( Table 2, Table S2). After heat treatment, there were more DEGs with 296, 128, 153 genes at 1, 6, and 12 h than those from other stages, respectively ( Table 2). The highest number of DEGs occurred at heat treatment for 1 h. After 1 h of heat stress, 119 and 177 genes were up-regulated and down-regulated in heat-tolerant variety "XK, " respectively ( Table 2). Some of the DEGs were shared between the two consecutive stages (Figure 3A). A comparison of the number of DEGs between 1 and 6 h of heat stress showed that 41 DEGs were shared at the two stages ( Figure 3A). The number of DEGs simultaneously expressed among 1, 6, and 12 h of heat stress was higher than that expressed simultaneously in the other three successive stages (Figures 3B,C). Moreover, one gene (Bra027596, NAC069) was differentially expressed throughout all the time periods under heat stress ( Figure 3D).

Functional Annotation of the DEGs
To functionally annotate the DEGs of heat-sensitive "GHA" and heat-tolerant "XK" varieties, we aligned all of the DEGs against the Gene Ontology (GO) (Ashburner et al., 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) database. For GO analysis, we annotated genes to three  Figure 4 by WEGO application (Ye et al., 2006). In the CC category, the terms of "cell" (GO:0005623) and "cell part" (GO:0044464) were the enriched components. "Binding" (GO:0005488) and "catalytic activity" (GO:0003824) were the top two MF terms. The most enriched components of the BP category were "cellular process" (GO:0009987), "metabolic process" (GO:0008152) and "response to stimulus" (GO:0009628).
During the initial period of heat stress (CK, 1 and 6 h), most DEGs were enriched in the BP category of "response to stimulus" (Figures 5A,B, Table 3, and Tables S3, S4). After heat stress for 12 h, many DEGs were involved in "programmed cell death" or "innate immune response" process ( Figure 6A, Tables S3, S4). Then, DEGs of "ribosome biogenesis" was enriched under heat stress for 18 h (Figure 6B, Tables S3, S4).
KEGG analysis suggested that 294 networks were involved in heat response, including energy metabolism, carbon metabolism, starch and sucrose metabolism, photosynthesis and others ( Table S5). The top 10 metabolic pathways possibly regulated by heat stress are presented in Table 4.

Validation of RNA-seq Data by Quantitative Real-Time PCR
To verify the RNA-seq expression data, we selected nine genes displaying diverse expression profiles during the heat treatment for real-time RT-PCR analysis. The expression levels of five genes were only confirmed at one or two heat treatment stages. Two genes of Bra010598 and Bra038432 were up-regualted in "GHA" at the early heat treatment (CK and 1 h) (Figures 5C,D). However, HSP70 (Bra034104), Bra007144 (rpl18) and Bra028816 (RGP, glycoprotein family protein) were highly expressed in "XK" at 12 or 18 h heat treatment ( Figure 5D and Figures 6C,D). These results confirmed the high reliability of the RNA-seq data obtained in the present study. The other four DEGs were further validated for all the heat treatment stages, including NAC005 (Bra010895) and NAC069 (Bra027596), Oxalate oxidase (Bra002316), HSPB27 (Bra030036). As expected, the expression of oxalate oxidase and HSPB27 increased during the heat stress (Figures 7A-C).

Cytochrome P450 Genes Associated with Heat Stress
In higher plants, cytochrome P450s play crucial role in biotic and abiotic stress. Thirteen cytochrome P450 genes were differentially expressed during the heat stress, and most of them were up-regulated in the "GHA" in all the treatments ( Figure 7A, Table S4). Five genes (Bra011759, Bra009148, Bra010598, Bra011821, and Bra034941) were highly induced at the early stages of heat stress (CK, 1 and 6 h), while three genes (Bra009148, Bra010598, Bra011821) were down-regulated in "XK" under 1 and 12 h heat stress, respectively ( Figure 7A). However, two other genes (Bra035844 and Bra001078) were highly expressed at 18 h of heat stress ( Figure 7A). The  cytochrome P450 gene CYP96A2 (Bra033626) was especially lowly expressed in "GHA" under CK and 12 h heat stress ( Figure 7A, Table S4).

Genes Related to Photosynthesis
Photosynthesis is particularly sensitive to heat stress and photosystem II (PSII) with its oxygen-evolving complex (OEC) is the stress-sensitive sites (Allakhverdiev et al., 2008). Two PSII subunit PSBTN (Bra023909) and PSBX (Bra035533) were highly expressed in the heat-tolerant variety "XK" (Figure 7A, Table S4). Moreover, one photosystem I subunit D-1 (Bra036240) was induced in "XK." Another gene (Bra034200) encoding photosynthetic electron transfer Cytb6/f was also slightly up-regulated in "XK" (Figure 7A, Table S4). The results showed that heat stress affect the photosynthesis in non-heading Chinese cabbage.

Differentially Expression of Transcription Factors (TF) during the Heat Stress
HSPs are known to contribute to thermotolerance. We also found two HSPs, Bra034104 (HSP70-1), and Bra030036 (HSP), that were highly expressed in the heat-tolerant variety "XK" during the initial period of heat stress (Figures 7A,C). Some members of the WRKY gene family were reported to be affect by high temperature (Qiu et al., 2004). In this study, one WRKY gene (Bra015372, WRKY7) was highly expressed in "XK" at 1 h of heat treatment. However, two other WRKY genes (Bra017561, WRKY8, and Bra006178, WRKY75) were up-regulated in "GHA" for nearly all the treatments ( Figure 7A, Table S2). The asterisks indicate significant differences between "GHA" and "XK," as determined by Student's t-test (**P < 0.01). Bra028816, RGP (hydroxyproline-rich glycoprotein family protein). The asterisks indicate significant differences between "GHA" and "XK," as determined by Student's t-test (**P < 0.01). NAC proteins are plant-specific TFs which have been shown to be related to plant development and to abiotic and/or biotic stress responses. In this study, one NAC TF (Bra010895, NAC05) was highly induced in "GHA" at the 6 and 18 h of heat treatments, while the other four NAC TFs were up-regulated in "XK" in all the treatments (Figures 7A,D, Table S4). Furthermore, the NAC069 gene was highly expressed in the "XK" at all the stages of heat treatment ( Figure 7E).

DISCUSSION
High temperature influences plant development and can reduce crop yield. To mitigate the effects of heat stress, it is critical to contrive plants that can withstand environmental challenges. In the present study, we investigated the genes that were responsive to high temperature at different stages for the leaves of non-heading Chinese cabbage using RNA-seq. Compared with the transcriptome of a heat-sensitive variety "GHA" and a heat tolerant variety "XK, " potential candidate genes involved in the high temperature response were identified in non-heading Chinese cabbage.
The differential expression analysis of RNA-seq data in "GHA" and "XK" revealed that more DEGs were identified in "GHA" than that of "XK"at 1 h of heat treatment ( Figure 1B, Table 2, and Table S2), indicating that "GHA" is more sensitive to heat stress at the early time. However, at the 12 h of heat treatment, many DEGs were up-regulated or specifically expressed in "XK" ( Table 2). That may be because many genes involved in heattolerance were expressed at this stage in "XK, " such as immune response genes and so on.
Cluster analysis showed that many DEGs were up-regulated in heat tolerant variety "XK" and most of them were involved in "oxidation reduction" or "regulation" process (Figure 2). For example, these genes included oxalate oxidase (Bra002316), HSPB27 (Bra030036), cytochrome P450 (Bra035844, Bra001078, and Bra000845) and NAC069 (Bra027596). These results suggested that oxidative damage was induced during the heat stress and more genes were highly expressed for heat tolerance in non-heading Chinese cabbage.
GO classification analysis showed that the majority of DEGs between "GHA" and "XK" were significantly overrepresented in "response to stimulus, " "secondary metabolic process, " "metabolic process, " "immune response, " "death, " and "ribosome biogenesis" (Figure 4, Table S3). For the DEGs in the four different heat treatment stages, it suggested that the process of responsive to heat stress were divided into three steps. Firstly, GO enrichment analysis revealed that most of the DEGs at the early heat stress stages (1 h) were involved in "response to stimulus" (Figure 5B). These results showed that a number of genes were differentially expressed during the early stage for immediately responding to heat stress. However, at the next time (after 12 h), many genes associated with programmed cell death (PCD) and other immune responsive genes were increased. That may be because the heat stress lead to increased cell death, starting with the PCD and immune system (Figure 6A). At the third step, most of the cell need to be repaired, and many genes involved in ribosome biogenesis were highly expressed at 18 h of heat treatment ( Figure 6B).
Cytochromes P450 (P450) is a family of heme-thiolate enzymes involved in the oxidative metabolism of a variety of endogenous and exogenous lipophilic compounds. In this study, many Cyp genes were up-regulated in the heat-sensitive "GHA, " which indicated that much reactive oxygen species (ROS) were generated during the response to heat stress at the early stage (CK, 1 and 6 h; Figure 7A, Table S4). ROS can directly modify cellular macromolecules, leading to a variety of toxic effects, including lipid peroxidation, protein dysfunction, oxidative stress and cell death. GO enrichment analysis also confirmed that a number of DEGs involved in PCD were expressed the 12 h of heat treatment ( Figure 6A). DEGs associated with ribosome biogenesis were also found to be expressed at the 18 h of heat treatment ( Figure 6B). In rice, transcriptome changes of high night temperature at the early milky stage reported that cytochrome P450 unigene was highly up-regulated in a heat-sensitive rice line (Liao et al., 2015). The results indicated that a higher expression level of the P450 genes would be an index of damage for heat stress.
As signaling molecules, ROS plays an important role in the regulation of abiotic stimuli in plants. ROS include superoxide and hyperoxide, which both contribute to the oxidant state, leading to PCD (Yoshinaga et al., 2005). Oxalate oxidase is involved in the burst of H 2 O 2 which can induce PCD in the aleurone layer during seedling development (Fath et al., 2000;Lane, 2002). In this study, Mn superoxide dismutase (SOD) (Bra034154) and oxalate oxidase (Bra002316) were upregulated in "XK" after heat treatment (Figure 7, Table S2). These accumulations of ROS responsive to the heat stress induced PCD associated genes at 12 h of heat treatment ( Figure 6A).
Heat stress results in the misfolding of newly synthesized proteins and the denaturation of existing proteins. HSP are responsible for protein folding, and can assist in protein refolding under stress conditions (Wang et al., 2004;Song et al., 2014a). In rice, the interplay between HSP101 and HSA32 can affect thermotolerance in seedlings (Lin et al., 2014). In Arabidopsis and carnation, the expression of heat shock transcription factors (Hsfs) HsfA3 in response to heat stress, has been shown to be dependent on the DREB2A transcription factor (dehydration-responsive element binding protein 2A; Schramm et al., 2008;Wan et al., 2015). In Chinese cabbage NHCC, Hsfs was also showed high connection with DREB2 (Song et al., 2016). In the present study, many HSP genes including HSP70, HSP101 and HSPB27 were induced during the heat stress ( Figure 7C, Table S2). Other studies demonstrated that HSP70s play roles in response to heat stress and hydrogen peroxide can enhance ABA-dependent expression of HSP70 to tolerate heat stress (Li et al., 2014;Zhang et al., 2015;Suzuki et al., 2016). Therefore, the up-regulated HSP genes play important role in heat tolerance in non-heading Chinese cabbage.
Many transcription factors are reported to be involved in multiple abiotic stresses including high temperature and cold stress (Bhardwaj et al., 2015;Yang et al., 2015). In this study, we also identified three differentially expressed transcription factors in response to heat stress, including HSF, NAC, and WRKY. The WRKY gene family has been reported to play an important role in plant stress responses. In A. thaliana, WRKY25, WRKY26, and WRKY33, were involved in the regulation of resistance to heat stress (Li et al., 2009;Chen et al., 2012). In our study, two WRKY genes were up-regulated in "GHA" for nearly all the treatments and one WRKY gene was highly expressed in "XK" at 1 h of heat treatment ( Figure 7A, Table S4). WRKY33 was also reported to be up-regulated in Chinese cabbage Chiifu and may play a role in defense signaling . These results suggested that WRKY participated in the response to heat stress and played different roles in the damage patterns between the heat-tolerant and heat-sensitive varieties. Previous studies have reported the enhancement of heat stress tolerance in rice by genetic engineering, such as the overexpression of the Arabidopsis molecular chaperone HSP101 or the rice transcription factor OsWRKY11 in transgenic rice. Overexpression of these genes decreased the necrosis of leaves after heat stress treatment in rice plants (Katiyar-Agarwal et al., 2003;Wu et al., 2009).
Most of the NAC transcription factors were also higher expressed in "XK" than in "GHA" (Figure 7A, Table S4). Previous studies suggested that many NAC TFs had been used to improve stress tolerance in crop plants by genetic engineering (Shao et al., 2015). In Arabidopsis, a NAC transcription factor, NAC078 was reported to be responsive to heat stress (Morishita et al., 2009). Over-expression of NAC2 was also reported to increase abiotic stress resistance in rice (Hu et al., 2008). In this study, many NAC TFs were highly expressed in the heat-tolerant variety "XK, " especially, NAC005 and NAC069 (Figures 7A,D,E). Up-regulation of the NAC transcription factors in "XK" suggested that NACs may participate in heat tolerance in Chinese cabbage. In Arabidopsis, NAC069 (AT4G01550) had been reported to be up-regulated after osmotic stress and 50 µM ABA treatment and in NaCl treated roots (Jiang and Deyholos, 2006;Mishra et al., 2014). In Chinese cabbage, Bra019599 (NAC062/NTL6) was greatly increased upon HS in two inbred lines, Chiifu and Kenshin, suggesting their functions in both cold and HT stress . Therefore, further transgenic research of NAC069 TF in the non-heading Chinese cabbage will elucidate its function in heat tolerance.

CONCLUSION
In this study, we compared the transcriptome of heat-sensitive and heat-tolerant non-heading Chinese cabbage varieties in response to heat stress using RNA-seq. Approximately 625 genes were differentially expressed between the two varieties "GHA" and "XK." The responsive genes can be divided into three phases according to increasing time periods of heat treatment: response to stimulus, programmed cell death, and ribosome biogenesis. Futhermore, NAC069 was up-regulated in heattolerant "XK" at all the heat treatment stages, indicating its roles in heat tolerance. The candidate genes will provide genetic resources for improving the heat-tolerance characteristics in non-heading Chinese cabbage. These results could be also used for future studies of the molecular mechanism of heat stress in plants.

AUTHOR CONTRIBUTIONS
Designed the experiments: GZ, AW, and JH. Performed the experiments: AW, JH, XH, and XL, Analyzed the data: JH, XL, and ZY. Wrote the paper: AW, JH, and GZ. All authors read and approved the final manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016. 00939 Figure S1 | All the expression patterns of the 625 DEGs between "GHA" and "XK" at different stages of heat treatment by clustering analysis in non-heading Chinese cabbage.
Table S1 | Primers used in the study for quantitative RT-PCR. Table S2 | Differentially expressed genes between "GHA" and "XK" at different heat treatment stages.