Transcriptome analysis reveals heat stress-responsive genes in octoploid strawberry seedlings and expression pattern analysis of FaHSF gene family under heat stress

TL;DR

This study identifies genes in octoploid strawberry seedlings that respond to heat stress, offering insights for breeding heat-resistant varieties.

Contribution

The study identifies and validates heat stress-responsive genes, particularly FaHSF family members, in octoploid strawberries.

Findings

Heat stress increased proline content and antioxidant enzyme activities in strawberry seedlings.

11,526 heat stress-responsive genes were identified through transcriptomic analysis.

FaHSFA2/A3/A7/B1/B2 genes were significantly up-regulated and validated by RT-qPCR.

Abstract

Extreme high temperatures severely affect the cultivation of octoploid strawberry seedlings. However, the molecular mechanisms of heat stress response in strawberry seedlings remain unclear. The results demonstrated that when strawberry seedlings were subjected to 40 °C treatment for 6 hours, the proline content and activities of superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD) all exhibited significant increases compared to the control group maintained at 25 °C. Consequently, strawberry seedlings exposed to either 25 °C or 40 °C for 6h and 24h were selected for subsequent transcriptomic analysis. A total of 11,526 heat stress-responsive genes were identified. Multiple metabolic pathways associated with heat stress were uncovered by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. Sixty-six FaHSF genes were identified in octoploid strawberries, among which FaHSFA2/A3/A7/B1/B2 were significantly up-regulated under heat stress according to transcriptomic analysis, and were further validated by quantitative reverse transcription polymerase chain reaction (RT-qPCR) assay, confirming the reliability of the transcriptome data generated in this study. Collectively, these results facilitate the identification of key heat stress resistance genes in octoploid strawberry and offer valuable insights into the mechanisms of heat tolerance, providing candidate targets for the breeding of heat-resistant varieties.

Figures

Physiological and phenotypic responses of octoploid strawberry seedlings under both control and heat stress conditions.(A–D) Quantitative analysis of proline content, SOD activity, CAT activity, and POD activity at distinct time points post-heat stress. The x-axis denotes the time intervals following heat stress: 0 h, 1 h, 3 h, 6 h, 12 h, and 24 h, respectively. (E) Representative phenotypes of octoploid strawberry seedlings exposed to control conditions (25 °C) and heat stress (40 °C). Enlarged views of selected regions (outlined by boxes in the top row) are displayed in the bottom row for clarity. Damaged root tips are indicated by black arrowheads. Data are presented as means ± SD from three independent replicates. An asterisk (*) denotes statistically significant differences.

DEGs analysis in response to heat stress.(A) PCA analysis among different time points under heat stress. (B) Venn diagram showing overlap between control samples and those subjected to heat stress at 6 h and 24 h. (C–D) Volcano plots depicting the distribution of DEGs under heat stress at 6 h (C) and 24 h (D).

GO analysis showing the top 20 significantly enriched biological processes of DEGs under heat stress at 6 h (A) and 24 h (B).GeneRatio indicates the proportion of DEGs annotated to each KEGG pathway relative to the total number of DEGs. Dot color reflects the q-value, and dot size corresponds to the number of DEGs mapped to the reference pathways. Green and light blue arrows denote the principal biological process terms responsive to heat stress at 6 h and 24 h, respectively. Dark blue arrows highlight newly identified biological process terms among the top 20 significantly enriched categories.

KEGG pathway analysis with the top 20 most enriched KEGG pathways of DEGs under heat stress at 6 h (A) and 24 h (B).The annotation details in this figure are consistent with those described for Fig. 3.

Statistics of transcription factors among DEGs.Transcription factors highlighted in red font denote those exhibiting substantial changes in expression patterns between 6 h and 24 h of heat stress.

Phylogenetic relationships among HSF proteins A. thaliana, F. vesca, and F. ananassa.All 14 HSF subgroups were clearly separated into distinct clades. HSFA, HSFB and HSFC subfamilies were represented by a unique color in the phylogenetic tree.

Predicted conserved motifs and protein domains in the HSFs.(A) Phylogenetic tree displaying HSFs from A. thaliana, F. vesca, and F. ananas, with predicted conserved motifs and protein domains shown adjacent to the respective proteins. (B) Identification of conserved protein motifs within the HSF family using the MEME program, where each motif is indicated by a distinct color. (C) HSF protein domain architecture predicted via the NCBI CDD, with protein length designated by the scale bar at the bottom.

Expression patterns of FaHSFs in octoploid strawberry seedlings under heat stress.(A) Heatmap visualizing the expression profiles of 66 FaHSFs in response to heat stress at 6 h and 24 h. Clustering analysis results for gene expression levels, illustrated using different color scales, are displayed in the upper right corner. The data were extracted from transcriptome datasets and analyzed using TBtools. Pink arrows indicate FaHSFs upregulated after 6 h of heat stress, while red arrows denote FaHSFs that are upregulated at both 6 h and 24 h. (B) RT-qPCR validation of the expression levels of selected upregulated FaHSFs. Data are presented as means ± SD from three biological replicates. Distinct letters above the bars indicate statistically significant differences according to one-way analysis of variance (ANOVA) followed by Duncan’s multiple comparison test (P ≤ 0.05).

Full text

Introduction

