Genome-Wide Identification and Characterization of Heat-Shock Transcription Factors in Rubber Tree

Heat-shock transcription factors (Hsfs) play a pivotal role in the response of plants to various stresses. The present study aimed to characterize the Hsf genes in the rubber tree, a primary global source of natural rubber. In this study, 30 Hsf genes were identified in the rubber tree using genome-wide analysis. They possessed a structurally conserved DNA-binding domain and an oligomerization domain. On the basis of the length of the insert region between HR-A and HR-B in the oligomerization domain, the 30 members were clustered into three classes, Classes A (18), B (10), and C (2). Members within the same class shared highly conserved gene structures and protein motifs. The background expression levels of 11 genes in cold-tolerant rubber-tree clone 93-14 were significantly higher than those in cold-sensitive rubber-tree clone Reken501, while four genes exhibited inverse expression patterns. Upon cold stress, 20 genes were significantly upregulated in 93-114. Of the upregulated genes, HbHsfA2b, HbHsfA3a, and HbHsfA7a were also significantly upregulated in three other cold-tolerant rubber-tree clones at one or more time intervals upon cold stress. Their nuclear localization was verified, and the protein–protein interaction network was predicted. This study provides a basis for dissecting Hsf function in the enhanced cold tolerance of the rubber tree.


Introduction
Heat-shock transcription factors (Hsfs) allow for abiotic-stress signal transduction via the transcriptional activation of stress-related genes [1]. The highly conserved DNA-binding domain (DBD), which is responsible for activating the expression of heat-shock protein (HSP) genes by binding to heat-shock elements (HSEs), is located in the N-terminal of all Hsfs. It is followed by an oligomerization domain (OD) that is composed of two hydrophobic heptad repeats (HR-A and HR-B) and connected to the DBD by a flexible linker of variable length (15-80 amino acid residues). On the basis of peculiarities of the HR-A/B region and the phylogenetic relationship, plant Hsfs are divided into three main classes: A, B, and C [1][2][3]. In addition, certain Hsfs possess a nuclear localization signal (NLS), nuclear export signal (NES), and C-terminal activation (AHA) domain [2,4].
The rubber tree (Hevea brasiliensis Muell. Arg.) is a primary global source of natural rubber (NR). As a species indigenous to the Amazon Basin in South America, cold stress adversely affects its survival and growth. To meet the demand for NR, rubber-tree planting has been continuously expanded and is no longer limited to areas without cold stress. For example, rubber-tree plantations in China frequently suffer from sudden cold waves. Decades of selective breeding have resulted in a gradual improvement in the cold tolerance of rubber trees. Several cold-tolerant rubber-tree clones, such as 93-114, Zhanshi327-13, GT1, and Guiyan73-165, have been selected. Hsfs may be involved in enhanced cold tolerance. However, information regarding Hsfs is limited [16,17] despite four versions of the rubber-tree genome having been published [18][19][20][21].
In the present study, 30 HbHsf genes were identified from the rubber-tree genome published by Tang et al. [21] and isolated from rubber-tree clone 93-114. We then analyzed evolutionary relationships, gene structure, and protein domains, and performed expression analysis using real-time quantitative RT-PCR in cold-tolerant rubber-tree clone 93-114 and cold-sensitive rubber-tree clone Reken501. The expression of several Hsfs with differential expression between 93-114 and Reken501 was further analyzed in three other cold-tolerant and cold-sensitive rubber-tree clones. As a result, three HbHsfs were suggested to be involved in the enhanced cold tolerance of rubber-tree clones. Their nuclear localization was verified, and the protein-protein interaction network was predicted. The present work provides a basis for dissecting Hsf function in the enhanced cold tolerance of rubber tree.

Plant Materials and Growth Condition
Plantlets of 4 cold-tolerant rubber-tree clones (93-114, Zhanshi327-13, GT1, and Guiyan73-165) and 4 cold-sensitive rubber-tree clones (Reken501, Haiken1, Reken514, and Reken515) budded on rootstocks were grown on the Experimental Farm of the Chinese Academy of Tropical Agricultural Sciences (CATAS) on Hainan Island. When plantlets had developed 1 extension unit and leaves were completely matured in the field, they were transferred into a growth chamber (PGR15, CONVIRON, Canada) with light for 16 h at 125 µmol·m −2 ·s −1 , temperature of 28 • C, and relative humidity of 80% for 2 days. Thereafter, plantlets were transferred into another chamber under the same conditions at 4 • C instead of 28 • C.

