Genome-wide analysis of the heat shock transcription factor family reveals saline-alkali stress responses in Xanthoceras sorbifolium

HSF genes identification

The genomes of Arabidopsis thaliana and rice (Oryza sativa) were downloaded from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html) (Nover et al., 2001; Lamesch et al., 2012). The genomic data of X. sorbifolium (Bi et al., 2019), Dimocarpus longan (Lin et al., 2017), and https://www.ncbi.nlm.nih.gov/genome/83651 (Yang et al., 2019) were collected from the Giga database (http://gigadb.org/). The candidate genes of the HSF family were extracted from the X. sorbifolium, D. longan, and https://www.ncbi.nlm.nih.gov/genome/83651 genome databases using the software HMMER (Johnson, Eddy & Portugaly, 2010) based on a Hidden Markov Model (HMM) of the DBD-domain profile (PF00447), obtained from the Pfam database (https://pfam-legacy.xfam.org/). The e-value cut-off was 0.01. To confirm the presence of the DBD domain in candidate XsHSF family members, all of the obtained protein sequences were analyzed by the online CD-search software (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). The putative XsHSF genes were named depending on their genomic location. The subcellular location of HSF proteins was predicted by the Cell PLoc2 software (Chou & Shen, 2010).

Chromosomal location and gene duplication analyses

