De novo RNA sequencing transcriptome of Rhododendron obtusum identified the early heat response genes involved in the transcriptional regulation of photosynthesis

Rhododendron spp. is an important ornamental species that is widely cultivated for landscape worldwide. Heat stress is a major obstacle for its cultivation in south China. Previous studies on rhododendron principally focused on its physiological and biochemical processes, which are involved in a series of stress tolerance. However, molecular or genetic properties of rhododendron’s response to heat stress are still poorly understood. The phenotype and chlorophyll fluorescence kinetics parameters of four rhododendron cultivars were compared under normal or heat stress conditions, and a cultivar with highest heat tolerance, “Yanzhimi” (R. obtusum) was selected for transcriptome sequencing. A total of 325,429,240 high quality reads were obtained and assembled into 395,561 transcripts and 92,463 unigenes. Functional annotation showed that 38,724 unigenes had sequence similarity to known genes in at least one of the proteins or nucleotide databases used in this study. These 38,724 unigenes were categorized into 51 functional groups based on Gene Ontology classification and were blasted to 24 known cluster of orthologous groups. A total of 973 identified unigenes belonged to 57 transcription factor families, including the stress-related HSF, DREB, ZNF, and NAC genes. Photosynthesis was significantly enriched in the Kyoto Encyclopedia of Genes and Genomes pathway, and the changed expression pattern was illustrated. The key pathways and signaling components that contribute to heat tolerance in rhododendron were revealed. These results provide a potentially valuable resource that can be used for heat-tolerance breeding.


Introduction
Rhododendron (Rhododendron spp. L.) is one of the most popular woody flowers worldwide; more than 800 rhododendron species are diploid with 2n = 2x = 26 and with a monoploid genome size ranging from 0.67 pg to 0.83 pg [1]. In temperate and tropical zones, a greater risk of heat stress occurs in the introduction and cultivation of rhododendron plants. Four rhododendron species showed the highest absolute and relative net photosynthetic rates in PLOS  sequence information, especially for the functional genes involved in the formation of importance traits, remains a major obstacle that severely restricts the future application of molecular-assisted selection (MAS) or GM breeding. Considering that many species of rhododendron are important ornamental plants, a better understanding of the heat tolerance of different species and provenances and the specific characteristics and molecular mechanisms that influence heat tolerance are necessary for the MAS and the improvement of rhododendron breeding in different adverse environments. The objectives of this study were as follows: 1) to assess the heat tolerance ability of the commonly used rhododendron varieties in South China and select a suitable material for the executive RNA sequencing; and 2) to achieve the expression profiles and unigene sequences, especially for the TF regulator involved in the stress response, by using the RNA sequencing transcriptome, which assists in further gene cloning and for conducting functional research on molecular mechanisms underlying different heat response variations in different species. Among the 92,463 unigenes in transcriptome, 38,724 unigenes with annotation were identified. The results from this study provide abundant genomic information for future research and clarify the molecular mechanisms governing the heat responsive process involved in the transcriptional regulation of photosynthesis.

Plant materials and growth conditions
Samples of rhododendron "Yanzhimi" (R. obtusum), "Xiajinpao" (R. hybridum), "Jinpao" (R. hybridum), and "Yanchun" (R. hybridum) were used as experimental materials. Two-year-old self-rooting plants were transferred and grown in the greenhouse for 3 months in Wuhan Academy of Agricultural Sciences and Technology. Control plants were then grown in a growth chamber under continuous light (20-40 μm sec -1 m -2 ) at 25˚C. For heat treatments, the plants were transferred to 40˚C under the same light source. Fresh leaves were collected after 4 and 24 h of treatment. Tissue from 15 randomly selected mature leaves were pooled to minimize biological variability for each sample and frozen in liquid nitrogen. Samples of leaves were separated and frozen in liquid nitrogen for further analysis, and three biological replicates were set for each time point.

