Genome-wide characterization of tea plant (Camellia sinensis) Hsf transcription factor family and role of CsHsfA2 in heat tolerance

Background Heat stress factors (Hsfs) play vital roles in signal transduction pathways operating in responses to environmental stresses. However, Hsf gene family has not been thoroughly explored in tea plant (Camellia sinensis L.). Results In this study, we identified 25 CsHsf genes in C. sinensis that were separated by phylogenetic analysis into three sub-families (i.e., A, B, and C). Gene structures, conserved domains and motifs analyses indicated that the CsHsf members in each class were relatively conserved. Various cis-acting elements involved in plant growth regulation, hormone responses, stress responses, and light responses were located in the promoter regions of CsHsfs. Furthermore, degradome sequencing analysis revealed that 7 CsHsfs could be targeted by 9 miRNAs. The expression pattern of each CsHsf gene was significantly different in eight tissues. Many CsHsfs were differentially regulated by drought, salt, and heat stresses, as well as exogenous abscisic acid (ABA) and Ca2+. In addition, CsHsfA2 was located in the nucleus. Heterologous expression of CsHsfA2 improved thermotolerance in transgenic yeast, suggesting its potential role in the regulation of heat stress response. Conclusions A comprehensive genome-wide analysis of Hsf in C. sinensis present the global identification and functional prediction of CsHsfs. Most of them were implicated in a complex gene regulatory network controlling various abiotic stress responses and signal transduction pathways in tea plants. Additionally, heterologous expression of CsHsfA2 increased thermotolerance of transgenic yeast. These findings provide new insights into the functional divergence of CsHsfs and a basis for further research on CsHsfs functions.

, Triticum aestivum [15], Salix suchowensis [16], and five model angiosperms (Arabidopsis thaliana, Cucumis sativus, Oryza sativa, Populus trichocarpa, and Vitis vinifera) [17,18], no systematic identification of Hsf family is available in tea plant. The elucidation of CsHsfs function and regulation in tea plant will provide foundation for further functional research of Hsf genes in evergreen woody crops.
Hsf transcription factors play vital roles in the regulation of plant response to abiotic stresses, such as salinity, drought, osmotic stress, cold, and HS [9,19,20]. For example, HsfA2, a typical representative of plant Hsfs, was a HS-inducible gene and regulated a subset of down-stream stress response genes [19]. Overexpression of HsfA2 from A. thaliana, Zea mays, Lilium longiflorum, or O. sativa, conferred heat tolerance in transgenic Arabidopsis [21][22][23][24]. GmHsfA1 overexpression also enhanced thermotolerance of transgenic soybean under HS [25]. The drought and heat stress induced expression of HsfA3 in Arabidopsis depends on the dehydration-responsive element-binding protein 2A (DREB2A) [26,27]. Although O. sativa HsfB2b expression was strongly induced by heat, salt, abscisic acid (ABA) and polyethylene glycol (PEG) treatments but not cold stress, OsHsfB2b overexpression increased the drought and salt sensitivity in rice [28]. However, Hordeum vulgare HsfB2c was proposed to be positively regulating the drought stress tolerance in barley [29], and CmHSFA4 and AtHSFA7b positively regulated salt stress tolerance in chrysanthemum and Arabidopsis, respectively [30,31]. Moreover, OsHsfC1b and T. aestivum HsfC2a overexpression improved salt tolerance and thermotolerance in transgenic O. sativa and wheat, respectively [32,33]. In addition to its role in stress responses, it is interesting to note that overexpression of sunflower HsfA9 improved seed longevity in transgenic tobacco [34]. Thus, the unique functions of specific Hsf from different plant species remain to be elucidated.
Tea plants [Camellia sinensis (L.) O. Kuntze] are one of the most important commercial perennial evergreen leaf crops used for the production of non-alcoholic beverages in china and worldwide [35]. Most of tea plant species are diploids with a genome size of thirty chromosomes [36]. Being sessile organisms, tea plants have evolved a series of complex mechanisms to cope with fluctuating environmental stresses, such as extreme temperatures [37], drought [38,39], soil acidification [40,41], and heavy metals [42,43]. Long-term hot summers hinder tea plants growth and development, seriously affecting the yield and quality of spring tea in the next year [44]. Hence, it is necessary to investigate the thermotolerance mechanisms of tea plants, which may be useful for HS-resistant cultivation and breeding in the future. We previously have attempted to discover the key HS-responsive genes using a suppression subtractive hybridization approach [45]. However, due to the technical limitations of suppression subtractive hybridization and unavailability of tea plant genome at that time, only 12 differentially expressed heat-responsive genes were identified, including one Hsf gene [45]. Although 16 CsHsfs have been identified from C. sinensis RNA-seq data [46], little information on the genome-wide identification of CsHsf genes is available. Given the potential value of Hsfs in resistance to HS, it is necessary to identify Hsf genes in tea plant genome. In this study, we initiated the characterization of C. sinensis Hsf gene family and the expression patterns of putative CsHsfs in responses to a variety of stress treatments. We also examined the heat resistance function of CsHsfA2 using a yeast heterologous expression system. These findings provide new insights into the functional divergence of CsHsf genes and a basis for genetic engineering breeding of plant species.

