Transcriptome Sequencing of Two Kentucky Bluegrass ( Poa pratensis L . ) Genotypes in Response to Heat Stress

Kentucky bluegrass (Poa pratensis L.) (KBG) is a major cool-season turfgrass. However, as its complex genetic background and production modes, limited genomic and transcriptomic information of KBG was known so far. In this study, a transcriptome study between wild type material Ninglan (summer stress sensitive) and cultivar material KBG03 (summer stress tolerant) was conducted, under optimal (25 °C) and high (40 °C) temperatures. A total of 81.42 Gb clean reads were generated and de novo assembled into 110,784 unigenes with an average length of 1,105 bp. About 50% KBG unigenes were similar to the non-redundant (NR) database. Among the NR BLASTx top hits, 27.47% unigenes were matched to Brachypodium distachyon. After heat stress, a massive amount of unigenes showed significantly differential expression in both genotypes. After 2h heat stress, more specially up-regulated differentially expressed unigenes (DEGs) and less down-regulated DEGs were detected in Ninglan than in KBG03. At 24h, the expression of 1710 and 730 unigenes were increased and decreased uniquely in Ninglan, and 1361 up-regulated DEGs and 757 down-regulated DEGs were just found in KBG03. More heat shock proteins (HSPs) and heat shock transcription factors (HSFs) DEGs were also identified at 2h than 24h in both genotypes. In addition, by Gene Ontology (GO) enrichment analysis, three core terms (“protein folding”, “response to heat”, and “response to hydrogen peroxide”) of biological process (BP) ontology were found in both genotypes under different heat stress time. The DEGs shared in both genotypes might be related to the basic mechanism of thermal response in KBG.


Introduction
Kentucky bluegrass (Poa pratensis L.) (KBG) is a perennial cool-season turfgrass that is extensively distributed in temperate regions of the world.Due to its softness, medium-to fine-leaf texture, high shoot density, good spring green-up rate and outstanding recuperative capacity (Huff, 2003), it is widely used on golf courses, athletic fields, road sides and home lawns (Puyang et al., 2015).Within KBG, chromosome numbers range from 24 to 124, and the base chromosome number is 7.A complex series of polyploidy and aneuploidy often exists among different strains (Huff, 2003).Their seeds could be produced by both asexually (apomixis) and sexually.Furthermore, the levels of apomixis among different genotypes vary considerably (Meyer, 1982).
As global warming intensifies, the threat on plants by high temperature has become more and more serious.It is considered as heat stress that temperature rises beyond a threshold (10-15 °C above ambient) and it may cause irreversible damage to plant growth and development (Wahid et al., 2007).The extent to heat damage on plants is a complex issue, depending on many factors (intensity, duration and rate of increase in temperature, when and where the high temperature occurs etc.) (Sung et al., 2003;Weis and Berry, 1988).As a complex trait, heat tolerance of plants is controlled by multiple genes.Comparing with warm-season grass, high temperature is more detrimental to cool-season grass.Heat injuries would eventually lead to growth inhibition, ion flux reduction, toxic compounds and reactive oxygen species (ROS) accumulation (Howarth, 2005).
Although Illumina sequencing technology was widely used to detect the DEGs related to heat stress in various plants, but it was limited applied in turfgrass.Besides RNAseq, suppression subtractive hybridization (SSH) and gene chips were adopted to analyze the difference of transcriptome.For example, between two fescue (Festuca sp.) genotypes with different heat tolerance, the transcripts of the heat tolerance-related genes were identified by SSH (Zhang et al., 2005).Using gene chips, the DEGs in response to long-term heat stress treatment were analyzed in switchgrass, and some core genes were selected to identify the heat sensitive plants (Li et al., 2013).Using the Brachypodium genome array, Priest et al. (2014) presented an overall differential expression analysis of Brachypodium transcriptome after chilling, high-salinity, drought, and heat stress.
As a non-model species, we know little genomic and transcriptomic information about KBG, so RNA-seq is a valuable tool for studying the response genes under heat stress.Several transcripome sequencing researches on KBG were reported in recent years, which were focus on the cold stress, salt stress, dwarf and wax deposition (Gan et al., 2016;Ni et al., 2016;Bushman et al., 2016;Zhang et al., 2016).This is the first report of transcriptome profiling of KBG in response to heat stress.In the present study, we analyzed the transcriptome of two KBG genotypes under normal and heat stress conditions using RNA-seq.The objective of this study was to analyze the differences between diverse genotypes on a transcriptome scale, and to reveal the DEGs responded to heat stress.It could provide us some useful information to understand the genes and genetic pathways relating to heat response in KBG.

