Transcriptomic Responses of Dove Tree ( Davidia involucrata Baill.) to Heat Stress at the Seedling Stage

: The dove tree ( Davidia involucrata Baill.), a tertiary relic species, is adapted to cool climates. With the progression of global warming, high-temperature stress has become the primary environmental factor restricting geographic distribution, ex situ conservation, and landscape application for D. involucrata resources. However, the detailed molecular events underlying D. involucrata responses to heat stress are poorly understood. Here, we conducted RNA-Seq-based gene expression proﬁling in D. involucrata seedlings during the time course of a 42 ◦ C heat treatment (0, 1, 6, and 12 h). After de novo assembly, we obtained 138,923 unigenes, of which 69,743 were annotated in public databases. Furthermore, 19,532, 20,497 and 27,716 di ﬀ erentially expressed genes (DEGs) were identiﬁed after 1 h (HS1), 6 h (HS6), and 12 h (HS12) of heat treatment in comparison to 0 h (HS0), respectively. Based on a KEGG enrichment analysis, the two pathways “protein processing in endoplasmic reticulum” and “plant hormone signal transduction” are hypothesized to play vital roles during heat response in D. involucrata , and their potential interactions during heat stress are also discussed. In addition, 32 genes encoding putative heat shock transcription factors (Hsfs) were found to be associated with the response of D. involucrata to heat stress. Finally, the expression patterns of eight heat-responsive genes derived from qRT-PCR were in agreement with their transcript level alterations, as determined by a transcriptome analysis. Taken together, our transcriptomic data provide the ﬁrst comprehensive transcriptional proﬁle a ﬀ ected by heat stress in D. involucrata , which will facilitate further studies on the improvement of heat tolerance in this rare and endangered species. These results suggest that most of the genes engaged in protein processing were induced during heat stress to promote the e ﬃ ciency of protein processing for stress response.


