The genomic, coding, and protein sequences for AcHsfs and AeHsfs were obtained from the Kiwifruit Genome Database (http://kiwifruitgenome.org/) (Yue et al., 2020). The protein sequences for AtHsfs were collected from The Arabidopsis Information Resources Database (TAIR, https://www.arabidopsis.org/) (Lamesch et al., 2012), while the protein sequences for OsHsfs and SlyHsfs were retrieved from the Heat Stress Transcription Factors (HEATSTER, http://www.cibiv.at/services/hsf/) database (Berz et al., 2019). To identify homologous sequences in both kiwifruit genomes (Ac and Ae), the protein sequences of AtHsfs, OsHsfs, and SlyHsfs were used as queries in the BLASTp search. Blast hits with a score of ≥ 100 and an e-value of ≤ 1×e^-10 were considered candidate Hsf proteins from the Ac and Ae genomes. The candidate genes from AcHsfs and AeHsfs were then identified using the HMMER 3.0 (https://www.ebi.ac.uk/Tools/hmmer/) based on the Hidden Markov Model (HMM) of the HSF_DNA-binding protein (PF00447) (Potter et al., 2018). To determine the conserved domains in kiwifruit Hsfs, the Conserved Domain Database (CDD, https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml), Pfam and Simple Modular Architecture Research Tool (SMART, http://smart.embl.de/) were employed (Lu et al., 2019; Letunic et al., 2021; Mistry et al., 2021). The motifs and domains were analyzed based on the HEATSTER database (Berz et al., 2019). The protein sequences containing the HSF_DNA-binding protein were subjected to subsequent analyses.

The length, molecular weight, grand average of hydropathicity (GRAVY), and theoretical isoelectric point (pI) of the kiwifruit Hsf proteins were calculated using the ExPASy server (http://web.expasy.org/protparam/) online tool (Duvaud et al., 2021).

The distribution pattern of intron and exon in kiwifruit Hsf genes was analyzed using the Gene Structure Display Server (GSDS 2.0, https://gsds.cbi.pku.edu.cn/), based on the coding and genomic sequences. In addition, the MEME tool (https://meme-suit.org/tools/meme/) was employed to identify a maximum of 12 conserved motifs in AcHsfs and AeHsfs (Bailey et al., 2009).

Multiple sequence alignments were conducted for AtHsf, OsHsf, SlyHsf, AcHsf, and AeHsf proteins using Clustal X with default parameters (Larkin et al., 2007). The phylogenetic tree was constructed using the neighbor-joining (NJ) method with 1,000 bootstrap replicates through the MEGA X software (Kumar et al., 2018).

Information about the length of a chromosome and the location of genes on the chromosomes was extracted from a GFF file of kiwifruit genomes using an in-house Perl script (Morris et al., 2008). To visualize the genes on chromosomes, the MG2C v2 (Map Gene2chromosome) (http://mg2c.iask.in/mg2c_v2.0) online software was utilized. The identification of duplication events in kiwifruit Hsfs was performed using the MCScan X program with default parameters (Wang et al., 2012). The TBtools software was employed to calculate the Ka (synonymous) and Ks (non-synonymous) substitution rates in duplicated gene pairs of kiwifruit Hsfs (Chen et al., 2020). The syntenic blocks within and between the Ae and Ac genomes were generated using the MCScan X software with default parameters (Wang et al., 2012), and the duplicated gene pairs of kiwifruit Hsfs were visualized using TBtools (Chen et al., 2020).

Four transcriptome datasets (PRJNA691387, PRJCA003268, PRJNA277383, and PRJNA514344) were retrieved from the NCBI (https://www.ncbi.nlm.nih.gov/) and NGDC (https://ngdc.cncb.ac.cn/). These datasets comprised samples from various plant tissues, different developmental stages of fruit, fruit treated with ethylene, and fruit treated with low temperature ( Table S5 ). The kiwifruit genome ‘Ac’ was used as a reference for the reanalysis of the transcriptome datasets (Pilkington et al., 2018). The clean reads were aligned using HISAT2 (v2.0.1) (Kim et al., 2019), and STRINGTIE (v2.1.5) was used to assemble and quantify the reads (Pertea et al., 2015).

The fresh leaf samples of Ae ‘Maohua 1’ and Ac ‘Donghong’ were collected from Kiwifruit Plant Resource Nursery at Lushan Botanical Garden, Chinese Academy of Sciences, Nanchang, China. A portion of the samples was used for gene cloning, while the remaining samples were utilized for callus induction. The calluses were subjected to light and dark conditions for further investigations. In-vitro kiwifruit seedlings were grown and exposed to temperature intervals of 25°C, 28°C, 32°C, 37°C, and 42°C for 3 h in a growth chamber with 55% relative humidity. Each treatment was replicated thrice to achieve uniform results. Tobacco and A. thaliana plants were cultivated under growth conditions including a 16/8 h light/dark photoperiod, a temperature of 26°C, and a humidity level of 55%.

The Hipure Plant RNA Mini Kit (Magen, Shanghai, China) was utilized to extract total RNA from the samples, and the EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (Transgen, Beijing, China) was employed to synthesize cDNA by following the manufacturer’s protocol. The reaction mixture for RT-qPCR was prepared according to the instructions of the MonAmpTM ChemoHS qPCR mix kit as follows: each 20 µl reaction mixture contained 1 µl of template cDNA, 1 µl of the qPCR mix (MonAmpTM ChemoHS), and 0.5 µl of each primer. The reaction conditions comprised an initial denaturation at 95°C for 5 min, followed by 45 cycles of denaturation at 95°C for 10 sec, annealing at 60°C for 20 sec, and extension at 72°C for 20 sec. The kiwifruit actin gene was used as a reference for data normalization. The relative expression was calculated by the 2^-△△Ct method (Livak and Schmittgen, 2001). The kiwifruit AcActin gene (Acc08081) and AeActin gene (DTZ79_07g10460) were selected as the reference genes (Huang et al., 2020; Abid et al., 2022). All the primer pairs used for RT-qPCR can be found in Table S4 .

The coding sequences (CDS) of AcHsfA2a and AcHsfA7b were cloned by using specific pairs of primers AcA2a/AcA7b-EGFP F/R and EGFP-AcA2a/AcA7b F/R from cDNA template of Ac ‘Donghong’ kiwifruit ( Table S4 ). The amplification was performed using the following thermal cycling conditions: initial denaturation at 98°C for 30 s, followed by 35 cycles of denaturation at 98°C for 10 s, annealing at 57°C for 5 s, and extension at 72°C for 10 s. The resulting CDS, without stop codon, were inserted into an empty pGreen vector to generate AcHsf-eGFP and eGFP-AcHsf constructs. An empty pGreen vector was used as a positive control. The constructed plasmids were then introduced into the A. tumefaciens strain EHA105. Finally, the bacterial culture was introduced into tobaccos and A. thaliana plants (with nuclear localization marker 35s-nls-linker-mkate) for performing sub-cellular localization (Yoo et al., 2007; Zhao et al., 2017). The transformed leaves and protoplasts were visualized under a laser scanning confocal microscope Olympus IX83 (Olympus, Tokyo, Japan).

To study the activity of the candidate gene promoter, it was cloned into the pGreenII-0800-LUC vector. The resulting recombinant product was then transformed into the A. tumefaciens strain EHA105. A blank pGreenII-0800-LUC vector was considered as a control. The bacterial culture was injected into the tobacco leaves by using a needleless syringe. The transformed plants were kept in the dark for 24 h and then transferred to normal light conditions for an additional 36 h. Finally, the plants were treated with different temperature intervals: 25°C, 37°C, and 42°C, each for 3 h. The LUC (Firefly luciferase) and REN (Renilla luciferase) values were measured using the Dual-Luciferase Reporter Assay System kit (Promega, Madison, WI, USA).

All statistical analyses were carried out on GraphPad Prism 9 software. The significance levels of data were checked by performing one-way ANOVA with the Dunnett posttest. Mean differences between groups were then examined using Tukey’s test, and mean differences were considered significant at p ≤ 0.05.

We identified a total of 36 and 41 putative Hsf family members from the genomes of Ac (hereafter referred to as AcHsf) and Ae (hereafter referred to as AeHsf), respectively ( Figure S1 and Table S1 ). Notably, the kiwifruit species exhibited a higher abundance of Hsf genes compared to most reported plants ( Table S3 ). Our findings demonstrated that all identified putative Hsf proteins in Ac and Ae possessed the conserved HSF_DNA-binding domain ( Figure S1 ). Interestingly, we observed that several Hsf proteins in Ac and Ae also harbored other conserved domains ( Figure S1 ), indicating their potential roles in defining the functions of specific genes. The CDS length of AcHsfs ranged from 177 bp (AcHsfA4a) to 2190 bp (AcHsfB5a), while that of AeHsfs ranged from 621 bp (AeHsfA4d) to 1536 bp (AeHsfA1c) ( Table S1 ). In terms of protein length, AcHsf proteins varied from 59 to 730 amino acids (aa), and AeHsf proteins ranged from 207 to 512 aa ( Table S1 ). Moreover, the predicted molecular weight of AcHsf proteins ranged from 6772.54 to 83845.65 kDa, and for AeHsf proteins, it ranged from 24202.95 to 56935.26 kDa ( Table S1 ). Theoretical pI calculations revealed that the pI values of AcHsf proteins varied from 4.56 to 9.03, while for AeHsf proteins, the range was 4.64 to 8.96 ( Table S1 ). Subcellular localization prediction indicated that both AcHsf and AeHsf proteins localize in the nucleus of plant cells ( Table S1 ).

The 36 AcHsf genes were distributed randomly across 20 chromosomes in the Ac genome ( Figure 1A ). Our results revealed that chromosome 20 included the greatest number of Hsf genes, with five genes present. This was followed by chromosomes 3, 12, 18, and 24, which contained three genes each. Chromosomes 2, 17, 25, and 27 harbored two genes each, while the remaining chromosomes contained only one AcHsf gene ( Figures 1A, C ). Similarly, the 41 AeHsf genes exhibited an uneven distribution across 20 chromosomes in the Ae genome. Chromosomes 02 and 12 contained the highest number of AeHsf genes, with five genes each. Chromosomes 20 and 24 carried four AeHsf genes, while chromosomes 03 and 18 contained three AeHsf genes. Chromosomes 17, 25, and 27 contained two AeHsf genes, and the remaining chromosomes had a single AeHsf gene ( Figures 1B, C ).

To explore the phylogenetic relationships and evolutionary patterns of Hsf genes in kiwifruit, we constructed a phylogenetic tree using the neighbor-joining (NJ) method. The tree included a total of 150 Hsf proteins, consisting of 36 AcHsfs from Ac, 41 AeHsfs from Ae, 25 OsHsfs from Oryza sativa, 27 SlyHsfs from Solanum lycopersicum, and 21 AtHsfs from Arabidopsis. Based on homology with AtHsfs, OsHsfs, and SlyHsfs, the kiwifruit Hsfs were classified into three main groups, namely groups A, B, and C ( Figure 2 ). Group A was further subdivided into nine subgroups (A1-A9), while group B was divided into five subgroups (B1-B5) ( Figure 2 ). Among the Hsfs from kiwifruit, 21 out of 36 AcHsfs and 26 out of 41 AeHsfs belonged to group A ( Figures 2 , S2 ). The subgroups within group A exhibited varying numbers of kiwifruit Hsf genes, ranging from four to eight. Subgroup A4 had the most kiwifruit Hsf genes, with four AcHsfs and five AeHsfs ( Figures 2 , S2 ). In group B, subgroup B4 contained the highest number of kiwifruit Hsf genes, with four AcHsfs and four AeHsfs ( Figures 2 , S2 ). Lastly, subgroup C1 consisted of two AcHsfs and two AeHsfs genes ( Figures 2 , S2 ).

The conserved domains and motif architectures of kiwifruit Hsfs were predicted by Pfam and MEME, respectively ( Figures 3A, B ). A total of 12 conserved motifs, labeled as motif 1 to motif 12, were identified for the kiwifruit Hsfs ( Figure 3A ). The number of motifs within group A varied from 1 to 11, within group B it ranged from 4 to 9, and within group C it varied from 6 to 7 ( Figure 3A ). The identification of certain motifs in kiwifruit Hsfs suggested the presence of a specific occurrence pattern. For instance, all Hsf genes in both kiwifruit species contained motif 1, while motifs 2, 3, 4, 6, and 9 were found in a majority of the Hsf genes ( Figure 3A ). We further confirmed that motifs 1, 2, 3, and 6 constitute a conserved HSF_DNA-binding domain by using the Pfam and CDD databases ( Figure S4 ) (Lu et al., 2019; Mistry et al., 2021). Interestingly, motif 5 was exclusive to group A Hsfs, while motif 10 was particularly present in group B, except for subgroup B3 ( Figure 3A ). This implies that the presence or absence of specific motifs may play a role in the functional diversification of kiwifruit Hsf genes within particular groups.

The exon-intron structure of genes is an essential feature that provides insights into their evolutionary history, functional diversification, and classification (Zhu et al., 2009). The AeHsf genes possessed a higher number of exons compared to the AcHsf genes ( Figures 3C , S3A ). The exon numbers for AcHsf and AeHsf genes varied from two to eight, with AeHsfA3c from subgroup A3 having the highest number of exons (eight exons) ( Figures 3C , S3B ). This indicates that there are distinct gene structures between AcHsfs and AeHsfs ( Figure 3C ). The number of introns in a gene can also play a role in regulating gene function through alternative splicing of transcripts (Shang et al., 2017; Laloum et al., 2018). Most AcHsf genes (34 out of 36 AcHsfs) had only one intron, whereas AeHsfs exhibited a wide range of intron numbers, ranging from 1 to 7, suggesting the functional diversification observed in AeHsf genes may be influenced by alternative splicing events ( Figure 3C ).

The addition or deletion of genes is the primary evolutionary source causing the expansion or contraction of a gene family (Magadum et al., 2013). There were 11 duplicated gene pairs in the Ac genome and 14 duplicated gene pairs in the Ae genome ( Figure 4A and Table S2 ). The fact that all duplicated gene pairs experienced whole-genome duplication (WGD) events suggests that WGD accounted for the expansion of the kiwifruit Hsf gene family ( Table S2 ) (Zwaenepoel and Van De Peer, 2019).

We estimated the selection pressure experienced by duplicated gene pairs as follows: Ka/Ks < 1 indicates purifying selection, Ka/Ks > 1 denotes positive selection, and Ka/Ks = 1 stands for neutral selection (Hurst, 2002; Zhang et al., 2006). In the present study, the selection pressure ranged from 0.1 to 0.5 for AcHsfs and 0.1 to 0.7 for AeHsfs ( Figure 4B and Table S2 ). Our results suggested that duplicated gene pairs in kiwifruit had experienced purifying selection. The divergence time for Hsf paralog gene pairs in Ac and Ae ranged from 19.61 to 45.31 and 19.72 to 73.45 MYA (million years ago), respectively ( Figure 4C ). Additionally, the divergence time for ortholog gene pairs between the Ac and Ae genomes ranged from 3.26 to 75.93 MYA ( Figure 4C ). Notably, the divergence time for Hsf ortholog gene pairs between the Ac and Ae genomes was higher than that of Hsf paralog gene pairs in Ac ( Figure 4C ).

We utilized four transcriptome datasets to evaluate the expression patterns of AcHsfs. The expression bias of AcHsfs was estimated across eight different plant tissues (flower buds, flowers, fruit T1: no ethylene production, fruit T2: autocatalytic ethylene production, leaf sink, leaves, roots, and shoots) ( Figure 5A ). Our analysis revealed highly tissue-specific expression profiles for kiwifruit Hsfs, including AcHsfA1d, AcHsfA5b, AcHsfA7a, AcHsfA7b, AcHsfA8b, AcHsfB2a, and AcHsfB2c. Specifically, the expression levels of AcHsfA7a were relatively higher in flowers, fruits T1, and roots compared to other genes tested ( Figure 5A ). Furthermore, we investigated the expression of AcHsfs at different fruit developmental stages and found that several genes, namely AcHsfA1d, AcHsfA2a, AcHsfA7a, AcHsfA8b, AcHsfA9b, and AcHsfB2c, exhibited higher expression levels during the post-pollination stage ( Figures 5B , S6 ). On the contrary, genes like AcHsfB1a, AcHsfB1c, and AcHsfB2a showed higher expression levels during the post-harvest stage ( Figures 5B , S6 ). To further understand the potential functions of AcHsfs in fruit development and response to ethylene treatment, we explored the expression profile under different fruit developmental stages with and without ethylene treatment ( Figure 5C ). Interestingly, the results displayed that the expression levels of multiple AcHsfs, including AcHsfA1c, AcHsfA5a, AcHsfA5b, AcHsfA7a, AcHsfA8a, AcHsfA8b, and AcHsfA9b, were down-regulated after ethylene treatment. Notably, ethylene treatment severely inhibited the expression of AcHsfA9b, while AcHsfB1a showed the opposite response ( Figure 5C ). This suggests that these AcHsfs might be involved in the regulation of fruit development and response to ethylene signaling. Additionally, the AcHsfs presented different expression profiles under low temperatures ( Figure S5 ). AcHsfA5b and AcHsfA9b showed an increased expression with temperature increment. On the contrary, some AcHsfs (AcHsfA3a, AcHsfB2c, and AcHsfC1b) showed increased expression under low temperatures ( Figure S5 ).

We selected five AcHsfs and five AeHsfs from kiwifruit plants treated with different temperature intervals (25°C, 28°C, 35°C, 37°C, and 42°C) for 3 h for RT-qPCR analysis ( Figure 6 ) (Wang et al., 2022). The plant samples used in this analysis were collected from the Ac cultivar ‘Donghong’ (DH) and the Ae cultivar ‘Maohua no.1’ (MH) ( Figure 6 ). Our results illustrated that a temperature of 42°Csignificantly enhanced the expression levels of all selected candidate genes, except for AcHsfB1a, indicating that these genes are primarily responsive to high temperatures ( Figure 6 ). Interestingly, the expression levels of two AcHsf genes (AcHsfA2a and AcHsfA7b) in DH and the AeHsfA2a in MH were significantly different at different temperature intervals ( Figure 6 ). AcHsfA2a and AeHsfA2a displayed an upward trend with higher temperatures, whereas the expression levels of AcHsfA7b decreased after 37°C, suggesting that different kiwifruit cultivars exhibit distinct responses to varying temperatures ( Figure 6 ).

Based on the differential expression patterns observed in RT-qPCR results, we selected two candidate genes, AcHsfA2a and AcHsfA7b, for sub-cellular localization analysis. In-silico subcellular localization analysis predicted their presence in the nucleus of plant cells. Experimental analysis in this study confirmed that both AcHsfA2a and AcHsfA7b were indeed localized in the nucleus of tobacco leaves and Arabidopsis protoplasts, regardless of whether the Hsf proteins were connected to the N-terminus or C-terminus of GFP ( Figure 7 ). However, it is worth noting that a slight GFP-AcHsfA2a signal was also observed in the cytoplasm of Arabidopsis protoplasts ( Figure 7B ). In short, both GFP-AcHsfA2a/AcHsfA7b and AcHsfA2a/AcHsfA7b -GFP were primarily located in the nucleus.

To further determine the response of AcHsfA2a to different temperatures, we conducted a transient expression assay using the pGreenII 0800-LUC-pro-AcA2a reporter vector in tobacco leaves. The results demonstrated a significant difference in the expression of AcHsfA2a compared to the control. As the temperature increased, there was a notable increase in the ratio of LUC/REN of pro-A2a -LUC ( Figure 8 ), which suggests that heat stress has the ability to promote the transcription of AcHsfA2a.