Plant materials and treatments
Two KBG materials (KBG03 and Ninglan) were used in the study.KBG03 is a new breeding line, with the characteristics of narrow leaves and light green color, which was bred by School of Agriculture and Biology, Shanghai Jiaotong University, Shanghai, China.As a wild type material, Ninglan was collected from Ningxia Province of China, which was provided by Ningxia University.Contrasting to KBG03, Ninglan exhibits deep green and wide leaves.Through field tests for several years, KBG03 appeared obvious summer stress tolerance in Shanghai and surrounding areas in summer, but Ninglan did not (data not shown).The two genotype materials were preserved and reproduced by tillering in the field.
Eight individual plants (clones) of each material were transplanted from field to plastic pots (10×10×10cm 3 ) filled with horticultural medium [peat:vermiculite = 7:3 As an effective research method, RNA sequencing (RNA-seq) was widely applied in turfgrass and forage grass transcriptomic studies.By Illumina high-throughput sequencing, three Brachypodium sylvaticum populations from Spain, Greece and Corvallis, USA were characterized (Fox et al., 2013).The first transcriptome analysis of two buffalograss (Bouteloua dactyloides) cultivars (Wachholtz et al., 2013) and four different guinea grass (Panicum maximum Jacq.) genotypes (Toledo-Silva et al., 2013) were also performed.De novo transcriptome assemblies from multiple tissues of two perennial ryegrass (Lolium perenne L.) genotypes (Farrell et al., 2014), four Lolium-Festuca complex species (Czaban et al., 2015), two Hemarthria materials (Huang et al., 2016) were subsequently reported.Furthermore, many de novo transcriptome analyses were also used to identify the differentially expressed genes (DEGs) responsible for certain specific trait within contrary types of materials.In switchgrass (Panicum virgatum L.), the buds transcriptome of two types of tillering mutants (hightillering and low-tillering) was analyzed (Xu et al., 2015).The mining of genes relating to anthocyanin biosynthesis associated with tissue-specific pigmentation was performed in two zoysiagrass (Zoysia japonica Steud.)cultivars with purple and green spikes and stolons, respectively (Ahn et al., 2015).Dinkins et al. (2017) constructed a tall fescue (Lolium arundinaceum) transcriptome and compared gene expression profiles for endophyte-symbiotic vs. endophytefree clones in different tissues.In KBG, two dwarf mutants were analyzed for the DEGs compared with the wild type by the comprehensive transcriptomic analysis (Gan et al., 2016), and the genes related to cuticular wax deposition between leaf non-elongation zone and the emerged blade zone were also identified (Ni et al., 2016).
In the field of research on biotic and abiotic stress, more and more applications of whole RNA-seq were reported in turfgrass.For instance, perennial ryegrass transcriptome changes between resistance and susceptible genotypes with three different treatments were analyzed during early infection by snow mould (Microdochium nivale) (Kovi et al., 2016).In order to identify defense-related genes, Amaradasa and Amundsen (2016) compared transcriptome of two leaf spot susceptible and two resistant buffalograss lines with inoculated and uninoculated of Curvularia inaequalis.de novo assembly of Sclerotinia homoeocarpa L. and creeping bentgrass (Agrostis stolonifera L.) transcriptomes were done by 454 sequencing technology.Several DEGs were identified from either the host or the pathogen during the infection (Orshinsky et al., 2012).By de novo sequencing assembly of barnyardgrass (Echinochloa crus-galli), potential herbicide resistance genes from susceptible and resistant biotypes with and without herbicides treatment were identified (Yang et al., 2013).The similar research on herbicide (paraquat) resistance was performed in goosegrass (Eleusine indica L.) (An et al., 2014).A KBG transcriptome study between salt tolerant and susceptible accessions was conducted under control and salt treatments (Bushman et al., 2016).In order to improve our current understanding of the cold/freezing response and obtain the potential resistance genes, comparisons of transcriptome profiles (v:v)], and maintained in greenhouse under natural conditions for about two weeks.Then these clones were transferred to illumination incubator (HGZ-300, Shanghai Huitai Equipment manufacturing Co., Ltd.) with constant 25 °C and full light (about 400 µmol m -2 s -1 ), and irrigated every other day.After one week, the temperature was increased rapidly to 40 °C (heat stress temperature) in 15 minutes.Equal leaf samples of each pot were harvested, and they were mixed by genotypes.In all, six sequencing samples were obtained from the two genotypes at 0h, 2h and 24h heat stress treatment points, counted as 03-0h (KBG03 at 0h), 03-2h (KBG03 at 2h), 03-24h (KBG03 at 24h), NL-0h (NingLan at 0h), NL-2h (NingLan at 2h) and NL-24h (NingLan at 24h).These samples were immediately frozen in liquid nitrogen for RNA extraction.Two untreated samples (03-0h and NL-0h) were used as controls.