Introduction
Davidia involucrata Baill., the sole species placed in the genus Davidia in the family Davidiaceae, is a well-known ornamental plant called the "dove tree." It is an endangered and rare relict tree species endemic to China with prominent scientific, horticultural, ecological, and medicinal values [1][2][3][4]. Owing to the highly strict ecotope demand, the geographic distribution of D. involucrata populations is rare and scattered at high altitudes in southwestern and south-central China. D. involucrata is best adapted to cool climates, and high temperature is generally detrimental to its growth and development. Notably, it has been estimated that this species will be greatly vulnerable under future rapid climate warming [5]. For Davidia resources conservation, research has been undertaken on its introduction and ex situ cultivation since 1979 in China [6]. However, such research is not going smoothly, since high-temperature stress in summer is the major limiting factor for the efficient propagation of D. involucrata outside its natural habitats [7]. Previous studies of D. involucrata in relation to heat stress have mainly focused on physiological or biochemical alterations [8,9]. Presently, little information is available concerning the molecular regulatory mechanisms of D. involucrata when responding to heat stress.
High temperature can cause plant growth retardation and even death, which could become a critical issue in the near future due to global warming. Climate predictions suggest that, toward the end of the 21st century, the global average temperature is projected to be elevated by 1-4 • C [10]. Hence, a fundamental understanding of heat response mechanisms in plants is of vital importance for the effective and efficient conservation of rare and endangered species. In general, there are four main heat stress regimes: Direct applied heat stress (DAHS), which is used to examine basal thermotolerance; pre-induction heat stress (PIHS), which is used for acquired thermotolerance; gradually increasing heat stress (GIHS), which is used for determining acquired thermotolerance by imitating the natural condition; and mild chronic heat stress (MCHS), which is used to determine mild heat thermotolerance [11]. Plants respond to high-temperature stress through modulating the expression levels of numerous genes, thus changing physiological, cellular, biochemical, and metabolic processes. For instance, high temperatures typically induce the gene expression of heat shock proteins (Hsps), which act as molecular chaperones primarily engaged in the maintenance and/or restoration of protein homeostasis [12][13][14]. Based on their molecular weight, plant Hsps are classified into six major families: Hsp100 (Clp), Hsp90, Hsp70 (DnaK), Hsp60 (chaperonin and GroEL), Hsp40 (DnaJ) and small HSP (sHsp) [15]. Their functions have been thoroughly documented in heat response, during which Hsps are mainly transcriptionally controlled by heat shock transcription factors (Hsfs) [16,17]. Hsfs generally serve as the terminal components of a stress-triggered signal transduction chain, regulating the expression of downstream targets including HSP genes [18][19][20]. Following the initial characterization of three tomato Hsf genes [21][22][23], a lot of knowledge regarding plant Hsf function has been gained and comes predominantly from analyses of Hsfs in tomato (Solanum lycopersicum) and Arabidopsis thaliana [19,24]. Approximately 27 and 21 different Hsf members were identified in tomato and Arabidopsis, respectively [17,25]. These Hsfs have been differentiated into three classes (A, B and C), which are mainly based on the structural features of their oligomerization domains [16][17][18][19]. Each Hsf class is further divided into various sub-classes. Notably, the platform known as HEATSTER provides different webtools and datasets for the identification and classification of Hsfs in plants [17,26]. Apart from Hsps and Hsfs, many other regulators have been implicated in heat stress response, such as abscisic acid (ABA), auxin, salicylic acid (SA), ethylene, multiprotein bridging factor 1c (MBF1c), and FK506-binding proteins (FKBPs) [19,[27][28][29]. There have been a lot of breakthroughs in the study of plant heat stress response; however, many functional genes involved in thermotolerance remain to be discovered and identified.
Transcriptomic studies have been widely adopted to dissect the molecular mechanisms by which plants adapt to adverse environmental conditions. In recent years, a broad spectrum of heat-responsive genes and/or pathways have been identified from various plant species by transcriptomic approaches. For instance, the de novo transcriptome sequencing and gene expression profiling of spinach (Spinacia oleracea) leaves under heat stress indicated that calcium signaling molecules and certain transcription factors are required for early heat response [30]. The comprehensive transcriptomic analysis of Korean fir (Abies koreana) revealed many Hsps and transcription factors as being differentially expressed upon heat treatment [31]. In addition to the previously described heat stress regulators, the transcriptional profiling of perennial ryegrass (Lolium perenne) under heat stress also identified four genes primarily involved in C 4 carbon fixation as well as one target of rapamycin (TOR) gene [32]. In maize (Zea mays), a transcriptome analysis highlights the crucial roles of the protein processing in the endoplasmic reticulum pathway in heat stress response [33]. It has been suggested that heat tolerances of reproductive and vegetative tissues are distinct and do not always correlate; it has also been suggested that reproductive tissues are highly sensitive to heat stress and are directly associated with fruit and seed production [11]. Therefore, transcriptomic studies have also been applied to unravel molecular responses to heat stress in reproductive tissues at specific developmental stages. For example, an RNA-Seq analysis in rice floral organs revealed that the expressions of five pectinase genes are highly responsive to heat stress at the booting stage [34]. Furthermore, the analysis of the transcriptome and proteome during development and heat stress response of tomato pollen demonstrated that, with the exception of Hsfs and Hsps, protein levels under heat stress are modulated by the interplay between protein synthesis and degradation than by changes of transcript levels [35]. In D. involucrata, the mining of genes associated with seed abortion was enabled through a transcriptome analysis [36]. However, the transcriptomic response of D. involucrata to heat stress has not yet been characterized.
In the present study, heat stress-induced transcriptome changes in D. involucrata were investigated, and potential responsive genes and pathways were identified by RNA-sequencing (RNA-Seq). To our knowledge, this work presents the first transcriptome profiling of D. involucrata responding to high-temperature stress. These findings improve our understanding of the adaptation mechanisms of D. involucrata to environmental stress and will serve as valuable genetic resources to inform further studies regarding the improvement of heat tolerance in D. involucrata.

Sample Preparation and Heat Treatment
After stratification in moist sand for more than one year, D. involucrata seeds were sown in the field of Deyang nursery, Deyang city, Sichuan Province. Subsequently, one-year-old D. involucrata seedlings obtained from Deyang nursery were cultivated for one additional year in pots in the greenhouse of China West Normal University. Two-year-old seedlings at uniform stages of development were placed in growth chambers for a one-week adaptation with the following settings: Relative humidity 75%-85%, 22 • C/16 • C (day/night), and a 14-h photoperiod (lights from 7:00 AM to 9:00 PM) with 400 µmol s −1 m −2 light intensity. One week later, healthy mature leaves of approximately the same size were collected at 0, 1, 6, and 12 h after the onset of heat treatment (42 • C, starting at 8:00 AM). Three biological replicates were set for each time point. All harvested samples were immediately frozen in liquid nitrogen and subsequently stored in a −80 • C freezer until RNA extraction.

RNA Extraction, cDNA Library Preparation and Transcriptome Sequencing
Total RNA was separately extracted from each sample using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA quality was analyzed with 1.0% agarose gel electrophoresis and using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). The cDNA libraries were generated from these RNA samples for Illumina sequencing according to the standard Illumina protocol. The library preparations were then subjected to sequencing on an Illumina HiSeq 4000 platform, and 150-bp paired-end reads were yielded. The sequencing raw data are available in NCBI Short Read Archive (SRA) database (accession: SRP186742).

