Transcriptome Analysis Identifies Key Genes Involved in Response and Recovery to High Heat Stress Induced by Fire in Schima superba

Fire-resistant tree species play a crucial role in forest fire prevention, utilizing several physiological and molecular mechanisms to respond to extreme heat stress. Many transcription factors (TFs) and genes are known to be involved in the regulatory network of heat stress response in plants. However, their roles in response to high temperatures induced by fire remain less understood. In this study, we investigated Schima superba, a fire-resistant tree, to elucidate these mechanisms. Leaves of S. superba seedlings were exposed to fire stimulation for 10 s, 30 s, and 1 min, followed by a 24-h recovery period. Fifteen transcriptomes were assembled to identify key molecular and biological pathways affected by high temperatures. Differentially expressed genes (DEGs) analysis revealed essential candidate genes and TFs involved in the heat stress response, including members of the ethylene-responsive factors, WRKY, MYB, bHLH, and Nin-like families. Genes related to heat shock proteins/factors, lipid metabolism, antioxidant enzymes, dehydration responses, and hormone signal transduction were differentially expressed after heat stress and recovery, underscoring their roles in cellular process and recovery after heat stress. This study advances our understanding of plant response and defense strategies against extreme abiotic stresses.


Introduction
Global climate change has led to an increase in the frequency and intensity of forest fires, posing significant threats to ecosystems and human safety [1].The occurrence of wildfires not only results in the loss of biodiversity and land degradation but also has adverse effects on atmospheric environments [2].Fire-resistant tree species are commonly used as fire isolation belts and play a crucial role in mitigating the impact of wildfires on ecosystems and human communities.These species have developed mechanisms to endure and recover from the intense heat associated with fire-prone environments [3,4].Understanding the physiological and molecular basis of heat stress induced by high temperature in fire-resistant trees is essential for developing effective fire management strategies and selecting appropriate species for reforestation in fire-prone areas.
Fire-resistant tree species have a range of physiological and molecular adaptations that enable them to survive and recover from heat stress induced by high temperature [5].Thick bark acts as a heat shield, insulating the underlying tissues from extreme temperatures and preventing damage to the cambium layer.Additionally, some species have evolved strategies to regulate leaf temperature through transpiration and leaf orientation, reducing the risk of overheating and dehydration during fires [6,7].Other adaptations include the accumulation of heat shock proteins (HSPs) and antioxidants, which play crucial roles in protecting cellular structures and maintaining homeostasis under heat stress conditions [8,9].On the molecular level, these species activate a suite of genes involved in heat stress response pathways, including those coding for HSPs, which are regulated by heat stress transcription factors (HSFs) [10,11].Transcriptomic analyses have shown that genes encoding HSPs, chaperones, and antioxidant enzymes are up-regulated in response to heat stress [12,13].A variety of transcription factors (TFs) and signaling molecules regulate the expression of these heat stress response genes [14,15], providing valuable insights into the molecular mechanisms underlying heat tolerance in fire-resistant trees.
Schima superba Gardner & Champ, an evergreen broad-leaved tree in the Theaceae family, is commercially valued for its timber and fire protection [16,17].Known for its strong fire resistance, S. superba is widely used as a fire break in tropical and subtropical regions in China.Previous studies have explored its leaf structural characteristics in capturing smoke particles following forest fire [18].In the aspect of heat release rate (HRR) and total heat released (THR), S. superba has demonstrated superior performance, making it more suitable for use as a fuel break in southern China [19].Recent advances in understanding plant heat tolerance have highlighted various TFs and functional genes that regulate the physiological and molecular mechanisms activated under heat stress [20,21].The advent of RNA sequencing has facilitated the identification of numerous genes involved in the response to extreme high temperature in fire-resistant tree species, such as the heat stress response genes observed in Michelia macclurei Dandy following fire stimulation [22].Despite these advances, comprehensive insights into the mechanism underlying the response to extreme heat induced by high temperature remain limited.
This study aims to identify key genes involved in the response and recovery processes of S. superba under extreme heat stress, using transcriptome analysis of leaves subjected to heat stress and subsequently recover.The findings will contribute to understanding the molecular regulatory networks that underpin high temperature tolerance mechanisms in S. superba, providing a basis for developing new fire prevention strategies to enhance ecosystem resilience to extreme temperature conditions.

Plant Materials and Treatments
The seedlings of S. superba were planted in Nanjing Botanical Garden Mem.Sun Yat-Sen (Nanjing, China).Given the potential for this species to encounter high temperatures due to fire, we subjected the seedlings to heat stress treatment using the flame of an alcohol lamp, following the method outlined in our previous publication [22].To capture the dynamic changes in gene expression, fresh leaves were collected from the periphery of the flame after heat treatments of 10 s, 30 s, and 1 min.Untreated leaves served as the control group (CK), and leaves collected after a 24 h recovery period following the 1 min heat treatment were designated as the recovery group.The experimental groups were labeled CK, S1, S2, S3, or S4, corresponding to the control group, 10 s, 30 s, 1 min treatment groups, and the recovery group, respectively.All collected leaves were immediately frozen in liquid nitrogen and stored at −80 • C until RNA extraction.