RNA isolation and RNA sequencing
The total RNA of each sample was isolated using the RNeasy Plant Mini Kit (Qiagen, USA) according to the manufacturer's instructions.The quality of each RNA samples was evaluated using electrophoresis in 1% agarose gels, and quantified by a Nanodrop 2000 machine (Thermo Scientific, Delaware, USA).The six cDNA libraries' preparation and Illumina HiSeq TM 2500 paired-end sequencing were conducted by OE Biotech Co. Ltd. (Shanghai, China).The raw sequence data are available at NCBI under Bio-project number PRJNA383017 and SRA accession SRP104019.

Sequence assembly and annotation
The raw cDNA reads from RNA-Seq were filtered to obtain clean reads for data analysis.During pretreatment trimming, the raw reads with ambiguous 'N' nucleotides and low-quality sequences (quality threshold of 20, length threshold of 35bp) were discarded.Then the clean reads from all libraries were further de novo assembled using Trinity software (version: trinityrnaseq_r20131110) (Grabherr et al., 2011).The assembled unigenes were aligned against public databases, including NCBI nonredundant (NR), SwissProt, Gene Ontology (GO), clusters of orthologous groups for eukaryotic complete genomes (KOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG), with a significance threshold of e < 1e-5.

Identification of differentially expressed unigenes
We calculated the expression levels of the unigenes by reads per kilobase of exon model per million mapped reads (FPKM) method (Trapnell et al., 2010).Hierarchical clustering analysis for common unigenes from the six samples was constructed using the software Cluster 3.0.A fold-change > 2 and false discovery rate (FDR) < 0.05 were used as standards to determine the significant differences in unigene expression between any two samples.DEGs were subjected to GO enrichment analysis, using the hypergeometric distribution test.
Validation of RNA-Seq by Quantitative real time PCR (qRT-PCR) Ten DEGs were chosen for validation using qRT-PCR.Total RNA (1ug) of each sample was reverse-transcribed to cDNA using SuperScript II reverse transcriptase (Invitrogen, China).Primers (Table S1) were designed using NCBI Primer-BLAST Tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi) and synthesized by Sangon Biological Engineering Technology and Service Co. Ltd, Shanghai.The P. pratensis β-tubulin gene was used as an internal reference gene (Albertini et al., 2004).qRT-PCR was performed on the CFX96 Real-Time PCR Detection System (Bio-Rad) using the SYBR Premix Ex TaqII (Takara, China).The thermal cycling conditions were as follows: 95 °C for 3 min, followed by 40 cycles of 95 °C 15s, 58 °C 30s, and 72 °C 20s.Relative expression levels of the unigenes were calculated by the 2 -ΔΔCt method.All amplifications were performed with three replicates and the data were analyzed using Bio-Rad software.