Transcriptome Assembly and Functional Annotation
Raw reads/data in the FASTQ format were firstly processed using in-house perl scripts. In this step, clean reads/data were generated by filtering reads with adapter, reads in which unknown bases represented more than 10% as well as reads containing more than 50% of bases that had a Q-value ≤20 from raw data. At the same time, the Q20, Q30, GC-content and sequence duplication level of the clean data were calculated. All the downstream analyses were based on clean data with high quality. Since the genome information of D. involucrata remains unavailable, the resulting clean reads/data were used to de novo assemble the transcriptome with Trinity software [37]. The longest transcript of each gene was taken accordingly as the unigene. The putative functions of D. involucrata unigenes were annotated by searching against different public databases, including the protein family (Pfam), NCBI non-redundant protein sequences (Nr), Swiss-Prot, NCBI nucleotide sequences (Nt), Kyoto Encyclopedia of Genes and Genomes ortholog (KO), EuKaryotic Orthologous Groups (KOG), and Gene Ontology (GO) database.

Differential Expression Analysis
The read counts of each unigene in a sample were obtained through mapping clean reads back onto the assembled D. involucrata transcriptome by RSEM v1.2.15 [38], with the bowtie2 parameter set at mismatch 0. Differential expression analyses among distinct conditions/groups was implemented with the DESeq2 based on the read counts from the normalized RNA-Seq data [39]. DESeq2 provides statistical routines for determining differential expression in digital gene expression data using a model based on a negative binomial distribution [39]. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. The D. involucrata genes with an adjusted p-value of less than 0.05 detected by DESeq2 were deemed to be differentially expressed genes (DEGs). A GO functional enrichment analysis of DEGs was carried out using the GOseq R package based on a Wallenius non-central hyper-geometric distribution [40]. For the KEGG pathways analysis, the KOBAS software was utilized to determine the statistical enrichment of the DEGs [41].

qRT-PCR Analysis
To verify the reliability of our transcriptome sequencing data, eight heat-responsive D. involucrata genes were selected for a quantitative real-time PCR (qRT-PCR) analysis. qRT-PCR was conducted with gene-specific primers (Table S1). DiActin was used as an internal control for normalization [36]. The reactions were carried out using StepOnePlus Real-Time PCR system (Applied Biosystems, Foster, CA, USA) with the SYBR Premix Ex Taq Kit (TaKaRa, Dalian, China). Each reaction was performed in triplicate. The melt curve was generated to validate the specificity of amplification. The relative gene expression was calculated based on the 2 −∆∆Ct method [42].

Illumina Sequencing and Transcriptome Assembly
To unravel the molecular mechanisms underlying D. involucrata responses to heat stress, we prepared twelve cDNA libraries from plant samples grown for 0 h (HS0), 1 h (HS1), 6 h (HS6) and 12 h (HS12) at 42 • C. A total of 621,415,654 raw reads were yielded from these libraries (Table S2). After quality trimming, 599,448,362 clean reads were generated. Subsequently, the resulting clean reads were assembled to generate 405,798 transcripts (Table 1). In addition, we further clustered these transcripts into 138,923 unigenes, ranging from 201 to 17,633 bp in length ( Table 1). The mean length and N50 value of the unigenes were 941 bp and 1493 bp, respectively (Table 1).