RNA Extraction, Library Construction, and Transcriptome Sequencing
The total RNA was extracted using the Plant RNA Kit (R6827, Omega Bio-Tek, Norcross, GA, USA).RNA quality was assessed with a NanoDrop One Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and a Qubit 3.0 Fluorometer (Life Technologies, Carlsbad, CA, USA).The qualified RNA was used for cDNA library preparation with the MGIEasy RNA Library Prep Kit for BGI ® .First, mRNA was enriched by oligo (dT) magnetic beads and then the mRNA fragment was obtained by adding a fragmentation.The fragmented mRNA was used to synthesize first-strand and second-strand cDNA by PCR.The double-stranded cDNA then underwent end repair, adding an A tail to the 3 ′ end, sequencing adapters, purification, PCR amplification, and product circularization to construct the library.The fragment size and concentration of the library were detected using an Agilent 2100 Bioanalyzer, and the high-throughput platform of DNBSEQ-T7 was used for sequencing.

Transcriptome Assembly and Unigene Annotation
To ensure the reliability of the analysis, the raw sequencing data was processed by filtering out the low-quality reads, including those with more than five N bases, low-quality base numbers reaching 50% of reads, adapter-contaminated reads, and repetitive sequences caused by PCR amplification.The quality control of filtered data was performed using FastQC (version 0.11.9, default parameters).The clean reads were used for transcript assembly by Trinity (version 2.11.0, parameter -min_kmer_cov 2) [23].The cd-hit (version v4.8.1, parameters -c 0.85, -aS 0.85) was then used for transcripts clustering to obtain unigenes.BUSCO (Benchmarking Universal Single-Copy Orthologs) was used to assess the unigene integrity [24].TransDecoder software (version 5.5.0 parameters -m 50, -single_best_only) was then used for coding region prediction of unigenes (Haas et al. 2016).Seven databases were used for unigene annotation through sequence and motif similarity searches, which were Nr, Pfam [25], Uniprot [26], KEGG [27], GO (Gene Ontology) [28], KOG/COG [29], PATHWAY [30].Transcription factors were identified by the plant database PlantTFDB (version 5.0) [31].

Gene Expression Quantification and Differentially Expressed Genes Analysis
The number of reads mapped to each unigene for each sample were obtained using RSEM [32] and converted to FPKM (Fragments Per Kilobase per Million bases).FPKM represents the average number of fragments per kilobase length of a gene per million fragments and is the most commonly used method for estimating gene expression levels.Paired-end reads from the same fragment were counted as one fragment, and the expression quantification of the unigene was then calculated.To detect the repeatability between samples, the Pearson correlation coefficient was calculated based on the FPKM expression level of unigenes in each pair of samples.Differential expression analysis was performed using DESeq2 software (version: 1.26.0), with a screening threshold of p value < 0.05 and |log2FoldChange| > 1 [33].The differentially expressed unigenes were subjected to volcano plot analysis, cluster analysis, and enrichment analysis.Enrichment analysis included Gene Ontology and KEGG Pathway using the clusterProfiler software (version: 3.14.3)[34].The heatmaps of the DEGs were generated by TBtools-II [35].

Quantitative Real-Time PCR (qRT-PCR) Validation
To validate the transcriptome data, 16 DEGs involved in heat stress were randomly selected for a quantitative real-time PCR (qRT-PCR) experiment.The experiment included primer design, cDNA synthesis, qRT-PCR, and data analysis.The gene primers (Table S1) were designed by Primer3 Release 2.3.4 and synthesized by General Biosystems Co., Ltd.(Anhui, China).According to the previous study of reference genes in S. superba [36], SsuACT (actin) was chosen as the reference gene.First-strand cDNA was synthesized from 0.5 µg isolated RNA using the PrimeScriptTM RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China).SYBR ® Premix Ex Taq™ II (Tli RNaseH Plus) (TaKaRa, RR820A, Dalian, China) was used for qRT-PCR according to the protocol.The LightCycler 480 (Roche Molecular Biochemicals, Mannheim, Germany) was used for qRT-PCR experiment.Three biological and technical replicates were set for the quantification of each unigene.The relative expression level was calculated using the 2 −∆∆CT method [37].The standard errors of deviation were calculated using the STDEVA function in Excel.