Stress Treatments
To evaluate stress intensity, plantlets of cold-tolerant rubber-tree clone 93-114 and cold-sensitive rubber-tree clone Reken501 were treated for 2 days at 4 • C, and then recultured for 5 days at 28 • C. The length of withered stems was surveyed once a day. The extent of cold tolerance was represented by the ratio of withered length to the total length of the stem (%). Each variety included 5 plantlets.
The same stress intensity as above was performed for gene-expression analysis in 3 cold-tolerant rubber-tree clones and 3 cold-sensitive rubber-tree clones, in addition to 93-114 and Reken501. Each clone included 120 plantlets and remained at 28 • C for 2 days. Thereafter, plantlets were divided into 2 groups. Each group included 60 plantlets. One group was treated at 4 • C, while the other remained at 28 • C as a control. Bark samples were collected at 0, 4, 8, and 24 h. At each time interval, bark-tissue samples were collected from the stem of 15 plantlets. Bark tissue from every 5 plantlets was mixed as 1 sample. This was done in triplicate.

Real-Time Quantitative PCR Analysis
Total RNA was extracted using a plant RNA extraction kit (TIANGEN, Beijing, China). Approximately 1 µg of RNA was used for reverse transcription, and cDNA was synthesized using a RevertAid™ First Strand cDNA Synthesis Kit (Thermo Scientific Inc., Waltham, MA, USA). Primer Premier 5.0 (Premier Biosoft International, Palo Alto, CA, USA) was used to design gene-specific primers for cloning the full-length cDNA of each HbHsf gene (Supplementary Table S1). qRT-PCR experiments were carried out in triplicate with the CFX384 real-time PCR system (Bio-Rad Laboratories, Inc., v. Bio-Oxford, USA) using an SYBR Prime Script RT-PCR Kit (TaKaRa, Dalian). Relative gene expression was determined using Actin7a as an endogenous control gene [22]. Relative expression levels were calculated using the 2 −∆∆CT method. Primers for qRT-PCR are listed in Supplementary  Table S2.