Results
Identification of Hsf gene family in C. sinensis A total of 25 putative Hsf genes were identified from C. sinensis cultivar 'Shuchazao' genome. Members of the CsHsf gene family were subdivided into classes A, B, and C according to differences in the length of the flexible linkers between the A and B parts of the HR-A/B regions, including 15 CsHsfA genes, 9 CsHsfB genes, and 1 CsHsfC genes ( Table 1). The length of the CsHsfs coding sequences ranged from 624 bp (CsHsfA5a) to 1998 bp (CsHsfA4b) with 207-665 amino acid residues. The molecular masses of CsHsfs varied greatly, ranging from 23.56 kDa (CsHsfA5a) to 74.98 kDa (CsHsfA4b), the computed theoretical isoelectric points ranged from 4.71 (CsHsfA9b) to 10.01 (CsHsfB3c), and the grand average of hydropathicity (GRAVY) was between − 1.01 (CsHsfB3c) and − 0.34 (CsHsfA9c). Additionally, subcellular localization predictions indicated that most of CsHsf proteins were predicted to be located in the nucleus, while CsHsfA5a and CsHsfA8 were targeted to cytoplasm and chloroplast, respectively. The detailed information is listed in Table 1.

Identification of conserved motifs and domains in CsHsfs
A total of 25 conserved motifs were identified from CsHsfs using MEME [49]. Motifs 1 and 3 were found in all of the CsHsfs, while motif 13 only existed in class A CsHsfs (Additional file 1: Figure S1). In addition, the number of motifs in the CsHsfBs and CsHsfCs was less than those in the CsHsfAs, especially in CsHsfB3c, which harbored only two motifs.
To better understand the structural characteristics of the CsHsf family, six conserved domains including DND-binding domain (DBD), oligomerization domain (HR-A/B), nuclear localization signal (NLS), nuclear export signal (NES), activator motifs (AHA), and repressor domain (RD) were identified using SMART [50], Pfam [51], NLStradamus [52] and NetNES [53] (Table 2). DBD is the most conserved domain, comprised of three α-helices and four β-sheets in the form of α1-β1-β2-α2-α3-β3-β4 (Fig. 1). However, no intact α3, β3, and β4 were detected in the DBD domain of CsHsfA5a, CsHsfA9c, CsHsfA9d, and CsHsfB4b, which may result in their shorter sequences compared to other CsHsfs. HR-A/ B, which is characterized by a coiled-coil structure (coilcoil structure), was also present in most of the CsHsfs except for CsHsfA5a, CsHsfA5b, CsHsfA8, CsHsfA9c, CsHsfB3c, CsHsfB4a, and CsHsfB4b (Table 2). Moreover, NES and NLS domains were detected in most of the CsHsfs, which are vital for shuttling CsHsfs between the nucleus and cytoplasm. In addition, AHA and RD domains were specific to each class. AHA motifs were only detected in class A CsHsfs (except for CsHsfA5a and CsHsfA5b), and four proteins in subclasses A2, A3, and A4 had two AHA motifs. The tetrapeptide motif LFGV, as the core of the RD, was identified in all CsHsfB members except for CsHsfB3c. Interestingly, only one DBD domain was identified in CsHsfA5a and CsHsfB3c.