Statistics of the Transcriptome Sequencing in S. superba
A total of 15 transcriptome libraries of S. superba were constructed in this study.In each library, more than six Gb clean bases and approximately 43 million clean reads were obtained (Table S2).The Q20 rate (the percentage of base with Phred value ≥20) was more than 97.979% and Q30 was more than 94.226%, with the GC content 44.214-45.374%(Table S2).After assembly, a total of 336,090 transcripts and 160,449 unigenes (Table S3) were obtained, with N50 lengths of 1682 bp and 1357 bp, respectively (Figure S1A,B).BUSCO (Benchmarking Universal Single-Copy Orthologs) analysis showed that the complete BUSCOs of transcripts and unigenes reached 97.4% and 82.6% respectively (Figure S1C,D).

Unigene Annotation and Transcription Factors (TFs) Identification
Unigenes were annotated across seven major databases, with 34,308 (21.38% of all) unigenes annotated in at least one database.The largest amount of unigenes were annotated in the Nr database (21.19%), followed by Uniprot (21.00%),GO (15.02%),KEGG (12.61%),Pfam (10.37%), and KOG (0.12%) (Table 1).GO annotation consisted of three categories: biological process (BP), cellular component (CC), and molecular function (MF), and the top 20 terms of each category are shown in Figure 1.Most unigenes (6638) were annotated in terms of the integral component of membrane that belongs to the cellular component.Three thousand four hundred twenty-six unigens were annotated in ATP binding term (molecular function category).In the biological process, the number of annotated unigenes was much smaller than in the other categories, where the most unigenes were annotated in terms of defense response (Figure 1).In KEGG pathway annotation, unigenes were classified into 22 terms across five categories (Figure 2).The organismal systems category contained only one term related to environmental adaptation, encompassing 2225 annotated unigenes.The metabolism category contained 11 terms and had the highest number of unigenes, with 4247 annotated in terms of global and overview maps and 1370 in carbohydrate metabolism.The genetic information processing category consisted of five terms, with 800 unigenes annotated in terms of folding, sorting, and degradation.In the category of environmental information processing, 1473 unigenes were annotated for signal transduction and 175 for membrane transport.In terms of transport and catabolism, 572 unigenes were annotated.In the KOG database, orthologous unigenes were classified into 22 groups.The top groups, based on the number of annotated unigenes, were general function prediction only, signal transduction mechanisms, and posttranslational modification, protein turnover, chaperones (Figure 3).Species homology comparison in the Nr database showed that 73.31% of S. superba unigenes were aligned to Camellia sinensis, followed by C. sinensis var.sinensis, accounting for 7.83% (Figure S2A).A total of 1504 TFs were identified and classified into 52 families (Table S4), with the top 20 families shown in Figure S2B.The majority of TFs were families of MYB and MYB-related AP2 and ERF, accounting for 10.509%, 10.008%, 9.59%, and 9.34% respectively (Figure S2B).In KEGG pathway annotation, unigenes were classified into 22 terms acros egories (Figure 2).The organismal systems category contained only one term environmental adaptation, encompassing 2,225 annotated unigenes.The metab egory contained 11 terms and had the highest number of unigenes, with 4,247 a in terms of global and overview maps and 1,370 in carbohydrate metabolism.Th information processing category consisted of five terms, with 800 unigenes ann terms of folding, sorting, and degradation.In the category of environmental inf processing, 1,473 unigenes were annotated for signal transduction and 175 for m  S4), with the top 20 families shown in Figure S2B.The majority of TFs were families of MYB and MYB-related AP2 and ERF, accounting for 10.509%, 10.008%, 9.59%, and 9.34% respectively (Figure S2B).52 families (Table S4), with the top 20 families shown in Figure S2B.The majority of TFs were families of MYB and MYB-related AP2 and ERF, accounting for 10.509%, 10.008%, 9.59%, and 9.34% respectively (Figure S2B).

Unigenes Quantification and Differentially Expressed Genes (DEGs) Analysis
The results showed that the correlation coefficient in major sample comparisons was above 0.9 (Figure 4A).Using the criteria for identifying differentially expressed genes (DEGs), the number of DEGs was identified between the two groups (Figure 4B).Comparisons between the groups revealed that, relative to the CK group, most unigenes were upand down-regulated in the S4 group.In the comparison between S3 and S4, the number of DEGs was similar to that in the comparison between CK and S4, but significantly higher than in the comparison between CK and S1, S1 and S2, and S2 and S3 (Figure 4B).Detailed information on DEGs is provided in Table S5.