Bioinformatics Analysis
Hsf genes were predicted from the whole genome sequence of Hevea brasiliensis (https://www. ncbi.nlm.nih.gov/genome/?term=rubber%20tree%20genome) [21]. All the nucleotide and amino acid sequences of predicted HbHsf genes were used as queries to perform BLAST searches in the public NCBI database (http://www.ncbi.nlm.nih.gov/). If 2 or more protein sequences at the same gene locus overlapped, only the longest sequence was used. Multiple sequence alignment of Hsf proteins from the rubber tree was carried out using the DNAMAN program (Lynnon Biosoft, Vaudreuil, QC, Canada). Full-length amino acid sequences of Hsf proteins from the rubber tree (HbHsf), Arabidopsis (AtHsf), and Populus trichocarpa (PtHsf) were used to generate a phylogenetic tree through MEGA5.10 [12]. Firstly, the alignment was built and aligned by ClustalW, then analyzed using the unrooted neighbor-joining method with pairwise deletion. Bootstrap value analysis was performed with 1000 replications to assess tree reliability. Deduced amino acid sequences were analyzed using the Compute pI/MW Tool (http://www.expasy.org/tools/pi_tool.html) for computation of the theoretical isoelectric point and protein molecular weight. Exon-intron structures of HbHsf genes were visualized using online program GSDS 2.0 (http://gsds.cbi.pku.edu.cn/), which aligned the respective coding sequences with corresponding full-length sequences. The protein domains of the HbHsf proteins were identified using MEME online tools (http://meme.nbcr.net/meme/).

Subcellular-Localization Analysis
For subcellular-localization assays, the entire coding sequences of HbHsfA2b, HbHsfA3a, and HbHsfA7a with no stop codons were amplified with specific primers (Supplementary Table S3). HbHsfA2b, HbHsfA3a, and HbHsfA7a were fused with GFP under the control of the 35S promoter in the pCAMBIA1302 vector. HbHsfA2b::GFP, HbHsfA3a::GFP, and HbHsfA7a::GFP fusion constructs, and empty GFP vectors were introduced into living onion epidermal cells by Agrobacterium (GV3101)-mediated transformation [23]. The expression of fusion constructs was measured after 24 h incubation in the dark by a confocal laser scanning microscope (Zeiss LSM 800, Jena Germany). DAPI was used to stain the nucleus.

Predicted Protein-Protein Interaction-Network Construction
The predicted protein-protein interaction network of HbHsfA2b, HbHsfA3a, and HbHsfA7a was generated by STRINGV10.0 (University of Zurich, Zurich, Switzerlant) online software, which visualizes known and predicted protein-protein interactions. Parameters were set as follows: Meaning of network edges-confidence (line thickness indicates strength of data support); active interaction sources-check all; minimum required interaction score-high confidence (0.700); and max number of interactors to show-no more than 10 interactors [24].

Statistical Analysis
Data were shown as the mean ± standard deviation and tested for statistical significance with a paired Student's t-test. p < 0.05 was selected as the point of minimal statistical significance in all analysis.

Isolation and Identification of Hsf Genes in Rubber Tree
A total of 30 Hsf genes were predicted in the sequenced genome of rubber-tree clone CATAS 7-33-97 through the BLAST search in NCBI. They contained apparently complete Hsf-type DNA-binding domains and oligomerization domains. Using predicted HbHsfs coding sequences, 30 homologous genes were isolated from rubber-tree clone 93-114. Identified HbHsfs were designated on the basis of names of their presumptive Arabidopsis and P. trichocarpa orthologs (Figures 1 and 2A). ORF sequences of isolated HbHsfs shared identity with the corresponding HbHsfs from rubber-tree clone CATAS 7-33-97. HbHsf genes contained one or two introns. Introns ranged from 66 (HbHsfB4d) to 5347 bp (HbHsfA1b) in length. HbHsfs ranged from 216 (HbHsfB2a) to 566 (HbHsfA3a) amino acid residues in length, and predicted isoelectric points (pI) varied from 4.73 (HbHsfA3a) to 9.36 (HbHsfB2a), and molecular weights (MW) from 24 Table 4). 153  Exon-intron organization analysis of 30 HbHsf genes revealed that 28 genes possessed one 157 intron, while two genes (HbHsfB4d and HbHsfA9b) had two introns. Introns were inserted into the 158 highly conserved DNA-binding domain-coding regions, and most of the introns were phase-zero 159 introns, except for phase-one introns HbHsfA9b and HbHsfA2b ( Figure 1A). 160  subgroups of Arabidopsis Hsfs and P. trichocarpa Hsfs. All HbHsfs were divided into three classes 171 (A, B, and C), which were characterized by the insertion length between HR-A and HR-B regions. 172 In comparison with C, insertions in B and A contained seven and 21 amino acid residues, 173 respectively ( Figure 2B). Class A included 18 HbHsf members, Class B had 10 members, and Class 174 C only contained two members ( Figure 2).

Expression Patterns of Hbhsfs in Cold-Tolerant and Cold-Sensitive Rubber-Tree Clones in Response to 182
Cold Stress 183 The difference in the cold-tolerance phenotype between cold-tolerant rubber-tree clone 93-114 184 and cold-sensitive rubber-tree clone Reken501 was remarkable when plantlets were treated at 4 °C 185 for two days, and then recultured at 28 °C for five days. The leaves and stems of Reken501 began to 186 shrivel after being recultured for one day. Stems of most Reken501 plantlets had withered after 187 being recultured for five days, while stems of a few 93-114 plantlets exhibited withering ( Figure 3A).