Identification of DEGs
Subsequently, the gene expression profiles were compared between different treatment groups to identify significant DEGs (Padj < 0.05) based on the criteria of twofold differential regulation [|log2fold change)| > 1] (Figure 3a). There were a total of 37,169 DEGs between the three heat-treated samples and the control samples. The comparisons between the HS1 and HS0, HS6 and HS0, and HS12 and HS0 resulted in 19,532, 20,497 and 27,716 DEGs, respectively (Figure 3b). Among the DEGs, 10,557, 11,607 and 14,138 were found to be up-regulated at HS1, HS6 and HS12, respectively, and 8975, 8890 and 13,578 were found to be down-regulated at HS1, HS6 and HS12, respectively ( Figure 3b). Furthermore, 9537 D. involucrata genes displayed significantly altered expression levels under all three heat-treatment conditions (Figure 3c). A full list of assembled unigenes and their transcript abundance changes under heat stress have been provided in Table S4.

Identification of DEGs
Subsequently, the gene expression profiles were compared between different treatment groups to identify significant DEGs (Padj < 0.05) based on the criteria of twofold differential regulation [|log 2 fold change)| > 1] (Figure 3a). There were a total of 37,169 DEGs between the three heat-treated samples and the control samples. The comparisons between the HS1 and HS0, HS6 and HS0, and HS12 and HS0 resulted in 19,532, 20,497 and 27,716 DEGs, respectively (Figure 3b). Among the DEGs, 10,557, 11,607 and 14,138 were found to be up-regulated at HS1, HS6 and HS12, respectively, and 8975, 8890 and 13,578 were found to be down-regulated at HS1, HS6 and HS12, respectively ( Figure 3b). Furthermore, 9537 D. involucrata genes displayed significantly altered expression levels under all three heat-treatment conditions (Figure 3c). A full list of assembled unigenes and their transcript abundance changes under heat stress have been provided in Table S4.

GO and KEGG Enrichment Analyses of Heat-Responsive Genes
To obtain more insight into the biological functions of heat-responsive genes in D. involucrata, the DEGs were subjected to a GO term enrichment analysis. In general, all the GO-annotated DEGs-8471 (43.37%) at HS1, 8820 (43.03%) at HS6, and 12,426 (44.83%) at HS12-were enriched in 79 terms (Table S5). Moreover, "binding" and "catalytic activity" were found to be enriched by as many as 4658 (54.99%) at HS1 and 5524 (44.46%) at HS12, respectively (Table S5). Notably, the GO term "protein phosphorylation" was significantly enriched at all time points (Table S5), indicating the important role of protein phosphorylation in the D. involucrata response to heat stress.

GO and KEGG Enrichment Analyses of Heat-Responsive Genes
To obtain more insight into the biological functions of heat-responsive genes in D. involucrata, the DEGs were subjected to a GO term enrichment analysis. In general, all the GO-annotated DEGs-8471 (43.37%) at HS1, 8820 (43.03%) at HS6, and 12,426 (44.83%) at HS12-were enriched in 79 terms (Table S5). Moreover, "binding" and "catalytic activity" were found to be enriched by as many as 4658 (54.99%) at HS1 and 5524 (44.46%) at HS12, respectively (Table S5). Notably, the GO term "protein phosphorylation" was significantly enriched at all time points (Table S5), indicating the important role of protein phosphorylation in the D. involucrata response to heat stress.

Heat Shock Transcription Factors Responding to Heat Stress
The transcriptional responses of Hsf genes to high-temperature stress were then examined in D. involucrata. Hsfs with an altered expression in at least one of heat-treated conditions have been

Heat Shock Transcription Factors Responding to Heat Stress
The transcriptional responses of Hsf genes to high-temperature stress were then examined in D. involucrata. Hsfs with an altered expression in at least one of heat-treated conditions have been

Verification of D. involucrata DEGs by qRT-PCR
To prove the reliability of our RNA-Seq results, eight candidate genes associated with heat stress response were selected for qRT-PCR. These genes included HSFA2b

Discussion
The augmented prevalence of high temperature worldwide is of global concern, as it dramatically affects plant growth, productivity, and distribution [44][45][46]. Though recent years have seen tremendous progress in unraveling the complex mechanisms underlying heat stress response in many different plants, more investigation is still required, particularly for uncharacterized species. The rare and endangered plant D. involucrata favors cold and humid climatic conditions in the mountains of southwestern and south-central China [5]. Therefore, high temperature is the

