Alternative splicing diversifies the heat response and evolutionary strategy of conserved Heat Shock Protein 90 in bread wheat (Triticum aestivum

Background Crops are challenged by the increasing high temperature. Heat shock protein 90 (HSP90), a molecular chaperone, plays a critical role in the heat response in plants. However, the comprehensive response profiles of HSP90s are largely unknown. Results Using the newly released wheat reference sequence, we identified 18 HSP90s that were highly conserved and presented as triplets. Transcriptome analysis showed that the expression and transcriptional regulation of HSP90s displayed conserved trends among the three HSP90 homologs within a triplet. Based on our previous data generated by combining full-length single-molecule sequencing and Illumina short read sequencing, 126 isoforms of HSP90 were identified and each HSP90 comprised one to three major isoforms. Intriguingly, unlike the conserved transcriptional regulation of HSP90s , different alternative splicing events were characterized among major isoforms generated by HSP90s in triplets under heat stress. Furthermore, many heat-responsive alternative events were also characterized, especially in the early heat response when transcriptional regulation had not been initiated, and the alternative splicing regulation between HSP90s in the triplets was distinctly different. Conclusions This study shows that alternative splicing diversified the heat response of the conserved HSP90 gene family and highlighted the evolutionary divergence of polyploidy wheat.

Plants have evolved complex systems to cope with heat stress. Heat shock proteins (HSPs), which function as molecular chaperones, have well-established roles in the heat response [7]. Ninety kDa HSPs (HSP90s) function not only in protein folding, degradation and transportation, as characteristic of chaperones, but also in signaling transduction and protein kinase activity regulation, in an ATP-dependent manner [8,9]. Under heat stress, the 70 kDa HSPs/HSP90s complex disassociates and releases the master heat response regulator, heat shock transcription factor A1s (HSFA1s), resulting in the activation of the plant heat-responsive transcription cascade [10,11]. More recently, HSP90s have been proven to have the ability to stabilize the circadian clock ZEITLUPE and auxin cofactor Fbox protein to maintain plant growth under heat stress [12,13]. Under normal conditions, the downregulated expression of Arabidopsis HSP90s results in abnormal growth and development, such as embryo defects [14]. HSP90s are also considered to exert capacitor buffering effects of genetic perturbation, such as genetic variations and mutations, on plant morphology and phenotype [15,16].
In wheat, nine cytosolic HSP90s have been reported, of which three HSP90AAs are highly expressed in the reproductive organs and are necessary for seedling growth, while six HSP90ABs are constitutively expressed and are essential for disease resistance [21]. However, the detailed landscape of heat response and functional redundancy and diversity of HSP90s in hexaploid wheat is limited, although the upregulated expression of some wheat HSP90s under heat stress has been reported [22,23].
With post-transcriptional regulation, many heat response genes, including several types of transcription factors and HSPs, change their alternative splicing (AS) patterns and produce new isoforms, expanding the diversity of proteome and regulation modes in the heat response [24][25][26]. For example, under normal conditions, HSFA2 mainly encodes a truncated isoform without transcriptional activation activity, while under heat stress, the intact and transcriptionally active isoform is largely expressed to induce the expression of heat response genes [27,28]. Interestingly, a new, small truncated form is induced and this results in the autoregulation of HSFA2 under severe heat stress [29]. For droughtresponsive element binding factors, the inactive nonfunctional isoform is prevalent under normal conditions, whereas the active functional isoform dramatically accumulates under heat stress and thus induces the expression of downstream genes, leading to the heat stress response [30]. Whether and how HSP90s respond to heat stress with posttranscriptional regulation remain to be determined.
In this study, we performed a genome-wide identification of the HSP90s in bread wheat using the newly released wheat reference sequence (IWGSC RefSeq v1.0) [31] and investigated the transcriptional and alternative splicing reprogramming of HSP90s using the dynamic and intensive heat response transcriptomes reported in our previous study [32]. The 18 HSP90s were highly conserved in sequence and expression pattern within a subfamily, while the major isoforms of homeologous HSP90s within a subfamily appeared to undergo different AS events and the corresponding isoforms exhibited diversified heat response patterns. Our findings indicated that AS regulation diversifies the heat response of conserved HSP90s and the advantages of the polyploidy nature of wheat in environmental adaptation.

Results
HSP90s are highly conserved within triplets As a result of sequence searches and domain confirmation, 18 HSP90s were identified in the wheat genome. The identified HSP90s were only present in homeologous groups 2, 5 and 7, and members in each subfamily were evenly distributed among the A, B and D subgenomes (regarded as a triplet), exhibiting a conserved copy number during wheat polyploidization. The 18 HSP90s were classified into HSP90AAs (three, triplet AA), HSP90ABs (three, triplet AB-5 and three, triplet AB-7), HSP90C1s (three, triplet C1), HSP90C2s (three, triplet C2) and HSP90Bs (three, triplet B) based on the phylogenetic tree ( Fig. 1, Table 1). The characterized cytosolic HSP90s were identical to the previously reported cytosolic HSP90 sequences [21], suggesting the high reliability of our characterization procedure.  Table S1. Identities between each HSP90 proteins. Table S2. Detail information for MEME motifs. Table S3. Isoform type classification based on FPKM value. Additional file 2: Figure S1. (A) Protein sequence motifs of HSP90s, motif analysis was performed by MEME (http://meme-suite.org/). (B) Gene structure analysis of HSP90s was carried out at Gene Structure Display Server 2.0 (http://gsds.cbi.pku.edu.cn/). Additional file 3: Figure S2. Expression profiles each HSP90 homeolog in different triplet under various conditions. Log2-transformed fold-change values were used to construct a heat map. The samples are indicated by a character that refers to the tissue, followed by the heat stress duration. G, grain. L, flag leaf. "0" is regarded as a control. The left panel shows the phylogenetic tree. Additional file 4: Figure S3. Expression abundance and splicing modes of major isoforms among HSP90 homologs in triplets AB-5 and − 7. Expression abundance of the major isoforms of HSP90s in triplet AB-5 (A) and AB-7 (C). Expression is represented as the log2-transformed value (FPKM + 1). Splicing patterns of the major isoforms of HSP90s in triplet AB-5 (B) and AB-7 (D). Solid boxes represent exons, lines represent introns, arrows indicate the splicing sites that differ from the longest intact coding sequence. Additional file 5: Figure S4. Expression abundance and splicing modes of major isoforms among HSP90 homologs in triplets B, C1 and C2. Expression abundance of the major isoforms of HSP90s in triplet B (A), C1 (C), and C2 (E). Expression is represented as the log2-transformed value (FPKM + 1). Splicing patterns of the major isoforms of HSP90s in triplet B (B), C1 (D), and C2 (F). Solid boxes represent exons, lines represent introns, arrows indicate the splicing sites that differ from the longest intact coding sequence.
Sequence analysis showed that the protein sequence identities of three HSP90 homologs in a triplet were above 96% (Additional file 1: Table S1), indicating that their sequences were highly conserved. Accordingly, the protein sequence motifs (Additional file: 2, Figure   S1A, Additional file 1: Table S2) and gene structures (Additional file 2: Figure S1B) were also highly conserved within a triplet. In detail, the HSP90AAs and HSP90ABs contained three exons and the HSP90C1s, HSP90C2s and HSP90Bs contained 19, 20 and 15 exons, respectively.
Conserved heat response pattern among HSP90 homeologs within a triplet Using the previously determined dynamic and intensive heat response transcriptomes of filling grain and flag leaves of wheat [32], we investigated the heat response patterns of wheat HSP90s (Fig. 2). Consistent with reports that HSP90AA is highly heat inducible in Arabidopsis and rice, the wheat HSP90AAs were remarkably induced (fold change > 2 and FDR-adjusted P-value < 0.01) after 10 min and 30 min heat treatment in leaves and grains, respectively. However, the wheat HSP90ABs were also highly heat inducible, whereas they are regarded as constitutively expressed genes in Arabidopsis and rice [17,18]. It is also worth noting that HSP90C2s and HSP90AB-5 did not respond to heat stress in grains, but did respond to heat stress in leaves. Our results demonstrated that the heat response patterns of HSP90s were highly dynamic and more intensive studies should be considered to uncover the array of heat responses of HSP90s.
We next compared the expression profiles of wheat HSP90s and surprisingly found a similar trend in the response of the three HSP90 homologs within a triplet (Fig. 2, Additional file 3: Figure S2). For example, all three of the HSP90 homologs in triplet AA were lowly expressed under normal conditions and 5 min heat treatment, but were sharply upregulated at later heat treatment time points. In conclusion, consistent with the high level of sequence conservation, the expression and heat stress response profiles of HSP90 homeologs were also conserved within a triplet, implying that no significant functional divergence of HSP90 homologs occurred during wheat polyploidization.  Table 1). The isoform number did not correlate with the exon number for each gene in this gene family. Secondly, we detected the expression abundances of these isoforms and found them to vary remarkably among the isoforms generated by the same gene. Intriguingly, some isoforms responded to heat stress with transcriptional regulation according to the following criteria: fold change > 2 and FDR-adjusted P-value < 0.01 (Fig. 3), while the related genes that produced these isoforms were not significant heat response genes under these criteria, for example, the 5B-and 5D-homeologs in triplet AB and the 5D-homeolog in triplet C2. Thus, these transcriptionally heat-responsive isoforms extended our understanding of the transcriptional regulation of HSP90s and further revealed the complexity of the heat stress response for this gene family.
Next, to characterize the predominant isoforms of each gene, we introduced the isoform expression percentage (IEP), which was calculated as the expression abundance ratio of one isoform to all isoforms produced by the same gene. An isoform with an average IEP of more than 30% across all of the time points in a tissue was regarded as a major isoform, an isoform with an IEP less than 5% in all time points was regarded as a rare isoform and all other isoforms were classified as minor isoforms. This analysis led to the classification of isoforms into major (30), minor (44) and rare (52) isoforms (Additional file 1: Table S3).
For 18 HSP90s, HSP90AB ( TraesCS5B01G258900) produced three major isoforms and 12 HSP90s generated two major isoforms. Interestingly, among the two or three major isoforms, one was already annotated in IWGSC RefSeq v1.0 and contained the longest complete coding region comprising the HATPase domain and HSP90 domain. However, another major isoform was newly discovered from our hybrid sequencing data and this produced truncated peptides that possessed only the HSP90 domain. Furthermore, expression analysis showed that the major isoforms generated by the same genes had comparable expression levels and response patterns (Fig. 4, Additional file 4: Figure S3, Additional file 5: Figure S4), making it intriguing as to what roles these truncated peptides play in the heat stress response.
Varied number and splicing modes of major isoforms generated by HSP90 homologs within triplets The above transcription analysis showed that the expression abundance and response patterns of HSP90 homeologs within triplets were conserved. However, regarding the number of major isoforms produced by each HSP90 homolog within a triplet, only one or two homologs produced two major isoforms, with the exception of the homeologs of triplet AA and triplet AB-5, which produced two major isoforms for each homolog (Fig. 4, Additional file 4: Figure S3, Additional file 5: Figure S4). Furthermore, the newly identified Differential AS responses of HSP90 homologs within a triplet By comparing the highly expressed isoform sets and their expression changes between the control and heat stress treatment samples, we identified 11 differentially spliced HSP90s under heat stress (i.e. HSP90s that respond to heat stress with AS regulation) for at least one time point sample (Fig. 5A) using the criterion defined in our previous study [32].
Interestingly, some HSP90s that did not respond to heat stress with transcriptional regulation, did respond to heat stress with AS regulation by generating new isoforms or changing the expression level of highly expressed isoforms, extending our understanding of the heat stress response of HSP90s and their functional conservation. Significantly, for the three HSP90 homeologs within a triplet, only one or two responded to heat stress with AS regulation, demonstrating the differential responses at the AS regulation level and suggesting diverged evolution in the heat stress response among the three HSP90 homologs.
Furthermore, using the qualitative and quantitative isoforms, we investigated the isoform with the highest abundance in each sample (Fig. 5B). The results showed the highest abundance isoform of HSP90s that generated more than one major isoform and produced intact HSP90 proteins or truncated peptides that altered among heat stress treatment time points. Interestingly, some minor isoforms defined by our criterion switched to become the highest abundance isoforms, highlighting the important role of AS regulation.
It was also worth noting that the highest abundance isoforms generated by three HSP90 homeologs at specific time points were also different, providing further evidence of the differential responses and diverged evolution. In conclusion, inconsistent with the conserved sequences and transcriptional regulation, AS regulation and the most highly abundant isoforms differed among the three HSP90 homeologs within a triplet, diversifying their heat stress response.

Discussion
HSP90s play vital roles in plant growth and the stress response [7,12,13] However, in our analysis, the expression trends and fold changes of HSP90s were not significantly distinguishable between each other in a triplet. Thus, it seems like the three HSP90 homeologs in a triplet are conserved during wheat polyploidization and have not undergone subfunctionalization or neofunctionalization for heat responses, which was consistent with the highly conserved sequences and motifs.
However, in our subsequent analysis, we found that the number and expression of isoforms generated by each HSP90 homeolog were distinct, and further investigations revealed that the AS modes of the major isoforms in the HSP90 triplet were not rigidly conserved, suggesting that the differentiation of HSP90 homeologs occurs at the AS level.
In allotetraploid Brassica napus, about 20% of the assayed AS events are stress-or organspecific [39] and different splicing patterns have also been widely characterized in duplicated or homoeologous genes in other analyses [40,41]. The differentiation in AS that results in distinct transcripts may be responsible for the diversified functions of homoeologous genes [42]. Thus, divergent AS patterns and differential AS responses may contribute to the functional divergence and differential evolution of HSP90 homeologs, changing our understanding of the conservation of HSP90s in terms of function and expression profile.
AS is an efficient strategy to enhance plant thermotolerance [34,43]. Under heat stress, the percentage of genes undergoing different AS events ranges from about 25-70% in wheat, Physcomitrella patens and grapes [24,25,44], and these genes include those encoding a set of HSPs and heat shock factors. Unlike the functional characterization of HSFA2 isoforms, the identification and analysis of the isoforms of HSPs in heat stress are limited. In this study, wheat HSP90s were found to produce many novel isoforms in the grain-filling stage under heat stress. Contrary to a recent report that the expression of the abnormal isoforms were generally lower than those of the full-length isoforms of Lipoxygenase members in the tea plant in response to pathogen infection and low temperature [45], the expression levels of the truncated major isoforms and their fulllength counterparts were found to be comparable in the current study. Notably, in the AB subfamilies, many truncated major isoforms were preferentially expressed in grains, suggesting the possible tissue-specific function of these genes.
About 40% of the differentially spliced genes are also found to be regulated at the transcriptional level, inferring the vital role of the cooperation of AS and transcriptional regulation in the wheat heat response [24]. In the present study, a much higher percentage was found. All of the 18 HSP90s were transcriptionally regulated and 11 of these were AS regulated. Therefore, transcriptional regulation and AS regulation cooperated in the modulation of HSP90s in response to heat stress, revealing the complexity of HSP90 regulation.

Conclusions
We performed a genome-wide analysis of the HSP90s in wheat and presented the comprehensive expression patterns of these genes and their isoforms over a heat stress time series. Our findings identified a large number of novel isoforms and demonstrated that AS diversifies the heat response of conserved wheat HSP90 genes, highlighting the important role of AS in conferring an evolutionary advantage on polyploid wheat in terms of environmental adaptation.
The hexaploid wheat high confidence protein sequences were downloaded from the URGI (https://wheat-urgi.versailles.inra.fr/). A blastp search was performed using the HSP90 sequences of Arabidopsis and rice as queries and the following parameters: an e-value lower than 1 e˗5 and an identity score above 50%. Using the HMM profiles of the HSP90 domain (PF00183) downloaded from the Pfam database, the Hmmsearch engine in the HMMER3.0 program was used to search the local protein database with a threshold of 1 e˗5 .
Then, the blastp and HMMER results were combined and redundancy was removed. All of the obtained HSP90 candidate sequences were subjected to searches using the

Phylogenetic relationship analysis
Multiple sequence alignment was performed by the CLUSTAL X2 program using the longest full-length proteins with a complete open reading frame. To investigate the evolutionary relationship of HSP90, a neighbor-joining tree with 1000 bootstraps was constructed using MEGA 7.0 software [46]. Sequences from Arabidopsis and rice were used as markers to group the identified HSP90s into different subfamilies.

Transcriptional regulation and AS analysis of TaHSP90s
The RNA-seq data generated by Illumina (150 bp paired-end sequencing by Illumina's HiSeq X Ten platform) and PacBio (PacBio RS II platform) sequencing, and data processing were performed as described in a previous study [32]. The abundances of genes and isoforms were calculated by fragments per kilobase of exon per million fragments mapped (FPKM) values and imported into HemI software [47] to generate a heat map using log2-

List of abbreviations
AS, alternative splicing.
FPKM, fragments per kilobase of exon per million fragments mapped.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The data used in this article could be found at NCBI SRA database under accession no.