Genome-wide identification, classification, and expression analysis of the HSF gene family in pineapple (Ananas comosus)

Identification and Characterization of AcHSF genes in pineapple

The amino acid sequences of pineapple HSFs were retrieved using the HSF-type DBD domain (Pfam: PF00447) as a query in Phytozome JGI. A total of 30 pineapple HSFs were obtained from JGI. Only 22 AcHSFs left, after filtering out the non-typical HSF-type DBD and repeated sequences or canonical coiled-coil structures by SMART online tool (Letunic, Doerks & Bork, 2012). The information of Chromosome localization, CDS, and AA length for AcHSFs was obtained from JGI Phytozome v12.1. The biophysical properties of coding AcHSFs were calculated using the Expasy ProtParam tool. The subcellular localization of AcHSFs was analyzed using BUSCA (http://busca.biocomp.unibo.it/).

Chromosome Localization phylogenetic relationships

The information of all pineapple HSFs’ chromosome localization site was acquired from Phytozome v12.1, including chromosome length, chromosome location, and gene start site. The MapChart v2.0 (https://mapchart.net/) was adopted to map the chromosomal location. The phylogenetic relationship of different HSF proteins was explored, and the phylogenetic tree was created using the AcHSF amino acid sequences and the other three species, Arabidopsis thaliana, A. thaliana, Oryza sativa (O. sativa), and Populus. trichocarpa (P. trichocarpa) by the MEGAX with a bootstrap value of 1000. The HSF gene in pineapples was referred to as AcHSF genes and classified according to HSFs in the phylogenetic tree classes A, B, and C.

Genetic structure and cis-acting elements

The gene structures of AcHSF, including exons, introns, and UTR were displayed by the GSDS online tool (Guo et al., 2007). The promoter sequence of AcHSFs found in Phytozome is located 2kb upstream of the translation initiation site. These sequences were analyzed using a plant cis-acting element database New PLACE (Higo et al., 1999), to identify cis-elements necessary for gene expression, development, and hormone signaling under abiotic stress.

Conserved domains and motifs analysis of AcHSFs

Clustal X 2.0 and DNAMAN software was used to align and edit the DBD domain and HR-A/B regions (OD). Using the cNLS online tool, NLS domains were predicted and the NetNES 1.1 online server identified NES domains in the AcHSFs. MEME server (http://meme-suite.org/) was applied to define the conserved motifs of AcHSFs following the parameters: the number of repetitions = any, the maximum number of motifs = 10, minimum width ≥10, maximum width ≤50, and motifs with an E-value <0.01 were retained.

Expression patterns analysis

The transcriptomic data generated from different organs and developmental stages of pineapple have been described previously (Ming et al., 2015; Wang et al., 2020). Briefly, the different organs include 3 stages of petal tissues, 4 stages of sepal tissues, 6 stages of stamen tissues, 7 stages of ovule tissues, 7 stages of gynoecium tissues, and root, mixed flower, leaf, and 6 stages of fruit were used to generate heatmap using the pheatmap package of R software.

Stress treatments

One-month-old uniform tissue-cultured seedlings in a rooting medium were used for the stress treatment analyses (Priyadarshani et al., 2018). The stress treatments of pineapple seeding as follows: cold (4 ° C), heat (45 ° C), and ABA (100 mM). Leaves were collected from three independent plants at 12 h, 24 h and 48 h after treatment and immediately stored in liquid nitrogen before RNA extraction.Untreated pineapple seedlings of the same size were used as controls.

RNA extraction and qRT-PCR

Total RNA of pineapple leaf tissues was obtained using the RNeasy Plant Mini Kit (50) (Qiagen). For cDNA synthesis, a total of 1 µg RNA per sample was used with cDNA Synthesis SuperMix (Transgen, Beijing, China). qPCR was conducted by the Bio-Rad CFX manager machine using TransStart^® Top Green qPCR SuperMix (Transgen, Beijing, China). Pineapple Actin2 is used as the reference gene for qPCR (Wang et al., 2020). There are total of three biological replicates for each sample, and the results are shown as the mean ± standard deviations.

Genome-wide identification of HSF genes in pineapple

HSF-type DBD domains (Pfam: PF00447) amino acid sequences were submitted into Ananas comosus v3 Phytozome database v12.1. A total of 30 putative pineapple HSF sequences were acquired. And then, checked by the SMART online tool and Pfam database, 1 pineapple HSF sequence was rejected due to the absence of typical HSF-DBD domains, and 7 HSF sequences were abandoned due to the absence of coiled-coil structure. As a result, 22 non-redundant pineapple HSFs were identified (Table 1). The comprehensive information of these 22 AcHSF genes, including gene name, gene ID, CDS and protein length, isoelectric points, molecular weights, predicted subcellular location, and other features are presented in Table 1. The gene with the longest amino acid length is ACHSF-A5, which contains 601 amino acids and also has the largest molecular weight of 65.81 KDa; and the gene with the shortest amino acid length was AcHSF-A9b, which contains 129 amino acids and also has the smallest molecular weight of 13.77 KDa. Prediction of protein isoelectric points (pI) can aid in the purification and isolation of proteins. The predicted isoelectric points (pI) of AcHSFs ranged from 4.68 (AcHSF-B2b) to 9.63 (AcHSF-B4c). Detailed information on other parameters has been given in Table 1.

According to the detailed gene information, 19 AcHSF genes were mapped to the 11 pineapple chromosomes and 3 AcHSF genes located in the scaffold (Table 1). The number of pineapple HSF genes for each chromosome varied significantly, and there is no discernible pattern in the location of these genes on chromosomes. For example, three AcHSF genes were located in chromosome 5, whereas only one was present in chromosomes 2, 6, 17, and 18 respectively (Fig. 1).

Phylogenetic analysis of AcHSFs gene family

A phylogenetic analysis of 31 Populus trichocarpa HSFs, 25 rice HSFs, and 21 Arabidopsis HSFs was performed to classify the phylogenetic relationships (Guo et al., 2016), together with those of AcHSFs by generating a neighbor-joining phylogenetic tree. The HSFs were grouped into three clusters, A, B, and C, according to the difference between the amino acid sequences of the DBD domain, the HR-A/B region, and the linker between them (Guo et al., 2016; Scharf et al., 2012). Class A consisted of 10 subclusters, designated A1 to A10. Class B contained B1 to B4, and Class C comprises C1 and C2 sub-clusters. In pineapple (Ananas comosus), according to their phylogenetic relationship, 12 AcHSFs out of 22 proteins belong to class A, followed by seven AcHSFs belonging to class B, and three copies of class C (Fig. 2). As a monocot, the pineapple was more similar to rice, rather than the dicot Arabidopsis and Populus trichocarpa. However, none of the AcHSFs were found in the subclass A8 and B3, which was reported to only exist in the monocots (Li et al., 2014). It is strange that the pineapple and rice subclass A7 HSFs showed higher similarity to A2 rather than the Arabidopsis and Populus trichocarpa subclass A7, and the AtHSF-A6a also shows abnormal clustering (Fig. 2).

Gene structures and cis -acting elements analysis of AcHSFs

Intron, exon, 5′ UTR, and 3′ UTR structures were analyzed using Gene Structure Display Server (GSDS) v2.0.0 to reveal the gene structural features of AcHSFs. The number of exons for AcHSFs ranged from 1 to 5 (Fig. 3), while only in the longest AcHSF-A9a (genomic sequence 34,598 bp.) was found 5 exons. The 5′ UTR, and 3′ UTR sequence of the AcHSF genes are incomplete, 8 out of 22 AcHSFs (AcHSF-A3, AcHSF-A7b, AcHSF-B1c, AcHSF-B2a, AcHSF-B2b, AcHSF-B4c, AcHSF-C1a, and AcHSF-C1b) do not have 5′ UTR and 3′ UTR sequences, 3 AcHSFs (AcHSF-A1b, AcHSF-A9a, and AcHSF-B4b) have only 5′UTR sequences, 6 AcHSFs (AcHSF-A4, AcHSF-A5, AcHSF-A7, AcHSF-B1a, AcHSF-B1b, and AcHSF-C2) have only 3′UTR sequences, and 5 AcHSFs (AcHSF-A1a, AcHSF-A1c, AcHSF-A2, AcHSF-A6, and AcHSF-A9) both have 5′ UTR and 3′ UTR sequences.

Stress responses elements of the dehydration-responsive element (DRE), ABA-responsive element (ABRE), low-temperature responsive element (LTRE), MYB, MYC, and WRKY elements has been reported to play important roles in drought, salt, cold, ABA, and GA responses (Chai et al., 2020; Li et al., 2012; Zhang et al., 2009). The 2kb sequences upstream of AcHSFs gene were selected for analysis. The cis-acting elements analysis of AcHSFs promoter demonstrated that every pineapple HSF contains at least 2 MYB, MYC, and WRKY elements, except for AcHSF-B4c (Table 2). But for the AcHSF-B4c, only 110bp promoter sequence can be found in the upstream area, among the 110bp promoter sequence, the main core component of the promoter TATA-box and CAAT-box, the light regulatory element (RYREPEATBNNAPA), and root-specific expression related elements (ROOTMOTIFTAPOX1) can be found. We also detected the ABRE, DRE, and LTRE in the AcHSFs promoter area. The result showed that the AcHSF-A1a and AcHSF-A4 lacked ABRE, the AcHSF-A6 and AcHSF-C1b lacked DRE, the AcHSF-C1b lacked LTRE, and the AcHSF-A9a, AcHSF-B4b, and AcHSF-B4c did not have these three stress response elements (Table 2). The cis-element studies in the promoters indicate that HSFs are highly related to the response to stress.

Conserved domains and motifs of pineapple HSFs

The modular structure of the HSFs contains 5 typical conserved domains: DBD, OD, NLS, NES, and AHA domains from N to C-terminal (Table 3). The most conserved DBD domain composed of approximately 100 amino acids, containing three α-helices and a four-stranded antiparallel β-sheet (α1-β1-β2-α2-α3-β3-β4) (Fig. 4A). In addition to the DBD domain, the HR-A/B next to the DBD domain is also important and plays a crucial role in HSF-HSF interaction (Scharf et al., 2012). Besides, HR-A/B also presents in all AcHSFs (Table 3, Fig. 4B). According to the previous studies, HSFs were artificially divided into A, B, and C classes by the distinction between the HR-A and HR-B motifs (Cheng et al., 2015; Giesguth et al., 2015; Singh et al., 2012). In general, the variable length of the flexible linker between parts A and B of the HR-A/B motif of classes A and C HSFs is approximately 15 to 80 amino acids, while the HR-A/B region is tightly connected without the embedded sequence in the middle in class B members. But strangely, the insert lengths between the HR-A and HR-B have almost no difference in pineapple HSFs (Fig. 4B). And the length of the total HR-A/B domain is about 42 amino acids almost the same in pineapple classes A, B, and C HSFs, while the length of classes A and C HSFs is about 50 amino acids and 29 amino acids of class B HSFs in Arabidopsis, rice and soybean (Chauhan et al., 2011; Guo et al., 2008; Jin, Gho & Jung, 2013; Li et al., 2014).

The nuclear localization signals (NLS) and nuclear export signals (NES) are necessary for proteins to import and export the nucleus. The intracellular distribution of HSFs varies dynamically between the nucleus and the cytoplasm, depending on nuclear import and export balance. (Heerklotz et al., 2001; Scharf et al., 1998). After detecting, almost all the HSFs contained NLS sequences rich in basic amino acid residues (K/R), except for AcHSF-B2a, AcHSF-B2b and AcHSF-B4b. However, a total of 8 AcHSFs did not find the NES motifs. As reported in other plants, the transcription activator AHA motif was only located in class A AcHSFs, but the difference is AcHSF-A3 lacks the AHA motif (Table 3).

In addition to the typical conserved domains of HSF, we also detected the putative motifs by Multiple Em for Motif Elicitation (MEME). A total of 10 different motifs were identified in AcHSFs with lengths ranging from 20 to 50 aa (Fig. 5). The motif composition of the same group members is similar, but there are great differences among different group members. The conserved motifs in HSFs indicated that all AcHSFs contained motif 1, motif 2, except for AcHSF-A9a and AcHSF-B4c lack of motif 1. Motif 3 only exists in class A and C HSFs, not in class B. However, motif 7 only present in class A HSFs, and motif 5 only presents in class B HSFs. Additionally, some motifs were only discovered in a certain subfamily of AcHSFs, for example, motif 9 was present in the B1 subclass (Fig. 5).

Expression analysis of AcHSFs in different tissues

Gene expression profiles are related to their functions (Su et al., 2017). To better understand the functions of 22 pineapple AcHSF genes, the tissue-specific expression patterns were detected by 36 different tissues transcriptome sequencing, including flower (mixed stage), leaf, root, fruit S1, S2, S3, S4, S5, and S7, Se1-4, Petal 1-3, Ov 1-7, St 1-6, and Gy 1-7 from Wang et al. (2020).

The results showed some genes are highly expressed in certain tissues, while others are expressed gradually with the development of tissues (Fig. 6). For example, AcHSF-A1c and AcHSF-A7b have high expression levels in 7 fruit tissues, the expression of AcHSF-A9a gradually increased in petal development and have the highest expression value in the P3 development stage. The AcHSF-B4b and AcHSF-B4c are highly expressed in the 7 ovule development stages, which illustrate their important roles in the pineapple ovule development process. We also found that some genes showed tissue-specific expression patterns, such as the AcHSF-B2a was mainly expressed in the fruit S7 stage, AcHSF-A2 and AcHSF-A6 are highly expressed in leaf and flower tissues. In addition, the expression profiles of the genes in the same class are significantly different. For instance, three members of AcHSF-A1 have different expression patterns in all detected tissues and development stages.

Expression profiles of AcHSFs response to various stresses

To extend our understanding of AcHSFs in response to stresses, we performed qRT-PCR to investigate the expression patterns of 6 randomly selected AcHSF genes (AcHSF-A1a, AcHSF-A2a, AcHSF-A9a, AcHSF-B2a, AcHSF-B4a, and AcHSF-C1a) in heat, cold and ABA stresses. The results illustrated that almost all of the selected AcHSFs showed similar expression patterns under the same stress conditions.

Cold stress drastically affects plant growth and development, and leads to a significant reduction in crop yield (Cai et al., 2015); therefore, plants must respond quickly to cold stress. As shown in the results, under the cold stress treatment, the expression of all the 6 AcHSFs increased rapidly from 0 h to 24 h and then reduced at 48 h (Fig. 7A). These may indicate that AcHSFs are commonly up-regulated within a short-timer by cold stresses, and then the expression is down-regulated rapidly. Heat shock transcription factors play crucial roles in response to heat shock induction. The result showed that the expression of 6 AcHSFs continues to increase from 0 h to 48 h in heat stress treatment (Fig. 7B). After ABA treatment, the expression of most selected AcHSF genes increased from 0 h to 12 h, and then decreased after 12 h, while the expression of AcHSF-A2a continued to increase (Fig. 7C). This result implies that the expressions of AcHSFs were suppressed under the longtime ABA treatment and might play crucial roles in different stress (cold and heat, etc.) response pathways.