Discussion
The augmented prevalence of high temperature worldwide is of global concern, as it dramatically affects plant growth, productivity, and distribution [44][45][46]. Though recent years have seen tremendous progress in unraveling the complex mechanisms underlying heat stress response in many different plants, more investigation is still required, particularly for uncharacterized species. The rare and endangered plant D. involucrata favors cold and humid climatic conditions in the mountains of southwestern and south-central China [5]. Therefore, high temperature is the crucial environmental factor that restrained its natural population distribution, ex situ conservation, as well as landscape application. At present, information on the large-scale regulation of gene expression profiles in response to heat stress is still lacking for D. involucrata. In this study, RNA-Seq, an invaluable approach allowing for the comprehensive characterization of transcriptomic events occurring in certain conditions, was utilized to uncover the heat response mechanism in D. involucrata. After de novo assembly, we generated 138,923 unigenes, of which 69,743 unigenes (50.2%) were annotated in public databases. Nevertheless, more complete transcriptome data of D. involucrata could be further obtained by integrating full-length transcriptome analyses in the future. Moreover, 19,532, 20,497 and 27,716 DEGs were identified after 1 h, 6 h, and 12 h of heat treatment in comparison to 0 h, respectively, and a functional enrichment analysis was subsequently performed for these DEGs. Taken together, our study gives the first insights into a global transcriptome picture of D. involucrata under heat stress, thus opening the door for adaptation research at the molecular level in this rare and endangered plant.

Heat Stress Treatment Used for Transcriptome Analysis in D. involucrata
Heat stress may result in either eustress or distress, both of which are determined by the thermic threshold of the organism, its state of acclimation and adaptation, as well as stress severity and duration [11,47]. Depending on the duration and strength of heat stress, plants can activate transient or long-term responses for protection and recovery which engage massive transcriptional and translational alterations. By contrast, extremely high temperature or long-term heat stress often lead to an irreversible disturbance of homeostasis, which can result in cell death [48]. Apart from basal thermotolerance, plants also employ conditionally expressed acquired thermotolerance at the cellular and molecular levels to withstand acute heat stress [11,49]. Acquired thermotolerance can be experimentally achieved through subjecting plants to a moderate level of heat stress (known as "priming"). Alternatively, priming can also be induced by a slow and continuous increase in temperature [11,49]. Heat stress is initially perceived by the leaf, which subsequently brings about physiological and biochemical alterations at the whole plant level [50,51]. Moreover, researchers have expended a lot of effort on studying heat response mechanisms in leaf tissues, including analyses of heat-induced leaf transcriptome changes [20,30,32,33,43]. Therefore, healthy mature leaves of approximately the same size were used for the transcriptome analysis of heat stress response in D. involucrata. Since photosystems in the leaves of D. involucrata seedlings were affected during a heat treatment of around 42 • C for several hours and for a longer period [9], a time-course experiment was performed at 42 • C in this study. After a one-week adaptation in growth chambers, D. involucrata seedlings were subjected to the 42 • C heat treatment, and leaf tissues were collected at time points 0, 1, 6, and 12 h after the start of heat treatment for transcriptome sequencing. Similarly, time-course gene expression profiling was conducted in carnation (Dianthus caryophyllus L.) leaves in response to the 42 • C heat treatment within 12 h [43].