Phylogenetic analysis of CsHsfs
To evaluate the evolutionary relationships among the Hsf family, a phylogenetic tree was generated using the Hsf sequences of 25 proteins from C. sinensis, 21 from A. thaliana, and 30 from P. trichocarpa (Fig. 2). According to this tree, the CsHsfs could be divided into three main classes (i.e., A, B, and C). Class A included 15 members from nine subclasses (i.e., A1-A9), Class B was divided into subclasses B1, B2, B3, and B4, while class C contains only one member (i.e., CsHsfC1). Interestingly, two P. trichocarpa Hsfs (i.e., PtHsfA1c and PtHsfB2b) were not clustered with their corresponding subclasses, but were closer to the HsfA3 subclass. In addition, compared with A. thaliana, P. trichocarpa Hsfs were closer to C. sinensis Hsf proteins.   N.D.

Gene structures and cis-elements analyses
The structures of CsHsf genes were analyzed by comparing their cDNA sequences and genomic DNA sequences. Generally, most of the CsHsf genes contained one or two introns, while CsHsfA4b and CsHsfB1 were comprised of five and eight introns, respectively (Fig. 3).
To further investigate the potential regulatory networks of the CsHsf family genes, the cis-elements in the 2 kb upstream sequences of the translation initiating site of 25 CsHsf genes were analyzed using PlantCARE [54] (Additional file 2: Table S1). A total of 37 types of cis-elements were discovered, including 4 plant growth   (Fig. 4). Among the plant growth-related cis-acting elements, 4 CsHsfs had O2-site and CAT-box elements, 5 CsHsfs contained GCN4 motif, and 2 CsHsfs had circadian related elements. In the category of hormone responsiveness, a total of 19 abscisic acid responsive elements (ABRE), 17 MeJA-responsive elements (CGTCA-motif and TGACG-motif), 8 salicylic acid-responsive elements (TCAelement), 5 auxin-responsive elements (TGA-element), and 17 gibberellin-responsive elements (GARE-motif, TATCbox, and P-box) were detected in the promoter regions of 8, 8, 7, 4, and 11 CsHsfs, respectively. Moreover, many stress-responsive cis-elements were detected, including ARE (anaerobic induction element), LTR (low-temperature responsiveness), TC-rich repeats (defense and stress responsiveness), GC-motif (anoxic specific inducibility), and MBS (drought-inducibility). The light-responsive elements accounted for the largest part of all cis-elements, especially the Box 4 element, with a total of 48 distributed in 25 CsHsf promoter regions. In addition, G-box (23), AE-box (11), and GT1-motif (10) are involved in light responsiveness. Overall, these results indicated the involvement of CsHsfs in responses to hormone treatments, low temperature, and drought stresses in tea plants. Interestingly, no HSE element was detected in these CsHsf promoter regions. Hence, it remains unknown whether the expresssion of CsHsf genes was regulated by heat stress conditions.

MiRNA target sites prediction
Tea plant degradome libraries (unpublished data) were used to predict target transcript candidates of CsHsfs. As shown in Additional file 3: Table S2, 7 CsHsfs are predicted to be targeted by 9 miRNAs. CsHsfB3c has 4 target sites, CsHsfA4a and CsHsfA4b have two target sites, while the other four CsHsfs (i.e., CsHsfA1a, CsHsfA1b, CsHsfB1, and CsHsfA9d) have only one target site. Hence, we inferred that a single miRNA can regulate multiple CsHsfs and a single CsHsf can be regulated by multiple miRNAs.

Expression profiles of CsHsf genes in different tissues
To examine the expression patterns of CsHsf genes among eight tissues, a heat map was drawn based on the transcriptome data downloaded from tea plant genome database [55]. The expression pattern of each CsHsf gene was significantly different in eight tissues ( Fig. 5; Additional file 4: Table S3). CsHsfA3 was highly expressed in old leaf (values > 2), CsHsfB2b was highly expressed in fruit and young leaf, and CsHsfC1 was highly expressed in root. In comparison, CsHsfA6, CsHsfB3c, and CsHsfB4a were only expressed in some specific tissues. However, the expression of CsHsfA9b, CsHsfA9c, CsHsfA9d was undetectable in all tissues.