188
The stem-wither ratio was 82.77% for Reken501 and 6.06% for 93-114 ( Figure 3B). By using this Exon-intron organization analysis of 30 HbHsf genes revealed that 28 genes possessed one intron, while two genes (HbHsfB4d and HbHsfA9b) had two introns. Introns were inserted into the highly conserved DNA-binding domain-coding regions, and most of the introns were phase-zero introns, except for phase-one introns HbHsfA9b and HbHsfA2b ( Figure 1A).
Ten motifs were detected in HbHsf proteins using MEME software ( Figure 1B). Prediction by the InterProscan database annotated Motifs 1, 2, and 5 as DBDs, Motif 3 as an HR-A/B domain, Motif 4 as an HR-A/B motif, Motif 7 as an AHA motif, Motif 8 as an NES motif, and Motif 9 as an NES motif. Motifs 6 and 10 were not annotated in the database. Motifs 1 and 2 were present in all of the HbHsfs members. Motif 5 was found in most of the HbHsf members except for HbHsfC1a, HbHsfC1b, and HbHsfB4b. Motif 4 was only present in HbHsfBs, while Motif 3 was absent from HbHsfBs except for HbHsfB3a. Motif 7 was only found in eight HbHsfAs. Motif 8 was present in seven HbHsfAs, and Motif 9 was absent from HbHsfBs except for HbHsfB4b.
Thirty HbHsfs, 21 Arabidopsis Hsfs (AtHsfs), and 28 P. trichocarpa Hsfs (PtHsf) were used to construct a phylogenetic tree (Figure 2A). Thirty HbHsf proteins were grouped into different subgroups of Arabidopsis Hsfs and P. trichocarpa Hsfs. All HbHsfs were divided into three classes (A, B, and C), which were characterized by the insertion length between HR-A and HR-B regions. In comparison with C, insertions in B and A contained seven and 21 amino acid residues, respectively ( Figure 2B). Class A included 18 HbHsf members, Class B had 10 members, and Class C only contained two members ( Figure 2).

Expression Patterns of Hbhsfs in Cold-Tolerant and Cold-Sensitive Rubber-Tree Clones in Response to Cold Stress
The difference in the cold-tolerance phenotype between cold-tolerant rubber-tree clone 93-114 and cold-sensitive rubber-tree clone Reken501 was remarkable when plantlets were treated at 4 • C for two days, and then recultured at 28 • C for five days. The leaves and stems of Reken501 began to shrivel after being recultured for one day. Stems of most Reken501 plantlets had withered after being recultured for five days, while stems of a few 93-114 plantlets exhibited withering ( Figure 3A). The stem-wither ratio was 82.77% for Reken501 and 6.06% for 93-114 ( Figure 3B). By using this experimental system, expression patterns of 30 HbHsfs in the two clones in response to cold stress were analyzed (Figure 4). Background expression levels of 11 genes in the bark tissue of 93-114 were significantly higher than those of Reken501, while the reverse was true for HbHsfA9a, HbHsfA9b, and HbHsfB2c. Upon cold stress, 20 of 30 HbHsfs were upregulated at 4 h and/or 8 h in 93-114. By contrast, only four genes were upregulated at corresponding time intervals in Reken501. Three genes-HbHsfA2b, HbHsfA4a, and HbHsfB4a-were significantly upregulated in 93-114, while HbHsfA9b was upregulated in Reken501 during cold stress. Expression patterns of HbHsfC1a and HbHsfC1b were similar between 93-114 and Reken501. The expression of all four HbHsfB4 members was higher in 93-114 than in Reken501.  HbHsfs that were upregulated in both 93-114 and Reken501, and those that were significantly upregulated at two or more time intervals upon cold stress in 93-114, were selected for further analysis of their expression patterns in three other cold-tolerant and cold-sensitive rubber-tree clones. The same experimental system was used. They included HbHsfA1b, HbHsfA2b, HbHsfA3a, HbHsfA4a, HbHsfA7a, and HbHsfB4a (FIguer 5). Of the six genes, HbHsfA2b, HbHsfA3a, and HbHsfA7a were significantly upregulated at one or more time intervals in all three cold-tolerant rubber-tree clones in comparison with the three cold-sensitive rubber-tree clones. Background expression of HbHsfA2b was significantly upregulated in the three cold-tolerant rubber-tree clones, but only in the most cold-tolerant rubber-tree clone (GT1) during cold stress ( Figure 5B). Background expression of HbHsfA3a was significantly upregulated and then further upregulated during cold stress in cold-tolerant rubber-tree clones Zhanshi327-13 and GT1. It was significantly upregulated in Guiyan73-165 at 24 h during cold stress in comparison with the other cold-sensitive rubber-tree clones ( Figure 5C). Background expression of HbHsfA7a was significantly upregulated in the three cold-tolerant rubber-tree clones. It was also upregulated during cold stress, except for at 4 h, where there was no difference between cold-tolerant rubber-tree clone Guiyan73-165 and the three cold-sensitive rubber-tree clones ( Figure 5E). during cold stress in cold-tolerant rubber-tree clones Zhanshi327-13 and GT1. It was significantly 216 upregulated in Guiyan73-165 at 24 h during cold stress in comparison with the other cold-sensitive 217 rubber-tree clones ( Figure 5C). Background expression of HbHsfA7a was significantly upregulated 218 in the three cold-tolerant rubber-tree clones. It was also upregulated during cold stress, except for at 219 4 h, where there was no difference between cold-tolerant rubber-tree clone Guiyan73-165 and the 220 three cold-sensitive rubber-tree clones ( Figure 5E).   Figure 6A). The protein-protein interaction 235 network of HbHsfA2b, HbHsfA3a, and HbHsfA7a was predicted. Ten proteins were involved in the 236 interaction with the three HbHsfAs. They were DREB2A, ROF1, HSP101, HSP90.1, HSF1, MBF1C, 237 HSP70, HSF3, HSP70b, and HSPB ( Figure 6B). 238  HbHsfA2b (B), HbHsfA3a (C), HbHsfA4a (D), HbHsfA7a (E) and HbHsfB4a (F) genes in three cold-sensitive rubber-tree clones (Haiken1, Reken514, and Reken515), and three cold-tolerant rubber-tree clones (GT1, Zhanshi327-13, and Guiyan73-165) under cold treatment. ** and *: very significant difference (p < 0.01) and significant difference (p < 0.05), respectively.