Heat Stress Affected Expression of Genes Associated with Protein Processing in the ER
In plants, sensitive cellular processes including protein folding in the endoplasmic reticulum (ER) can be disrupted by environmental stresses. The limited capacity of the ER to process proteins leads to massive accumulation of unfolded or misfolded proteins, which is the greatest risk for plant development and survival during heat stress [52,53]. Hence, protein processing in the ER is a very important pathway under heat stress. In this study, a KEGG pathway analysis revealed that "protein processing in endoplasmic reticulum" was significantly enriched under all three heat-treatment conditions ( Figure 4 and Table S6). There were 555 DEGs at HS1, with 14 of them down-regulated; there were 390 DEGs at HS6, with 21 of them down-regulated; and there were 424 DEGs at HS12, with 46 of them down-regulated ( Figure 5, Figure S1, and Table S6). These results suggest that most of the genes engaged in protein processing were induced during heat stress to promote the efficiency of protein processing for stress response.
In higher plants, heat stress is known to redirect protein synthesis via a reduction in the synthesis of normal proteins concomitantly with a remarkable enhancement in the transcription and translation of Hsps [54]. Mounting evidence suggests that transcript abundance of Hsps along with chaperones is implicated in heat stress adaptation [30][31][32][33]43]. Hsps serve as molecular chaperones that are responsible for the folding, assembly, translocation and degradation of protein molecules, and they play paramount roles in combating stresses through re-establishing normal protein conformation and, thereby, cellular homeostasis [13,44]. In protein processing in the ER pathway, we observed that a lot of DEGs encoding Hsps (i.e., HSP90, HSP70, HSP40 and sHSP) were substantially up-regulated upon heat treatment. These results imply that inductions of Hsps facilitate protein processing in the ER and are important for acclimating to high-temperature stress in D. involucrata.
To eliminate misfolded proteins in the ER and maintain homeostasis, an ER-associated degradation (ERAD) process is usually activated [55]. Through ERAD, misfolded proteins are directed for degradation via the ubiquitin-proteasome system. DERLIN is the core component of the retrotranslocation machinery of the plant ERAD pathway, and it was reported to be induced responding to ER stress [56,57]. Consistently, we observed an up-regulation of DERLIN at HS1 and HS6 ( Figure 5 and Figure S1). Ubc6 in yeast is a ubiquitin-conjugating enzyme that localizes to the ER membrane and is part of the Doa10 complex in the ERAD pathway [58]. In this study, the plant homolog of Ubc6, UBC32, was induced at all three heat treatment time points ( Figure 5 and Figure S1). Likewise, the transcriptional levels of UBC32 in Arabidopsis were found to be elevated by various abiotic stresses and ER stress [55].
CHIP in Arabidopsis acts as the chaperone-associated E3 ubiquitin ligase and has been implicated in protein degradation and stress responses [59]. Importantly, CHIP is induced under stress conditions including high temperature [60]. In accord with this, we found CHIP was up-regulated in D. involucrata at the treatment time points HS1 and HS6 ( Figure 5 and Figure S1). It has been demonstrated that disruption in CHIP resulted in elevated sensitivity to abiotic stresses [59]. However, transgenic plants overexpressing CHIP are hypersensitive to temperature stress [60]. The exact functions of CHIP in D. involucrata during heat stress need to be dissected in the future.

Influence of Heat Stress on Plant Hormone Signal Transduction
Phytohormones are regarded as master regulators of plant growth, development, and stress responses [61][62][63]. Our KEGG pathway analysis of DEGs revealed that "plant hormone signal transduction" was significantly altered under heat stress ( Figure 4 and Table S6). More specifically, all hormone regulatory pathways were affected at three heat treatment time points, except for the salicylic acid pathway at HS6 ( Figure S2). Though abscisic acid (ABA) is the best-characterized stress hormone, the roles of other hormones such as auxin, cytokinin (CK), ethylene (ET), jasmonic acid (JA) and brassinosteroid (BR) in plant stress responses are beginning to be better understood [64,65]. Furthermore, emerging evidence indicates that crosstalk between different hormonal signaling pathways leads to synergetic or antagonistic interactions, which modulate stress adaptations [66][67][68]. Based on our transcriptome data, it is tempting to assume that multiple phytohormones and their signaling crosstalk are engaged in the heat stress response in D. involucrata.
Many studies have revealed that phytohormone signaling and protein processing in the ER interact coordinately and cooperatively in stress responses. For instance, Arabidopsis UBC32 was demonstrated to be engaged in salt stress responses through the regulation of BR signaling by the ERAD pathway [55]. In maize, the crosstalk between the ER quality control (ERQC) and ABA signaling engages the transcription factor ZmbZIP17 [69]. Fascinatingly, ER stress was shown to negatively regulate auxin signaling in Arabidopsis [70]. Furthermore, the ER-based auxin homoeostasis is crucial for unfolded protein response (UPR) activation [70]. In our study, auxin signaling might have been the primary hormone regulatory pathway in the heat stress response in D. involucrata based on the maximum amount of DEGs (i.e., 82 DEGs) assigned to the auxin pathway (Table S7). Moreover, majority of the DEGs (86.3% at HS1, 91.5% at HS6, 90.1% at HS12) in the auxin pathway were suppressed during heat stress (Table  S7). It has been previously demonstrated that hormone synthesis, transport, and response for most part of the phytohormone pathways are conserved among different plant species [71]. In particular, the regulatory regime of auxin signaling can be transferred to all phototrophic eukaryotes containing co-orthologues except for the unicellular green algae Chlamydomonas reinhardtii [71]. Therefore, it is likely that the UPR in D. involucrata leads to changes in the auxin signaling pathway to cope with heat-induced ER stress.