Unigenes Quantification and Differentially Expressed Genes (DEGs) Analysis
The results showed that the correlation coefficient in major sample comparisons was above 0.9 (Figure 4A).Using the criteria for identifying differentially expressed genes (DEGs), the number of DEGs was identified between the two groups (Figure 4B).Comparisons between the groups revealed that, relative to the CK group, most unigenes were up-and down-regulated in the S4 group.In the comparison between S3 and S4, the number of DEGs was similar to that in the comparison between CK and S4, but significantly higher than in the comparison between CK and S1, S1 and S2, and S2 and S3 (Figure 4B).Detailed information on DEGs is provided in Table S5.To investigate the specific functional categories of DEGs, GO and KEGG pathway enrichment were conducted.In the GO enrichment analysis, the top 20 terms in each category are shown in Figures 5 and 6.Compared to the CK group, when the plant suffered a short high-temperature stimulation, the largest amount of unigenes changed in the cellular component category, such as cellular anatomical entity, membrane and intrinsic component of membrane (Figure 5A).Compared to short heat stress (10 s), when the leaves suffered 30 s stress, DEGs were significantly enriched in term of oxidoreductase activity (Figure 5B).As the heat stress duration increased to 1 min, a growing number of DEGs were enriched in terms of membrane, intrinsic component of membrane, and integral component of membrane (Figure 6A).During the recovery process, the dominant DEGs were enriched in the cellular component and molecular function, including the terms of cellular anatomical entity, catalytic activity, and membrane (Figure 6B).To investigate the specific functional categories of DEGs, GO and KEGG pathway enrichment were conducted.In the GO enrichment analysis, the top 20 terms in each category are shown in Figures 5 and 6.Compared to the CK group, when the plant suffered a short high-temperature stimulation, the largest amount of unigenes changed in the cellular component category, such as cellular anatomical entity, membrane and intrinsic component of membrane (Figure 5A).Compared to short heat stress (10 s), when the leaves suffered 30 s stress, DEGs were significantly enriched in term of oxidoreductase activity (Figure 5B).As the heat stress duration increased to 1 min, a growing number of DEGs were enriched in terms of membrane, intrinsic component of membrane, and integral component of membrane (Figure 6A).During the recovery process, the dominant DEGs were enriched in the cellular component and molecular function, including the terms of cellular anatomical entity, catalytic activity, and membrane (Figure 6B).
KEGG enrichment revealed the top 20 terms across four comparisons (Figure 7).The metabolism class was the predominant term that enriched in four comparisons.Within this class, the term of metabolism pathways was enriched in three comparisons and had the most DEGs, followed by the term of biosynthesis of other secondary metabolites (Figure 7A,B,D).The second major class was the environmental information processing, which contained two terms of MAPK signaling pathway and plant hormone signal transduction.These terms were enriched in three of the comparisons, excluding S1_vs_S2 (Figure 7A,C,D).The cellular processes class contained two terms of motor proteins and phagosome, with both enriched in the comparison of S3_vs_S4, while only the motor proteins term was enriched in comparison of S2_vs_S3.The plant-pathogen interaction term, part of the organismal systems class, was only enriched in the comparison of S2_vs_S3, and had 196 DEGs, representing the most unigenes.The term circadian rhythm was enriched in the comparison of S3_vs_S4 (Figure 7).