Measurement of chlorophyll fluorescence parameters
The third fully expanded leaf attached to the plant was subjected to a pulse−amplitude modulation fluorometer (PAM-2500, Walz, Germany) to detect the chlorophyll fluorescence parameters. After 20 min of dark adaptation, the minimum fluorescence level (Fo) was determined under low-intensity measuring light. Maximum fluorescence level (Fm) was measured after 0.5 s saturating pulse at 4,000 μmol m −2 s −1 . Steady-state fluorescence level (Fs) was obtained after 20 min of actinic light (234 μmol m −2 s −1 ) adaptation. Light-adapted maximum fluorescence level (Fm') was measured with a second saturating pulse (0.5 s, 4,000μmol m −2 s −1 ). The actinic light was then closed, and the light-adapted minimum fluorescence level (Fo') was determined using a far-red light for 5 s. These parameters were used to obtain four identification indexes, as follows: Fv/Fm = (Fm−Fo)/Fm, Fv'/Fm' = (Fm'-Fo')/Fm' [21]. The parameters Fv/Fm and Fv'/Fm' represent the maximum photochemical quantum yield of PS II and the excitation energy capture efficiency of PSII reaction center, respectively. The survival rate of 30 leaves was observed and statistically calculated for each sample with three biological and technical repeats.

RNA isolation and construction of mRNA-Seq library
Total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. The RNA quality was evaluated using Bioanalyzer 2100 (Agilent Technologies, CA, USA), and the A260/A280 ratios of all samples varied from 2.0 to 2.1. Poly (A) + mRNA was purified with oligo(dT) beads. The mRNA was randomly segmented into small fragments by the divalent cations (Fragmentation Buffer, Illumina, Hayward, CA) at 94˚C for 5 min, which were then used as templates to synthesize the first-strand cDNA using random hexamer primers. The second-strand cDNA was synthesized using RNAseH and DNA polymerase I. Short cDNA fragments were purified with QiaQuick PCR extraction kit. The cDNA fragments were connected with sequencing adapters according to Illumina's protocol (San Diego, CA USA). After agarose gel electrophoresis, the target fragments of 300-500 bp were selected for PCR amplification to create the final cDNA library. All of these raw data have been deposited and released in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/) under bioproject accession PRJNA398170.

