Global Gene-Expression Analysis to Identify Differentially Expressed Genes Critical for the Heat Stress Response in Brassica rapa

Genome-wide dissection of the heat stress response (HSR) is necessary to overcome problems in crop production caused by global warming. To identify HSR genes, we profiled gene expression in two Chinese cabbage inbred lines with different thermotolerances, Chiifu and Kenshin. Many genes exhibited >2-fold changes in expression upon exposure to 0.5– 4 h at 45°C (high temperature, HT): 5.2% (2,142 genes) in Chiifu and 3.7% (1,535 genes) in Kenshin. The most enriched GO (Gene Ontology) items included ‘response to heat’, ‘response to reactive oxygen species (ROS)’, ‘response to temperature stimulus’, ‘response to abiotic stimulus’, and ‘MAPKKK cascade’. In both lines, the genes most highly induced by HT encoded small heat shock proteins (Hsps) and heat shock factor (Hsf)-like proteins such as HsfB2A (Bra029292), whereas high-molecular weight Hsps were constitutively expressed. Other upstream HSR components were also up-regulated: ROS-scavenging genes like glutathione peroxidase 2 (BrGPX2, Bra022853), protein kinases, and phosphatases. Among heat stress (HS) marker genes in Arabidopsis, only exportin 1A (XPO1A) (Bra008580, Bra006382) can be applied to B. rapa for basal thermotolerance (BT) and short-term acquired thermotolerance (SAT) gene. CYP707A3 (Bra025083, Bra021965), which is involved in the dehydration response in Arabidopsis, was associated with membrane leakage in both lines following HS. Although many transcription factors (TF) genes, including DREB2A (Bra005852), were involved in HS tolerance in both lines, Bra024224 (MYB41) and Bra021735 (a bZIP/AIR1 [Anthocyanin-Impaired-Response-1]) were specific to Kenshin. Several candidate TFs involved in thermotolerance were confirmed as HSR genes by real-time PCR, and these assignments were further supported by promoter analysis. Although some of our findings are similar to those obtained using other plant species, clear differences in Brassica rapa reveal a distinct HSR in this species. Our data could also provide a springboard for developing molecular markers of HS and for engineering HS tolerant B. rapa.

Introduction responses have been performed in B. rapa in recent years, but these have primarily focused on cold, salt, and drought stresses [24][25][26][27][28][29][30]. Previously, the expression pattern of HSR genes in Chinese cabbage had been analyzed using 24K microarray, which does not include enough genes to cover the entire B. rapa genome [31]. In this study, we used version 3 microarrays (Br135K) to analyze gene expression in the Chiifu and Kenshin lines. We identified differentially expressed genes (DEGs), specifically expressed genes (SEGs), HSR genes, HS marker genes, and membrane leakage-related (MLR) genes, and we discuss these genes in the context of HS.

Plant materials and heat treatment
Seeds of two Chinese cabbage (Brassica rapa ssp. pekinensis) inbred lines, Chiifu and Kenshin, were obtained from the Korea Brassica rapa Genome Resource Bank, and the plants were grown for 4 weeks in a growth chamber at 22°C under a 16 h light/8 h dark photoperiod with a photon flux density of 140 μmol m -2 s -1 . For heat shock treatments, leaf discs (1 cm in diameter) were incubated at 45°C for 0.5, 1, 2, 3, or 4 h by floating on a water bath. For control samples, the leaf discs were incubated at 22°C for 0.5h by floating on distilled water. Then, leaf discs were blotted, frozen immediately in liquid nitrogen, and stored at -70°C.

Electrolyte leakage test
Immediately after treatment, electrolyte leakage from HS-treated and control plants was measured as previously described [32,33] with some modifications [28]. Briefly, five leaf discs were randomly selected from mixed discs after excision from fully expanded leaves of six plants (10 discs/plant and mixed), and placed in a glass tube with 10 ml distilled water. The samples were then incubated on an orbital shaker at 150 rpm for 30 min at room temperature, after which the initial conductivity (I) was measured using a CON110 conductivity meter (Oakton Ins. USA). The leaf discs were then kept in a boiling water bath for 10 min, after which they were cooled to room temperature, and the final conductivity (F) was measured. The relative electrolyte leakage was calculated using the formula I/F ×100.