Results and Discussion
Sequencing and assembling A total of 675,628,156 raw reads (82.47 Gb) were produced from six samples.After quality assessment and data filtering, 667,118,738 clean reads (81.42 Gb) were selected and de novo assembled into 110,784 unigenes (≥ 300 bp) for further analysis (Table 1).Among them, 72,720 unigenes (66%) were ≥ 500 bp and 39433 (36%) were ≥ 1000 bp.The N50 of these unigenes was 1694 bp, the average length 1105.18bp, and the maximum length was 17,946 bp.Meanwhile, we individually assembled the data according to the genotypes.Then we obtained 333,570,076 and 333,548,662 clean reads from KBG03 and Ninglan, respectively.In KBG03, 127,500 unigenes (≥ 300 bp) were assembled, with an average length of 889.15 bp and an N50 of 1302 bp (Table 1).In Ninglan, 162,678 unigenes (≥ 300 bp) were assembled, with an average length of 815.55 bp and an N50 of 1101 bp.Evidently more reads were assembled, less and longer unigenes were obtained.In addition, the result showed more number and total length of unigenes was detected in Ninglan than KBG03.

Functional annotation and classification of assembled unigenes
The assembled unigenes from all samples were queried against five public databases (Table 2).More than half unigenes (55,781 unigenes, 50.35%) were similar to proteins in the NR database.A total of 35,425 unigenes (31.98%), 30,108 unigenes (27.18) and 29,210 unigenes (26.37%) were similar to SwissProt, KOG, and GO databases, respectively.Only 8.31% of the total unigenes were annotated in KEGG database.

Analysis of DEGs between different heat stress time and genotypes
Based on the transcript abundance of all common unigenes from six samples, linkage hierarchical clustering analysis was conducted (Fig. 3).Depending on the diversity of the expression profiles from different genotypes and heat stress time, we identified the relationship of the six samples at transcriptional level.Cluster analysis showed that the expression profiles of the same stress time were grouped together and these of 03-0h and NL-0h samples had the closest relationship.The result suggested the transcriptional differences were mainly attributed to heat stress rather than genetic background.Heat map further showed that the expression profiles of 0h and 24h samples were more similar than these of 0h and 2h samples.It suggested that the expression profiles of short heat treatment had larger variety than long time treatment.
Further analysis of DEGs showed it was far more upregulated unigenes than down-regulated unigenes under heat stress (Fig. 4).After 2h heat stress, more special upregulated differentially expressed unigenes (DEGs) and less down-regulated DEGs were detected in Ninglan (up: 2446; down: 379) than in KBG03 (up: 1514; down: 410).Meanwhile, there were 1364 up-regulated DEGs and 66 down-regulated DEGs were shared between them.At 24h, the expression of 1710 and 730 unigenes were increased and decreased specially in Ninglan, and 1361 up-regulated DEGs and 757 down-regulated DEGs were uniquely detected in KBG03.Another 696 and 341 DEGs were shared between the two genotypes with increased and decreased expression, respectively.Overall, more DEGs and up-regulated DEGs were identified at 2h than 24h in KBG.It was because the expression levels of some high-level expression unigenes (at 2h) were recovered to the normal levels (at 24h), such as heat shock protein (HSPs) unigenes and heat shock transcription factors (HSFs) unigenes.During the whole heat stress process, the down-regulated unigene numbers of the two genotypes were always similar, while more up-regulated unigenes were detected in Ninglan than KBG03 significantly.The Venn diagram (Fig. 4B) also showed 1433 and 1042 shared DEGs were detected in both materials at 2h and 24h, respectively.The common upregulated unigenes at 2h (1364) was nearly twice the number of 24h (696), while the common down-regulated unigenes at 24h (341) were much more than those at 2h (66).
In the study, along with the heat stress proceeding, less DEGs and up-regulated DEGs were detected in both genotypes.However, more down-regulated DEGs were found at 24h.Similar result was reported in carnation (Dianthus caryophyllus L.) in response to heat stress, the number of DEGs reached a peak at 2h and reduced at 12h.However, along with the stress increased, more decreased transcript abundance genes were identified only in heattolerant tall fescue materials (Hu et al., 2014).The difference might be result from different resistance of the materials (and may also involve species, heat stress time, developmental stages, etc.).