Expression profiles of CsHsf genes in responses to drought and salt treatments
To examine the roles of CsHsf genes in responses to drought and salt stresses, the transcriptome data of these genes were analyzed. Genes with the changes in transcription level greater than 1.5-fold were considered to be significantly regulated. After exposure to drought and salt stresses, most CsHsfs were up-regulated, whereas CsHsfA1a and CsHsfB4a were down-regulated ( Fig. 6; Additional file 5: Table S4). CsHsfB3a was significantly up-regulated by drought stress (Fig. 6a), and CsHsfA3, CsHsfA7, CsHsfB3a, CsHsfB3b, CsHsfB3c, and CsHsfB4b were significantly up-regulated by salt stress (Fig. 6b). Interestingly, the expression level of CsHsfA5a, CsHsfA7, CsHsfA9a, CsHsfB2c, and CsHsfB4b showed the opposite trends under two abiotic stresses, indicating that they may have different roles in responses to drought and salt stresses in tea plants.

Expression profiles of CsHsf genes with exogenous calcium treatments
To explore whether calcium ions (Ca 2+ ) are involved in Hsf-mediating response to heat stress, we analyzed expression profiles of CsHsf genes in tea plants foliar-sprayed with exogenous Ca 2+ . According to RNA sequencing data, the fold changes of CsHsfs greater than 1.3-fold were considered to be significantly regulated. When exposure to exogenous Ca 2+ , the transcript abundance of CsHsfA4a was significantly up-regulated; and the transcription levels of the other four CsHsfs (i.e., CsHsfA4b, CsHsfB1, CsHsfB2b, CsHsfB2c) were slightly increased ( Fig. 8a; Additional file 7: Table S6). Interestingly, the expression of the remaining CsHsfs was downregulated. These results suggested that Ca 2+ appeared to be involved in Hsf-mediating response to heat stress in tea plants. To further confirm the reliability of the RNAseq results, four up-regulated and four down-regulated CsHsf genes were tested by qRT-PCR using the same tea plant variety (Fig. 8b; Additional file 8: Table S7). The results showed that the expression of six selected CsHsfs were well correlated with RNA-seq data, while the expression of CsHsfA2 and CsHsfC1 showed the opposite trends with RNA-seq results.  [47]. HemI software [56] was used to generate the heatmap. The color bar represents the normalized transcript per million (TPM) values (log10-transformed fold-changes). Red and blue colors indicate up-and down-regulated genes and the gray represents no expression. Detailed TPM values are listed in Additional file 4: Table S3 CsHsfA2 localizes to the nucleus in onion epidermal cells CsHsfA2 up-regulated strongly by heat stress was selected for subcellular localization analysis. As shown in Fig. 9, the control GFP signal was uniformly distributed throughout the cytosol and nucleus in onion epidermal cells, whereas the diffuse CsHsfA2-GFP and GFP-CsHsfA2 signals were only detected in the nucleus. Thus, CsHsfA2 is a nuclear protein, possibly serving as a transcription factor.

Heterologous expression of CsHsfA2 confers thermotolerance in transgenic yeast
Because of its dominant role in thermotolerance [9] and strongly up-regulated expression, we constructed a yeast expression vector and transformed it into yeast cell to evaluate the possible roles of CsHsfA2 in response to heat stress. Under normal temperature conditions (30°C), there were no obvious differences in yeast cells expressing CsHsfA2 compared with the empty vector cells (i.e., pPIC3.5 K) (Fig. 10). However, when exposed to heat stress, the growth of yeast cells expressing CsHsfA2 was better than the control cells. These results suggested that heterologous expression of CsHsfA2 improved thermotolerance in transgenic yeast.

