Deep evaluation of the evolutionary history of the Heat Shock Factor (HSF) gene family and its expansion pattern in seed plants

Identification of HSFs and phylogenetic analysis

Here, we sampled 23 species representing three main taxa for pteridophytes, gymnosperms and basal angiosperms including seven genomes and 17 transcriptomes, (Table S1). Most of the transcriptome data was obtained from the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov) (Data from Ran et al., 2018). Multiple databases were screened for the genome assemblies including: ConGenIE (http://congenie.org/), GigaDB (http://gigadb.org/dataset/100209), Dryad (https://datadryad.org/stash/dataset/), WaterlilyPond (http://eplant.njau.edu.cn/waterlily/), FernBase (https://www.fernbase.org/), and the Liriodendron chinense database (http://120.78.193.56:8000/). To increase the reliability of the data, we analyzed both the genomes and transcriptomes of Ginkgo biloba in this study. The methods of our RNA-seq dataset analysis were drawn from a study by Ran et al. (2018). We used the predicted proteome of each genome as a query to search for HSF-type DBD domains (HSF_DNA_bind_PF00447) from Pfam-A.hmm (Pfam release 32.0) using PfamScan software (https://www.ebi.ac.uk/Tools/pfa/pfamscan/), which were considered as candidate genes. We then extracted the amino acid sequences of the HSFs. We also downloaded 537 HSF sequences extracted from 23 plant species (Table S2) representing seven main taxa in the Heatster database (http://www.cibiv.at/services/hsf) and used them in BLAST searches for analyzed species to further identify candidate HSF proteins. For those candidate sequences, we examined the facticity of DBDs and ODs using the SMART 7 software (Letunic, Doerks & Bork, 2012) (http://smart.embl-heidelberg.de/) and the HEATSTER website (https://applbio.biologie.uni-frankfurt.de/hsf/heatster/). The candidate proteins without an integrated DBD domain or HR-A/B domain were removed.

For the phylogenetic reconstruction in this study, we used MUSCLE (http://www.drive5.com/muscle) to conduct the alignment of the candidate genes. Phylogenetic trees were generated using both the NJ and the ML methods. The NJ tree was constructed by TreeBeST (version 1.9.2, http://treesoft.sourceforge.net/treebest.shtml, parameters: -t mm –b 100). Approximately maximum-likelihood (ML) phylogenetic trees were constructed using FastTree (version 2.1.11, http://www.microbesonline.org/fasttree/treecmp.html, with default parameters). Then, phylogenetic analyses were conducted using RaxML version 8.0.19 (Stamatakis, 2014) with 100 bootstraps, the PROTGAMMAAUTO model, and maximum likelihood reconstruction using rapid hill-climbing and rapid bootstrap analyses (-f ad). Phylogenetic trees were examined and manipulated with Evolview v2 (He et al., 2016). We classified the HSF subfamilies using both the phylogenetic tree and the annotation from the HEATSTER website. Some results from the HEATSTER website were inconsistent with the phylogenetic tree, so we performed follow-up checks to confirm the subfamily classification. The final results were based on the appearance of domain characteristic motifs.

We used all HSFA, HSFB, and HFSC genes identified for phylogenetic tree reconstruction in order to better understand the evolutionary relationship within subfamilies and for in depth phylogenetic analyses of the HSFB clade and HSFA-HSFC clade. To understand the complicated evolutionary relationship of the HSFA2, HSFA6, HSFA7, and HSFA9 clades of subfamilies and the HSFC clade, we extracted those two group genes for phylogenetic tree reconstruction, with Chlamydomonas reinhardtii used as an outgroup. Our methods for protein sequence alignments and phylogenetic analyses followed the same steps as previously outlined.

Synteny analysis and molecular dating analyses

We used MCScanX (Wang et al., 2012) to detect the gene replication events and included a total of 21 plant genomes in a synteny analysis covering green algae, mosses, ferns, gymnosperms, basal angiosperms and angiosperms (Table S3). We analyzed all protein models from these genomes for all possible intra- and inter-species genome-wide comparisons and downloaded all genome annotation and corresponding protein sequences for those species. Homologous genes are classified as either orthologous in different species if they are separated by a speciation event, or paralogous in the same species if they are separated by a gene duplication event. We identified the paralogous and orthologous genes in or between those genomes through synteny detection using MCScanX with default parameters (minimum match size for a collinear block = five genes, max gaps allowed = 25 genes). The output files from all the intra- and inter-species comparisons were integrated into a single file named “Total_Synteny_Blocks”, including the headers “Block_Index”, “Locus_1”, “Locus_2”, and “Block_Score”, which served as the database file. We performed the all-against-all protein sequence comparisons necessary for MCScanX using DIAMOND v 0.8.25 (Buchfink, Xie & Huson, 2015). The gene list containing all candidate HSF genes was queried against the “Total_Synteny_Blocks” file. We used these results to identify whether or not HSF genes exist in a syntenic block. We chose eight representative species for gymnosperms (Gnetum montanum, Ginkgo biloba), basal angiosperms (L. chinense, Amborella trichopoda), monocots (Oryza sativa, Zea mays), and eudicots (A. thaliana, Solanum lycopersicum) to do a synteny analysis between species on close taxa. The methods and procedures used were the same as those previously outlined.

The HSFC1-C2 subfamily genes and the HSFA2 and HSFA9 subfamily genes were extracted from the database and used to estimate divergence time. We calibrated a relaxed molecular clock on the node and found the divergence time of monocots and eudicots to be between 140 Mya (a minimum age) and 200 Mya (a maximum age) (represented by the divergence of A. thaliana and O. sativa, Gensel & Andrews, 1984). We performed a Bayesian dating analysis in the Markov chain Monte Carlo (MCMC) tree (Yang, 2007) using an approximate likelihood calculation for the branch lengths, an auto-correlated model of among-lineage rate variation, the GTR substitution model, and a uniform prior on the relative node times. We used Markov chain Monte Carlo sampling to estimate posterior distributions of node ages, with samples drawn every two steps over 200,000 steps following a burn-in of 10,000 steps. We could then trace the gene duplication time based on the resulting gene divergence times.

The phylogeny and evolution of HSFs in land plants

A total of 670 HSF sequences from 44 species were used for the phylogenetic analysis (Tables S1 and S2). We identified 287 new candidate HSF sequences from 24 species, with 228 of those divided into known subfamilies (A1-A9, B1-B5, C1-C2) (Table 1) on the HEATSTER website. Across the comprehensive samples studied, the number of HSF gene subfamilies identified varied greatly, ranging from two in chlorophyta to 30 in angiosperms. The unrooted phylogenetic tree inferred from amino acid sequences was well resovled to three main clades: HSFA, HSFB and HSFC (Fig. 1, Figs. S1–S4). The newly identified HSF genes were re-confirmed on a phylogenic tree. Most subfamilies of clades (A3, A4, A5, A8, A9, B2, B3, B5, C1, C2) were accordingly recovered, while the relationships between these clades were weakly supported. The HSF subfamilies displayed a strong diversification in structure, composition and function (Scharf et al., 2012; Guo et al., 2016; Wang et al., 2018), thus, significant genetic differentiation between clades, especially for HSFA and HSFB were likely resulted from the unstable topology observed. The HSFA group was found in all sampled taxa, while the HSFB group was absent in chlorophyta, and the HSFC group was only present in the angiosperms.

The HSFA group contains major regulators in the HS response of plants (Wang et al., 2018), and, as a result of diversification during plant evolution, displayed variations in different taxa. Interestingly, the A4-A9 subfamily clades are only occurred in angiosperms and A9 genes are only identified in Eudicots. Some subfamilies clustered as a branch, such as A3, A4, A5, and A9, while others were clustered as several branches (Figs. S1 and S5). HSFA1 is a master regulator which cannot be replaced by any other HSF (Scharf et al., 2012) and probably be the most ancient HSFA group. Although all the HSFA1 genes with HSFA8 in angiosperms clustered as a clade, most of the HSFA1 genes from pteridophytes and gymnosperms were dispersed into several clades. The deep divergence of HSFA1 in pteridophytes and gymnosperms indicates the early diversification of HSFA1 before the radiation of all seed plants. Meanwhile, HSFA2, HSFA6, HSFA7, and HSFA9 were blended into a complex clade, and HSFA9 formed a monophyletic group, but others remain unclear. We also noticed that the HSFA2 gene and the HSFA6 gene clustered together with very little genetic difference in some angiosperm species such as O. sativa, Phoenix dactylifera, Citrullus lanatus, as the HSFA6 and HSFA7 did in C. lanatus. The relationship between the HSFA4 clade and the HSFA5 clade was closer in the tree, with two HSFA5 genes sneaked into the HSFA4 clade. It has previously been suggested a close relationship between HSFA3 and HSFC, however, due to the increased number of ferns and gymnosperms, one HSFA1 clade of gymnosperms, rather than HSFA3, was clustered with HSFC.

HSFC only displayed the pattern as the angiosperms clade clustered to one clade of HSFA1 in gymnosperms (Figs. S1, S2, S3 and S5). It is assumed that a duplication event occurred in the ancestral angiosperms which could have contributed to the rise of HSFC. HSFC1 is a common gene subfamily and varies in gene numbers between monocots and eudicots (Table 1, Fig. S5); there are usually two members in most monocots and only one member in eudicots. These results indicate that HSFC experienced steady expansion during the evolution of monocots, and may be involved in important developmental pathways (Wang et al., 2018). Notably, HSFC2 was only present in monocots, but HSFC1 was present in all angiosperm species except for A. trichopoda. In monocots, HSFC1 and HSFC2 clustered together with strong support. HSFC1-HSFC2 clade of monocots group to HSFC1 of eudicots, based with HSFC1 of basal angiosperms. This suggests that HSFC2 is the result of recent duplication that occurred early in the divergence between monocots and eudicots.

Contrary to a previous study (Wang et al., 2018), the results of this study suggest that the HSFB subfamily (HSFB1-HSFB5) is moderately supported as a monophyletic group (Figs. S1 and S6). HSFB1, HSFB2, and HSFB4 have been widely observed across land plants, while both HSFB3 and HSFB5 are only present in the eudicots and basal angiosperms. Although HSFB5, unlike other subfamily members of the HSFB group, has a conserved tetrapeptide LFGV in the C-terminal domain, it is closely related to HSFB3 (Fig. 1). Additionally, the number of HSFB1 genes in gymnosperms is far more than that in angiosperms, with the common number reduced from 3 or 4 in gymnosperms to 1 or 2 in angiosperms (Table 1). In particular, the number of HSFB1 genes in conifers (Picea abies, Pinus taeda, Picea glauca) is significantly increased than that of other seed plants. Multiple copies of the HSFB1 gene in P. abies, P. taeda, and P. glauca clustered and formed a strongly supported monophyletic group. This result indicates that the evolution of these three conifers probably involved both polyploidy and repetitive element activity (Drewry, 1988; Ahuja, 2005; Li et al., 2015). The multi-copy genes may be attributed to two whole genome duplication (WGD) events in the ancestry of major conifer clades (Li et al., 2015). Though many angiosperm lineages have experienced additional rounds of genome duplication (Soltis, Visger & Soltis, 2014; Jiao et al., 2014; Li et al., 2015), there is no obvious proliferation in the member. This result is consistent with the speculation that WGD in angiosperms did not give rise to a remarkable expansion of HSFB1 genes. The HSFB1 genes in angiosperms, gymnosperms and pteridophytes were independently found in different branches, which suggests that HSFB1 is an ancient group which diverged during the evolutionary history of different taxa. HSFB1 genes in gymnosperms experienced several expansions including ancient duplication, while HSFB1 genes in angiosperms rarely retained duplication except for a few recent duplicates. All HSFB2 genes in gymnosperms and angiosperms clustered as a group, respectively. We were unable to trace out a remarkable expansion in gymnosperms, but more than two genes in angiosperms were assumed to be the result of recent duplication. In some species, such as Selaginella moellendorffii, we observed that some genes identified as different subfamilies, such as HSFB1 and HSFB4, and have genetic similarities to highly supported clades. The complicated relationship of these two subfamilies may be a result of recent duplication events. In this study, the HSFB3 and HSFB5 subfamilies were only present in eudicots and basal angiosperms. This is likely the result of duplication events occurring in ancestral angiosperms with subsequent loss of paralogue genes in the monocots.

Gene duplication analysis

To examine the expansion patterns and genetic divergences of the HSF family, a synteny analysis was performed to identify gene duplication events across 21 species (Table S3). We also conducted a synteny analysis between different species on the closely related taxa.

Gene duplication events were identified in 11 species including pteridophytes, basal angiosperms, monocots and eudicots (Table 2). In green algae, moss, and gymnosperms, we did not detect any HSF genes in synteny blocks. In S. moellendorffii, the only non-seed plant analyzed, we identified one pair of duplication genes. These two genes, ‘SelmoHSFB1b’ and ‘SelmoHSFB4,’ belong to different subclasses of the HSF gene family which were observed as being syntenic to each other. We speculate that these genes may be derived from a duplication event and have evolved with differences at the gene sequence level. In L. chinense, the only basal angiosperm analyzed, we identified five pairs of duplication genes with four of those five pairs from the same gene subclass (HSFA2, HSFB1, HSFB2, HSFC1) and the remaining pair from a different gene subclass (HSFA4-HSFA5). Gene duplication events were detected in all sampled eudicot and monocot species. In five eudicots (A. thaliana, Populus trichocarpa, Prunus persica, S. lycopersicum, Mimulus guttatus), we identified 29 pairs of duplication genes out of which 33 pairs belonged to the same gene subclasses (HSFA1, HSFA4, HSFA5, HSFA6, HSFA8, HSFB2, HSFB3, HSFB4, HSFB5) and four pairs belonged to different gene subclasses (HSFA2-HSFA9, HSFA6-HSFA7). In four monocots (O. sativa, Sorghum bicolor, Z. mays, Brachypodium distachyon), we also identified 29 pairs of duplication genes out of which 33 pairs belonged to the same gene subclasses (HSFA1, HSFA2, HSFA4, HSFA6, HSFB1, HSFB2, HSFB4, HSFC1, HSFC2) and four pairs belonged to different gene subclasses (HSFA2-HSFA6, HSFB1-HSFB2, HSFB2-HSFB4). In general, all HSF gene subfamilies except HSFA3 showed the signature of gene duplication. These results also demonstrated that gene pairs from different subfamilies, such as HSFA2-HSFA6, HSFA2-HSFA9, HSFA4-HSFA5, HSFA6-HSFA7, HSFB1-HSFB4, HSFB1-HSFB2, and HSFB2-HSFB4, were paralogous gene pairs.

Beyond that, synteny analysis among different species identified the orthologous genes in different taxa (Table 3). In detail, only HSFA1 genes from different sources were found as orhologous genes between two gymonosperms (G. montanum, G.biloba). As a result of the analysis of gymnosperms (G. biloba) and basal angiosperms (L. chinense), HSFA1, HSFA4, and HSFA5 were detected as orthologous genes. Though the analysis among basal angiosperms (A. trichopoda, L. chinense) and eudicots (S. lycopersicum, A. thaliana) found several orthologous genes, such as HSFA6-HSFA7, HSFA4-HSFA5, HSFA2-HSFA9, and HSFB2-HSFB5, among basal angiosperms and monocots (O. sativa, Z. mays), we only identified out HSFA2-HSFA6 and HSFA2-HSFA7 as orthologous genes. Interestingly, the analysis of eudicots-monocots reveal a consistent pattern as basal angiosperms-monocots with HSFA1-HSFA5 and HSFA2-HSFA7 being identified as orthologous genes in monocots and HSFA6-HSFA7 as orthologous genes in eudicots.

Our results indicate that gene duplication in HSF genes has been a frequent event during the evolution of plants, significantly contributing to their expansion and functional diversification (Fig. 2). Our results also suggest that HSFA4 and HSFA5 have a close genetic relationship, the origin of which may be related to the ancient duplication of HSFA1. It is possible that HSFA6 and HSFA7 originated from gene duplication, most probably derived from HSFA2. HSFA9 was proven to be derived from HSFA2 after the divergence of ancestral angiosperms. HSFB1 is considered to be the most ancient among the HSFB genes, and we predict that HSFB2 and HSFB4 derived from HSFB1 considering the close relationship between them.

Divergence time analysis

The estimated divergence dates of HSFA2 and HSFA9 in eudicots are indicated in Fig. 3. The divergence time of those two gene subfamilies in this study ranges from 131 Mya to 155.2 Mya, which is within the period of Late Jurassic to Lower Cretacous. The estimated split time of the HSFC2 clade and HSFC1 in monocots is indicated in Fig. 4, and ranges from 125 Mya to 190.4 Mya, which is within the Jurassic and Lower Cretaceous periods. The time of the occurence of these gene duplications are consistent with the origin of the most recent common acestors of all living angiosperm, which likely be around 140–250 Mya (Magallón et al., 2015; Foster et al., 2017). Although uncertainty remains for other characters, our reconstruction of the differentiation time scale between gene subfamilies allows us to propose a new plausible scenario for the early diversification of angiosperms at genomic level. The origin and rapid diversification of angiosperms represent one of the most intriguing topics in evolutionary biology (Sauquet & Magallón, 2018), and the evolution research of this gene family (such as the origin, expansion and loss of genes) provides an unprecedented opportunity to explore remarkable long-standing questions that may hold important clues toward understanding present-day biodiversity and adaption to different environments.