Construction of Br135K chips
The Br135K microarray (Brapa_V3_microarray, 3 0 -Tiling microarray) is a high-density DNA array prepared by NimbleGen (http://www.nimblegen.com/) using Maskless Array Synthesizer (MAS) technology as described in Jung et al. [28]. Labeling, data processing, and background correction were performed as described previously [28]. To assess the reproducibility of the microarray analysis, we repeated the experiment using independently prepared total RNA from two or three biological replicates. The complete raw microarray data have been deposited in the Omics database of NABIC (http://nabic.rda.go.kr) with accession numbers NC-0023-000001-NC-0023-000024.

Gene chip data analysis
Genes with adjusted P-value or false discovery rate (FDR) below 0.05 were collected and further selected for those genes with expression greater than or less than in at least one treatment. Multivariate statistical tests such as principal component analysis and multidimensional scaling were performed with Acuity 3.1 (Molecular Devices, USA). Clustering analysis was carried out using MultiExperiment Viewer version 4.9 (MeV4.9, http://www.tm4.org/mev.html). To obtain insights regarding the putative biological functions and biochemical pathways of DEGs, we carried out enrichment analyses by searching Gene Ontology (GO) [34], agriGO [35], and the Kyoto Encyclopedia of Genes and Genomes [36].

RNA extraction, RT-PCR, and qRT-PCR analysis
Total RNA was treated with RNase-free DNase I (Promega, USA) to remove genomic DNA contamination. First-strand cDNA synthesis was carried out using the Ace-α kit (Toyobo, Osaka, Japan). RT-PCR was performed using 25 ng of cDNA from plants exposed to HS treatments. The gene-specific primers for the stress-responsive genes are listed in S1 Table. Reactions were performed in 0.2 mL PCR tubes containing 10 pmol of each primer, 150 μM of each dNTP, 1.2 U of Taq polymerase, and 1X Taq polymerase buffer; double-distilled water was added to a total volume of 20 μL. The PCR cycle consisted of pre-denaturation at 94°C for 5 min, followed by 25 cycles of denaturation at 94°C for 30 s, annealing at 58°C for 30 s, and extension at 72°C for 45 s, followed by an additional extension step for 5 min at 72°C. PCR products were electrophoresed on 1.5% agarose gels. For qRT-PCR, pre-denaturation was performed at 95°C for 30 seconds, followed by a 40-cycle three-step reaction (5 seconds at 95°C, 20 seconds at 60°C, and 15 seconds at 72°C) with the primer sets shown in S2 Table. Analysis of enriched conserved cis-elements HSR genes were subjected to the Short Time-series Expression Miner (STEM) [37], and 1,000 bp upstream sequences of DEG start codons were analyzed for the presence of plant cis-acting elements (PLACE) [38]. The occurrences of these motifs were compared with their frequencies among all promoters represented in the 135K microarray. A P-value was then calculated for each motif and profile combination, based on the hypergeometric distribution [39]. Motifs with P-values below 10 -4 were considered to be significantly enriched.

Membrane leakage test
Along with survival percentage, lipid peroxidation, and peroxidase (POD) activity, electrolyte leakage was considered to be a good indicator of thermotolerance in Brassica species [40]. Therefore, we measured leakage during HT exposure of two inbred Chinese cabbage lines, Chiifu and Kenshin. As shown in Fig 1, the two lines exhibited similar electrolyte leakage after 1 h exposure to 45°C, but a significant difference was observed after 2 h. On the basis of this result, we concluded that Kenshin is more resistant to > 2 h exposure to HT at 45°C, exactly the reverse of the two strains' relative tolerance to freezing temperatures [28]. The major parameters used to phenotype thermotolerance of crop plants include viability, pollen development and fertility, photosynthetic rate, and germination [41]. However, breeders and farmers of the leafy vegetable B. rapa pay much more attention to leaf defense against pathogens, growth, and the heading process, complicating the choice of good marker for thermotolerance. Our results indicate that electrolyte leakage could be a useful marker for temperature sensitivity in B. rapa.

Microarray analysis to identify HSR genes
To identify HSR genes, we identified DEGs under HS from among 41,173 genes represented on Br135K microarrays (S3 Table). After removing genes with PI (Probe Intensity) values below 500 for all samples (usually it's hard to detect the gene expression level using RT-PCR under 25 cycles condition), DEGs were defined as genes with > 2-fold changes in at least one time point. A large number of genes, 4,008 for Chiifu and 2,611 for Kenshin, were up-regulated (Fig 2A and 2B, S4 and S5 Tables). Among these, 729 and 654 genes in Chiifu and Kenshin, respectively, were up-regulated > 2-fold during HS treatments of 0.5-4 h, suggesting that about 20% of HS up-regulated genes (18.2% in Chiifu; 25% in Kenshin) were expressed more strongly during HSR (Fig 2A and 2B, S4 and S5 Tables). In addition, many genes were downregulated by HS (4,164 in Chiifu and 3,158 in Kenshin) (Fig 2C and 2D, S6 and S7 Tables). Among these, 1,413 and 881 genes in Chiifu and Kenshin, respectively, maintained > 2-fold down-regulation during HS treatments of 0.5-4 h; thus, HS suppressed about 30% of HSdown-regulated genes at all-time points used (Fig 2C and 2D, S6 and S7 Tables). The proportion of DEGs under HS was similar to that obtained in an analysis of the tomato transcriptome (2,203 genes, 9.6%) [13], but higher than the 2% reported by other studies [42][43][44].
Next, we functionally analyzed HSR genes by GO enrichment using agriGO [35], based on information about Arabidopsis homologs. All Arabidopsis counterparts of our microarray probes were used as background references, and significantly represented GO items were defined as FDR values below 0.05. At a significance of < 0.05 FDR, 63 and 106 GO items were significantly enriched among the up-regulated genes in Chiifu and Kenshin, respectively (S8 Table). On the other hand, 126 and 145 GO items in Chiifu and Kenshin, respectively, were significantly enriched among the down-regulated genes (S9 Table). The most enriched GO items among the up-regulated genes in both genotypes were associated with HS: response to heat (GO:0009408), response to ROS (GO:0000302), response to temperature stimulus (GO:0009266), response to abiotic stimulus (GO:0009628), MAPKKK cascade (GO:0000165), etc. (S8 Table). Some GO terms exhibited genotype-specific up-regulation: for Chiifu, response to water deprivation (GO:0009414), and for Kenshin, calcium-mediated signaling (GO:0019722), protein targeting to membrane (GO:0006612), and respiratory burst (GO:0045730) (S8 Table).
To determine the association of HSR genes with specific pathways, we placed the Arabidopsis counterparts of the identified DEGs on KEGG pathway maps. This analysis identified 126 pathways (S10 Table), including plant hormone signal transduction, starch and sucrose metabolism, plant-pathogen interaction, glutathione metabolism, and so on. Most genes related to the RNA degradation pathway were up-regulated by HS in both Chiifu and Kenshin. Most of genes related to cyanoamino acid metabolism and the N-glycan biosynthesis pathways were down-regulated by HS in both Chiifu and Kenshin. We identified no pathway specific for the HS response in Chiifu or Kenshin.

Intrinsic transcriptome differences between Chiifu and Kenshin prior to HT treatment
Because previous reports suggested that plant HS tolerance was largely due to constitutive expression of many genes prior to stress treatment [45,46], we compared the intrinsic differences between two Chinese cabbage inbred lines. Intrinsic genes were selected by fold change on the basis of PI values. A large number of genes, 522 genes for Chiifu and 205 genes for Kenshin (S11 Table), exhibited a > 4-fold change compared to the other, but the functions of most genes exhibiting high levels of expression could not be associated with HS. Many unknown genes (including orphan genes) were in this category, suggesting that they play roles in the HSR. To analyze the expression patterns of SEGs in each inbred line, we selected some of the top-ranking genes from S11 Table and subjected them to qRT-PCR analysis (Fig 3). None of these genes have yet been functionally characterized, but expression of all of them except Bra026915 was reduced upon HS treatment. To isolate the candidate genes which may work for the heat tolerant in Kenshin, we isolated genes with intrinsically expressed in Kenshin and  showed up-regulation in Keshin after heat treatment at the same time. In the end, 14 genes, including one Hsf gene (Bra029292, HSFB2A), one methionine sulfoxide reductase (Bra019187, PMSR4), three orphan genes (Bra039796, Bra016685, and Bra028606), were identified.
To obtain functional clues about the genes in this category, we carried out GO annotation on the basis of the Arabidopsis homolog data (Table 1). Genotypic SEGs included those encoding kinases, TFs, and factors involved in response to stress, response to abiotic or biotic stimulus, transferase activity, hydrolase activity, and signal transduction. Noteworthy categories include cellular respiration, unsaturated fatty acid biosynthesis, and Hsf genes, which were over-represented in Kenshin relative to Chiifu. For instance, 5 Hsf genes, Bra013409, Bra029292 (HSFB2A), Bra029291 (HSFB2A), Bra004272 (HSFA8), and Bra000557 (HSFA2) were highly intrinsic expressed in Kenshin. Otherwise, only one Hsf gene, Bra012828 (HSFA7A), showed highly intrinsic expressed in Chiifu. We also carried out KEGG pathway analysis to reveal differences between the two genotypes (S12 Table). In total, 100 and 86 Table 1. GO annotation of intrinsically expressed genes. One gene can be associated with two or more GO terms. NA, no Arabidopsis homologs.

Cellular respiration 8 2
Unsaturated fatty acid biosynthetic process 6 2 Heat shock transcription factor 5 1 pathways were represented by the up-regulated genes in Chiifu and Kenshin, respectively. Metabolic pathways were most prominent in both genotypes, suggesting that a metabolomics approach will be required in future studies.

Analysis of some HSR genes
As a first step toward identifying putative regulatory genes for thermotolerance, we identified several genes exhibiting very high levels of expression or specific expression changes, and confirmed them by RT-PCR (Fig 4, S1 Fig, and S13 Table). The top ten genes induced by HS in both genotypes included genes encoding small Hsps, chaperones, kinases, and unknown proteins. Many of these are well-known genes involved in thermotolerance in plants [3,19,47]. Among these, Bra002538 (HSP18.2) exhibited the strongest induction in both genotypes, > 410-fold at 0.5 h and > 630-fold at 1 h, and maintained its level thereafter. An unknown gene (Bra033343, AT1G02700) exhibited a peculiar pattern of expression upon HS treatments: in Chiifu, 1,213-fold induction at 0.5 h and 217-fold induction at 1 h, but in Kenshin, 1,573-fold at 0.5 h, 2,462-fold at 1 h, and 176-fold at 2 h. These patterns are similar to Hsp90 and Hsp70 expression, as expected for the HSR. Several genes were specifically induced in either genotype. Among them, Kenshin SEGs (Bra028923, Bra029764, Bra014525, Bra027662, and Bra034336) could be good targets for thermotolerance studies in Chinese cabbage because Kenshin has been traditionally used as a breeding stock to develop heat-tolerant plants [22,23].

Hsfs and Hsps
Many Hsfs and Hsps play important roles in plant heat tolerance. In the Br135K microarray, 38 Hsf-and 206 Hsp-encoding genes (Hsp100/ClpB, Hsp90, Hsp70, Hsp60, Hsp40, and sHsp) are represented (S14 Table). Among them, 49 genes (including 11 Hsf genes) exhibited very low levels of expression (PI values < 500 of in all samples). Upon 1 h HS treatment, however, 10 and 13 Hsf genes were up-regulated > 2-fold in Chiifu and Kenshin, respectively. Among these, nine Hsf genes were up-regulated in both genotypes:  [13]. Some Chinese cabbage Hsfs (HsfB1, HsfB2A) exhibited expression profiles similar to those observed in tomato, whereas others did not, indicating the presence of distinct regulatory mechanisms in different plant species. The Arabidopsis genome contains 21 Hsf genes that can be categorized into three classes (A, B, and C) [10,48]. HsfA1a is a master regulator of the HSR [10,12,13,49]. However, given the low expression level of HsfA1 following HS, our data imply that other Hsfs could serve as the major TFs in B. rapa. HsfA4a and HsfA8 are sensors of reactive oxygen species (ROS), which are produced as a secondary stress signal during the HSR [50,51]. In particular, HsfA4 acts as a sensor of the H 2 O 2 signal [15], whereas HsfA5 acts as a negative regulator of this pathway [16]. In our microarrays, both genes were highly expressed. Although the expression of HsfA4 increased slightly upon HS, the expression of HsfA5 decreased slightly (S14 Table). All these results suggest that heat tolerance in B. rapa might be primarily regulated at the protein level, with only minor regulation at the transcriptional level. Alternatively, heat tolerance might involve key regulators that are different from those Arabidopsis, as proposed in Wang et al. [52].
Like Hsfs, Hsps also play important roles in heat tolerance in plants, especially in protein folding, assembly, translocation, and degradation [17]. In our microarrays (S14 Table), most Hsp70, Hsp90, and Hsp101 genes were constitutively and highly expressed in all samples. However, small Hsps encoded by genes such as Hsp18.2, Hsp20, Hsp22, and Hsp23.6 exhibited strong induction (15-656-fold) by HS treatment in both genotypes. We confirmed the expression of some genes by RT-PCR (S2 Fig). These results are similar to those in tomato, in which the majority of HS-induced chaperone genes belong to the sHsp family [13].
In Arabidopsis, HS-associated 32 kD protein (HSA32) (AT4G21320), an Hsp, is essential for acquired thermotolerance during long-term recovery after acclimation treatment [53] and in rice [54]. B. rapa does not contain ortholog of this gene, providing further evidence for the existence of different HSRs among plants. Recently, Wang et al. [52] reported that different key genes and regulatory mechanisms are involved in abiotic stress responses in B. rapa and Arabidopsis.

Upstream components of the HSR, including ROS
Upstream events of the HSR include perception and signal transduction of HS. Jia et al. [55] demonstrated that Hsf activity and Hsp production involved in thermotolerance are a result of cross-talk among H 2 O 2 , nitric oxide (NO), Ca 2+ channels, and calmodulin (CaM). In particular, CaM activates Hsfs, protein kinase, phosphatase, and cyclic nucleotide-gated ion channels. Among these, ROS including H 2 O 2 are considered to be the first signaling components produced by HS [15,56,57]. The ROS signal can then be recognized by histidine kinases and Hsfs (HsfA4 and HsfA8) [19,50,51]. Hsf sensing activates multiple TFs (Zat family, WRKY, and MBF1c) through the MAPK signal pathway [19].
Regarding the expression of ROS-scavenging genes upon initiation of the HSR, two possibilities exist: such genes could be induced in order to reduce levels of ROS produced by HS, or they could be suppressed in order to maintain the levels of ROS involved in downstream signaling. Among 170 PODs on the Br135K microarray, only seven were differentially expressed upon HS treatment (S15 Table). Three genes were up-regulated in both genotypes: ascorbate peroxidase 1 (APX1) (Bra031598), Fe-superoxide dismutase 3 (FSD3) (Bra026503), and Fesuperoxide dismutase 2 (FSD2) (Bra029190). Previous work revealed APX1 as a central component of the ROS gene network [58] and a key player in stress responses [59], whereas FSD2 and FSD3 confer protection against oxidative stress [60,61]. Glutathione peroxidase 2 (GPX2) (AT2G31570; Bra022853), which exhibits Kenshin-specific up-regulation, protects against stress through activation of cytosolic superoxide dismutase [61]. Thus, BrGPX2 is a promising candidate thermotolerance-related gene in Kenshin. In addition, one peroxidase superfamily protein (Bra016930), which is also likely to be involved in ROS signaling, was down-regulated upon HS in both genotypes.
Because calcium signaling can engage in cross-talk with ROS signaling to trigger thermotolerance, we investigated the expression profiles of several calcium-related genes represented on the Br135K chip: 24 calcium channels, 39 CaM, 31 calcium-dependent protein kinases (CDPKs), five CaM-binding protein phosphatases, and 33 cyclic nucleotide-gated ion channels (CNGCs) genes (S16 Table). Most calcium channel genes were highly expressed at all-time points, suggesting that they are regulated at the protein level, if at all. In contrast to the nine Arabidopsis CaM genes up-regulated by HS treatment [62], 14 BrCaM genes (with the exception of CaM8) exhibited high constitutive expression under most conditions. These data further support the idea that the HSR is differentially regulated among species. Expression of CDPK and CNGC genes exhibited no significant changes after HS treatment, with the exception of Bra001676, a homolog of Arabidopsis salt-induced CNG20 [63], which had Kenshinspecific expression. Taken together, these data suggest that regulation of calcium-related proteins in B. rapa HSR is exerted at the protein level.
Although most phosphatase genes were highly expressed in all samples (S18 Table), three were differentially expressed: two were up-regulated genes upon HS (AT3G51470; Bra036796, AT4G29690; Bra011129), and one was expressed specifically in Kenshin (probable apyrase 5, Bra019669). The functions of these genes have not yet been determined. As with the kinases, the expression patterns suggest that most phosphatases function at the post-transcriptional level.

Transcription factors (TFs)
TFs regulate many genes involved in plant growth and development at the transcriptional level. The Br135K microarray includes 3,354 TFs (8% of the total of 41,174 genes), of over 43 different classes ( Table 2, S19 Table). Among them, 109 TFs exhibited significant changes in their expression levels following HS treatment ( Table 2). Apart from the heat-induced Hsfs, the majority of these TFs belong to the Integrase-type, NAC, homeodomain, HB, bHLH, and DREB families. Members of these TF families are involved in the HSR and thermotolerance in different plant species [45,64,65], indicating the existence of complex transcriptional regulatory networks beyond the direct regulatory control of Hsfs.
According to a recent study, the ER membrane-associated basic leucine zipper (bZIP) transcription factor bZIP28 is activated by ER stress resulting from the accumulation of misfolded or unfolded proteins following HS; this factor stimulates expression of stress response genes [66]. Two BrbZIP28 genes were highly expressed in all Chinese cabbage samples; expression of one of them (Bra034147) increased slightly upon HS treatment, whereas the other gene (Bra023224) exhibited no significant change.
Several TFs were selected and subjected to qRT-PCR analysis. Three genes (Bra022602, Bra039069, Bra019599) were induced in both lines at similar levels upon HS, and three and two genes exhibited prominent expression in Chiifu (Bra003801, Bra005852, Bra027501) and Kenshin (Bra021735, Bra024224), respectively (Fig 5). Due to the large standard deviation resulting from normalization, the qRT-PCR data of several genes (MYB41, MYB95, MYB82, and bHLH92) were not consistent with the microarray values. Except for Bra022602 (MYB82), expression of most genes was highly increased by HS treatment. With the exception of three genes (CCH-type, MYB82, and MYB95) whose functions have not been identified, these genes are related to abiotic stress responses in Arabidopsis. Two genes with prominent expression in  Chiifu, Bra005852 (DREB2A: dehydration-responsive element binding protein 2A) and Bra027501 (bHLH92), are responsive to cold stress [28] suggesting that they serve another role in HSR. In Arabidopsis, bHLH92 functions in response to NaCl, dehydration, and cold [67], but its homolog, Bra027501, would be expected to function in HS of B. rapa, based on its expression change. Arabidopsis DREB2A controls the levels of HsfA3, Hsp18.1-CI, Hsp26.5-MII, and Hsp70 during heat and drought stress [68][69][70], implying that it responds to both heat and cold stresses. B. rapa contains three DREB2A genes: two of these (Bra005852, Bra019162) were expressed at low basal levels under normal growth conditions and were greatly up-regulated by HS, whereas the other (Bra009112) was expressed at a relatively high level that increased further upon HS treatment (S19 Table). These results indicate that expression of Hsps in B. rapa might be also induced by DREB2A. Two genes exhibiting a similar response to HS treatment in both inbred lines, Bra003801 (MYB95) and Bra019599 (NAC062/NTL6), are predicted to function in both cold and HT, even though only the cold response has been characterized in Arabidopsis [71,72]. In Arabidopsis, NAC062/NTL6 processing is stimulated by cold and ABA to activate the expression of defense response genes, but in B. rapa, its processing is reduced by HT [71]. Given that the expression of Bra019599 was greatly increased upon HS in both lines, it is possible that BrNAC62 functions in both cold and HT stress.
Two genes with a notable response to HS treatment in Kenshin, Bra024224 (MYB41) and Bra021735 (a bZIP/AIR1 [Anthocyanin-Impaired-Response-1]), may function in HS tolerance. Arabidopsis MYB41 activates suberin synthesis under abiotic stress conditions [73], implying a relationship to HS. Arabidopsis AIR1 regulates various steps in the flavonoid and anthocyanin accumulation pathway [74]. Its expression is regulated by salt stress, but not by HT or drought. It is possible that BrAIR1 could function under HS conditions, particularly for Kenshin. Therefore, BrAIR1 could be a useful molecular marker for HT tolerance in B. rapa.
HSR genes and proteins can be divided into two groups, signaling components (protein kinases and TFs) and functional genes (HSPs and catalase) [75,76]. Zat7 [77], Zat10 [15], and Zat12 [78] respond to HS. In particular, Zat12 is necessary for the expression of APX, Zat7, and WARKY 25. B. rapa contains Zat1, 6, 10, and 12, but not Zat7. The expression levels of Zat12 (four alleles) and Zat10 (two alleles) were very high and slightly increased by HS. Also, WRKY25 expression was very high. These results suggest the possible involvement of TFs through the MAPK signaling pathway.

HS marker genes
According to a previous report [41], Arabidopsis plants exhibit four classes of tolerance to HS regimes: basal thermotolerance (BT), short-term acquired thermotolerance (SAT), long-term acquired thermotolerance (LAT), and thermotolerance to moderately high temperature (TMHT). Because genes representing each regime have been identified in Arabidopsis, we examined expression of the B. rapa genes that correspond to the Arabidopsis genes involved in the tolerance phenotypes (S20 Table, Fig 6). Expression of most genes did not differ between Chiifu and Kenshin, but only expression of XPO1A (exportin 1A) (Bra008580, Bra006382), which is involved in BT and TMHT in Arabidopsis, seemed to be correlated with thermotolerance as well as membrane leakage. LAT-related genes could be up-regulated (induced or stimulated) even after 30 min HT treatments, but no difference was observed between the two genotypes. The Hsa32 (HS-associated 32 kD: AT4G21320) homolog, which is also a marker of LAT, was absent from the B. rapa gene database. These results suggest that (1) HSR markers differ between Arabidopsis [41] and B. rapa, and (2) BT and LAT from B. rapa are controlled by different genes.

Membrane leakage-related (MLR) genes
Even though HSR genes were immediately induced by HT exposure, membrane leakage of B. rapa remained constant until 1 h HT exposure (Fig 1). Regarding electrolyte leakage, Kenshin is less sensitive to prolonged HT exposure than Chiifu, implying that the HT responses of these lines are negatively correlated with leakage. Because electrolyte leakage reflects the traits of these two inbred lines, analysis of the MLR genes is important for an understanding of the B. rapa HT-response. MLR genes can be defined as those whose expression patterns correlate with the membrane leakage pattern either positively or negatively, i.e., expression of genes should not exhibit a significant change from no treatment to 1 h HS treatment, but should start to change upon 2 h HS treatment (over 2-fold, up or down) and increase thereafter. As shown in S21 Table, expression levels of 15 and 65 genes were changed over 2-fold in Chiifu and Kenshin, respectively. Among them, four genes specifically expressed in Kenshin were selected as putative MLR genes: two up-regulated genes [Bra032994 (NA), Bra025083 (cytochrome P450; CYP707A3)], and two down-regulated genes [Bra006613 (FAD/NAD(P)-binding oxidoreductase) and Bra016265 (protein kinase superfamily protein)] (Fig 7A and 7B). In particular, the Bra025083 expression pattern was well matched to the membrane leakage pattern. B. rapa includes three CYP707A3 genes, and two of them, Bra025083 and Bra021965, exactly followed the pattern of membrane leakage (Fig 7C). The reported functions of CYP707A3, a major ABA 8'-hydroxylase, includes the dehydration response [79,80], stomata closure [81], and salt stress [80]. Our data suggests another function for this gene in HS conditions.

Orphan genes
Orphan genes, which are protein-coding genes unique to a species, are widespread across all organisms [82]. Estimates of the percentage of orphan genes in various species range from 1% to 71%, with 5-15% being fairly typical [83]. Orphan genes are associated with abiotic stresses, including HS [83][84][85]. About 13% of genes on the Br135K microarray were orphan genes (5,349 out of 41,173 genes) (S22 Table). Among them, 34 were induced or up-regulated by HS treatment, and 16 and 13 genes were induced in a Chiifu-or Kenshin-specific manner (S22 Table, Fig 8). Bra026318 (15-37-fold induction) and Bra024533 (25-34-fold induction) exhibited the strongest induction in Chiifu and Kenshin, respectively. In addition, 54 and 14 genes were specific to Chiifu and Kenshin, respectively, suggesting that they were associated with specific traits of the two genotypes, ncluding their HS responses. Several genes that were up-regulated in both inbred lines were confirmed by qRT-PCR (Fig 9). In contrast to the PI value changes observed in the microarray experiment, the induction revealed by qRT-PCR was less obvious, even though the expression levels clearly increased upon HS. Analysis of conserved cis-elements enriched for HSR genes Because the full genome sequence has already been reported for Chiifu (http://brassicadb.org/ brad/), we performed an analysis of cis-elements of the HSR genes in order to extend our knowledge of HSR mechanisms in B. rapa. Although the genome sequence of Kenshin has not yet been published, results obtained using the Chiifu genome sequence should be sufficient to provide information about Kenshin because expression of most TFs exhibited similar patterns in both inbred lines (Fig 5). We performed our clustering analysis with STEM [37] using HSR genes that were > 2-fold DEGs (8,160 genes) in Chiifu. Fifty profiles were obtained from STEM, and nine of them significant overrepresentation (p < 0.01) (S23 Table, Fig 10). We then subjected the 1,000 bp upstream sequences of the genes in the nine significantly overrepresented profiles to PLACE database for motif scanning [38]. The occurrences of these motifs were compared with their frequencies among all promoters represented in the 135K microarray. A Pvalue was then calculated for each motif and profile combination, based on the hypergeometric distribution [39]. Motifs with P-values below 10 -4 were considered to be significantly enriched significantly enriched motifs (SEMs) (S24 Table). Among the SEMs in profiles P35, P42, P47, and P19, we identified several motifs related to the HSR: CCAATBOX1 (CCAAT) in P35 is required for heat shock promoter activity [86,87]; and ACGTABREMOTIFA2OSEM (ACGTGKC) in P42, ABRELATERD1 (ACGTG) in P47, and two calcium-responsive motifs in P47, ABRERATCAL (MACGYGB) [88], and CGCGBOXAT (VCGCGB) could be responsible for the HS response. These three motifs are also known as ABREs binding by a bZIP TF responding to multiple stresses, such as drought and cold [88][89][90][91][92][93]. Our motif and transcriptome analyses suggest that ABREs are also necessary for the HSR. SEMs found in down-regulated genes were also related to the HSR: the DRECRTCOREAT (RCCGAC) motif, represented in P1 and P3, is an AP2 TF-binding site involved in cold and ABA-responsive processes [94]; and SITEIIATCYTC (TGGGCY), in P1 and P13, is a TCP-domain TF-binding site important for mitochondrial oxidative phosphorylation. Furthermore, we subjected genes containing the same cis-elements to agriGO for GO enrichment analysis. Processes including response to heat, heat acclimation, MAPKKK cascade, calmodulin binding, and regulation of oxygen and ROS metabolic process were represented by up-regulated genes containing shared cis-elements (S24 Table). In addition, shoot system development, respiratory burst during defense response, and  very-long-chain fatty acid metabolic processes were represented by down-regulated genes containing shared cis-elements (S24 Table). All of this information will be useful because combinatorial analysis of promoters, co-expression networks, and transcriptome can provide further insight into gene functions and signaling pathways [39,95].

Conclusion
Genome-wide analysis of the transcriptome following heat treatment in two inbred lines of Chinese cabbage provided important information regarding the HSR. First, even though Kenshin was more resistant to membrane leakage than Chiifu, the number of DEGs was higher in Chiifu. Second, expression of small Hsps was highly induced by HS treatment, whereas highmolecular weight Hsps exhibited constitutively high expression. Third, although expression of several upstream genes of the HSR was induced by HS treatment, most genes associated with the signaling pathway were constitutively expressed, suggesting that they are regulated at the protein level. Fourth, besides well-known TFs, many TFs and orphan genes seemed to be related to HSR in B. rapa. Fifth, most of the B. rapa HSR is likely to use mechanisms identified in other plants, but the critical or core genes could be different. To strengthen HSR mechanism in Chinese cabbage, it will be essential to perform combinatorial analyses of transcriptomes, proteomes, metabolomes, and small RNAs.   Table. List of GO terms significantly enriched among genes down-regulated by HS in Chiifu and Kenshin. Yellow shade indicates the FDR value below 0.05 (significantly represented items). (XLSX) S10 Table. KEGG pathway maps of HSR genes in Chiifu and Kenshin. Gene No. represents Arabidopsis counterparts involved in each pathways. (XLSX) S11 Table. Intrinsic genes expressed in Chiifu or Kenshin under normal growth conditions. Expression levels of all genes in this list changed by over 4-fold in both genotypes (XLSX) S12 Table. KEGG pathway maps of intrinsic genes. Gene No. represents Arabidopsis homologs involved in each pathway. (DOCX) S13  Table. Expression profiles of membrane leakage-related (MLR) genes. MLR genes were defined as those whose expression matters correlated negatively or positively with the membrane leakage pattern, i.e., genes should exhibit no significant change from no treatment to 1 h HS treatment, but should change after 2 h HS treatment (over 2-fold, up or down). (XLSX) S23