The chromosomal location information of HSF genes was obtained from X. sorbifolium genome database, and the chromosomal location map was produced by TBtools (Chen et al., 2018). Tandom duplicated genes were identified and paralogs were separated by five or fewer genes in a 100 kb region, as proposed by Tu et al. (2010). Identification of segmentally duplicated genes was conducted based on the plant genome duplication database (Lee et al., 2017; Lee et al., 2013). The non-synonymous substitution (Ka) and synonymous substitution (Ks) values of duplicated pairs were calculated on PAL2NAL (http://www.bork.embl.de/pal2nal) (Suyama, Torrents & Bork, 2006) after alignment of protein sequences. The duplication time of duplicated genes was calculated by T = Ks/2λ × 10⁻⁶ Mya (λ = 1.5 ×10-8 for dicots) (Blanc & Wolfe, 2004).

Gene characteristics, structure, and cis-elements in promoter

Physical parameters of putative XsHSF proteins, including the protein length, molecular weight (kDa), theoretical isoelectric point (pI), and grand average of hydropathicity (GRAVY), were calculated by the ExPASy tool (http://www.expasy.org/tools/protparam.html). The exon-intron structures in each XsHSF genes were displayed by alignment of the CDS with DNA sequences using Gene structure display server (GSDS) program (http://gsds.cbi.pku.edu.cn/). The orthologous gene of XsHSF in Arabidopsis were annotated using the Blast2GO (Conesa et al., 2005). The online database PLACE (https://www.dna.affrc.go.jp/PLACE/?action=newplace) was employed to investigate putative cis-regulatory elements in the sequences about 1,500 bp upstream of the transcription start site (ATG).

Protein domains and conserved motif analyses

The online MEME program was used to display the motif structure of HSF proteins (http://meme-suite.org/tools/meme) on the default parameters except for the optimum motif widths with 6 to 200 residues. SMART (http://smart.embl-heidelberg.de) and CD-search (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) web servers were used for structural motif annotation.

Sequence alignment and phylogenetic analysis

The Clustal X program was used to conduct the multiple sequence alignment of the full-length HSF protein sequences from X. sorbifolium, D. longan, A. yangbiense, Populus trichocarpa, A. thaliana, and O. sativa (Thompson et al., 1997). The Minimum Evolution (ME) phylogenetic tree was constructed by MEGA 6.0 software under the Jones-Taylor-Thornton (JTT) amino acid and Gamma Distributed (G) substitution model. Bootstrap values from 500 replicates were indicated at each node. Finally, the phylogenetic tree was beautified using online software iTOL https://itol.embl.de/itol.cgi.

MicroRNA target site prediction

The PsRNATarget tool (https://www.zhaolab.org/psRNATarget/) was used to predict microRNA target sites in each full-length XsHSFs sequence according to published microRNA sequences from Citrus sinensis, belonging to the Sapindales (Song et al., 2012).

Preparation of plant materials

X. sorbifolium was cultivated in the greenhouse in Qingdao Argriculture University. To explore the expression of XsHSFs, seedlings that had been germinating for 30 days were treated with water, 200 mM NaHCO ₃, or 200 mM NaCl. After 0.5, 6 h, and 7 days of treatments, the roots and leaves of seedlings with different treatments were sampled and immediately frozen in liquid nitrogen for total RNA extraction and real-time RT-PCR assays. All samples were obtained in three biological replicates. All samples were collected and frozen in liquid nitrogen immediately, then stored in −80 °C for RNA extraction.

RNA extraction and qRT-PCR

Total RNA was extracted from mature leaves and young roots, separately, on different treatment using the RNA extraction Kit (Tiangen, Beijing, China). A total of 1 µg of total RNA was reverse-transcribed to first-strand cDNA by the PrimerScript RT Reagent Kit (TaKaRa Biotechnology, Dalian, China). Subsequently, qRT-PCR was performed using fluorescent dye SYBR Green I (TaKaRa, Dalian, China) according to the instructions. GADPH gene was selected as reference genes. Three biological replicates were carried out and expression levels were calculated by the 2^−ΔΔCt method (Schmittgen & Livak, 2008). In addition, 500 ng total RNA was used to synthesis miRNA first-strand cDNA by the miRNA 1st strand cDNA synthesis kit (by stem-loop) (Vazyme, China, Nanjing). The 5.8s gene was used as reference gene for miRNA qRT-PCR. Gene primer sequences were listed in Table S1.

Identification and characteristic analysis

A total of 21 gene sequences that are part of the HSF family were identified in the X. sorbifolium genome using the PF00447 Hidden Markov Model (HMM) sequence as a query (E < 1e⁻⁵). These sequences were named XsHSF1-XsHSF21 based on their chromosomal positions from top to bottom successively (Table 1). The number of HSF genes in X. sorbifolium was similar to Arabidopsis (21), D. longan (20), rice (25), Poplar (27), whereas it was 0.50-fold less than that in Soybean (52) and 1.4-fold greater than that in Acer yangbiense (15) (Nover et al., 2001; Baniwal et al., 2004; Wang, Zhang & Shou, 2009; Li et al., 2014). Gene lengths were between 796 to 6,687 bp. The corresponding protein sequences were submitted to CD-search (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) web servers to confirm the presence of the DBD domain. The XsHSF proteins ranged in length from 211 and 525 amino acids (aa) with an average of 371.4 aa and the isoelectric points (pI) ranged from 4.75 to 9.03. According to the prediction of the protein subcellular localization by Cell-PLoc 2.0 package, all HSF proteins were located in nucleus (Table 1).

Chromosomal location and gene duplication

The 20 XsHSF genes were distributed on 11 chromosomes (Chr) and no HSF genes were found on Chr2, Chr9, Chr14, and Chr15 (Fig. 1). The XsHSF21 gene was the only 1 found on an unmapped scaffold. Gene duplication events, leading to gene functional diversity and plant adaptations, are mainly divided into three categories: segmental duplication, tandem duplication, and transposition (Kong et al., 2007). We identified 1 pair (XsHSF1 and XsHSF4) in the segmental duplication category and 3 genes (XsHSF2, XsHSF9, and XsHSF13) that might have arisen from transposition in the X. sorbifolium genome, accounting for 23.81% of the genes in the family (Data S1). The Ka/Ks ratios of gene pairs were calculated to investigate the gene divergence trend. A Ka/Ks ratio <1 indicates purifying selection. None of the XsHSF were randomly distributed.

Phylogenetic analysis and conservative domain alignment

To reveal the classification of HSF proteins from X. sorbifolium and their evolutionary relationships, full-length HSF proteins in X. sorbifolium, D. longan, Acer yangbiense, Populus trichocarpa, O. sativa, and A. thalian were selected to construct a phylogenetic tree (Fig. 2). The full-length sequences of all HSF proteins identified in this study are presented in Data S2. The HSF proteins of X. sorbifolium are allocated into five major groups (Group α, β, γ, δ, and ɛ). Group α included eight genes (XsHSF5, XsHSF6, XsHSF14, XsHSF8, XsHSF16, XsHSF17, XsHSF18, and XsHSF19), group β included four genes (XsHSF3, XsHSF7, XsHSF10, and XsHSF21), group γ and δ both have four members, and group ɛ only included XsHSF12.

Gene structure and promoter Cis-acting element analysis of XsHSFs

The intron/exon arrangements (Fig. 3, left) illustrated the structural similarities of the various XsHSF genes. Most of the XsHSF genes have one intron in their CDS region except for XsHSF1 with two introns. The cis-acting elements of the promoter play a leading role in regulating the expression of genes involved in phytohormone and environmental response. The 2,000 bp upstream sequences from ATG of XsHSF genes were extracted from X. sorbifolium genome and predicted by Plant CARE to identify cis-acting elements. A total of 21 cis-acting elements were identified in the different XsHSF promoters (Fig. 3 right, Table S2). Among these predicted cis-acting elements, most were related to light responsiveness, including ACE, G-Box, GT1, 3-AF1 binding site, Sp1, MRE, Box 4, LAMP-element, TCCC-motif, and AE-box. Moreover, a large number of cis-elements related to phytohormones were present on XsHSF promoters. For instance, the TGA-element and AuxRR-motif related to auxin responsiveness, the P-box, GARE-motif, and TATC-box related to gibberellin responsiveness, as well as ABRE, which were present in almost all XsHSF promoters except those of XsHSF13, XsHSF15, and XsHSF20, were activated in response to abscisic acid, implying that XsHSF genes could be induced by phytohormone signals. Various elements, including LTR, MBS, ARE, and TC-rich repeats, that are activated in response to stress responsiveness, were found in XsHSF promoters, suggesting that XsHSFs may participate in the response to environmental stimulation. Among these, ARE was found in all XsHSF promoters except for XsHSF2, XsHSF13, and XsHSF18. In addition, the RY-element, involved in seed-specific regulation, O2-sit, involved in the regulation of zein metabolism, CAT-box, involved in meristem development, and the GCN4-motif, involved in endosperm development, were found in several XsHSF promoters, indicating their essential function in plant development.

Conserved domains and motif analysis

The alignment of the two conservative domains (DBD and OD) is shown in Fig. 4. DBD and OD domains were present in all XsHSF proteins. And XsHSF transcription factors were divided into three classes according to the number of amino acids between HR-A and HR-B in OD domain, XsHSF A subfamily includes members in group β, γ, and δ, and XsHSF B and C include the members in group α and ɛ, respectively. The conserved motifs, analyzed by MEME online tool, are shown in Fig. 5. Each putative motif was annotated using Pfam and SMART softwares. The motifs 1 and 2 in the DBD domain as well as motif 3 in the OD domain were relatively conserved in all XsHSF proteins. However, motif 4 encode NES and are extant only in A class. Motif 6 was present only in the HSF A and B classes but not in the C class. Motif 8 and motif 9 were present in the A and C classes while motif 5 was seen in B class.

Expression pattern analysis of XsHSF genes during salt and alkaline stress

To investigate the roles of XsHSFs in response to salt and alkali stress, qRT-PCR analysis was conducted in leaves (Fig. 6) and roots (Fig. 7) and raw data were showed in Data S3. Diverse expression profiles indicated that under normal growth conditions all XsHSFs were detected in leaves and roots but that their expression patterns significantly varied. In leaves, the expression of XsHSF11 and XsHSF15 was induced by NaCl treatment, while the XsHSF1, XsHSF9, and XsHSF7 expression levels were significantly decreased (Figs. 6A, 6E, 6G, 6I, 6K). The trend of XsHSF3, XsHSF8 and XsHSF14 genes expression was first up and then down, while the expression of XsHSF5 showed the opposite trend on NaCl stress. Upon NaHCO₃ treatment, except of XsHSF9 and XsHSF10, the expression levels of most genes, including XsHSF1, XsHSF2, XsHSF3, XsHSF5, XsHSF7, XsHSF11, XsHSF14, and XsHSF15, were significantly decreased in leaves (Fig. 6). In roots, the expression trend of several genes (XsHSF1, XsHSF2, XsHSF3, XsHSF10, XsHSF11 and XsHSF14 gene) under NaCl treatment was the same as that under NaHCO₃ treatment, XsHSF1 and XsHSF10 were up-regulated, XsHSF11 and XsHSF14 were down-regulated (Figs. 7A, 7B, 7C, 7H, 7J). What is noteworthy is that the expression of XsHSF10 was up-regulated more than 80 times after 0.5 h treatment with NaHCO3 solution. Moreover, XsHSF7, XsHSF8, and XsHSF11 were up-regulated on NaCl treatment, while their expression displayed first up and then down in 7d under NaHCO₃ treatment (Figs. 7E, 7F, 7G).

Prediction of microRNA target sites in XsHSFs

The prediction of microRNA target sites uncovered that the activity of several XsHSF proteins was regulated by microRNAs and that some genes might be regulated by mutiple microRNAs (Data S4). For example, the HSFA group gene XsHSF2 might be regulated by miR156, miR394, miR395, miR408, miR7129, and miR854, while the miR156, miR171, and miR390 target sites were found in the HSF B group XsHSF19 mRNA (Table 2). Due to XsHSF3 and XsHSF9 have the target sites of miR164 and miR172, respectively. qRT-PCR was performed to determine the expression levels of miR172 and miR164 in the leaves and roots treated with NaCl (Fig. 8). With the increase of NaCl concentration, the expression levels of miR172 firstly increased and then decreased in leaves, and decreased and then increased in roots, which was opposite to the expression level of XsHSF9 gene, and the correlation coefficient in leaves and root were 0.8791 and 0.6774, respectively. Similarly, the miR172 expression trends in leaves and roots following NaCl concentration were also opposite to those of XsHSF3.

Interaction network of HSF proteins between X. sorbifolium and A. thaliana

The protein–protein interactions of HSF proteins between X. sorbifolium and A. thaliana were analyzed using the online software STRING (http://string-db.org) (Fig. 8). The orthologous protein was matched with the highest bitscore, and the lines with different color represented different types of evidence for the interaction. The XsHSF11 and XsHSF13 proteins had highly similar amino sequences to AtHSFA2, and XsHSF9, XsHSF15, XsHSF4, XsHSF12, XsHSF20 had similar sequences to AtHSF3, AtHSFA8, AtHSF1, AtHSFC1, AtHSFA3, respectively (Table 1, Fig. 9). In Arabidopsis, AtHSFA8, AtHSFA3, AtHSFA2, AtHSFA4A, and AtHSF1C were reported to regulate heat stress (Giesguth et al., 2015; Friedrich et al., 2021). AtHSFA3 binds AtHSFA2 and other proteins form heteromeric complexes to acquire thermotolerance. In addition to heat stress, the HSFA4A gene was also reported to confer salt tolerance and is regulated by the MAP kinases MPK3 and MPK6 (Nover et al., 2001). AtHSF1a is the orthologous gene of XsHSF4 in Arabidopsis, and overexpressing the AtHSF1a gene can promote the tolerance to diverse stressors, including high pH, by promoting inducible HSP expression (Qian et al., 2014). Those results showed that XsHSF genes might regulate the diverse stresses in X. sorbifolium.

Characteristics of the HSF family in X. sorbifolium

The HSF transcriptional factors play a major role in plant development, biotic, and abiotic stress. The number of HSF genes is species-specific, for example, 21, 25, 27, 52, and 56 HSF genes were identified in Arabidopsis, rice (Guo et al., 2008), poplar (Wang et al., 2012), soybean (Li et al., 2014), and wheat (Ye et al., 2020), respectively. But the number of HSF family members did not correlate with the genome size of the species, for example, the number of HSF genes was similar in Arabidopsis, rice, and X. sorbifolia, but the genome sizes of Arabidopsis, rice, popular and X. sorbifolia are vastly different, at 135 Mb (Goodman, Ecker & Dean, 1995), 466 Mb (Yu et al., 2002), and 420 Mb, respectively. Previous studies showed that the number of plant HSF genes may be related to selective pressure during plant evolution in addition to genome replication events (Klaus-Dieter et al., 2012). Compared with soybean (52) and wheat (56), our result indicated a decrease of HSF members in X. sorbifolium, which may be caused by a lack of whole-genome duplication events in X. sorbifolium (Bi et al., 2019). In the present study, genes of the HSF family were also identified in X. sorbifolia, https://www.ncbi.nlm.nih.gov/genome/83651 and https://www.ncbi.nlm.nih.gov/genome/13043, three Sapindaceae trees. No recent whole-genome-wide duplication event was detected in Sapindaceae (Lin et al., 2017), but the number of HSF family members in these species ranged from 15 to 21, indicating that HSF genes in X. sorbifolia, https://www.ncbi.nlm.nih.gov/genome/83651, and https://www.ncbi.nlm.nih.gov/genome/13043 experienced different selection pressure.

Based on the phylogenetic analysis and the result of the multiple sequences alignment, 21 XsHSFs were further clustered into three classes, similar to HSF gene division in other plants (Busch, Wunderlich & Schöffl, 2005; Guo et al., 2008; Ma et al., 2014; Ma et al., 2015; Hu et al., 2015). The DBD domains of XsHSF s were highly conserved, suggesting that the DBD domains were founded before gene functional divergences occurred. Moreover, the number of amino acid residues inserted between HR-A and HR-B in OD domain was different among HSFA, B, and C classes, which is in agreement with the results of previous studies (Fig. 4). Several distinct motifs were observed in the three classes, suggested that three class genes underwent different selective constraints and function divergence.

Response of XsHSF genes to saline-alkali stress in X. sorbifolium

The HSPs play an important role in response to various biological and abiotic stresses by promoting cell repair, protein folding and degradation, and signal transduction. Moreover, the HSFs are the principal regulatory factors that directly activate the expression of HSPs. The expression of HSFs can be regulated by abiotic stresses, including heat, salt, drought, and cold. But the mechanism of HSF gene response to saline-alkali stress is less studied. Various studies have indicated that different HSF family members have varied functions in regulating plant salinity tolerance. Several HSF genes can enhance the salt tolerance in Arabidopsis (Ogawa, Yamaguchi & Nishiuchi, 2007), maize (Jiang et al., 2018), and Populus (Shen et al., 2015). AtHsfA6a acts as a transcriptional activator of stress-responsive genes and their overexpression can be induced by exogenous ABA, NaCl, and drought stress (Hwang et al., 2014). AtHSFA7b regulates its target genes, containing E-box-like and heat shock element motifs, to trigger serial physiological changes to improve salt tolerance, including adjusting osmotic potential to reduce water loss rate, removing reactive oxygen species (ROS), and maintaining cellular ion homeostasis (Zang et al., 2019). AtHSF4A mediates tolerance to salt and oxidative stress by interacting with mitogen-activated protein kinases MPK3 and MPK6, leading to transcriptional activation of HSP17.6A (Pérez-Salamó et al., 2014). The expression of maize ZmHsf04 can increase salt tolerance in Arabidopsis (Jiang et al., 2018). In addition, ssome HSFs are functioning as a negative regulator, such as ZmHsf08 in maize (Wang et al., 2021) and OsHsfB2b in rice (Xiang et al., 2013).

The salt tolerant response of XsHSF genes in X. Sorbifolium is still unclear. In this study, to explore the role of XsHSFs in regulating saline-alkali stress in X. sorbifolium, the seedlings were treated with neutral salt NaCl and alkaline salt NaHCO₃. The expression trends of different XsHSF gene varied in seedlings treated with a neutral salt and an alkaline salt. The expression profile of 11 genes with low sequence similarity was detected in leaves and roots following salt or alkali treatment. In all, four class A members (XsHSF10, XsHSF11, XsHSF15) and class B members (XsHSF5, XsHSF8, XsHSF14) exhibited high expression levels in leaves on salt treatment, suggesting those genes mainly participated in salt tolerance in leaf (Fig. 6). This is in accord with that the orthologous genes of XsHSF10 and XsHSF11, AtHSFA2 and AtHSFA1D, were essential positive regulator for salt tolerance (Ogawa, Yamaguchi & Nishiuchi, 2007; Liu, Liao & Charng, 2011). On alkaline stress, the expression of most XsHSFs were supressed in leaves, including XsHSF1, XsHSF2, XsHSF3, XsHSF5, XsHSF7, XsHSF11, and XsHSF15, but only two genes (XsHSF5 and XsHSF15) were also down-regulated in both leaves and roots under alkali stress, suggesting that XsHSF5 and XsHSF15 may function as negative regulators in alkali tolerance (Figs. 6 and 7). Most XsHSFs expression levels can be induced by alkali stress in less than 12 h, but were suppressed after seven days, indicating that these genes may respond to alkali stress in a short time, but the long-term high concentration of alkali inhibited their expression. Moreover, the functions of XsHSF genes in X. Sorbifolium are worthy of further investigation.