Differentially Expressed Genes Response to Strong Heat Stress in S. superba
To explore the key genes involved in the response and recovery process of high heat stimulation by fire in S. superba, differentially expressed key TFs and genes were identified across five sample groups.Eight main classes of TFs were summarized that expressed differentially in response to high temperature.Among these, ethylene-responsive transcription factor was predominant, with most genes up-regulated in the S4 group.These genes exhibited various expression patterns in the five sample groups.For example, some genes showed a gradual up-regulating trend from CK to S4, including unigene23512 (AP2/ERF), unigene15678 (ERF CRF4-like), unigene12438 (ERF RAP2-4 like), unigene34925 (ERF017-like), and unigene33602 (ERF5-like).Conversely, ERF113-like (unigene 57401) and two other ERFs (unigene49074 and unigene97535) had decreased expression with the increase of high-temperature stimulation time, but had significantly higher expression levels in plants that recovered for 24 h compared to the CK group.Three ERF ABR1-like unigenes did not show significant expression changes after hyperthermia stimulation compared with the CK group but were significantly up-regulated after 24 h of recovery.Two AP2like unigenes (unigenes6988 and unigene10710) showed the highest expression level in the S1 group and lowest expression in recovery plants (Figure 8A).After high-temperature stimulation recovery, WRKY family unigenes were up-regulated compared to the CK and S1-S3 groups (Figure 8C).The Nin-like family showed a similar expression pattern to the WRKY family, with the highest expression in recovery plants (Figure 8D).Unlike the WRKY and Nin-like families, the bHLH, MYB, TCP, bZIP and Trihelix families exhibited diverse expression patterns in CK and S4 groups (Figure 8E-I).Beyond the eight main TFs, some other TFs also showed high or low expression in recovery plants.For example, the NAC83-like, KUA1-like and SRM1 families were up-regulated in S4, but the genes of HD-ZIP, ATHB-13, ATHB-14, HBI1-like, NF-YB3, NF-YA3, NF-YA4, and YABBY2 were down-regulated in S4 (Figure 8B).
genes showed a gradual up-regulating trend from CK to S4, including unigene23512 (AP2/ERF), unigene15678 (ERF CRF4-like), unigene12438 (ERF RAP2-4 like), unigene34925 (ERF017-like), and unigene33602 (ERF5-like).Conversely, ERF113-like (unigene 57401) and two other ERFs (unigene49074 and unigene97535) had decreased expression with the increase of high-temperature stimulation time, but had significantly higher expression levels in plants that recovered for 24 h compared to the CK group.Three ERF ABR1-like unigenes did not show significant expression changes after hyperthermia stimulation compared with the CK group but were significantly up-regulated after 24 h of recovery.Two AP2-like unigenes (unigenes6988 and unigene10710) showed the highest expression level in the S1 group and lowest expression in recovery plants (Figure 8A).After high-temperature stimulation recovery, WRKY family unigenes were up-regulated compared to the CK and S1-S3 groups (Figure 8C).The Nin-like family showed a similar expression pattern to the WRKY family, with the highest expression in recovery plants (Figure 8D).Unlike the WRKY and Nin-like families, the bHLH, MYB, TCP, bZIP and Trihelix families exhibited diverse expression patterns in CK and S4 groups (Figure 8E-I).Beyond the eight main TFs, some other TFs also showed high or low expression in recovery plants.For example, the NAC83-like, KUA1-like and SRM1 families were up-regulated in S4, but the genes of HD-ZIP, ATHB-13, ATHB-14, HBI1-like, NF-YB3, NF-YA3, NF-YA4, and YABBY2 were down-regulated in S4 (Figure 8B).In addition to transcription factors, many important functional genes are involved in the biological processes of heat response and recovery after stimulation.These include genes related to heat shock protein genes, antioxidant enzyme genes, dehydrationresponsive genes, and lipid metabolism-related genes.The heat shock protein genes were significantly up-regulated in S4, including 24 HSPs and five HSFs (Figure 9A).Dehydrationresponsive genes also exhibited their highest expression in S4, except for unigene18792 (senescence/dehydration-associated protein), which reached the highest expression level and decreased after 24 h of recovery.Regarding antioxidant enzyme genes, those related to superoxide dismutase, ascorbate peroxidase, and iron superoxide dismutase showed decreased expression after high-temperature stimulation and recovery.However, L-ascorbate peroxidase genes were up-regulated.Catalase genes were significantly down-regulated after a 10 s fire stimulation but later up-regulated to a level close to the CK group (Figure 9E).Lipid metabolism, a complex biological process, involves many genes that are either up-or down-regulated in response to high temperature.For example, the lipid-binding protein gene and lipid transfer protein genes decreased expression in S4, while genes related to choline kinase, glycoside hydrolase, and lipoxygenase were up-regulated (Figure 9C).In plant hormone signal transduction, ethylene insensitive genes were up-regulated after heat stress and recovery.Conversely, many genes related to auxin-responsive protein, auxin-induced protein, and auxin-binding protein were down-regulated in S4.Genes of gibberellin 2-oxidase had the highest expression in S4, whereas gibberellin 20-oxidase genes showed decreased expression in S4.Abscisic acid (ABA) 8 ′ -hydroxylase showed different expression patterns in S4 (Figure 9D).
decreased expression after high-temperature stimulation and recovery.However, L-ascor-bate peroxidase genes were up-regulated.Catalase genes were significantly down-regulated after a 10 s fire stimulation but later up-regulated to a level close to the CK group (Figure 9E).Lipid metabolism, a complex biological process, involves many genes that are either up-or down-regulated in response to high temperature.For example, the lipidbinding protein gene and lipid transfer protein genes decreased expression in S4, while genes related to choline kinase, glycoside hydrolase, and lipoxygenase were up-regulated (Figure 9C).In plant hormone signal transduction, ethylene insensitive genes were upregulated after heat stress and recovery.Conversely, many genes related to auxin-responsive protein, auxin-induced protein, and auxin-binding protein were down-regulated in S4.Genes of gibberellin 2-oxidase had the highest expression in S4, whereas gibberellin 20-oxidase genes showed decreased expression in S4.Abscisic acid (ABA) 8'-hydroxylase showed different expression patterns in S4 (Figure 9D).

qRT-PCR Validation of DEGs
To validate the accuracy of RNA-Seq, 16 DEGs related to heat stress response were randomly selected for qRT-PCR.The results showed a high degree of consistency between RNA-Seq and qRT-PCR (Figure 10).For instance, genes such as HSF30 (Unigene21477), protein early-responsive to dehydration (Unigene8723), ABC transporter (Unigene28303), ABA 8 ′ -hydroxylase (Unigene14772) and WRKY31 (Unigene51504) exhibited significant up-regulation in the group that had undergone 24 h of recovery after heat stress.On the contrary, genes related to auxin-binding protein (Unigene32819), MYB4 (Unigene15957), and gibberellin 20-oxidase (Unigene8330) had the lowest expression levels in the recovery group (Figure 10).
ABA 8′-hydroxylase (Unigene14772) and WRKY31 (Unigene51504) exhibited significant up-regulation in the group that had undergone 24 h of recovery after heat stress.On the contrary, genes related to auxin-binding protein (Unigene32819), MYB4 (Unigene15957), and gibberellin 20-oxidase (Unigene8330) had the lowest expression levels in the recovery group (Figure 10).