Discussion
Tea plant, a perennial evergreen woody crop, has to cope with various abiotic stress during its lifecycle [57,58].  6 Expression profiles of CsHsf genes in tea plants (C. sinensis cv. 'Tieguanyin') exposed to (a) drought (25% polyethylene glycol 6000) and (b) salt (200 mM NaCl) treatments. HemI software [56] was used to generate the heatmap. The heatmap of CsHsf expression based on the transcriptome data [47]. The color bar represents the normalized transcript per million (TPM) values (log10-transformed fold-changes). Red and blue colors indicate upand down-regulated genes and the gray represents no expression. Detailed TPM values are listed in Additional file 5: Table S4 Fig. 7 Expression profiles of CsHsf genes in tea plants (C. sinensis cv. 'Echa 10') exposed to (a) heat (38°C) and (b) ABA (50 μM) treatments. HemI software [56] was used to generate the heatmap. The heatmap of CsHsf expression based on the qRT-PCR data of three biological and technical replicates. The color bar represents the normalized expression level (log2-transformed fold-changes). Red and blue colors indicate up-and downregulated genes and the gray represents no expression. qRT-PCR data are listed in Additional file 6: Table S5 Previous studies have showed that Hsf family genes play vital roles in responses to abiotic stress, especially high temperature stress [59]. Hence, it is necessary to investigate the Hsf family in tea plant. C. sinensis Hsf genes family were firstly identified according to the RNA-seq data by Liu et al. [46]. However, due to the limitations of transcriptome data with no reference genome, only 16 CsHsfs were identified. Here, we took advantage of the high quality tea plant reference genome [60] to identify and characterize the CsHsfs family bioinformatically.
The numbers of Hsf family are diverse in different plant species, and there are 25 Hsfs in pepper, 29 in Chinese white pear, 30 in sesame, 26 in soybean and 78 in bread wheat [12,[61][62][63][64], respectively. In this study, we identified 25 CsHsfs from tea plant genome and classified them into A, B, and C subfamily on the basis of their structures and phylogenetic relationships among C. sinensis, A. thaliana, P. trichocarpa. The number of subclass HsfA9 members in tea plant enlarged with four members compared with only one member in A. thaliana and tomato (Table 1), which suggested the possibility of gene replication events during evolutionary process. Likely, the subclass HsfA9 also had four members in pepper [64]. However, the expansion reasons of the CsHsfA9 genes remain to be elucidated by further investigations. The theoretical isoelectric point of CsHsfB3c was 10.01 (Table 1), implying that it was a basic protein, and the other CsHsfs were acidic proteins,  [37] and HemI [56] were used to create the heat map. The color bar represents the normalized FPKM values (log2-transformed), and the black represents no expression detectable. Detailed FPKM values are listed in Additional file 7: Table S6. b Verification of the RNA-seq results of eight CsHsf genes by qRT-PCR analysis. Three biological replicates and three technical replicates were performed for the experiment. The qRT-PCR data was analyzed by ANOVA followed by Fisher's LSD multiple comparison tests and ** represents significant differences at P < 0.01 which indicated that they might play roles in different microenvironments [65]. The GRAVY results were all negative ( Table 1), indicating that they were all hydrophilic proteins, which was consistent with the results in potato [66], carnation [67], and Chinese cabbage [68].
The highly conserved DBD domain consists of about 100 amino acid residues among various plant species [10]. However, it is noteworthy that the DBD of CsHsfA5a, CsHsfA9c, CsHsfA9d, and CsHsfB4b was shorter than the other CsHsfs ( Fig. 1; Table 2), which may be attributed to the incomplete assembly of tea plant genome. AHA domain, which is specific to HsfA subgroup, is indispensable to activate the transcription of heat shock proteins (HSPs), but was not detected in CsHsfA5a and CsHsfA5b, suggesting that they might play their roles by binding to other HsfAs to form hetero oligomers [69]. Phylogenetic analysis showed that the CsHsfs could be divided into three main classes corresponding with those in A. thaliana and P. trichocarpa (Fig. 2), which was consistent with those reported previously [9,61,70]. The length of intron exhibited certain degrees of variation (Fig. 3), which was similar to that in other plants like potato [66] and carnation [67]. Additionally, some homologous genes had differences in intron numbers, intron length, and intron position, implying that their functions may be differentiated.
Promoter analysis indicated that the quantity and variety of cis-elements in each CsHsf gene were obviously different (Fig. 4), presumably suggesting that the transcription of CsHsfs may be differentially regulated by the combination of response elements. Moreover, no HSE element was detected in these CsHsf promoter regions, which implied that the expression of these heat-related CsHsf genes might not be directly induced by heat stress [66,67]. The exact mechanism of gene expression needs further research. Further analysis of tea plant degradome data showed that 7 CsHsfs were predicted to be targeted by 9 miRNAs (Additional file 3: Table S2), implying that CsHsfs could be regulated at post-transcriptional level and miRNAs and their targets were not in one-to-one Fig. 9 Subcellular localization of CsHsfA2 in onion epidermal cells. GFP indicates the vector control. CsHsfA2-GFP and GFP-CsHsfA2 indicate CsHsfA2 was fused with GFP at the C-terminus and N-terminus, respectively. The GFP signal was detected using a Zeiss LSM700 confocal laserscanning microscope (Carl Zeiss Inc., USA) at 488 nm. Bars = 50 μm correspondence [71], but this hypothesis needs to be experimentally validated.
The exploration of gene expression patterns may help in understanding their biological functions [72]. In this study, the expression profiles of each CsHsf gene in eight tissues or exposure to different stresses (i.e., drought, salt, heat, and exogenous ABA) were investigated. Several CsHsfs showed tissue-specific expression patterns, such as CsHsfA3 in old leaf, CsHsfB2b in fruit and young leaf, and CsHsfC1 in root (Fig. 5), suggesting that these CsHsfs might be involved in the development of various organs and tissues. However, the expression of CsHsfA9b, CsHsfA9c, CsHsfA9d was undetectable in all tissues, including apical bud and young leaf, which was confirmed by qRT-PCR alalysis in control samples (one bud and two leaves), suggesting that the three members may have functional redundancy in tea plant. Previous studies have showed that the expression of Hsfs can be regulated by various abiotic stress, especially heat stress [22,33,73]. AtHsfA2, a key heat-inducible gene, could also be induced by salt and osmotic stresses in A. thaliana [74]. Likely, LlHsfA2 expression could be induced by heat shock, but not by salt treatment in lily (Lilium longiflorum) [23]. Moreover, CmHsfA4 was highly induced by salt stress in chrysanthemum [30]. In this study, the transcription of CsHsfA2 was highly up-regulated by heat, salt, and drought stresses (Figs. 6, 7), but was down-regulated by exogenous ABA, suggesting its different roles in responses to various stresses. AtHsfA3 was involved in heat and oxidative stresses responses [26,75], while CsHsfA3 could be induced by heat, salt, drought, and exogenous ABA. CsHsfA7 was found to be up-regulated by heat shock stress, which was consistent with our previous findings [45]. Overexpression of OsHsfA7 enhanced salt and drought tolerance in transgenic rice [76], while AtHsfA7b positively regulated salt stress tolerance in A. thaliana [31]. Hence, the specific functions of CsHsfA7 in response to heat stress remain to be elucidated. TaHsfC2a-B played important roles in developing wheat grains via an ABA-dependent pathway in response to heat stress [33]. In addition to HS, CsHsf genes were also regulated by ABA (Fig. 7b). CsHsfA3, CsHsfA8, CsHsfB1, and CsHsfC1 were up-regulated by exogenous ABA, while CsHsfA2, CsHsfB2a, and CsHsfB2c were downregulated, suggesting that these CsHsfs played different roles in ABA-mediated regulatory pathway.
Ca 2+ is a ubiquitous secondary messenger, and plays vital roles in response to a variety of environmental stresses [77,78]. Here, exogenous Ca 2+ pre-treatment induced the expression of CsHsfA4a under heat shock stress (Fig. 8), implying its potential roles in response to heat stress by Ca 2+ signal pathway, but testing this hypothesis requires further research. Overexpression of AtHsfA4A resulted in tolerance to salt stress in A. thaliana [79]. Moreover, ectopic overexpression of BnHSFA4a in enhanced desiccation tolerance in A. thaliana seeds [80]. However, overexpression of sunflower HaHsfA4a alone did not confer thermotolerance in transgenic tobacco [81]. These results suggested the diverse roles of HsfA4s in responses to different abiotic stresses. Interestingly, the qRT-PCR results of CsHsfA2 and CsHsfC1 expression showed the opposite trends with RNA-seq results (Fig. 8), which was a normal phenomenon in other studies [46,61,64] and might be due to the different normalized method between FPKM and qRT-PCR.
Previous studies have demonstrated that the majority of Hsfs were localized to the nucleus [64]. In our study, CsHsfA5a and CsHsfA8 were predicted to be localized in the cytoplasm and chloroplast, respectively, while the other 23 CsHsfs were targeted to the nucleus (Table 1). To further confirm the subcellular localization prediction, CsHsfA2, CsHsfA5a, and CsHsfA8 were selected and transiently expressed in onion epidermal cells. Interestingly, CsHsfA2, CsHsfA5a, and CsHsfA8 fusion proteins were localized in the nucleus of onion epidermal cells ( Fig. 9; Additional file 9: Figure S2). Hence, the subcellular localization prediction results could not reflect the true locations of target proteins, and needed to be verified by experiments. It is worth noting that HsfA2 is a key heatresponsive gene in resistance to heat stress [9,22,82]. Since tea plant has no genetic transformation system available, Fig. 10 Heterologous expression of CsHsfA2 confers thermotolerance in transgenic yeast. pPIC3.5 K and CsHsfA2 indicate the yeast cells transformed with the empty vector and recombinant plasmid, respectively. P. pastoris cells were exposed to heat stress (45°C for 30 min), and then spotted onto YEPD medium plates (OD 600 = 0.1, 0.01, and 0.001). The control samples were incubated under normal temperature conditions (30°C) we heterologously expressed CsHsfA2 in eukaryotic model organism yeast to dissect the biological function of CsHsfA2 in response to heat stress. Thermotolerance assays indicated that heterologous expression of CsHsfA2 improved yeast resistance to high temperature (Fig. 10). To further elucidate the functions and regulation mechanisms of CsHsfA2 in tea plant through virus-induced gene silencing, as underway in our lab, we would be able to explain the specific molecular mechanisms of CsHsfA2 regulating tea plant heat response.

Conclusions
In this study, we identified and comparatively analyzed 25 full-length tea plant Hsfs in both coding sequences and gene-expression profiles, as well as their expression patterns in responses to abiotic stress. It is worth noting that Ca 2+ signal and ABA pathway seemed to be involved in the CsHsf-mediated heat response. Moreover, CsHsfA2 was located in the nucleus. Additionally, CsHsfA2 conferred thermotolerance when heterologous expressed in transgenic yeast. All these results provide useful information for elaborating the Hsf-mediated stress-response system in tea plant as well as other plant species.

Methods
Identification and sequence analysis of Hsf genes from C. sinensis The amino acid sequences of 21 A. thaliana Hsf genes were downloaded from National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm. nih.gov/) as queries to search against the C. sinensis var. sinensis genome [60]. Furthermore, all obtained CsHsf proteins were analyzed to detect DBD domains and coiled-coil structures by the SMART [50] and CCD programs (https://www.ncbi.nlm.nih.gov/cdd/). Finally, to verify the accuracy of these sequences, BLASTN similarity searches against the published data of C. sinensis were performed with a threshold E-value of less than 1.0E− 90.
The physicochemical parameters of CsHsf proteins were calculated using the ExPASy program (http://web. expasy.org/compute_pi/) [47] with default parameters. WoLF PSORT (https://wolfpsort.hgc.jp/) [48] was used to predict subcellular localizations of the CsHsf proteins. The typical functional structure domains were analyzed using SMART [50], Pfam (http://pfam.xfam.org/search) [51], NLStradamus [52] and NetNES 1.1 server [53]. Multiple alignments of CsHsf DNA-binding domain (DBD) were analyzed using MultAlin [83]. The MEME (http://meme-suite.org/tools/meme) [49] and WebLogo (http://weblogo.berkeley.edu/logo.cgi) [84] programs were used to analyze and visualize the CsHsf conserved motifs with optimum motif width ≥ 6 bp and ≤ 200 bp and maximum number of motifs = 25. Phylogenetic trees were constructed using the neighbour-joining method in MEGA (version 7.0) software and bootstrap test replicate was set to 1000 [85]. The structures of CsHsf genes were analyzed using the online Gene Structure Display Server [86]. Furthermore, the cis-regulatory elements in the 2000 bp promoter regions of CsHsf genes were analyzed in the PlantCARE program [54], and visualized the number of cis-elements in each CsHsf on the HemI software [56]. In addition, psRNATarget online tool [87] was used to predict the miRNAs targeting the CsHsf genes according to the degradome libraries constructed in our lab (unpublished data).

Expression profiles of CsHsf genes based on transcriptome data
To study the expression of CsHsf genes in eight tissues (i.e., apical bud, flower, old leaf, young leaf, fruit, mature leaf, stem, and root) and responses to drought and salt stresses, the C. sinensis expression data were downloaded from Tea Plant Information Archive (TPIA) [55]. The TPM (Transcript per million) values of each CsHsf gene were identified and log10 transformed.
The transcription of each CsHsf under heat stress in C. sinensis leaves treated with exogenous calcium was calculated using RNA sequencing data [37]. The expression profiles of each CsHsf were visualized using the HemI software [56].

Plant materials and treatments
One-year-old cutting seedlings of tea plants [Camellia sinensis (L.) O. Ktze. 'Echa No. 10'; an individual of 'Enshi-taizicha' group species] were cultivated in a growth chamber at Huazhong Agricultural University (Wuhan, China) with a photoperiod of 12 h light (24 ± 1°C)/12 h dark (20 ± 1°C) for 22 days before treatments. To simulate heat stress, the plantlets were placed in an illumination incubator (38°C). For ABA treatment, the plantlets were treated with 50 μM ABA as described by Wang et al. [58]. Young shoots (one bud and two leaves) were harvested at 0, 1, 2, and 12 h after each treatment, immediately immersed in liquid nitrogen, and stored at − 70°C prior to RNA extraction. Additionally, exogenous Ca 2+ pre-treatment (20 mM CaCl 2 ) was conducted as described previously [37], and the samples were collected after exposure to heat stress for 4 h. Three plants were pooled and taken as one biological replicate and three biological replicates were used.
Quantitative RT-PCR analysis of CsHsf genes in C. sinensis Total RNA was extracted using the Quick RNA Isolation Kit (Huayueyang, Beijing, China). Equal amounts of total RNA (1 μg) in all samples were treated with DNase I to eliminate genomic DNA contamination, and then used for cDNA synthesis using a TransScript® II All-in-One First-strand cDNA Synthesis SuperMix for qPCR (One- Step gDNA) Kit (TransGen, Beijing, China). The resulting cDNA was diluted 25-fold in distilled deionized water for qRT-PCR assay. The qRT-PCR assays were performed as described by Wang et al. [88] in a StepOne Plus™ Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). Gene-specific primers (Additional file 10: Table S8) were designed according to the Minimum Information for Publication of Quantitative Real-Time PCR Experiments guidelines [89]. C. sinensis β-actin (Genbank accession number HQ420251) served as the internal reference gene for qRT-PCR normalization analysis, and the relative expression levels of CsHsf genes were calculated using the 2 -ΔΔCT method [90]. All experiments were conducted with three biological and three technical replicates.

Subcellular localization of CsHsfA2, CsHsfA5a, and CsHsfA8
The coding regions of CsHsfA2, CsHsfA5a, and CsHsfA8 without stop codon was firstly fused to the plant expression vector pCAMBIA2300-C-GFP and pCAMBIA2300-N-GFP using a Seamless Assembly and Cloning Kit (Aidlab, Beijing, China), respectively. Then, the onion epidermal cells were transformed and detected as described by Wang et al. [58].

Thermotolerance analysis of transgenic yeast
The ORF of CsHsfA2 was inserted into the pPIC3.5 K yeast expression vector (Invitrogen, Carlsbad, CA), and then the recombinant plasmid was transformed into Pichia pastoris SMD1168 competent cells (Invitrogen) using the freezethaw method. The thermotolerance assays were conducted as described by Jiang et al. [91] with minor modifications. Briefly, P. pastoris cells harboring the recombinant plasmid (OD 600 = 1.0) or pPIC3.5 K were incubated in a water bath at 45°C for 30 min, and then the yeast cells were spotted onto YEPD medium plates after 10-fold dilutions. The photograph was taken after 3 d cultivation under normal temperature conditions (30°C).