Hsfs May Play a Pivotal Role During Heat Stress
Accumulating data indicate the existence of a heat stress-activated transcriptional network that modulates the expression of heat-responsive genes [53]. In such a network, Hsfs were documented to play a central and highly conserved role in eukaryotes [72]. In tomatoes, HsfA1a possesses a unique role as a master regulator for acquired thermotolerance [20,73]. The knockdown of HsfA1a expression in tomato practically eliminates heat stress-triggered synthesis of HsfA2, HsfB1, and molecular chaperones, suggesting no substitution of function by other Hsfs [73]. In Arabidopsis, the HsfA1 group comprises Hsf A1a, A1b, A1d and A1e. It has been demonstrated that HsfA1s in Arabidopsis act as major regulators in heat-responsive gene expression [74,75]. Transcriptomic studies using HsfA1s quadruple knockout mutants revealed that more than 65% of up-regulated genes during heat stress were HsfA1-dependent [74]. These induced genes that are under the control of HsfA1s include several Hsfs, such as HsfA2, HsfA7a, HsfB1, HsfB2a, and HsfB2b [74]. Consistent with these studies, we observed that most of the DEGs belonging to class A1a, A2, A7a, B1, B2a, and B2b were up-regulated upon heat treatment (Table S8), implying the possibility that HsfA1a-dependent transcriptional regulation networks also exist in D. involucrata. Other DEGs encoding Hsfs belong to class A3, A4, A6, B3, B4, and C1 (Table S8). Future studies may dissect the molecular functions of these Hsfs in D. involucrata during heat stress.

Conclusions
The endangered and rare tertiary relict tree D. involucrata is highly sensitive to high temperatures. To dissect the dynamic and complicated gene regulatory network in D. involucrata seedlings under heat stress, we carried out time-course gene expression profiling using Illumina HiSeq 4000 sequencing. Compared to the control (0 h), we identified 19,532, 20,497 and 27,716 DEGs after 1 h, 6 h, and 12 h of heat treatment, respectively. Furthermore, a KEGG pathway analysis of DEGs revealed that protein processing in the ER and plant hormone signaling may play crucial roles in D. involucrata responses to high-temperature stress. Additionally, the expression levels of 32 genes encoding putative Hsfs were significantly altered during heat stress. Our transcriptomic data will facilitate future studies on such heat-sensitive and endangered species and will provide the theoretical basis for conservation and utilization of these valuable plant resources.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/10/8/656/s1, Table S1: Primers used for the qRT-PCR; Table S2: Summary of RNA-Seq data from all samples; Table S3: Annotation of unigenes in different databases; Table S4: A full list of assembled unigenes and their annotation and expression levels during heat stress; Table S5: DEGs enriched in GO terms; Table S6: KEGG pathway analysis of DEGs; Table S7: Classification of DEGs associated with plant hormone signal transduction; Table S8: DEGs encoding heat shock transcription factors (Hsfs); Figure S1: Analysis of DEGs associated with protein processing in endoplasmic reticulum in two comparisons [HS6vsHS0 (a) and HS12vsHS0 (b)]; Figure S2: Analysis of DEGs associated with plant hormone signal transduction in three comparisons [HS1vsHS0 (a), HS6vsHS0 (b) and HS12vsHS0 (c)]; Figure S3: Correlation analysis between data of RNA-Seq (x-axis) and qRT-PCR (y-axis); File S1: cDNA sequences of 32 Hsfs that displayed differential expressions under heat stress; File S2: cDNA sequences of eight heat-responsive genes used for qRT-PCR.
Author Contributions: Q.L. and X.X. conceived and designed the research. Q.L. performed the experiments. Q.L., R.R.V., W.X. and X.X. analyzed the data. Q.L. wrote the manuscript with input from R.R.V. and X.X. All authors read and approved the manuscript.