Discussion
Plant responses and recovery to high heat stress involve many genes and biological processes, forming complex transcriptional regulatory networks [20,38].Transcriptome analysis in many plants, such as switchgrass [39], potato [40], and rice [41] has identified many genes related to high heat stress.S. superba, an important fire-resistant tree species, has previously been studied for its leaves' physiological and biochemical response to smoke by simulating forest fire [18].However, the molecular mechanisms underlying its response to high temperature remain underexplored.In this study, transcriptomic analysis of S. superba under fire-induced high temperature stress and subsequent recovery

Discussion
Plant responses and recovery to high heat stress involve many genes and biological processes, forming complex transcriptional regulatory networks [20,38].Transcriptome analysis in many plants, such as switchgrass [39], potato [40], and rice [41] has identified many genes related to high heat stress.S. superba, an important fire-resistant tree species, has previously been studied for its leaves' physiological and biochemical response to smoke by simulating forest fire [18].However, the molecular mechanisms underlying its response to high temperature remain underexplored.In this study, transcriptomic analysis of S. superba under fire-induced high temperature stress and subsequent recovery provides significant insights into the genetic and molecular mechanisms that enable this species to withstand and recover from such extreme temperature conditions.
GO and KEGG functional annotation of unigenes revealed the diverse biological processes, molecular functions, and cellular components involved in stress response and recovery.Differentially expressed genes exhibited significant changes between the control (CK) and recovery groups by heat stress in S. superba leaves.This suggests a sophisticated regulatory network that coordinates the plant's response to heat stress and subsequent recovery.HSPs, a main group of molecular chaperones, are highly conserved in plants and are essential for maintaining protein folding, stability, and cellular homeostasis under heat stress conditions [42][43][44][45].The identified HSPs encompass various families from many species, such as 21 HSFs encoding genes in Arabidopsis [46], 56 in wheat [47], etc.These HSPs perform diverse biological functions to prevent cellular toxicity under prolonged heat exposure.For example, overexpression of Arabidopsis HSFs can enhance tolerance to heat and other stresses [48].Additionally, 25 ZmHsf and 22 ZmHsp70 genes were identified as having potential effects on heat stress response [49], and 141 Hsf-like and Hsp-like genes were identified in eggplant, with most increasing after heat stress [50].In this study, upon exposure to fire and 24 h of recovery, 24 HSPs and five HSFs exhibited pronounced up-regulation after recovery, suggesting their pivotal role in protecting cellular proteins and reflecting the immediate adaptive response to mitigate the detrimental effects of thermal stress on cellular components.
The activation of oxidative metabolism has been shown to overlap with the transcription of HSFs and HSPs during heat stress, indicating its integral role in the heat response and enhanced heat tolerance of plants [51].Under high temperature stimulation, some specific antioxidant enzyme genes will be activated.Enzymes such as superoxide dismutase (SOD), peroxidase (POD), ascorbate peroxidase (APX), and catalase (CAT) typically respond more significantly.For instance, these enzymes were found to be increased in the heat-tolerant genotypes of Cucurbita moschata (Duchesne ex Lam.) Duchesne ex Poir.[52].Similarly, hightemperature treatment in mulberry cultivars led to elevated antioxidant enzyme activities, including SOD, CAT, POD, APX, and glutathione reductase (GR) [53].Over-expressed CgHSP70 plants showed higher POD activity and proline content in Chrysanthemum [54], suggesting that these enzymes could help clear reactive oxygen in plants and reduce oxidative damage after heat stress.In S. superba, CAT genes, which are essential for hydrogen peroxide breakdown, were initially down-regulated after heat stress, but increased in the longer heat stimulation and recovery groups.The genes of SOD and ascorbate peroxidase were down-regulated after extreme heat stress and decreased to a much lower expression than the CK group, while L-ascorbate peroxidase genes showed high expression in the recovery group.This indicates that CAT and L-ascorbate peroxidase genes facilitate the clearance of reactive oxygen after heat stress and reach a balance with SOD and POD to manage oxidative stress.Together, these enzymes are regulated to protect cells from ROS damage and facilitate recovery after exposure to extremely high temperatures.
Reactive oxygen species induced by heat stress can cause damage to lipids and proteins in plants, which affects the stability of plant cell membranes.This triggers changes in the expression of some related genes in lipid metabolism pathways to maintain membrane integrity and function [55,56].This study identified a series of lipid metabolism-related genes showing a complex response to high temperature.Genes related to lipid-binding protein, lipid transfer protein, phospholipase A, and fatty acid desaturase decreased after long heat stress and reached their lowest expression after 24 h of recovery.However, genes for choline kinase, glycoside hydrolase, and lipoxygenase were up-regulated after recovery from heat stress.In grass species, heat stress enhances the accumulation of phospholipids and glycolipids [57], suggesting that the ability of lipid binding and transfer decrease after heat stress.The up-regulated genes may participate in the biosynthesis of phospholipids and glycolipids during lipid metabolism, promoting tolerance and recovery potential from extreme heat stress in S. superba.During the recovery phase, gene expression towards the restoration of cellular functions is disrupted by heat stress.Genes associated with photosynthesis and carbon metabolism, such as calcium-binding related genes, transmembrane protein, and mitochondrial protein genes showed gradual recovery.This indicates the plant's ability to rebuild energy reserves and metabolic pathways essential for growth and development.This adaptive response likely supports the regeneration of damaged tissues and enhances overall resilience following high temperatures.
In addition, it has been proven that all major plant hormones play a critical role in response to heat stress, including auxin, abscisic acid (ABA), gibberellins, and ethylene [58].ABA is a key hormone involved in abiotic stress responses.In the context of heat stress, the ABA regulatory network includes several HSFs, Dehydration-Responsive Element Binding Protein 2A (DREB2A), bZIP, MYB, and NAC proteins.DREB2A induces the expression of HSFA3 under high temperatures, thereby activating many stress-related genes [14,59].In S. superba, dehydration-responsive genes closely related to ABA were significantly up-regulated after heat stress, and the expression patterns of ABA 8 ′ -hydroxylase, which is involved in ABA catabolism, varied during the recovery phase (S4).Differential expression of ABA-related genes indicates a complex regulatory mechanism.Increased ABA-related gene expression after heat stress can enhance stress tolerance by promoting stomatal closure and reducing water loss, while reduced levels during recovery can facilitate growth resumption.This dynamic regulation of ABA metabolism has been observed in other plants, such as Triticum aestivum L., where ABA levels are closely regulated during drought and heat stress to optimize survival and growth [60].Gibberellin-related genes also showed different expression patterns; gibberellin-2-oxidase genes had the highest expression in S4, while gibberellin 20-oxidase genes decreased in expression.In plants, high temperature also commonly suppresses GA and IAA biosynthesis [61,62].Most of the auxin-response and auxin-binding genes were down-regulated in S4, which may have a common function with decreased ABA genes to the response of high temperature in S. superba.
In plants, various transcription factors and functional genes play crucial roles in responding to environmental stresses.In this study, a total of 1504 TFs were identified and classified into 52 families, with ERF, WRKY and MYB being the most prominent.In other plants, families such as AP2/ERF, MYB, WRKY, and bZIP have been shown to play significant roles in heat stress responses, as seen in rice [63] and fire-resistant tree species M. macclurei [22].Ethylene responsive factors (ERFs) are known to mediate plant responses to abiotic stresses by regulating genes involved in ethylene signaling, stress tolerance, and recovery processes [64,65].In S. superba, most ERFs exhibited a gradually increasing pattern from CK to S4, peaking after 24 h of recovery, accompanied by increasing expression of ethylene-related genes.This indicates their involvement in stress response and recovery.The WRKY family showed up-regulation in the recovery phase compared to CK and S1-S3 groups, suggesting their role in restoring normal cellular functions post-stress.Nin-like (NLP) TFs exhibited a similar expression pattern to the WRKY family, with the highest expression in recovery plants but low expression in long heat stress of S3 group, indicating a potential collaborative role in stress recovery.In Capsicum annuum (L.), the expression of CaNLP1 increased after various abiotic stress [66], which is consistent with the expression in S. superba's response to heat stress.Other families including MYB, bHLH, bZIP, TCP, Trihelix, NAC and NF-Y exhibited diverse expression patterns in CK and S4 groups.These TFs control diverse biological processes in plants, such as development and abiotic stress responses [67,68].For example, the expression of TCP gene changes in response to abiotic stresses including heat, drought, and salt in ginger [69].Under heat stress treatment, PgMYB2, PgMYB9, PgMYB88 and PgMYB151 are differentially expressed in Pennisetum glaucum (L.) R. Br. [70].The diversity of these TFs in S. superba reflects their varied roles in different stages of heat response and recovery, potentially regulating a broad spectrum of stress-responsive genes.