The octoploid cultivated strawberry (Fragaria × ananassa) is a globally important fruit crop with substantial economic value, yet its cultivation and yield are highly susceptible to fluctuations resulting from climate change (Zhao et al., 2020; Arief et al., 2023). Strawberries generally thrive in temperate regions, elevated temperatures and extended daylight periods are known to suppress flowering, yet these conditions enhance stolon production, promote seedling growth, and subsequently shorten the cultivation period (Hardigan et al., 2018; Guo et al., 2021). More recently, the increasing frequency of extreme heat events attributed to global warming has further threatened the cultivation of strawberry seedlings (Hu et al., 2015). High temperatures restrict morphological development, elevate the risk of bacterial infection, and disrupt hormonal balance, among other detrimental effects (Matsubara et al., 2004; Sato et al., 2024; López, Denoyes & Bucher, 2024). Therefore, elucidating the gene expression networks that mediate the response of strawberry seedlings to heat stress is essential for breeding new strawberry cultivars with enhanced heat stress resistance.

Heat stress leads to the accumulation of misfolded proteins and excess reactive oxygen species (ROS) in plant cells (Berrios & Rentsch, 2022; Ullah et al., 2024). The heat shock response (HSR) facilitates protein refolding or degradation, and is one of the most common pathways enabling plants to survive heat stress (Ohama et al., 2017). The transcription factors (TFs) known as heat shock factors (HSFs) and their chaperones, heat shock proteins (HSPs), are central to the regulation of the HSR (Friedrich et al., 2021). HSF gene families are widely present in plants, for instance, the Arabidopsis genome encodes 21 HSF genes that display both redundancy and functional divergence (Fu et al., 2022). Under normal conditions, HSFs form inhibitory cytoplasmic complexes with HSPs (such as HSP70/HSP90) (Zhang et al., 2015; Ding, Shi & Yang, 2020). During heat stress, HSPs preferentially bind misfolded proteins, releasing HSFs (Liu, Liao & Charng, 2011). The liberated HSFs undergo trimerization and phosphorylation, becoming active and translocating into the nucleus, where they bind to heat shock elements (HSEs) in the promoter regions of HSP genes, thus promoting HSP expression (Chen et al., 2022; Fragkostefanakis, Schleiff & Scharf, 2025). Newly synthesized HSPs assist in refolding misfolded proteins, thereby reducing their cellular accumulation (Nagaya et al., 2010; Hirai et al., 2011). Subsequently, free HSPs can re-associate with HSFs, prompting the dissociation of HSF trimers into monomers and returning the transcription factors to an inactive state (Song et al., 2015; Figaj, 2025). This regulatory circuit forms a dynamic negative feedback loop, ensuring that the HSR is rapidly activated and promptly attenuated after heat stress subsides (Rao et al., 2022). Excess ROS generated during heat stress can induce protein oxidation, misfolding, and aggregation, thereby inflicting cellular damage. A crucial plant adaptation to heat stress involves scavenging excess ROS (Mittler et al., 2022). The expression of antioxidant genes, such as genes that translate into superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT), is upregulated, boosting the enzymatic antioxidant system that supports plant viability under heat stress conditions (Waszczak, Carmody & Kangasjarvi, 2018; Nadarajah, 2020). These findings highlight the pivotal role of transcriptional regulation in plant responses to heat stress.

Transcriptome sequencing is a powerful approach for revealing gene expression patterns, identifying functional genes, and elucidating regulatory mechanisms underpinning plant responses to heat stress (Song et al., 2024). Previous transcriptome analyses of octoploid cultivated strawberries at the flowering and fruiting stages were limited by the available genome assemblies in 2016, resulting in numerous unannotated genes and a lack of in-depth analysis of the metabolic pathways associated with these differentially expresses genes (DEGs) (Liao et al., 2016). The transcriptomic and metabolic changes of octoploid cultivated strawberry fruits during postharvest cooling and heat storage have been studied, revealing the endogenous metabolic responses of strawberries to temperature variation. In order to examine molecular evidences related to fruit discoloration of octoploid cultivated strawberries, transcriptome of strawberry fruit under high temperature have been analyzed (Zheng et al., 2022). To investigate the heritability of heat stress-induced epigenetic and transcriptomic changes, transcriptome analysis of woodland strawberry (F. vesca) after reproduction has been studied (López, Denoyes & Bucher, 2024).

Despite these advances, research on transcriptomic responses of strawberry seedlings to heat stress remains limited. In this study, physiological and phenotypic differences in ‘Benihoppe’ octoploid strawberry seedlings exposed to control and heat stress conditions were compared, followed by comprehensive transcriptome sequencing at 0 h, 6 h, and 24 h of heat stress. The DEGs were identified and their functions annotated using GO and KEGG databases. Additionally, HSF gene family members were identified in ‘Benihoppe’ and their expression patterns under heat stress were analyzed to validate the transcriptomic data. These results will facilitate the further discovery of key genes associated with heat stress resistance and provide a theoretical foundation for the breeding of heat-tolerant strawberry varieties.

Materials & Methods

Plant material and heat stress treatments

All octoploid strawberry seedlings (‘Benihoppe’) used in this study were provided by the Engineering and Technical Center for Modern Horticulture, Nanjing, China. The ‘Benihoppe’ plantlets were hydroponically cultured in Hoagland nutrient solution and grown in a controlled-environment incubator (25 °C, 80% relative humidity, 10 h dark/14 h light cycle) for six weeks. For heat stress treatment, subsets of seedlings were transferred to another incubator with the temperature raised to 40 °C while maintaining all other conditions. Seedlings were exposed for 1 h, 3 h, 6 h, 9 h, 12 h, or 24 h, with corresponding controls harvested at 25 °C for each time point; 0 h controls were used as the baseline.