GO enrichment of DEGs
GO enrichment analysis of DEGs was conducted in different compared combinations, and the terms of BP ontology were analyzed in the subsequent analyses.As a whole, more enriched DEGs were detected in Ninglan than KBG03 and more DEGs enriched at 2h than those at 24h.The result was similar with DEGs analysis in above.
Only two unigenes in GO: 0006457 appeared diverse trend of regulation: 1) comp148926_c0_seq2_1 had no significant expression change at 2h.Nevertheless, it was down-regulated significantly at 24h (< 5% in KBG03; < 8% in Ninglan).2) In KBG03, the expression of comp146338_c0_seq1_1 raised 5764 folds and 465 folds at 2h and 24h, respectively.But in Ninglan, it was downregulated to 2% at 2h, and up-regulated 75 times at 24h.Further analysis showed ten terms were the same in both materials at 2h (03-0h vs. 03-2h; NL-0h vs. NL-2h) and nine were the same at 24h (03-0h vs. 03-24h; NL-0h vs. NL-24h).It suggested these GO terms might be related to the basic thermal response mechanism in KBG.
Secondly, we analyzed the GO enrichment between the two genotypes at 0h, 2h and 24h heat stress points (Table S3).Among the three compared groups (03-0h vs. NL-0h; 03-2h vs. NL-2h; 03-24h vs. NL-24h), there were twelve common GO terms: "DNA integration", "DNA recombination", "RNA-dependent DNA replication", "DNA biosynthetic process", "mitochondrial genome maintenance", "telomere maintenance", etc. Fig. 3. Hierarchical clustering analysis of all common unigenes from six samples.Each column represents a different sample: 03-0h, 03-2h and 03-24h represent KBG03 under heat stress for 0h, 2h and 24h, respectively; NL-0h, NL-2h and NL-24h represent Ninglan under heat stress for 0h, 2h and 24h, respectively.Red: up-regulation; green: down-regulation; black: no change responded to heat stress in heat-resistant material.Also in tall fescue, they found the heat-tolerant sample had more DEGs relative to heat-sensitive one (Hu et al., 2014).These response genes might participate in heat stress regulatory mechanisms in heat-tolerant grass.In this study, more DEGs and HSPs were abundant in the wild genotype material Ninglan than self-breeding line KBG03.So we speculated that Ninglan accumulated more heat response genes through a long-term acclimation process, and it might be more resistant to heat stress compared with KBG03.However, as mentioned in methods, Ninglan exhibited weaker summer stress tolerance than KBG03, so high temperature might not be the main environmental stress factor in the summer for these years in Shanghai and surrounding areas.Further single-factor experiments (including thermal stress, submergence stress, drought stress, UV-B stress, etc.) would be conducted to verify the tolerance of the two materials.
HSP are generally grouped into five classes named for their molecular masses (Vierling 1991): HSP100, HSP90, HSP70, HSP60, and small HSPs (sHSPs).A total of 93 HSPs DEGs (no HSP100 and HSP60 classes) were detected under heat stress treatments (03-0h vs. 03-2h; 03-0h vs. 03-24h; NL-0h vs. NL-2h; NL-0h vs. NL-24h) (Table S4).In our study, all the HSPs were grouped into three families (Howarth 1991):42 unigenes were low molecular weight (LMW: 15-30 k Da) HSPs; 33 unigenes HSPs in DEGs HSPs play a central role in response to high temperature stress and thermotolerance acquisition in plants (Sun et al., 2014).HSPs would be synthesized upon multiple unfavorable environments, including heat, cold, salinity, UV-B light and anoxic stress (Driedonks et al., 2015).According to the annotation in NR and SwissProt databases, we analyzed the HSPs related unigenes in DEGs (Table 3).In KBG03, a total of 78 and 46 HSPs were separately found at 2h and 24h.Among these, there were 11 (2h) and 46 (24h) unigenes had no expression before heat stress.In Ninglan, 99 and 62 HSPs were separately detected at 2h and 24h, in which 16 (2h) and 8 (24h) unigenes had no expression at 0h.Overall, more HSPs were induced at 2h than 24h, and more HSPs were detected in Ninglan than KBG03.The data demonstrated, in response to heat stress, more HSPs showed significant differential expression in Ninglan than KBG03, and a large portion of them presented a strong response under short-term high temperature stress.
As molecular chaperones, HSPs could prevent protein aggregation and assist refolding of non-native protein by heat shock treatment (Xu et al., 2011).Plants could acquire thermotolerance by synthesizing a set of HSPs under high temperature (Howarth 1991).In the orchardgrass (Dactylis glomerata L.) (Huang et al., 2015), more unigenes were Among all types of TFs (Table 4), HSFs were the most (38) (03-0h vs. 03-2h: 24; 03-0h vs. 03-24h: 10; NL-0h vs. NL-2h: 31; NL-0h vs. NL-24h: 16) and they were all upregulated by high temperature treatment.HSFs were detected more largely at 2h than those was detected at 24h and more in Ninglan than in KBG03.In KBG03, 17 significantly up-regulated HSFs were only detected at 2h, and 3 were just detected at 24h (Table S4).In Ninglan, 17 differential expression HSFs were only detected at 2h, and 2 were just found at 24h.In both genotypes, more HSFs DEGs were detected at 2h than those at 24h.It indicated that HSFs were early response genes to heat stress.
Base on the presence of conserved DNA-binding domains plus the nearby HR-A/B region amino acid structure, HSFs could be classified into three classes: A, B and C (Nover et al., 2001).In tomato (Solanum lycopersicum), as a critical gene, LeHsfA1a appeared the master regulator function in heat stress response (Hahn et al., 2011), which formed a regulatory network with HsfA2 and HsfB1 (Mishra et al., 2002).In carnation, the expression levels of HsfA1, HsfA2 and HsfB1 were also increased under heat stress (Wan et al., 2015).In Arabidopsis, expression of rice OsHsfA2e enhanced the tolerance to environmental stresses (Yokotani et al., 2008).In Chrysanthemum nankingense (Asteraceae), HsfA2 and HsfA3 played an important role in the heat stress, but HsfA1a and HsfA2 were not detected (Sun et al., 2014).In this study, HsfA2, HsfA3, HsfA4, HsfA6, HsfA9, HsfB1, HsfB2, HsfC1 and HsfC2 were all detected, except for HsfA1 (Table S7).Among them, HsfA2a, HsfA3, HsfA9, HsfB2b and HsfB2c were detected in all compared groups, but HsfA2c and HsfA2d were only detected at 2h heat stress treatment.HsfB2c and HsfA2d hold the most number in all types of HSPs (seven and five were separately detected in Ninglan at 2h).
At the same heat stress time, more HSFs showed significant differential expression between the two genotypes at 24h (9) than at 2h (3), and no differentially expressed HSFs were detected at 0h (Tables 3 and 4).Furthermore, more HSFs had higher expression in Ninglan (8) than those in KBG03 (3) (Table S6).Zhang et al. (2005) found prominently higher expression of TFs were detected in heat tolerant plants, especially in early heat stress condition.As the analysis in HSPs DEGs, we surmised Ninglan maybe has stronger resistance than KBG03 again.Significantly, most HSFs DEGs were only detected at one stress time, except for comp136900_c0_seq9_1.It had no expression in Ninglan all the time.
In response to high temperature, thirteen HSPs DEGs were separately shared between the two genotypes at 2h and 24h, which were more than those detected at optimum temperature (4) (Table 3).Further analyses of these HSPs (Table S6), two DEGs (comp131554_c0_seq5_1, comp63245_c0_seq1_2) were only detected at 0h, and they were separately related to HSP81 and HSP70.The other two heat shock cognate protein 70 (HSC70: comp192598_c0_seq21_2, comp199473_c0_seq1_2) DEGs were detected at both 0h and 2h, and not found at 24h.Their expression was never detected in KBG03 all the time (0h, 2h and 24h).It was reported that many HSP70 not only induced by stress, but also play a vital role in housekeeping activities under normal conditions.Moreover, some HSC70 (HSP70 homologs) are constitutively expressed in eukaryotic cytosol (Xu et al., 2011).We found HSPs DEGs between the two materials under normal temperature were mainly HSP70 and HSC70 indeed.Six HSPs had differential expression between the two genotypes at both heat stress time and five were unique at 2h and 24h respectively.
Under normal condition, the expression of CL26Contig1 was very low in both genotypes.After heat stress, its expression was sharply increased in KBG03 but not in Ninglan.Combined with the future resistance identification, these special HSPs DEGs could be focused for further study.