Conclusions
In this study, 15 transcriptomes of S. superba subjected to heat stress and 24 h of recovery were assembled, revealing the primary molecular and biological pathways of response to high temperature in this fire-resistant tree species.Analysis of differentially expressed genes identified key candidate genes and transcription factors involved in the heat stress response in S. superba.Ethylene-responsive factors, WRKY, MYB, bHLH, and Nin-like families play a central role in regulating the networks of heat stress response and adaption, while other TF families contribute to different aspects of stress.Functional genes related to heat shock proteins/factors, lipid metabolism, antioxidant enzymes, dehydration response, and hormone signal transduction factors are crucial for maintaining cellular homeostasis during and after exposure to high temperatures.Changes in several hormonalrelated genes indicate more complexity and stress response mechanisms.These findings provide comprehensive insights into the molecular mechanisms governing the response and recovery of S. superba to high temperatures.The identified key genes not only enhance our understanding of plant defense and resilience strategies but also offer potential genes

Figure 2 .
Figure 2. Kyoto Encyclopedia of Genes and Genomes (KEGG) classification of S. superba.

Figure 2 .
Figure 2. Kyoto Encyclopedia of Genes and Genomes (KEGG) classification of S. superba.

Figure 2 .
Figure 2. Kyoto Encyclopedia of Genes and Genomes (KEGG) classification of S. superba.