Determination of proline content and SOD, CAT, POD activity

To proline content and SOD, CAT, POD activity measurement, 0.1 g strawberry seedlings grown under normal(25 °C) or 40 °C for 0 h, 1 h, 3 h, 6 h, 12 h and 24 h were rapidly ground in liquid nitrogen, respectively. Proline, CAT, and POD assay kits were used, the detailed method described by Li et al. (2023). SOD activity assessment involved homogenizing 0.1 g of liquid nitrogen-ground tissue in one mL of ice-cold extraction buffer, followed by centrifugation at 4 °C for 10 min at 8,000 g, with subsequent analysis of the supernatant using SOD assay kits maintained on ice.

All proline, SOD, CAT, and POD assay kits used were supplied by Suzhou Comin Biochemistry Co., Ltd. (Suzhou, China).

RNA extraction, cDNA library construction, and sequencing

For transcriptome analysis, five groups of strawberry seedlings grown under two distinct temperature regimes (25 °C for 0 h, 6 h, 24 h as controls and 42 °C for 6 h, 24 h as treatments) were sampled, each with three biological replicates. All experimental procedures were conducted within rigorously controlled artificial climate chambers to ensure experimental reproducibility and methodological rigor. Total RNA was isolated from samples using the Trizol-based Total RNA Extraction Reagent (Taraka, Shiga, Japan). RNA was quantified using a Nanodrop UV-Vis spectrophotometer (Thermo Fisher, Waltham, MA, USA), ensuring a minimum yield of > 7 µg with A260/280 ≥ 1.9 and A260/230 ≥ 1.7. cDNA libraries were constructed following Illumina standard protocols and sequenced by Tsingke Biotechnology Co., Ltd. (Beijing, China) on the Illumina HiSeq 2000 platform.

Transcriptome assembly, annotation, and identification of the DEGs

