Genome sequences were downloaded from published celery genome database (http://apiaceae.njau.edu.cn/celerydb and http://celerydb.bio2db.com/) for constructing the local BLAST (Li et al., 2022b). The hidden Markov model profile of the HSF domain (PF00447) from the Pfam database (http://pfam.xfam.org/) was downloaded and utilized to identify the HSF members in celery using the HMMER 3.0 program with default parameters. Candidate proteins were further submitted to NCBI-CDD (https://www.ncbi.nlm.nih.gov/cdd) to confirm the conserved domains. The physicochemical properties of the HSF genes identified from celery were analyzed using online ExPASy software (http://web.expasy.org/protparam). The subcellular localization was predicted with WOLF PSORT (https://www.genscript.com/tools/wolf-psort) and SoftberryProtComp 9.0 (http://www.softberry.com). The HSF sequences of Arabidopsis were downloaded from TAIR (https://www.arabidopsis.org/). ClustalX (v2.1) was used to align multiple sequences of all HSFs with default parameters, and phylogenetic tree was constructed by MEGA7.0 using neighbor-joining (NJ) method with 1000 bootstrap replicates.

The chromosomal location of AgHSFs was visualized through TBtools software (Chen et al., 2020). The exon/intron structure of AgHSFs was analyzed by online tool Gene Structure Displayer Server v2.0 (http://gsds.gao-lab.org/). Conserved motifs were analyzed using MEME Suite v5.4.1 with the parameters set as follows: the maximum number of motifs, 20; and the optimum width of each motif, between 10 and 50 residues.

The cis-regulatory elements in the 1.5 kb promoter regions of AgHSFs were predicted using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html) (Lescot et al., 2002). The STRING software was used to output the protein interaction value data between AgHSFs and other proteins in celery genome, and established protein interaction network. Results were visualized with Cytoscape software v3.7 (Shannon et al., 2003).

The heat-tolerant variety ‘Jinnan Shiqin’ and the non-heat-tolerant variety ‘Liuhe Huangxinqin’ were used as plant materials. Celery seeds were placed in an incubator at 20 °C in dark for germination, and then transferred them into growth chambers at 25 °C/20 °C (14 h light/10 h dark cycles). When the plants reached the stage of about 10 leaves, they were transferred into 38 °C for 24 h heat stress. The celery leave samples were collected after at 0 h (QC1), 4 h (QC2), 12 h (QC3), and 24 h (QC4) of treatment. Plant RNA samples were extracted using a total RNA kit (Tiangen, Beijing, China), and then converted into cDNA using PrimeScript RT Kit (TaKaRa, Dalian, China). The sequences of AgHSF genes were amplified for RT-PCR using EasyTaq DNA Polymerase Kit (TransGen Biotech, Beijing, China). The PCR product was ligated with pEASY-T1 vector and sequenced in Tsingke Biotechnology Co., Ltd. (Beijing, China).

The full-length of AgHSFa6-1 was amplified using specific primers ( Table S1 ), and inserted into green fluorescent protein (GFP)-fusion expression vector pA7. The expression vector 35S:AgHSFa6-1:GFP was constructed with two BamH I restriction enzyme sites. The empty vector containing 35S::GFP protein served as a control. The fusion expression and empty vector were transformed in tobacco leaves by DNA particle bombardment system (PDS-1000, Bio-Rad, USA). The tobacco was firstly disposed in dark for 24 h, fluorescence was then observed using a laser confocal microscope (LSM800, Zeiss, Germany).

The full-length open reading fragments of AgHSFa6-1 were amplified by PCR using specific primers, and inserted into the pYES2 vector ( Table S1 ). Then the recombinant vector was transformed into Saccharomyces cerevisiae strain INVSC1, and positive transformants were selected in the SD-Ura medium at 30 °C. For heat stress treatment, the positive yeasts were cultured on the SG-Ura medium for 3 d at 30 °C, 39 °C, 41 °C and 45 °C, respectively.

For the overexpression of AgHSFa6-1, the double-cauliflower mosaic virus 35S (CaMV 35S) promoter and NOS terminator from pSAT1-cEYFP-N1 were amplified and subcloned into the vector pCAMBIA1301. To obtain transgenic Arabidopsis, the amplification product was inserted downstream of the double CaMV 35S promoter into the modified pCAMBIA1301 vector between BamH I and Sac I sites. The construct contains the CaMV 35S promoter driving AgHSFa6-1 and the CaMV 35S promoter driving the Hyg gene for hygromycin resistance as a selectable marker. The pCAMBIA1301 vector also contained a GUS reporter gene following CaMV 35S promoter. The recombinant vector was introduced into the Agrobacterium tumefaciens strain GV3101 through electroporation. Then, it was transformed into Arabidopsis Columbia (Col-0) wild-type (WT) plants by floral-dip method (Zhang et al., 2006). The transgenic T₀ lines were selected using hygromycin on MS medium, and the hygromycin-resistant plants were selected out to obtain T₁ seeds for PCR assays. Then, consistently selected twice to harvest T₂ seeds of transgenic Arabidopsis, and the T₃ homozygous lines were used for further verification.

Seeds of Arabidopsis were placed on 1/2 MS medium to germination. The seeds were vernalized in a dark incubator at 4 °C for three days and then changed the growth condition at 25 °C/15 °C (16 h light/8 h dark cycles) and 70% relative humidity. Ten-day-old seedlings were subjected to heat stress for three days at 38 °C/25 °C (16 h light/8 h dark cycles), and then measured the primary root length. Two-month-old seedlings were subjected to 38 °C/25 °C heat stress for 24 h. They were sampled at 0 h, 4 h, 12 h, and 24 h for RNA extraction utilized for subsequent RT-qPCR analysis. Moreover, leaf-samples of 4 h were collected for stomata analysis and samples of 24 h were collected for physiological measurements. Then plants were domesticated at 42 °C for two days to compare their phenotypic difference.

The free proline content was measured according to the acidic-ninhydrin-based colorimetric method (Bates et al., 1973). The activities of superoxide dismutase (SOD) and peroxidase (POD) were measured using the nitrogen blue tetrazolium (NBT) and the guaiacol methods, respectively (Macadam et al., 1992; Cao et al., 2018). The concentration of malondialdehyde (MDA) in leaves was measured by the TBA colorimetric method (Li et al., 2020). For stomata analysis, the plant samples were immersed into FAA fix solution, dehydrated, embedded in paraffin, sectioned, and stained with toluidine blue. The width and length of the stomatal aperture were analyzed using ImageJ software v1.8.0.112 (Graphics software; NationalInstitutes of Health: Wayne Rasband, USA, 1997).

During above experiment, three independent replicates of each transgenic lines were performed. All of the data were analyzed using SPSS 26.0 software (SPSS Inc., Chicago, IL, USA), and the significance levels were based on P < 0.05 (*) and P < 0.01 (**).

The transcriptome data of celery under variable high temperature was obtained by this research group (Li et al., 2022a), and the transcript abundance of AgHSF was calculated in log2(FPKM) to generated cluster heatmaps.

RT-qPCR was performed using SYBR Premix Ex Taq (TaKaRa, Dalian, China). All the steps were followed the manufacturer’s instruction (CFx384TM Real-Time System, Bio-Rad, USA), and the expression level was calculated by the 2^-△△Ct method (Pfaffla, 2001). AgTUB and AtACT2 were served as reference genes for celery and Arabidopsis genes, respectively. All gene primer sequences used in this study were designed using Primer Premier software (version 6.0; Premier Biosoft International: Palo Alto, CA), and the sequences are listed in Table S1 .

A total of 29 AgHSF genes were identified and re-named from AgHSFa6-1 to AgHSFc1-1 ( Table S2 ). The results demonstrated that almost all the genes were localized in the nucleus. The polypeptides of AgHSF proteins ranged from 232 (AgHSFa9-2)-503 aa (AgHSFa1-3), with the molecular weight ranging from 26450.87 (AgHSFa9-2)-55586.82 Da (AgHSFa1-3). The aliphatic amino acids accounted for 18%, whereas the aromatic amino acids accounted for only 8%. The theoretical isoelectric point (pI) of proteins belongs to class A ranged between 4.83-6.25, indicating that they were acidic. However, except for B2, the pI values of class B proteins were more than 7, suggesting they were alkaline. All AgHSF proteins had negative GRAVY values, which means that all proteins were hydrophilic.

A phylogenetic tree was created between 29 celery and 20 Arabidopsis HSF proteins to determine their categorization and explore their evolutionary relationship ( Figure 1 ). The result revealed that the relationships of HSF genes between the two species were highly homologous. The AgHSFs in celery were classified into three classes and 14 subgroups. The class A had the most abundant proteins, which were clustered into nine subgroups, corresponding to HSFA1-HSFA9, and class B were clustered into four subgroups, corresponding to HSFB1-HSFB4. In contrast, the class C only had one subgroup for one gene, corresponding to HSFC1.

The evolution of gene family members was mainly determined by the variation of their conserved motifs and the diversity of gene structures (Hu et al., 2010). Results revealed that the longest intron was nearly 4 Kb (AgHSFa5-2), and all AgHSF genes had between two and four exons and one to three introns ( Figure 2A ). However, the types and numbers of motifs varied considerably ( Figure 2B ). A total of 20 motifs were detected in AgHSF proteins; the sequence of each motif is listed in the supplementary table ( Table S3 ). All 29 AgHSFs shared only one motif type (motif 4), and the number of motifs varied from two to 11. Similarities were more prevalent among the same groups or classes. AgHSFs belongs to Class A usually had more motif types, and except for AgHSFa9-2, they had motifs 1 and 3 in common, and members belongs to class B had motifs 1, 4, 5, and 8 in common. Moreover, motifs were always arranged in a particular order, such as motifs 2- 1- 4- 3 and 9. The highly conserved exon-intron distribution and motif divergence indicated that AgHSFs evolved according to a group, establishing the basis for the functional diversity of different classes.

To identify the localization of AgHSF genes, the TBtools software was used to visualize the location of these genes on different chromosomes in the celery genome ( Figure 3 ). Result demonstrated that except for Chr01, 28 AgHSF genes were unevenly distributed between Chr02 and Chr11, whereas AgHSFa8-1 was located at an unanchored contig. The Chr04 and Chr10 had the largest number of genes (four genes), followed by Chr02, 06, 09, and 11 with three genes each. Genes belongs to class A were located on each chromosome from Chr02 to 11, while the class B genes were only present on Chr02, 03, 04, 05, and 10, and the C class genes were only located on Chr10. Some AgHSF genes demonstrated high cross collinearity, clustering at similar sites on the same chromosome. For instance, AgHSFa3-1 and AgHSFa4-3 were localized at a similar site on the same chromosome, revealing that they might have sequence similarity.

Transcription factors can regulate gene expression by binding to cis-regulatory elements in response to different biotic or abiotic stress signals (Li et al., 2021). Various cis-regulatory elements were predicted in promoter sequences of all the AgHSF genes, and they were divided into three groups based on their functions: plant growth and development, phytohormone responsive, and abiotic and biotic stress ( Figure 4B ). Result revealed that some elements were widely distributed in promoter regions of AgHSFs, with light response elements (Box 4, G-box, and GT1-motif) being the most frequently distributed, followed by antioxidant response elements (ARE) and stress-response element (STRE). Furthermore, abscisic acid response elements (ABRE) were also detected in promoter regions of several AgHSFs ( Figure 4A ), indicating these genes can regulate the abscisic acid (ABA) secretion or synthesis process in response to heat shock.

To further investigate the regulatory mechanism of HSF genes in celery, a protein interaction network comprising the homologous proteins in Arabidopsis was constructed ( Figure 5 ). Except for AgHSFa9-2, 28 AgHSF proteins and 15 other proteins in the celery genome were present in this interaction network. Out of these, three HSP proteins constituted the central node and were connected with multiple AgHSF proteins. AgHSF proteins also exhibited significant relationships with DREB, PBS, MPK, and SQN proteins. HSF genes were proved to respond to heat shock by regulating the expression of HSP (Jacob et al., 2017). In our study, members of HSF family and HSP formed an expansive regulation network, indicating that these genes may play important roles in sensing and responding to heat stress.

According to the transcriptome data, all the 29 AgHSF genes were expressed after 0 h, 4h, 12 h, and 24 h of heat stress, referring to QC1, QC2, QC3, and QC4, respectively. Previous studies demonstrated that HSFA2 genes in Arabidopsis and tomato could be induced by heat stress, resulting in enhanced permeability and antioxidant capacity of plants (Meiri et al., 2010; Fragkostefanakis et al., 2016). Several AgHSFs, including AgHSFa2-2, AgHSFa2-3, and AgHSFa6-1, displayed a typical trend of initially increasing and declining. The expression levels of AgHSFa6-1 increased significantly at QC2 and reached their peak at this stage but decreased during QC3 and QC4 ( Figure 6A ). Despite this, their expression levels were significantly higher than any other AgHSF genes.

Nine highly expressed AgHSFs (AgHSFa2-1, AgHSFa2-2, AgHSFa2-3, AgHSFa6-1, AgHSFa6-2, AgHSFb1-1, AgHSFb1-2, AgHSFb2-1, AgHSFb2-2) were further detected by RT-qPCR to validate their expression levels in two celery cultivars (heat-tolerant celery variety ‘Jinnan Shiqin’ and the non-heat-tolerant celery variety ‘Liuhe Huangxinqin’) ( Figure 6B ). The expression levels of nine AgHSFs were all exhibited higher expression level in ‘Jinnan Shiqin’ after 4, 12, 24 h heat stress, especially AgHSFa6-1 exhibited a thousand times higher expression level than CK after 4 h heat stress. This result revealed that AgHSFa6-1 were involved in the response to heat stress, and the heat-tolerant celery variety could response to heat stress more rapidly to adapt to high temperature. Overall, similar expression trends in gene expression were detected from transcriptome data and RT-qPCR analyses.

Result demonstrated that the control was strongly expressed in the cytomembrane and nucleus, whereas 35S:AgHSFa6-1:GFP was detected only in the nucleus, indicating that the AgHSFa6-1 protein was targeted in the nucleus ( Figure 7 ). The result was consistent with previous subcellular localization predictions based on protein structures and demonstrated that AgHSFa6-1 could execute the function of a transcription factor.

The yeasts were cultured on SG-Ura medium at 30, 39, 41, and 45 °C for three days ( Figure 8A ). Result demonstrated that the control strain only grew at 30 °C, while the transgenic strains grew significantly at 39 and 41 °C. However, at the temperature of 45 °C, neither the control nor transgenic strains grew, suggesting that the AgHSFa6-1 gene responds to temperature stresses in yeast.

The transgenic Arabidopsis plants were obtained using hygromycin selection, and their validity was verified by PCR and GUS staining ( Figure S1 ). A total of 25 transgenic lines were obtained and designated as from OE-1 to OE-25. Two lines of transgenic lines (OE-17 and OE-19) and WT plants were selected to verify the AgHSFa6-1 function in response to heat stress. The morphological and root-length differences were measured using ten-days-old and two-month-old seedlings. As illustrated in Figures 8B, C , under normal growth conditions, OE-17, OE-19, and WT grew consistently, whereas after 24 h of continuous heat stress at 38/25 °C, the WT plants visibly withered to death, but the transgenic seedlings were still alive. Additionally, the growth of WT seedlings was stalled by heat stress, whereas OE-17 and OE-19 seedlings survived, and the average root length of transgenic plants was significantly longer than that of WT ( Figures 8C, G ).

Stomata are the main channels for gas exchange, and plants normally adjust the stomatal conductance to vary their temperature (Wang et al., 2022). When exposed to 4 h of heat stress, the conductance of WT plants was more open than that of transgenic plants, and both the aperture size and length of WT stomata were significantly larger than those of transgenic stomata ( Figures 8D–F, H, I ). These findings demonstrated that over-expression AgHSFa6-1 enhanced the heat-endurance ability of transgenic plants by regulating stomatal conductance and opening to control water loss.

Leaf samples from two-month-old Arabidopsis plants were used to measure all physiological parameters. The physiological indices of WT and transgenic plants under heat stress were determined, including proline, MDA, soluble protein, and antioxidant enzymes (SOD and POD). Results demonstrated that the content of proline in transgenic plants were significantly higher than in WT plants, while the MDA content displayed the opposite trend ( Figures 9A, B ). Proline is a signal molecule that regulates osmotic pressure, and MDA reflects the degree of membrane damage (Liu et al., 2016; Fu et al., 2018). Examining soluble protein content revealed that the average values of OE-17 and OE-19 were higher than those of WT plants, consistent with the proline level ( Figure 9C ). Studies indicated that heat stress could promote the secretion of ROS, and antioxidant enzymes like SOD and POD have ROS-scavenging mechanisms (Wang et al., 2003; Gechev and Petrov, 2020). SOD and POD in transgenic plants were significantly higher than in WT plants following a heat treatment ( Figures 9D, E ). These results suggest that transgenic plants are more tolerant to high-temperature stress than WT plants due to overexpression of AgHSFa6-1.

To understand the regulatory function of AgHSFa6-1 in wild-type and transgenic Arabidopsis plants under heat stress, the expressions of AgHSFa6-1 and seven stress-responsive genes in two lines of transgenic Arabidopsis were measured using RT-qPCR (Figures 9F–N ). Both AgHSFa6-1 and stress-responsive genes showed similar expression trend in OE-17 and OE-19. The expression levels of seven stress-responsive genes showed that AtHSP98.7 and AtAPX1 at all stages of heat stress were higher than those of CK. However, expression levels of BOB1, AtCPN60B, AtHSP70-1, AtADH2, and AtGOLS1 were insignificantly higher than WT during 0-4 h, but expression levels of AtCPN60B, AtHSP70-1, AtADH2 were even lower than WT after 12 h. The expression levels of AtHSP98.7 initially increased between 0 and 4 h, then decreased between 4 and 8 h, and peaked at 24 h. On the contrary, the expression levels of AtCPN60B and AtADH2 were initially decreased and then increased. These genes were significantly upregulated by high-temperature stress and reached their peak expression at 24 h. Moreover, their expression levels were much higher in the transgenic plants than in WT plants. Expression trend of AgHSFa6-1 in OE-17 and OE-19 was consistent with both transcriptome data and expression analysis in celery, and transcription and translation were relatively faster compared to the other seven genes. This revealed that AgHSFa6-1 improved the resistance of transgenic Arabidopsis to heat stress by positively regulating the expression of abiotic stress-sensitive genes.