Transcription factors (TFs) in DEGs
Besides HSPs, TFs are also an abundant type of stressrelated genes (Nover et al., 2001).As the earliest group of responding genes to biotic and abiotic stresses, they play essential roles in signal transduction pathways (Zhang et al., 2005).Some families of TFs (e.g., bHLH, MYB, WRKY, ERF, HSF) are reported to be regulated under different stress responses in plants (Sun et al., 2014;Wan et al., 2015).335  qRT-PCR validation of the DEGs To validate the expression profile data of the RNA-Seq results, ten common DEGs in all compared combinations (03-0h vs. 03-2h; 03-0h vs. 03-24h; NL-0h vs. NL-2h; NL-0h vs. NL-24h) were selected and examined by qRT-PCR (Fig. 5).According to the annotation in SwissProt database, the function of these unigenes were as follow: comp223146_c0_seq6_2 and comp148676_c0_seq4_1, encoding heat stress TFs; CL1177Contig1, comp232094_c0_seq1_2 and comp136236_c2_seq4_1, encoding HSPs; comp214191_c4_seq5_2 and comp230124_c0_seq1_2, encoding chaperone protein; comp150388_c1_seq13_1, encoding zinc finger protein; CL691Contig1, encoding floral homeotic protein APETALA 2; comp152630_c0_seq10_1, encoding peptidyl-prolyl cis-trans isomerase.qRT-PCR result showed almost all unigenes were up-regulated, except for CL691Contig1.In addition, at different heat stress treatment time, the change trends of these unigenes were consistent with RNA-Seq.So it indicated that the RNA-Seq data was accuracy and reliable.