High-quality clean reads were generated using fastp software (version: fastp-V0.20.1) through removal of adapter sequences, contaminated oligonucleotides, poly-N sequences, and low-quality reads (parameters: ∼fastp -i sample.R1.fastq -I sample.R2.fastq -o out.R1.fastq -p sample.R2.fastq -q 20 -l 100 -n 10). The clean data were subsequently aligned to the reference genome (https://www.rosaceae.org/Analysis/18085091) using HISAT2 (version v2.2.1) to generate BAM files (parameters: ∼hisat2-1 R1.fastq.gz-2 R2.fastq.gz-S sample_mapped.sam). Transcript assembly was performed using StringTie (version V2.0.4) with the Fragaria × ananassa ‘Benihoppe’ genome v1.0 as reference (parameters: ∼stringtie -p 4 -G genome.gtf -o output.gtf -l sample input.bam), followed by merging of initial assembly results from all samples. Transcript comparisons against reference annotations were conducted using gffcompare software (version: gffcompare−0.9.8.Linux_x86_64) to obtain final assembly annotations. FPKM quantification was performed using the ballgown package with parameters: ∼stringtie -e -B -p 4 -G merged.gtf -o samples.gtf samples.bam. All downstream analyses utilized these high-quality processed data. The raw transcriptome data have been deposited in the NCBI database under BioProject accession number PRJNA1369106.

Gene annotation for all assembled sequences was conducted using several public databases: Nr (NCBI non-redundant protein sequences), Pfam (Protein family), KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (manually curated protein sequence database), KO (KEGG Orthology), and GO (Gene Ontology).

DEGs analyses between experimental groups were conducted using the DESeq_2_ R package (version 1.26.0, the parameter is https://bioconductor.org/packages/release/bioc/html/DESeq2.html). Genes with a false discovery rate (FDR) ≤ 0.05 and —log_2_ (fold change)—≥ 2 were designated as differentially expressed genes (DEGs).

GO and KEGG analysis

GO enrichment analyses of DEGs were performed with the GOseq R package, utilizing the Wallenius non-central hypergeometric distribution (Young et al., 2010). Statistical enrichment of DEGs within KEGG pathways was conducted using the KOBAS software suite (Mao et al., 2005).

Identification of transcription factors

Transcription factor genes were identified based on the genomic annotation (GFF3) file for ‘Benihoppe’ octoploid strawberry, obtained from the GDR database (https://www.rosaceae.org/Analysis/18085091). Promoter sequences were extracted via the TBtools Fasta Extract or Filter (Quick) function (Chen et al., 2023), and TFs were identified using promoter sequence data analyzed through the PlantTFDB database.

Identification of FaHSF

21 Arabidopsis HSF protein sequences were retrieved from the TAIR database as query references to identify homologous HSFs in F. vesca and F. ananassa. Protein sequences for Fragaria vesca (Whole Genome v4.0.a2) and Fragaria × ananassa ‘Benihoppe’ (Genome v1.0) were obtained from the GDR database. BLAST analyses were conducted using TBtools’ “Blast Several Sequences to a Big Database” function. Structural integrity of the HSF domains was confirmed using the NCBI Conserved Domain Database. Phylogenetic analysis of identified HSFs from Arabidopsis, woodland strawberry, and octoploid strawberry was conducted using the neighbor-joining method in MEGA12 software based on the full protein sequence (Li et al., 2019). Chromosomal locations of FaHSFs were plotted according to the genomic annotation (GFF3) file of the ‘Benihoppe’ octoploid strawberry, and visualized using the TBtools Gene Location Visualize function. Coding sequence fasta and gff3 files of ‘Benihoppe’ octoploid strawberry were used to construct collinearity and gene duplication relationships via MCscanX (https://github.com/wyp1125/MCScanX/releases/tag/v1.0.0), with visualization performed using the TBtools One step MCDcanX function. Gene duplication type analysis was performed using DupGen_finder (https://github.com/qiao-xin/DupGen_finder). Conserved motifs and domains were characterized via the MEME database and NCBI CDD, Cis-regulatory elements within the promoter sequences of each HSF gene in F. ananassa were systematically analyzed using the PlantCARE software platform, with graphical visualization implemented using the TBtools “Gene Structure View” function.

Analysis of the expression pattern of FaHSFs

The expression patterns of all FaHSFs genes under heat stress at 6 h and 24 h were analyzed using RNA-seq data generated in this study, and clustered heatmaps were constructed using the TBtools “HeatMap” function.

Strawberry seedlings subjected to 40 °C heat treatment for 1 h, 3 h, 6 h, 12 h, and 24 h (with parallel controls at 25 °C for matched time points) were ground to a fine powder in liquid nitrogen, and 0.1 g tissue was used for total RNA isolation, and were isolated by using RNAprep Pure Plant Plus Kit (Tiangen, Beijing, China). The first strand cDNAs were synthesized by using PrimeScript™ RT reagent Kit (Takara, Dalian, China). Octoploid strawberry GAPDH gene used as an internal control, quantitative reverse transcription polymerase chain reaction (RT-qPCR) analysis of FaHSFs were conducted by using an Applied Biosystems 7500 real-time qPCR system with the One Step TB Green^®^ PrimeScript™ RT-PCR Kit II (Takara, Dalian, China) (Li et al., 2023).

Statistical analysis

All statistical analyses were performed using the One-way ANOVA algorithm implemented in IBM SPSS Statistics software (Armonk, NY, USA). The statistical power of this experimental design was calculated in RnaSeqPower (https://rodrigo-arcoverde.shinyapps.io/rnaseq_power_calc/), with parameters set to effect = 0.5 and alpha = 0.05, the calculated RNASeqPower was 0.9944105.

Results

Physiological and phenotypic responses to heat stress

To assess the heat resistance of octoploid strawberry seedlings, plants were exposed to 40 °C for a series of time intervals (0 h, 1 h, 3 h, 6 h, 12 h, and 24h), and their physiological parameters (proline content, and SOD, CAT, POD activities) were measured. Significant increases in proline content and SOD activity were observed after 3 h, whereas heightened CAT and POD activities were detected from 6 h onwards (Figs. 1A–1D). Phenotypic analysis revealed that no apparent changes occurred after 6 h of heat treatment, but marked root ulceration became evident in heat-stressed seedlings at 24 h, while such symptoms were absent in the controls (Fig. 1E). In summary, these data suggest that octoploid strawberry seedlings begin to sustain heat damage after 6 h of exposure to 40 °C, while damage is considerably exacerbated after 24 h.

Physiological and phenotypic responses of octoploid strawberry seedlings under both control and heat stress conditions.(A–D) Quantitative analysis of proline content, SOD activity, CAT activity, and POD activity at distinct time points post-heat stress. The x-axis denotes the time intervals following heat stress: 0 h, 1 h, 3 h, 6 h, 12 h, and 24 h, respectively. (E) Representative phenotypes of octoploid strawberry seedlings exposed to control conditions (25 °C) and heat stress (40 °C). Enlarged views of selected regions (outlined by boxes in the top row) are displayed in the bottom row for clarity. Damaged root tips are indicated by black arrowheads. Data are presented as means ± SD from three independent replicates. An asterisk () denotes statistically significant differences.*

Transcriptome and DEGs analyses

To uncover the transcriptional gene responses to heat stress, transcriptome profiling was performed on strawberry seedlings exposed to 25 °C (as control) and 40 °C for 0 h (control, C0), 6 h (control, C6 and heat stress, H6), and 24 h (control, C24 and heat stress, H24). A total of 15 RNA-seq libraries were constructed, with each library yielding approximately 20 million high-quality reads. After quality control, a total of 89.2 GB of clean data was generated, with guanine-cytosine (GC) content exceeding 45.96%, and Q20 and Q30 scores above 98.4% and 95.06%, respectively (Table S1, each experimental point per condition was biological replicate). Moreover, 4,575 novel genes were identified in this study (Table S2), representing a potential resource for future functional investigations. The gene expression distribution was consistent across all samples (Fig. S1A). Principal component analysis (PCA) demonstrated high correlation and repeatability within the same group (Fig. 2A). A lower pearson correlation coefficient (PCC) was observed between C0, C6, and C24 samples, while a high PCC was detected between all control and heat-treated samples (H6, H24) (Fig. S1B). The numbers of DEGs identified in comparisons C0 vs. C6 and C0 vs. C24 were 109 and 190 respectively (Fig. S2). These variations may stem from inevitable operational variations, therefore, subsequent analyses focused exclusively on comparisons between C6 vs. H6 and C24 vs. H24. These results confirm the reliability and accuracy of the transcriptome sequencing data for subsequent analyses. To further study the characteristics of gene expression and dissect the mechanism of heat-tolerance, we identified 11,526 DEGs (heat stress-responsive genes) by comparing gene expression levels between C6 and H6 and between C24 and H24 within each line (Tables S3 and S4). Samples with and without heat stress were clearly separated based on the expression heatmap (Fig. S3). Of 11,526 DEGs, 7,203 were identified in C6 vs. H6 and 7,748 in C24 vs. H24, with 3,425 shared (Fig. 2B). The difference in the number of DEGs in C6 vs. H6 and C24 vs. H24 was mainly due to the addition of 470 down-regulated genes in C24 vs. H24 (Figs. 2C and 2D). The difference of DEGs between C6 vs. H6 and C24 vs. H24 indicates that strawberry seedlings have different molecular mechanisms in response to continuous high temperature.

DEGs analysis in response to heat stress.(A) PCA analysis among different time points under heat stress. (B) Venn diagram showing overlap between control samples and those subjected to heat stress at 6 h and 24 h. (C–D) Volcano plots depicting the distribution of DEGs under heat stress at 6 h (C) and 24 h (D).

GO enrichment and KEGG pathway analysis of DEGs

To investigate the functional roles of DEGs in heat stress response, GO enrichment analysis focusing on biological processes was conducted for C6 vs. H6 and C24 vs. H24 comparisons. The DEGs were categorized into 520 and 507 GO terms, respectively (Tables S5 and S6). Based on adjusted P-values, the top 20 significantly enriched GO terms were selected, with primary processes pertinent to heat stress response being analyzed. Commonly enriched terms included “Response to stress” “Response to stimulus” “Phenylpropanoid biosynthetic process” and “Flavonoid metabolic process” in both C6 vs. H6 and C24 vs. H24 (Fig. 3). Notably, “L-phenylalanine metabolic process” was not among the top 20 in C24 vs. H24, whereas terms such as “Cellular detoxification” “Hydrogen peroxide metabolic process” and “Positive regulation of oxidoreductase activity” emerged within the top 20. Additionally, GO terms directly associated with heat response, such as “Response to heat” and “Response to red or far-red light” were also significantly enriched in both time points (Tables S5 and S6).

GO analysis showing the top 20 significantly enriched biological processes of DEGs under heat stress at 6 h (A) and 24 h (B).GeneRatio indicates the proportion of DEGs annotated to each KEGG pathway relative to the total number of DEGs. Dot color reflects the q-value, and dot size corresponds to the number of DEGs mapped to the reference pathways. Green and light blue arrows denote the principal biological process terms responsive to heat stress at 6 h and 24 h, respectively. Dark blue arrows highlight newly identified biological process terms among the top 20 significantly enriched categories.

KEGG pathway analysis with the top 20 most enriched KEGG pathways of DEGs under heat stress at 6 h (A) and 24 h (B).The annotation details in this figure are consistent with those described for Fig. 3.

To identify metabolic pathways activated by heat stress, KEGG enrichment analysis was performed. In the C6 vs. H6 group, DEGs were assigned to 132 KEGG pathways (Table S7). The top 20 most significantly enriched pathways included “Phenylpropanoid biosynthesis” “Glutathione metabolism” “Flavonoid biosynthesis”, “Carbon fixation in photosynthetic organisms” “Cutin, suberine and wax biosynthesis” and “Flavone and flavonol biosynthesis” (Fig. 4A). For the C24 vs. H24 group, 130 KEGG pathways were detected (Table S8), with the top 20 pathways showing two additional heat stress-related processes: “Fructose and mannose metabolism” and “Zeatin metabolism” (Fig. 4B). These results indicate that the molecular mechanisms of heat stress response differ between 6 h and 24 h treatment periods.

Transcription factors analysis of DEGs

TFs play essential regulatory roles in plant development and stress adaptation. Analysis of all annotated DEGs in this study identified 40 major TF families, with the eight most abundant being bHLH, WRKY, MYB-related, HSF, NAC, bZIP, ERF, and MYB (Fig. 5). In the whole genome, the eight TF families with the highest number of annotated genes included NAC, MYB-related, C2H2, bHLH, WRKY, MYB, DLSV, and FAR1 (Fig. S4). These results indicate that TF families such as bHLH, WRKY, MYB-related, HSF, bZIP, and ERF play prominent roles in the heat stress response of octoploid strawberries.

Statistics of transcription factors among DEGs.Transcription factors highlighted in red font denote those exhibiting substantial changes in expression patterns between 6 h and 24 h of heat stress.

Among the identified transcription factors, the bHLH, GRAS, MYB, MYB-related, WRKY, and GRF families showed higher member counts at 24 h post-heat treatment compared to 6 h, with GRF transcription factors exclusively appearing at 24 h, suggesting their predominant involvement during the late phase of heat stress response. In contrast, HSF, C3H, TALE, and Trihelix families displayed substantially larger quantities at 6 h relative to 24 h, indicating their primary role in the early response to heat stress (Fig. 5).

Identification and analysis of HSF in cultivated strawberries

To validate the accuracy of our transcriptomic data, we next analyzed the expression of HSF genes under heat stress. 21 AtHSF proteins from Arabidopsis were retrieved and used as queries for homology searches and conserved domain validation against the GDR database, resulting in the identification of 17 FvHSFs and 66 FaHSFs proteins in octoploid strawberry, with names assigned based on homology with AtHSF, all HSFs protein sequences were list in Data S1. To explore evolutionary relationships, an unrooted phylogenetic tree was constructed using HSF protein sequences from A. thaliana, F. vesca, and F. ananassa. These HSFs clustered into three subfamilies: HSFA, HSFB, and HSFC, with the HSFA subfamily being the largest and further subdivided into nine groups (Fig. 6). Additionally, the 66 FaHSF genes are randomly distributed across chromosomes, and the distribution among subgenomes is similar (the chromosome chr_2c may have undergone inversion during evolution) (Fig. S5). Collinearity analysis revealed 153 collinear gene pairs in octoploid strawberry (Fig. S6A, Table S9), and the all 66 FaHSFs originated from whole-genome duplication (WGD) events. Notably, diploid woodland strawberry and each subgenome of the octoploid strawberry exhibited collinear gene pairs (Fig. S6B, Table S10). These results indicate that HSF genes are highly conserved in plants.

Phylogenetic relationships among HSF proteins A. thaliana, F. vesca, and F. ananassa.All 14 HSF subgroups were clearly separated into distinct clades. HSFA, HSFB and HSFC subfamilies were represented by a unique color in the phylogenetic tree.

To further characterize FaHSFs, we analyzed HSF motif, domain conservation in A. thaliana, F. vesca, and F. ananassa (Fig. 7A). Nine conserved protein motifs were detected, with motifs 1, 2, 3, and 5 present in nearly all HSFs at consistent positions (Fig. 7B; see Fig. S7 for motif details). Each HSF in these three species contains the core HSF domain, highlighting their high degree of conservation (Fig. 7C). Furthermore, analysis revealed that the promoter region of FaHSF contains numerous cis-elements responsive to temperature and light stimuli, indicating that the FaHSF protein likely possesses the capacity to respond to heat stress (Fig. S8, detailed cis-element information is summarized in Table S11). Thus, these 66 FaHSFs were subsequently used for further expression analyses.

Predicted conserved motifs and protein domains in the HSFs.(A) Phylogenetic tree displaying HSFs from A. thaliana, F. vesca, and F. ananas, with predicted conserved motifs and protein domains shown adjacent to the respective proteins. (B) Identification of conserved protein motifs within the HSF family using the MEME program, where each motif is indicated by a distinct color. (C) HSF protein domain architecture predicted via the NCBI CDD, with protein length designated by the scale bar at the bottom.

FaHSF expression under heat stress

Analysis of the transcriptome data revealed that five out of fourteen HSF subgroups were induced by heat stress. Expression levels of FaHSFA2, FaHSFA3, and certain FaHSFB2 (FaHSFB2b) members were significantly upregulated after both 6 h and 24 h of heat stress, while FaHSFA7, FaHSFB1, and some FaHSFB2 (FaHSFB2a) members were upregulated solely after 6 h of heat stress (Fig. 8A). To validate these transcriptomic findings, RT-qPCR was performed on these five HSF subgroups (Figs. 8B–8F). Due to high sequence homology within subgroups, designed primers were specific to subgroups, with sequences provided in Table S12, and raw data for RT-qPCR were list in Table S11. RT-qPCR results indicated that expression of these subgroups was indeed induced by heat stress, with levels of FaHSFA7, FaHSFB1, and FaHSFB2 at 24 h being significantly lower than at 6 h and closely resembling those in controls (Fig. 8B). Overall, these data confirm the reliability of the transcriptome analysis and provide a valuable foundation for future studies on the molecular mechanisms of heat stress responses in strawberry.

Expression patterns of FaHSFs in octoploid strawberry seedlings under heat stress.(A) Heatmap visualizing the expression profiles of 66 FaHSFs in response to heat stress at 6 h and 24 h. Clustering analysis results for gene expression levels, illustrated using different color scales, are displayed in the upper right corner. The data were extracted from transcriptome datasets and analyzed using TBtools. Pink arrows indicate FaHSFs upregulated after 6 h of heat stress, while red arrows denote FaHSFs that are upregulated at both 6 h and 24 h. (B) RT-qPCR validation of the expression levels of selected upregulated FaHSFs. Data are presented as means ± SD from three biological replicates. Distinct letters above the bars indicate statistically significant differences according to one-way analysis of variance (ANOVA) followed by Duncan’s multiple comparison test (P ≤ 0.05).

Discussion

Physiological and phenotypic changes of ‘Benihoppe’ octoploid strawberry under heat stress

In this study, the ‘Benihoppe’ variety was selected as the model for investigating the transcriptomic response of octoploid strawberry seedlings to heat stress. ‘Benihoppe’ is an early-maturing cultivar from Japan, widely cultivated in China and extensively used for summer seedling cultivation, this rationale underpins the selection of seedling materials for investigating heat stress responses in the present study. The recent high-quality genome assembly of ‘Benihoppe’ makes it an optimal model for molecular research into heat stress physiological and genomic adaptation in octoploid cultivated strawberry (Song et al., 2024).

Proline is a non-proteinogenic amino acid that accumulates in plants in response to abiotic stress. Previous studies have shown that proline stabilizes cellular membranes and protects proteins from denaturation by maintaining hydration shell stability, thereby enhancing plant tolerance to environmental stress (Hayat et al., 2012). Our results revealed a significant increase in proline content in strawberry seedlings after 3 h of heat stress, indicating that cellular osmotic balance is altered early under heat stress and that elevated proline is required to sustain cell turgor, stabilize membranes, and reduce electrolyte leakage, thereby improving stress resistance (Fig. 1A). Basal ROS levels are essential for cellular homeostasis, but ROS accumulates in various cellular compartments under heat stress (Mittler et al., 2022; Molina-Moya et al., 2025). In our study, SOD activity increased significantly following three hours of heat exposure, reflecting the enzyme’s critical role in rapidly catalyzing O2^−^ to H_2_O_2_ and O_2_, and suggesting ROS accumulation in strawberry seedlings by this time point (Fig. 1B). POD and CAT activities were notably elevated after 6 h of heat stress; POD catalyzes the detoxification of H_2_O_2_ and other hydroperoxides, whereas CAT efficiently decomposes H_2_O_2_ into water and oxygen (Figs. 1C and 1D). These findings imply that physiological injury becomes apparent by 6 h of heat treatment. Examining the root morphology, root tips exhibited pronounced blackening and necrosis after 24 h of heat stress. Combined with physiological data, this indicates severe heat-induced damage in the seedlings after 24 h (Fig. 1E). These observations justified our selection of 6 h and 24 h time points for subsequent transcriptomic analysis.

Analysis of DEGs

Comparison of DEGs between C6 vs. H6 and C24 vs. H24 revealed 3,425 genes in common (Fig. 2B). This suggests that distinct molecular responses to heat stress occur at 6 h and 24 h in strawberry seedlings, as further evidenced by GO and KEGG enrichment analyses.

GO enrichment analysis identified “Phenylpropanoid biosynthetic process”, “Flavonoid metabolic process” and “L-phenylalanine metabolic process” among the top 20 significantly enriched GO terms in both C6 vs. H6 and C24 vs. H24 (Fig. 3). HSFs not only activate heat shock protein genes but also regulate key genes within the phenylpropanoid biosynthetic pathway. Phenylpropanoid compounds provide an initial chemical defense in plants against heat-induced oxidative damage by functioning as potent antioxidants that scavenge ROS (Fu et al., 2022; Molina-Moya et al., 2025). Flavonoid biosynthesis, a major branch of the phenylpropanoid pathway, is a crucial component of the heat stress response. The phenylpropanoid pathway is derived from phenylalanine via the shikimic acid pathway (Rahim, Zhang & Busatto, 2023). The synthesis of flavonoids is a key branch of the phenylpropanoid pathway in response to heat stress. The synthesis of phenylpropane compounds originates from phenylalanine and is achieved through the shikimic acid pathway, linking these three metabolic pathways to heat stress adaptation (Dong & Lin, 2021). In the C24 vs. H24 group, three newly significant GO terms “cellular detoxification”, “hydrogen peroxide metabolic process”, and “positive regulation of oxidoreductase activity” (Fig. 3B, dark blue arrows) are associated with redox reactions, indicating that ROS levels increase substantially after 24 h of heat stress, necessitating enhanced redox activity to protect cells.

Among the top 20 pathways enriched in KEGG analysis, several were closely associated with heat stress response in both C6 vs. H6 and C24 vs. H24, including “Phenylpropanoid biosynthesis” “Glutathione metabolism” “Flavonoid biosynthesis” “Carbon fixation in photosynthetic organisms”, and “Flavone and flavonol biosynthesis” (Fig. 4). Glutathione is a principal plant antioxidant involved in detoxification, regulation of cellular signaling, and maintenance of the redox state under heat stress (Dard et al., 2023). High temperatures negatively impact carbon fixation efficiency during photosynthesis, reducing adenosine triphosphate (ATP) and nicotinamide adenine dinucleotide phosphate (NADPH) production and leading to photoinhibition. Additionally, stomatal closure under heat stress diminishes transpiration rates, further limiting CO_2_ uptake and decreasing carbon fixation (Rydzy et al., 2023). Therefore, “Carbon fixation in photosynthetic organisms” is vital for heat stress adaptation. The “Cutin, suberine and wax biosynthesis” pathway supports plant survival in elevated temperatures by minimizing water loss, reflecting solar radiation, mitigating heat load, and offering both physical and antioxidant protection (Li et al., 2022). The presence of “Fructose and mannose metabolism” and “Zeatin metabolism” exclusively among C24 vs. H24 pathways indicate that carbohydrate metabolism (for osmotic regulation and energy provision) and cytokinin-related pathways become more prominent during prolonged heat stress.

FaHSF expression analysis to verify the accuracy of transcriptome data in this study

We examined the distribution of transcription factors among the differentially expressed genes and found that HSFs ranked fourth in abundance (Fig. 5). HSFs represent essential regulators of the heat shock response (HSR). Previous studies identified only 9 HSFs in cultivated strawberries (comprising six HSFA subfamily members and three HSFB subfamily members, with two HSF members subsequently validated functionally in Arabidopsis), which may be attributed to technological limitations in sequencing methodologies at the time, and the current study provides enhanced guidance for experimental design in this field (Liao et al., 2016). Utilizing the high-quality ‘Benihoppe’ gene annotation and Illumina HiSeq platform (Song et al., 2024), our study identified 66 FaHSFs genes. In certain subgroups, the number of homologous genes identified for F. vesca and F. ananassa from two A. thaliana homologous genes is one and four, respectively, as observed in the HSFA6 and HSFA7 subgroups (F. vesca, like A. thaliana, is diploid, whereas F. ananassa is octoploid). This differential homology accounts for the identification of 17 homologous genes in F. vesca and 66 homologous genes in F. ananassa using 21 Arabidopsis HSF proteins as references (Fig. 6). Motif and domain analyses revealed that all 66 FaHSFs contain conserved HSF domains, suggesting functional conservation across A. thaliana, F. vesca, and F. ananassa (Fig. 7). Transcriptome data indicated that only five subfamilies were upregulated under heat stress at 6 h and 24 h (Fig. 8A, marked by red and pink arrows), and these upregulated genes are evenly distributed across various subgenomes (Fig. S5), suggesting that the remaining FaHSFs may be transiently responsive to earlier stress or act at later time points. RT-qPCR validation confirmed the upregulation of these five subfamilies under heat stress, supporting the robustness of our transcriptome data (Figs. 8B–8F).

Conclusions

Extreme high temperatures severely affect the cultivation of octoploid strawberry seedlings. However, the molecular mechanisms of heat stress response in strawberry seedlings remain unclear. In this study, ‘Benihoppe’ octoploid strawberry seedlings were subjected to treatments at 25 °C (control) and 40 °C (heat stress) for 1 h, 3 h, 6 h, 12 h, and 24 h, followed by a comprehensive analysis of physiological parameters and phenotypic changes at each time point. Notably, both proline content and SOD activity exhibited significant increases at 3 h after exposure to heat stress, whereas CAT activity and POD activity were significantly elevated at 6 h. The root tips of strawberry seedlings began to display visible cellular damage at 6 h of heat stress and showed pronounced tissue decay by 24 h. Based on these observations, transcriptome analysis was performed for samples collected at 6 h and 24 h under heat stress. A total of 11,526 DEGs were identified. The subsequent GO and KEGG enrichment analyses of DEGs revealed the involvement of multiple metabolic pathways associated with the heat stress response. Further analysis of DEGs encoding transcription factors identified pronounced differential expression of HSF family genes at both 6 h and 24 h. Through transcriptome analysis, five FaHSF subfamily genes (FaHSFA2, FaHSFA3, FaHSFA7, FaHSFB1, FaHSFB2) were found to be upregulated under heat stress, and this upregulation was further validated by RT-qPCR assay, supporting the reliability of the transcriptome data. Collectively, these findings provide a theoretical framework for the investigation of heat stress responses in strawberry seedlings and present potential molecular targets for the molecular breeding of heat-tolerant strawberry cultivars.

Supplemental Information

10.7717/peerj.20932/supp-1Supplemental Information 1Transcriptome profile of octoploid strawberry seedlings under heat stress(A) Overall distribution of gene expression levels across all samples. (B) Pearson correlation coefficient values among biological replicates within all sample groups. Control groups are represented as CK: C0, C6, C24.

10.7717/peerj.20932/supp-2Supplemental Information 2Volcano plots depicting the distribution of DEGs in C0 vs. C6 (A) and C0 vs. C24 (D)

10.7717/peerj.20932/supp-3Supplemental Information 3Hierarchical clustering DEGs in control (C0, C6, C24) and heat stress (H6 and H24)Each column corresponds to a specific comparison group, while each row represents a distinct gene. In the heatmap, red coloration signifies genes with high expression, and blue signifies genes with low expression.

10.7717/peerj.20932/supp-4Supplemental Information 4The number of transcription factors in the whole octoploid strawberry ‘Benihoppe’ genome

10.7717/peerj.20932/supp-5Supplemental Information 5The distribution of FaHSF genes on the ‘ Benihoppe ’ octoploid strawberry chromosomesVertical colored bars represent the chromosomes of ‘ Benihoppe ’ octoploid strawberry. The gene and chromosome name are shown at the right and left of each chromosome, respectively. The scale bar on the left indicateschromosome length .

10.7717/peerj.20932/supp-6Supplemental Information 6Analysis of collinearity in the HSF gene family(A) Intra-genomic collinearity of the FaHSF genes. (B) Collinearity of the HSF gene family between* F. vesca* and F. ananassa.

10.7717/peerj.20932/supp-7Supplemental Information 7Detailed sequence of the predicted motifs

10.7717/peerj.20932/supp-8Supplemental Information 8Predicted cis-elements of FaHSF genesEach element is indicated by a specific color.

10.7717/peerj.20932/supp-9Supplemental Information 9Summary of transcriptome sequencing data of strawberry seedlings

10.7717/peerj.20932/supp-10Supplemental Information 10Novel genes identified in the transcriptome

10.7717/peerj.20932/supp-11Supplemental Information 11DEGs of C6 vs. H6 selected from RNA-seq data

10.7717/peerj.20932/supp-12Supplemental Information 12Table S4. DEGs of C24 vs. H24 selected from RNA-seq data

10.7717/peerj.20932/supp-13Supplemental Information 13Table S5 GO enrichment analysis of biological process in C6 vs. H6

10.7717/peerj.20932/supp-14Supplemental Information 14GO enrichment analysis of biological process in C24 vs. H24

10.7717/peerj.20932/supp-15Supplemental Information 15KEGG enrichment analysis of C6 vs. H6

10.7717/peerj.20932/supp-16Supplemental Information 16KEGG enrichment analysis of C24 vs. H24

10.7717/peerj.20932/supp-17Supplemental Information 17Collinearity gene pairs in F. ananassa

10.7717/peerj.20932/supp-18Supplemental Information 18Collinearity gene pairs between F. vesca and* F. ananassa*

10.7717/peerj.20932/supp-19Supplemental Information 19The cis-regulatory elements predicted of HSF genes in F. ananassa

10.7717/peerj.20932/supp-20Supplemental Information 20The list of primers used in this study

10.7717/peerj.20932/supp-21Supplemental Information 21Raw data for RT-qPCR

10.7717/peerj.20932/supp-22Supplemental Information 22MIQE checklist

10.7717/peerj.20932/supp-23Supplemental Information 23The sequence of HSFs in A. thaliana, F. vesca, and F.ananas