Figure 4 .
Figure 4. Pearson analysis of sample correlation (A) and the number of differentially expressed genes (DEGs) in different comparison groups (B).(A) The horizontal and vertical axes represent each sample respectively.The color depth represents the size of the correlation coefficient between the two samples.The closer the color is to red, the better the correlation between the two samples.The closer the correlation coefficient is to 1, the higher the similarity of the expression patterns between the samples.(B) The horizontal axis represents different comparisons, the vertical axis represents the number of unigenes, and red and green represent up-and down-regulated unigenes, respectively.

Figure 4 .
Figure 4. Pearson analysis of sample correlation (A) and the number of differentially expressed genes (DEGs) in different comparison groups (B).(A) The horizontal and vertical axes represent each sample respectively.The color depth represents the size of the correlation coefficient between the two samples.The closer the color is to red, the better the correlation between the two samples.The closer the correlation coefficient is to 1, the higher the similarity of the expression patterns between the samples.(B) The horizontal axis represents different comparisons, the vertical axis represents the number of unigenes, and red and green represent up-and down-regulated unigenes, respectively.

Genes 2024 , 19 Figure 5 .
Figure 5.The GO enrichment of differentially expressed genes in S. superba.(A,B) represent the comparisons of CK_vs_S1 and S1_vs_S2.

Figure 5 .
Figure 5.The GO enrichment of differentially expressed genes in S. superba.(A,B) represent the comparisons of CK_vs_S1 and S1_vs_S2.

Figure 6 .
Figure 6.The GO enrichment of differentially expressed genes in S. superba (A,B) represent the comparisons of S2_vs_S3 and S3_vs_S4.

Figure 6 .
Figure 6.The GO enrichment of differentially expressed genes in S. superba (A,B) represent the comparisons of S2_vs_S3 and S3_vs_S4.

Figure 7 .
Figure 7.The KEGG enrichment of differentially expressed genes in S. superba.(A-D) represents the different comparisons between samples.

Figure 7 .
Figure 7.The KEGG enrichment of differentially expressed genes in S. superba.(A-D) represents the different comparisons between samples.

Figure 8 .
Figure 8. Heatmap of the differentially expressed transcription factors involved in heat response and recovery in S. superba.Red means up-regulation, and white means down-regulation.(A), (C-I) represent different TF families, and (B) contains several TFs that changed expression levels in response to heat stress.

Figure 8 .
Figure 8. Heatmap of the differentially expressed transcription factors involved in heat response and recovery in S. superba.Red means up-regulation, and white means down-regulation.(A), (C-I) represent different TF families, and (B) contains several TFs that changed expression levels in response to heat stress.

Figure 9 .
Figure 9. Heatmap of the differentially expressed functional genes related to heat response and recovery in S. superba.(A) Heat shock protein genes; (B) Other genes; (C) Lipid metabolism-related

Figure 9 .
Figure 9. Heatmap of the differentially expressed functional genes related to heat response and recovery in S. superba.(A) Heat shock protein genes; (B) Other genes; (C) Lipid metabolism-related genes; (D) Hormone signaling; (E) Antioxidant enzyme genes; (F) Dehydration-responsive genes.Red means up-regulation, and blue means down-regulation.

Figure 10 .
Figure 10.qRT-PCR validation of 16 DEGs involved in heat stress response in S. superba.The xaxis represents the five sample groups, and the y-axis shows the FPKM by RNA-Seq and the relative quantitative expression level by qRT-PCR for each unigene.

Figure 10 .
Figure 10.qRT-PCR validation of 16 DEGs involved in heat stress response in S. superba.The x-axis represents the five sample groups, and the y-axis shows the FPKM by RNA-Seq and the relative quantitative expression level by qRT-PCR for each unigene.

Table 1 .
Statistics of unigene annotation in Schima superba.