Conclusions
This study represents the first RNA-seq analysis of KBG under heat stress treatments.A total of 81.42 Gb clean reads were de novo assembled into 110784 unigenes, and all unigenes were functionally annotated by NR, SwissProt, GO, KOG, and KEGG databases.Our data showed more DEGs, including HSPs and HSFs were detected under short heat treatment (2h) than long time treatment (24h), but less down-regulated DEGs were detected at 2h.Three core terms of BP ontology were found in both genotypes under different heat stress time.The DEGs shared in both genotypes might be related to the basic mechanism of thermal response in KBG and the unique DEGs in KBG03, or Ninglan, might be relevant to their difference in summer stress tolerance.

Fig. 1 .
Fig. 1.Distribution species of the NR BLASTx top hits

334Fig. 4 .
Fig. 4. Analysis of DEGs (fold change>2, P<0.05) in KBG in response to heat stress.(A) Numbers of the significantly regulated unigenes of KBG03 and NingLan at 2h and 24h heat stress points.(B) Venn diagrams analysis of the DEGs at the two heat stress time points.Green circles denote Ninglan, Red circles denote KBG03

Table 1 .
Summary of the unigenes of KBG species and two genotypes

Table 2 .
Annotation summary of all unigenes in public databases

Table 3 .
Statistics of the DEGs related to HSPs and HSFs DEGs function was predicted by the annotation in NR and Swiss-Prot databases.Total: combined the unigenes annotated in both databases.

Table 4 .
Statistics of several types of TFs in DEGs (by Swiss-Prot database)