Processing and mapping of Illumina reads
Nine cDNA samples from the control and from 4 and 24 h heat-treated plants were prepared for sequencing. The library was sequenced using Illumina HiSeq X Ten. Raw reads were cleaned by removing the adaptor sequences. Empty reads and reads with unknown sequences "N" were removed before data analysis. All the following analysis were built on getting a precise and unique trascript sequence. So in order to obtain the sequence, we separately subjected cleaned reads to de novo assembly using Trinity software (http://trinityrnaseq.sourceforge.net/ ) with default k mer size set to 25 bp [22]. The Trinity software carried out the assembly by three steps. Firstly, the clean reads were first assembled to produce longer fragments named as trinity. Secondly, the paired-end reads were used to calculate the distance and relation among these transcripts and to avoid obtaining trinities from the same transcript. Finally, the transcripts were connected until they cannot be extended on either end. These sequences were subjected to redundancy removal by using a sequence clustering software to obtain nonredundant transcripts, and the gene expressive quantity that based on the reads count were done at the same time by the software Corset [23].. These non-redundant transcripts were named unigenes. Bowtie was used to map the clean reads back to the assembled transcripts [24]. Transcripts with 70% mapping coverage were discarded from the final assembly.

Gene function annotations and classifications
Swiss-Prot, Cluster of Orthologous Groups (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases were used to annotated all the unigenes with local BLAST programs (E value < 1.0E -5 ). For the Swiss-Prot annotation results, the Blast2GO program was applied to classify the unigenes based on the Gene Ontology (GO) terms [25], and the WEGO tool was used to draw the GO tree [26]. The COG database was also used to predict the functional classification of the unigenes. All unigenes were matched with the KEGG pathway database by using the BLASTx to predict the pathway in which the unigenes participated (E value < 1.0E -5 ).

Identification of differentially expressed genes (DEGs)
Differentially expressed genes (DEGs) of different libraries were analyzed using the fragments per kb of exon region in a given gene per million mapped fragments (FPKM) method and were identified with edgeR package [27]. FPKM was used to determine the number of mapping reads for every unigene and to assess the unigene expression levels. Negative binomial distribution methods were applied to calculate the p values. Benjamini-Hochberg methods were chosen to adjust for the multiple tests. The unigenes were significantly and differentially expressed when they had a p-value < 0.0001 and a false discovery rate (FDR) < 0.05. GO functional enrichment and KEGG pathway enrichment analysis of DEGs were conducted via comparison with the whole-transcriptome background using the formula described in previous studies (Bonferroni-corrected p-value 0.05) [28]. Note that the database that used for KEGG and GO mapping is derived from kobas database (http://kobas.cbi.pku.edu.cn/), and the mapping of the differentially expressed gene is performed by kobas software [29].
Validation by real-time quantitative RT-PCR Thirteen unigenes were chosen to confirm their involvement in the response to high temperature using RT-qPCR method. Primer Premier 6.0 software was used to design gene-specific primers based on the selected unigene sequences (see S1 Table). RT-qPCR was conducted on an ABI Stepone plus System using Fast Start Universal SYBR Green Master Mix (Roche). The PCR amplifications were as follows: 95˚C for 3 min, followed by 40 cycles of 95˚C for 15 s, and then 60˚C for 30 s. A melting curve was obtained from 60˚C to 95˚C by increasing the temperature stepwise by 0.5˚C every 5 s. This step was performed to test the specificity of the amplified product. The expression levels of selected DEGs were normalized by comparing with two internal reference genes, HIS3 and GADPH [30]. Their relative expression levels were calculated via the 2−ΔΔCt method. All RT-qPCR experiments were repeated in three biological replications and in three technical replications.

Results
R. obtusum "Yanzhimi" shows better performance of high leaf survival rate and photosynthesis protection under heat stress After 3 days of treatment under 38˚C/35˚C, all the leaves of "Xiajinpao", "Jinpao", and "Yanchun" were withered yellow, and the plants died. In all gradients and high temperature treatments, "Yanzhimi" became pale and showed wilting symptoms after exposure to 38˚C/35˚C for 6 and 9 days, but grew well under other high temperature stress conditions (Fig 1).
The chlorophyll fluorescence kinetic parameters in the four species of rhododendrons were measured using the imaging-PAM system to estimate the effect of heat stress on photosynthesis and explore the relationship between heat tolerance and photosynthetic activities. Three parameters, including the maximum photochemical efficiency of photosystem II (Fv/Fm), efficiency of excitation capture of open PSII center (Fv'/Fm'), and initial fluorescence (Fo), were selected for further investigation.
Under 38˚C/35˚C treatment, the phenotype of rhododendron was severely damaged and showed diverse differences. Fig 2a shows that the Fo value increased at the 3-day time point and decreased after 6 and 9 days in the four cultivars. The 38˚C/35˚C temperature treatment more severely affected the "Jinpao" and "Yanchun" plants as compared with "Yanzhimi" after 3 days. The Fo value decreased rapidly to an undetectable amounts at the 6-day time point, thereby indicating that the leaf photosynthetic apparatus was irreversible damaged (Fig 2a).
Fv/Fm showed a decreasing trend in the four rhododendrons. Except for "Yanzhimi", the values drastically decreased after 3 days and were undetectable after 6 days. The Fv/Fm of "Yanzhimi" also decreased from 0.79 to 0.77, 0.718, and 0.66 after 3, 6, and 9 days, respectively (Fig 2b).  "Yanzhimi" showed a slower decrease in ChlF parameters and higher survival rate of leaves than the other three rhododendrons, which indicated than "Yanzhimi" had the strongest protective system that maintains PSII and the best heat tolerance.

Illumina sequencing and de novo assembly
Nine cDNA libraries were generated with mRNA from the three sample groups: control group (grown at 25˚C), group1 (heat treatment exposed at 40˚C for 4 h), and group2 (heat treatment exposed at 40˚C for 24 h) for the full expanded leaf. These cDNA libraries were then subjected to Illumina deep-sequencing. Clean reads with 4,349,881,212, 4,122,064,627, and 4,508,391,358 nucleotides were obtained from the control, group1, and group2, respectively ( Table 1). The clean reads were >114 Gb (total length), which was equivalent to~160 fold coverage of the genome of R. obtusum (about 0.7 Gb). All clean reads were de novo assembled with the Trinity method because R. obtusum does not have an appropriate reference genome sequence. Most of the reads could be mapped back to the assembled transcripts. A total of 325,429,240 clean reads were assembled into 395,561 transcripts and then clustered into 92,463 unigenes with a total unigene length of 97,052,187 bp. The length of the unigenes ranged from 201 bp to 15,666   (Fig 3a and 3b). For the GC content, most transcripts and unigenes were distributed within the range of 40%-50% (Fig 3c and 3d), which is identical to the GC content of each sample's clean data in Table 1.

Functional annotation of the assembled transcriptome
To predict and analyze the function of the assembled unigenes, we assessed the non-redundant sequences using a BLASTX search against the following databases: Swiss-Prot (a manually annotated and reviewed protein sequence database), protein family database (Pfam), Gene Ontology (GO), cluster of orthologous groups (COG), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Orthologous Groups of proteins (EGGnog  (Table 3). Overall, 40,518 (43.82%) unigenes were significantly matched to the known genes in the above-mentioned public databases (S2 Table). The use of GO terms for functional categorization was performed according to the Swiss-Prot annotation. The main GO terms included biological process (BP), cellular component (CC), and molecular function (MF). Based on the sequence homology, 92,463 unigenes were mainly categorized into 51 functional groups (Fig 4a). The GO analysis indicated that many identified genes were associated with various biological processes and molecular functions under heat stress. In the category of BP, the largest groups were associated with cellular and metabolic processes. About 23,106 unigenes were annotated in the metabolic process category, which allowed the identification of novel genes involved in the secondary metabolism pathways during heat acclimation. For the MF category, the unigenes with binding and catalytic activities comprised the largest groups. For CC, the top three largest categories were cell, cell part, and organelles (Fig 4a).
To obtain a better understanding of the biological functions of the unigenes, we used the annotated sequences to search in the KEGG database. Among the 92,463 unigenes, 27,514 had significant matches and were assigned to 121 KEGG pathways (Fig 4b, S3 Table). The most   strongly represented pathways for these genes were the metabolism pathways; 767 unigenes out of total 1910 metabolism-related unigenes were enriched under 24 h of heat treatment (Fig  4b, S3 Table). These annotations served as a basis for investigating the processes and pathways involved in the response to heat stress. A total of 23,514 unigenes were subdivided into 24 COG classifications to categorize the orthologous gene products (Fig 4c). Among these classifications, the cluster of "general function prediction only" (R group) represented the largest group, followed by the "posttranslational modification, protein turnover, chaperones" (O group). The two categories involving "cell motility" (N group) and "nuclear structure" (Y group) represented the smallest COG classifications.

DEGs among the heat stress treatment and control in rhododendron
Plants require an effective and quickly response to adapt to the heat stress and survive in high temperatures. This rapid reaction might predominantly rely on a gene network involving TFs, protein kinases, plant hormone signal transduction pathways, and ROS signaling. For the executive DEG analysis, all Illumina sequencing biological repeat samples were compared, and sample 0h-6 showed a relatively low correlation with the 0h-4 and 0h-7 samples. Hence, 0h-6 achieved an accurate DEG result (S1 Fig). Subsequent analysis showed that many rhododendron unigenes were significantly affected by the high temperature treatment. We identified 8,395 unigenes as DEGs, including 5,828 up-regulated unigenes and 2,557 down-regulated unigenes at 4 h; whereas 10994 unigenes were classified as DEGs, including 6,140 up-regulated unigenes and 4,854 down-regulated unigenes at 24 h. Control and heat stress samples were identified using a heat map (S2 Fig). HSPs are environmentally induced proteins that play a role in acquired stress tolerance and enable the plants to cope with heat and other environmental stresses. HSPs also operate and collaborate with other stress response mechanisms and functional components to reduce cellular damage. Analysis of the transcriptome profile in plants after heat treatment indicated that the HSP family plays a central role in the response to heat stress. HSPs families, including HSP100, HSP90, HSP70, HSP60, and sHSPs, are involved in folding and assembling of the protein, maintaining protein stabilization, activating the protein, and degrading the protein in many normal cellular processes and under stress conditions. In this study, 109 DEGs of the HSPs family were detected in the rhododendron transcriptome. These genes showed significant changes under heat treatment (S2 Table). The genes of high and middle molecular weight HSPs, including HSP70, HSP90, and HSP83, were similar to other de novo assembled transcriptomes under heat stress [31]. HSP90 (C18696) showed great differential expression under 4-h and 24-h heat treatments, which represent log2 ratios of 13.32 and 9.72, respectively.

Identification of TFs and target genes in R. obtusum "Yanzhimi"
TFs are important in gene transcriptional expression regulation network by affecting the downstream expression of many genes. TFs can guide many different biological processes, including secondary metabolism, cell division, and cell differentiation. Therefore, the identification of TFs and target genes will provide a better understanding of the regulatory networks of genes in the heat tolerance of R. obtusum.
A total of 973 protein sequences were obtained after ORF prediction and default filter screening. These sequences can be divided into 57 TF families, including the widely reported gene family involved in abiotic stress tolerance; for instance, several DREB/CBFs, MYBs, NACs, and HSFs were characterized in our transcriptome (Fig 5 and Table 4). Among these families, the 33 annotated HSFs can be subdivided into 3 subgroups including 22 HSFAs, 10HSFBs, and 1 HSFC. Their differential expression profiles are shown in S2 Table. The number of TFs in R. obtusum was lower than in the model plant. In Arabidopsis, 1922 known genes encode TFs that can be divided into 64 gene families. In rice, 2025 TF genes were annotated in PlnTFDB [32]. The expression levels of AP2/EREB, AP2/ERF, HSF, MYB WRKY, ZNF, and other gene family members were screened from the rhododendron transcripts by using log2 ratio> 2.0 and P-value <0.001 as the differential expression thresholds. We especially focused on the validation of HSF gene family members at 24-h heat treatment stage. AtDREB1D (C61022), AtDREB1C (C46630), AtDREB2C (C41140), AtDREB2A (C62227), and AtDREB2A (C62226) were reported by the DREB class of TFs. These genes were also reported as important TFs involved in cold, drought, high temperature, and high salt stress in A. thaliana. Fifty-two TFs were finally selected with up-regulated log2 ratio >2.0 and FPKM >10 filtering. These TFs are possible regulators that positively transcriptionally regulate the heat responsive genes in rhododendrons. The expression level of C62226 gene showed the highest up-regulation, and the log2 ratio was as high 9.11 at 24h. This dramatic response indicated that DREB plays an important role in rhododendron under heat stress. Moreover, HSFB2A (C75881), HSFA2 (C80422), HSFA3 (C63903, C63904, C63905, C63906), and HSFA6b (C52862) also show greatly up-regulated changing fold under 4-h or 24-h heat stress treatments, respectively, thereby showing the involvement of the HSF gene family in heat response. The expression of C80422, which hits AtNAC2 gene in Arabidopsis was high, and the log2 ratio was 9.52 at 24 h (Table 2).

Validation and expression pattern analysis
Twenty-two candidate DEGs were selected, and their expression profiles were compared between control and heat stress samples using quantitative real-time (qRT) PCR to further validate the reliability of the Illumina sequencing read analysis and to detect the expression of some candidate regulators. Among the chosen DEGs, dehydration-responsive element-binding protein 2A (C62226), dehydration-responsive element-binding protein 2C (C41140), dehydration-responsive element-binding protein 1D (C61022), heat stress TFA2 (C80422), heat stress TF A3 (C63903, C63904, C63905, and C63906), heat stress TFA6b (C52862), Zinc finger protein ZAT8/ZAT10/ZAT11 (C31931, C36935, C36936, and C36937), NAC domaincontaining protein NAC2/NAC29 (C3629, C81022, C45823, and C45824), and WRKY TFs (C86364) were related to the response to abiotic stresses. The expression patterns of all these  Table 4. The number of unigenes identified as transcription factors in rhododendron leaves. 24 h vs ck, Changed log2 ratio of Rhododendron exposed to 40˚C for 4 h. B. Changed log2 ratio of Rhododendron exposed to 40˚C for 24 h. DEGs obtained through qRT-PCR data were consistent with the high-throughput RNAsequencing (RNA-Seq) data, thereby confirming that the Illumina results were reliable (Fig 6). As shown in Fig 6, most of the selected TFs showed a significantly up-regulated expression pattern under 4 h and 24 h heat treatments. The highest expression levels of many DEGs were observed at 4 h as compared with 0 and 24 h, thereby indicating that these TFs responded to high temperature at the early stage of heat stress and subsequently accumulated to exhibit their regulatory effect.

Effect of heat stress on photosynthesis-related unigenes
Based on the KEGG pathway enrichment analysis, we found that the photosynthesis pathway was significantly changed under 4-h or 24-h heat treatment. In this pathway, 48 unigenes showed differential expressions out of 77 characterized unigenes in rhododendron transcriptome (Additional file 4). These DEGs covered all the sections of this pathway, including PSI, PSII, cytochrome b6/f complex, photosynthetic electron transport, and F-type ATPase. The activity of PSII is restricted by a variety of environmental stresses, and the ROS produced in this process activates the psbA mRNA [33]. Mattoo et al. [34] reported that both electron transport and the synthesis of ATP are important for the synthesis of D1 in Spirodela oligorrhiza. Studies on intact chloroplasts from spinach showed that the level of stromal ATP is correlated with the light-dependent synthesis of the D1 protein in intact chloroplasts [35]. D1 (psbA) and 5 ATPases (beta, alpha, gamma, delta, and b) were up-regulated in this process, thereby supporting the results of the previous study (Fig 7). Photosynthesis under heat stress might be hindered by the reduction of soluble proteins, Rubisco binding proteins, large-subunits (LS), small-subunits (SS) of Rubisco in darkness, and the increase in amounts of these proteins under light [36].

Discussion
Considering that Rhododendron consists of non-model and woody plants, current data on the molecular and genetic properties of species in this genus are severely insufficient. Only a few transcriptome databases were established for R. latoucheae Franch. and R. fortunei [37,38]. These databases focused on root nutrition adoption and development of new EST-SSR markers, respectively, but not on the abiotic stress. The lack of reference genomic and transcriptome data restricted the research on rhododendron, especially the molecular mechanism underlying its complex traits. As a feasible and economical technology for creating relatively comprehensive sequence data in a short time, transcriptome RNA sequencing technology has become popular in plant research [39][40][41][42][43]. In the present study, 92,463 unigenes were assembled. The obtained sequence data can serve as a basis for further studies on gene cloning, expression analysis, and SSR marker development. By comparing the morphological and physiological indexes in different cultivars, we screened a heat tolerant material to identify the related key genes that determine its traits. The up-regulated TFs screened in R. obtusum "Yanzhimi" transcriptome, including HSFs, DREBs, ZATs, and NACs, belonged to important gene families that participate in abiotic stress tolerance, and their homologs were functionally characterized in the model species. For instance, the constitutive expression of the C2H2-EAR zinc finger protein Zat10 in Arabidopsis elevated the expressions of reactive oxygen-defense transcripts [44], and the homologous gene named C65939 in rhododendron transcriptome was up-regulated. In our transcriptome, four HSFA3 homologous genes were identified and named as C63903, C63904, C63905, and C63906, and their expression levels were increased by 10,000-fold after heat treatment for 24 h. Considering that only one gene copy of HSFA3 was found in A. thaliana and rice, the four homologous unigenes detected in our transcriptome are greater in expression than in the model plants. This finding can be explained by several reasons. First, the horticultural variety R. obtusum "Yanzhimi" might be highly heterozygous in this locus, resulting in a larger allelic variation. Second, a draft reference genome of the rhododendron plants has not yet been published. Hence, the evolution of ancestral whole genome polyploidization has not yet been confirmed in rhodondendron but has already been established in other basal eudicots, such as Vitis vinifera and Populus [45,46]. Some HSFA3 paralogs were found. Third, the heat TF could form the truncated S-HsfA2, which is terminated prematurely by variable shear and binds to the promoter sequence [47]. Whether the rhododendron HSFA3 will pass through this mechanism, which further results in detecting multiple transcripts, still needs to be verified by subsequent experiments.
The homologous unigene AtHSFA2 (C80422) shows a significantly up-regulated expression under 24-h heat treatment (change fold > 16,000). Hence, we presumed that this gene functions positively in the heat tolerance of rhododendron. In Arabidopsis, AtHSFA2 directly regulates the chloroplast-localized HSP21 and other heat shock transcript proteins and repairs the thylakoid membrane damage to protect the chloroplasts. This phenomenon is consistent with our findings on its protective activity in the PSII reaction center. As another important member in HSFA subclade gene family, AtHsfA3 can bind directly with DRE cis-element by AtDREB2A/AtDREB2C, thereby playing an important role in plant cold and drought stress signaling by positively regulating the heat tolerance in Arabidopsis. This finding indicated that DREB-HSF-HSP pathway is conserved in plants, and HsfA3 serves as a converging node in the cross-stress transcriptional expression network [48].
In summary, the relatively low values of Fv/Fm and Fv'/Fm' under heat stress reveal that the electron transfer at acceptor side of PSII was blocked. Thus, the heat-induced damage of photosynthetic apparatus was relatively lower in R. obtusum "Yanzhimi" leaves (Fig 2). A series of heat responsive unigenes was characterized in our transcriptome. These unigenes include some important TFs that show significant fold changes under heat stress and are involved in the transcriptional regulation of photosynthetic structure and protective genes. These results are consistent with results of previous studies on grapevine and spinach [10,31]. Thus, this study provides the candidate unigenes and pathways involved in the heat tolerance, thereby deepening our understanding of the molecular and physiological processes. This study also provided a large amount of basic genomic data for further research.