The Hidden Markov Model (HMM) profile of the HSP20 domain (PF00011) was taken from Protein family database (Pfam 34.0; http://pfam.xfam.org/). The Cucurbita moschata Genome Database (CuGenDB, http://cucurbitgenomics.org/) was searched for this profile using BlastP methods with a cut-off E-value <10⁻⁵ (Yu et al., 2016). The Pfam, Simple Modular Architecture Research Tool (SMART, v9; http://smart.embl.de/smart/batch.pl) and Conserved Domain Database (CDD, v3.19; https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) confirmed the locations of the conserved HSP20 domain. After excluding the redundant sequence lacking the common HSP20 domain or with a molecular weight outside the range of 15–42 kDa (Ji et al., 2019), the final candidate HSP20 proteins of Cucurbita moschata were identified. Using the same method, the putative HSP20 members from Arabiodpsis thaliana and Oryza sativa were obtained from the TAIR database (https://www.arabidopsis.org/index.jsp) and the Rice Genome Annotation Project (http://rice.plantbiology.msu.edu/), respectively.

The coding sequences, genomic sequences, and amino acid sequences of the CmoHSP20 genes were obtained from the CuGenDB. The theoretical isoelectric point (pI) and molecular weight (MW) of each HSP20 protein was estimated using the pI/MW tool from ExPASy (v.3.0; http://web.expasy.org/protparam/), and the number of transmembrane regions per protein was calculated using TMHMM software (v.2.0; http://www.cbs.dtu.dk/services/TMHMM/). The structures of the HSP20 genes were visualized using the Gene Structure Display Server (GSDS 2.0; http://gsds.cbi.pku.edu.cn). The Multiple Em for Motif Elicitation (MEME) program (MEME Suite 5.3.3; http://meme-suite.org/tools/meme) was used to identify conserved motifs of HSP20 proteins, specifically recording the number of repetitions (any, maximum number of motifs-10, and the optimum motif widths set from 6 to 200 amino acid residues). The results were visualized using TBtools (Chen et al., 2020). The chromosomal position of each CmoHSP20 gene was acquired from the C. moschata genome browser at the CuGenDB and mapped using MapChart software (Voorrips, 2002).

The full amino acid sequences of 97 total HSP20 members from C. moschata (N = 33), Arabidopsis (N = 31) and Oryza sativa (N = 33) were aligned using CLUSTALW (Thompson et al., 1994) program (Supplementary Table S1). MEGA-X (Kumar et al., 2018) was used to construct a phylogenetic tree using the Maximum Likelihood (ML) method with bootstrap test of 1,000 times. The HSP20 proteins were classified into different groups according to the classifications and in silico subcellular localization of HSP20 proteins in O. sativa and Arabidopsis (Waters et al., 1996). The phylogenetic tree was visualized and enhanced using the EvolView online tool (Evolview v3; https://evolgenius.info/evolview-v2).

To investigate the degree of synteny, the homologous regions of the CmoHSP20 genes were first identified using Multiple Collinearity Scan toolkit (MCScanX) software (Wang et al., 2012). CmoHSP20 gene duplication events were determined based on whether the length of the shorter gene covered was equal or greater than 70% of the longer gene and if the similarity of the two aligned genes was equal or greater than 70% (Gu et al., 2002; Yang et al., 2008). Tandem and segmental duplications are reported to be the two main mechanisms underlying gene family expansion (Cannon et al., 2004). Genes located on the same chromosome fragment of less than 100 kb and separated by five or fewer genes were considered to be tandem duplicated genes (Wang et al., 2010). Genes found to be coparalogs located on duplicated chromosomal blocks were considered to be segmental duplicated genes (Wei et al., 2007). Tandem and segmental duplication events were visualized using Circos software (v0.69; Krzywinski et al., 2009). Ka/Ks values can be used to predict selection pressure for replicating genes. A Bio-linux system was used to calculate the nonsynonymous (Ka) and synonymous (Ks) nucleotide substitution parameters. If the ratio of Ka/Ks was greater than, equal to, or less than one, this indicated positive, neutral, and purifying selection, respectively (Hurst et al., 2002). Syntenic maps were generated using Circos software to display synteny relationships in the orthologous HSP20 genes obtained from pumpkin and the other focal species.

The tertiary structures of CmoHSP20 proteins were predicted using the online prediction tool SWISS-MODEL (https://swissmodel.expasy.org/) according to the default parameters. The cis-acting elements in the 1,500 bp upstream sequences of coding region of C. moschata HSP20 genes were retrieved from the C. moschata genome, and the types, numbers and functions of these elements were analyzed using PlantCARE software (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

Pumpkin cultivar Miben 2 (Cucurbita moschata cv. Miben 2) seeds were sprouting in an incubator under 37°C for 2 days under dark condition and then grown in a growth chamber under standard greenhouse conditions (light/dark cycle: 12 h/12 h at 25°C, 70% relative humidity). When seedlings sprouted three true leaves, uniform seedlings were transferred to a high temperature growth chamber (42°C). Root and leaf samples were collected after 0, 3, 6, 12, and 24 h under normal condition (25°C) and heat stress treatment (42°C) and frozen with liquid nitrogen for total RNA extraction and expression analysis. Each treatment was conducted with three independent replicates, and samples from five plants were collected for each replicate.

Total RNA was extracted from each sample using an RNAprep Pure Plant kit (Tiangen Biotech (Beijing) Co., Ltd., China, DP441) according to the manufacturer’s instructions. First strand cDNA was synthesized from 1 μg of RNA using the PrimeScript RT reagent kit (TaKaRa, Japan, RR047A). Quantitative RT-PCR (qRT-PCR) was performed to study the expression profiles of CmoHSP20s using gene-specific primers from TransStart Tip Green qPCR SuperMix (TransGen Biotech, AQ141-04) with an Applied Biosystems QuantStudio 3 & 5 system (ThermoFisher Scientific, United States). The housekeeping gene Actin was used as an internal control. All primers were designed to avoid the conserved region for specificity (Supplementary Table S2).

The relative expression levels of CmoHSP20 genes in roots and leaves were calculated by the 2^−△△CT method (Livak and Schmittgen, 2001). All data were calculated using the expression level under heat stress divided by that under normal condition at the same time points and presented as the means ± standard error (SE) of three replicates and differences were detected using the Student’s t-test. Asterisk (∗ or ∗∗) indicate a significant difference at p < 0.05 or 0.01, respectively.

Fourty-three HSP20 proteins were initially obtained by HMM search from pumpkin genome database. After removing the repetitive and incomplete sequences and the sequences with a molecular weight outside of the 15–42 kDa range, 33 sequences were confirmed as HSP20 genes and named based on their chromosomal locations. Sequences including genomic sequence, transcript sequence, CDS sequence and protein sequence of the gene family were shown in Supplementary Tables S3–S6. For information including gene name, gene identity, chromosomal location, open reading frame (ORF) length, amino acid (AA) number, molecular weight (MW) and isoelectric point (pI) of each CmoHSP20 gene, see Table 1. CmoHSP20 genes were distributed on 13 pumpkin chromosomes. The AA number of the CmoHSP20 proteins ranged from 121 (CmoHSP20-31) to 353 (CmoHSP20-21), with an average of 189 AAs. The MW of CmoHSP20s was between 14.29 kDa (CmoHSP20-31) and 39.57 kDa (CmoHSP20-21), with an average of 21.27 kDa, while the pI values of CmoHSP20 ranged from 4.42 (CmoHSP20-1) to 10.08 (CmoHSP20-20), with an average of 7.69. All but six proteins lacked transmembrane domains, suggesting that most CmoHSP20s were non-membrane protein.

To evaluate the evolutionary relationships of pumpkin HSP20 proteins, a ML phylogenetic tree was constructed. The pumpkin HSP20 proteins were divided into nine subfamilies: cytosol I (CI), CIII, CIV, CV, mitochondria I (MI), MII, endoplasmic reticulum (ER), plastid (P), and peroxisome (PX/Po). Cytoplasmic subfamilies (CI to CV) contained 20 members and constituted the largest clade, while no CmoHSP20s were found in subfamily CII. The M (MI and MII) contained two gene members while the ER, P, and PX/Po subfamilies each contained three HSP20s. Two CmoHSP20s (CmoHSP20-16 and CmoHSP20-25) could not be clustered into any subfamily and remained unclassified (Figure 1; Table 1).

To investigate the structural features of the HSP20 proteins, the conserved motifs were analyzed using MEME. A total of 10 distinct motifs, named motif 1 to motif 10, were detected (Supplementary Figure S1). The lengths of these conserved motifs varied from 8 (motif 8) to 49 (motif 7) AAs. The number of the conserved motifs for each HSP20 protein ranged from 1 to 9. While most CmoHSP20s had 3 to 7 conserved motifs, CmoHSP20-31 and CmoHSP20-33 contained only one motif each (motifs 1 and 2, respectively; Figure 2B). The ACD was formed from two conserved regions, CRI (β2, β3, β4, and β5) and CRII (β7, β8, and β9), which were separated by a hydrophilic domain β6-loop. Based on the Pfam and SMART analyses, the highly conserved ACD is formed from the full sequences of motif 2, 3, 1, and 5 (Figure 3A), in which the combined sequence of motif 2 and 3 contained the CRI of ACD, while the combined sequence of motif 1 and 5 contained the CRII of ACD (Figure 3B). These four motifs were detected in majority of the CmoHSP20 proteins (Figure 2B), with only four proteins (CmoHSP20-1, CmoHSP20-17, CmoHSP20-21, CmoHSP20-30) either lacking the β6-loop or containing variable sequences between β5 and β7.

To understand the evolutionary relationships of pumpkin HSP20 genes, the exon-intron structures were analyzed (Figure 2C). Among the HSP20s, 17 (52%) were intronless, 14 (42%) possessed one intron, and one gene (CmoHSP20-31) contained two introns and one gene (CmoHSP20-33) contained five introns. Gene structure analysis grouped genes with similar exon-intron patterns into the same clusters (Figures 2A,C).

In order to obtain a reasonable theoretical structure of the CmoHSP20s, protein homology modelling was performed using a SWISS-MODEL server. The predicted 3D structures are shown in Figure 4. Each CmoHSP20 protein was searched for template automatically in the software and then built model using the template (Supplementary Figure S2). The proteins can be divided into eight groups (group A to H) according to their structure similarity. Group A contained 22 proteins, with the most members. Group B, D, E and G all had only one protein. All CmoHSP20 proteins had β-turn, and the protein structures of the same template were basically similar, suggesting that the predicted results were credible.

A total of 33 CmoHSP20 genes were mapped on 13 chromosomes (Chr) and exhibited a non-uniform distribution (Figure 5). The maximum number of CmoHSP20 genes per Chr was six (Chr01 and Chr04) with only a single gene on Chrs 03, 09, 11 and 16.

According to the defined criteria, the analysis of gene duplication events showed that there were two pairs of tandem duplication genes (CmoHSP20-10/-11 and CmoHSP20-22/-23) and two pairs of segmentally duplicated genes (CmoHSP20-10/-28 and CmoHSP20-16/-25) in the pumpkin HSP20 gene family (Figure 5; Supplementary Table S7).

To understand the phylogenetic relationships of the HSP20 genes with other species, comparative synteny maps of four related genomes (C. moschata VS A. thaliana, C. moschata VS O. sativa, C. moschata VS Cucumis sativus, and C. moschata VS Cucumis melo) were created. Twelve CmoHSP20 genes displayed a syntenic relationship with genes in Arabidopsis, 19 with those in C. sativus, 18 with those in C. melo and only 1 syntenic relationship with O. sativa genes (Figure 6; Supplementary Table S8). The number of collinear gene pairs between pumpkin and other members of Cucurbitaceae (cucumber or melon) was greater than that between more distantly-related Arabidopsis or rice, with the number of collinear gene pair being lowest between pumpkin and rice.

The ratio of Ka/Ks for each orthologous CmoHSP20 gene pair was used to evaluative the type and strength of selective pressure. Both the segmental duplicated CmoHSP20 gene pairs possessed a Ka/Ks ratio <1 (purifying selection), with the higher in the CmoHSP20-16/CmoHSP20-25 pair (Ka/Ks value = 0.26).

To understand the role of cis-regulatory elements in CmoHSP20, cis-elements were identified in the 1.5 kb upstream sequence from the translation start site (ATG) of each CmoHSP20 (Figure 7; Supplementary Table S9). Four types of cis-elements, including hormone responsive, stress-responsive, plant development-related, and light responsive elements were identified. The largest number of cis-elements observed across the 33 CmoHSP20 genes was associated with light-related responsiveness, such as G-box, Box 4, AE-box and GT1-motif. The hormone-related cis-acting elements, including abscisic acid responsiveness (ABRE), auxin responsiveness (AuxRR-core and TGA-element), gibberellin responsive elements (GARE-motif, P-box, and TATC-box), MeJA-responsive (CGTCA-motif and TGACG-motif), and salicylic acid-responsive (TCA-element and SARE), were widely present in the promoter region. Among these elements, ABRE, CGTCA-motif and TGACG-motif accounted for the largest part of the hormone responsive category, while the SARE element was only found in the promoter region of CmoHSP20-21. The stress-related category cis-elements containing abiotic stress-related elements (LTR, TC-rich repeats and MBS) and biotic stress-related elements (WUN-motif), as well as plant development-related elements, including circadian control (circadian), differentiation of the palisade mesophyll cells (HD-Zip 1), endosperm expression (GCN_4 motif), meristem expression (CAT-box), flavonoid biosynthetic genes regulation (MBSI), seed-specific regulation (RY-element), and zein metabolism regulation (O2-site) were also identified. Only CmoHSP20-25 promoter region contained HD-Zip 1 element, and only CmoHSP20-26 promoter region contained RY-element element.

To investigate the expression changes of pumpkin HSP20 genes in response to heat stress, qRT-PCR was used to examine the transcript levels of the 33 CmoHSP20 genes. Overall, the relative expression level of CmoHSP20 genes fluctuated during the 24 h treatments under heat stress conditions (Figure 8). Compared to roots, the expression levels of almost all the genes in leaves were extremely up-regulated. For example, CmoHSP20-3, -5, -7, -9, -13, -14, -15, -17, -18, -22, -23, -26, -29, -30 and -32 were highly induced in leaves after short-term heat stress (42°C for 3 h), with the expression levels were more than 300-fold than under normal condition. The CmoHSP20-16, -24 and -25 genes were down-regulated under heat stress in both roots and leaves, while the expression level of CmoHSP20-31 did not change within neither roots nor leaves in response to heat stress.