Discussion
Plant Hsfs are terminal components of a signal-transduction chain mediating the expression of abiotic-stress-responsive genes [1]. The number of Hsf members is 18 in plums [15], 21 in Arabidopsis [5], 25 in rice [6], 25 in apples [8], and 32 in cassavas [25]; these plants are diplontic. The number of Hsf members increases in allopolyploids-there are 43 in triploid banana plants [13], 59 in tetraploid soybean plants [9], and 56 in sextaploid wheat plants [10]. On the basis of the version of the rubber-tree genome [21] and verification by RT-PCR, 30 HbHsf genes were isolated in the present study. The number of HbHsf genes falls in the range of diploid plants. Expansion of the Hsf gene family seems to be the result of the allopolyploid process. HbHsfs are structurally similar to Hsfs from other plants, in that DBD contains an intron insertion, HR-A/B domains are highly homologous to those of the corresponding Class A Hsf members from Arabidopsis and P. trichocarpa, and the number of inserted amino acids between the HR-A and HR-B domains is the same as that in Arabidopsis [5], rape [14], and cotton [26]. This structural conservation suggests that the functions of Hsf in different plants are conserved.
During evolution, gene-family expansion is implemented through tandem duplication and segmental duplication. Segmental-duplication genes with different functions are on different chromosomes or scaffolds. Tandem-duplication genes that have similar functions are located on the same chromosome, and the distance between the two genes is relatively close [27,28]. In P. trichocarpa, only gene pair PtHsfA6a/A6b is from tandem duplication, while other clusters are from segmental duplications [29]. In peanuts, cluster AiHsf7/8 is located on chromosome five, and the distance between the two genes is about 2 Kb. Gene pair AiHsf7/8 is from tandem duplication [27]. In the rubber tree, due to the lack of a high-density genetic map, identification of paralogs that are derived from Whole Genome Duplication (WGD) or segmental duplication has been largely restricted. On the basis of identification of tandem-duplication events in the rubber tree by Zou et al. [28], only gene pair HbHsfC1a/C1b was located on same scaffold 0795, and the distance between the two genes was about 2 Kb. Other clusters, such as HbHsfA1a/A1c, HbHsfA2a/A2b, HbHsfA3a/A3b, HbHsfA4a/A4b, HbHsfA5a/A5b, HbHsfA9a/A9b, HbHsfB3a/B3b, and HbHsfB4a/B4c were located on different scaffolds (Table S4). Therefore, gene pair HbHsfC1a/C1b was from tandem duplication, whereas other gene pairs were from segmental duplications. Compared to tandem-duplicated gene pair HbHsfC1a/C1b, some segmentally duplicated gene pairs had differential expression patterns in response to cold and/or heat stress (Figure 4 and Figure S1), suggesting these duplicated genes display functional diversity.
Although the rubber tree is tropical, the cold tolerance of cultivars has been gradually improved by selective breeding. Among 30 HbHsfs, HbHsfA2b, HbHsfA3a, and HbHsfA7a were significantly upregulated at one or more time intervals in four cold-tolerant rubber-tree clones in comparison with four cold-sensitive rubber-tree clones (Figures 4 and 5). Their corresponding homologs OsHsfA2 [6], MeHsfA3a [30], and OsHsfA7 [7] are upregulated upon cold stress. OsHsfA2 and OsHsfA7 genes are also induced by high temperature. In P. trichocarpa, transcription levels of PtHsfA2 and PtHsfA7a/b are increased under high-temperature treatment [29]. In rubber-tree clone 93-114, 20 of 30 HbHsfs were upregulated under cold stress (Figure 4). Of the 20 genes, the expression levels of 14 HbHsf genes, including HbHsfA2b, HbHsfA3a, and HbHsfA7a, were also upregulated upon heat stress (42 • C; Figure S1).
Hsfs act as key components in regulating the enhanced production of reactive-oxygen-species (ROS) scavengers and HSPs that facilitate plant stress tolerance [1]. OsHsfA2a/c/f, AtHsfA2a, and OsHsfA7 are induced by oxidative stress [6,7,31]. APX2, a stress-related gene in Arabidopsis, is the target gene of HsfA2 [2,32]. Ectopic expression of CtHsfA2b in Arabidopsis enhances the transcriptional activity of AtApx2 by directly binding to the HSE on the promoter of AtApx2 [33]. AtHsfA2 knockout lines and HsfA2 overexpression lines indicate that HsfA2 is one of the key regulators in protecting organelles against oxidative damage [32,34]. HbHsfA2b, HbHsfA3a, and HbHsfA7a may be involved in regulating the activity of ROS scavengers, because the activity of peroxidase, catalase, and superoxide dismutase in the leaves of cold-tolerant rubber-tree clone 93-114 was significantly higher than that in cold-sensitive rubber-tree clone Reken501 during cold stress [35]. In addition, available data showed that heat-shock proteins are associated with chilling tolerance in many plants and are transcriptionally regulated by Hsfs [36]. The interaction of Hsfs with other proteins is essential for activating their target stress-related genes, including HSPs and ROS scavengers [37]. Interactions of HbHsfA2b, HbHsfA3a, and HbHsfA7a with several HSPs, Hsfs, and DREB2A were predicted ( Figure 6B). In Arabidopsis, activation of sHSP by HsfA2 depends on the formation of an ROF1-HSP90.1-HsfA2 complex by interactions [38]. The predicted protein-protein interaction network showed that HbHsfA2b could interact with ROF1 and HSP90.1 ( Figure 6B). Moreover, some of the FvHsf proteins of Fragaria vesca were reported to localize to both the nucleus and cytosol, such as FvHsfA2a, FvHsfA3a, FvHsfA4a, FvHsfA5a, FvHsfB2b, and FvHsfC1a [12]. While HbHsfA2b was only detected in the nucleus ( Figure 6A), the expression of a number of sHSPs was significantly activated upon cold stress in cold-tolerant rubber-tree clone 93-114 in comparison with cold-sensitive rubber-tree clone Reken501 [16]. These data suggest that HbHsfA2b may be involved in activating some sHSPs in the rubber tree in a similar manner to HsfA2 in Arabidopsis. It would be worthwhile to elucidate the interaction between HbHsfA2b, HbHsfA3a, and HbHsfA7a, their interaction with proteins such as HSPs and DREBs, and the identity of their target genes.

Conclusions
A total of 30 HbHsf genes were identified in the rubber tree. They were structurally conserved and may function in a similar manner to corresponding members in other plants. Most of the HbHsfs were significantly upregulated in cold-tolerant rubber-tree clones upon cold stress. Specifically, the upregulation of HbHsfA2b, HbHsfA3a, and HbHsfA7a may contribute to the enhanced cold tolerance of the rubber tree. The present work provides a basis for dissecting HbHsf function in the enhanced cold tolerance of the rubber tree.

Conflicts of Interest:
The authors declare no conflict of interest.