The watermelon variety ‘Xi Nong 8’ was used in this study. Seedlings were cultured in a light incubator at day/night 28/22°C (16/8 h) for 30 days, and Hoagland nutrient solution was used for hydroponics. For tissue expression analysis, the roots, stems, true leaves and cotyledons of seedlings were sampled. The seedlings placed in the Hoagland nutrient solution containing 10% polyethylene glycol 6000 (PEG6000) were used for drought treatment, the seedlings under 4°C conditions were used for cold treatment, and the seedlings spraying 100 μM ABA (dissolved in 0.004% ethanol) for hormone treatment. The control seedlings for ABA treatment were sprayed with 0.004% ethanol, and the control for abiotic stresses were the hydroponic seedlings under normal conditions. The leaves were collected after treatment for 0, 0.5 (cold and drought stress), 1, 3, 6, 12, and 24 h. Three biological replicates were performed for each treatment, and samples for each replicate were collected from five seedlings.

The identification process of ClHSP70 family in watermelon was shown in Supplementary Figure S1. Briefly, the watermelon genome sequence was downloaded from EnsemblPlants database (EB, http://plants.ensembl.org/info/data/index.html). The Hidden Markov Model (HMM) profile of HSP70 domain (PF00012) was downloaded from Protein family (Pfam) database (Finn et al., 2014) and used as a BLAST query against the watermelon genome database (E-value<1 × 10⁻⁵). Meanwhile, HMMER-3.0 program was employed to identify watermelon ClHSP70s (Finn et al., 2011). Subsequently, SMART (http://smart.embl-heidelberg.de/) and Pfam (http://pfam.xfam.org/search) were used to identify the HSP70 domains in all putative ClHSP70s to verify whether they belonged to HSP70 family.

The chromosome distributions, amino acid (AA), genome sequences and gene structures of watermelon ClHSP70 members were obtained from EB database. In addition, EXPASY Proteomics Server (https://web.expasy.org/protparam/) was used to analyze the AAs, molecular mass, and theoretical isoelectric points (Lopes-Caiter et al., 2013). The cell localization prediction of ClHSP70s was performed using Wolf PSORT program (http://wolfpsort.org/).

The protein sequences of HSP70 family from various species, including A. thaliana, Solanum lycopersicum, Solanum tuberosum, Csativus sativus, Vitis vinifera, Zea mays, Oryza sativa and C. lanatus, were multiply aligned using the ClustalW program of MEGA11 (Tamura et al., 2021). The gene IDs of HSP70 family from different plants were shown in Supplementary Table S1. The phylogenetic tree was constructed based on the maximum likelihood (ML) method with 1000 bootstrap replications using TBtools (Chen et al., 2020).

Conserved motifs in ClHSP70 sequences were searched by MEME suite (https://meme-suite.org/meme/) (Bailey et al., 2009), with the following parameters: any number of repetitions; maximum number: 10; optimum widths: 6-200 AAs. SWISS-Model (https://swissmodel.expasy.org/) was used to model three-dimensional (3D) structures of ClHSP70s based on homology modeling (Waterhouse et al., 2018). Global model quality estimation (GMQE) was used to evaluate the model quality (Wang et al., 2023). To illustrate the exon/intron structures of all ClHSP70s, the cDNA sequences were aligned with their genomic sequences via TBtools.

The chromosomal localizations of ClHSP70s and the genome files of watermelon, rice and Arabidopsis were obtained from EB database. The chromosome distribution and gene duplication of ClHSP70 genes were analyzed and mapped using TBtools software, which was also used to visualize the collinearity map between different species. The duplication of ClHSP70s was defined according to the following criteria: a) the sequence alignment length covered more than 70% of the longer gene sequence; b) more than 70% similarity between aligned gene regions (Gu et al., 2002; Yang et al., 2008). Tandem duplicated genes were defined as two genes separated by no more than five genes in a 100 kb chromosomal region (Wang et al., 2010). Paralogous genes located on duplicated chromosomal blocks were considered as segmental duplication (Wei et al., 2007).

The ratios of nonsynonymous (Ka) and synonymous (Ks) (Ka/Ks) were calculated for the results of ClHSP70 duplications, and Ks was converted into the differentiation time (T) with the formula: T = Ks/(2 × 6.1×10⁻⁹) ×10⁻⁶ million years ago (Mya) (Lynch and Conery, 2000).

The promoter sequences of ClHSP70s (upstream sequence 2,000 bp) were identified by TBtools. The cis-elements were analyzed by PlantCARE online (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/). The interaction analysis of ClHSP70s was performed using the STRING online (http://cn.string-db.org/), and Arabidopsis was used as a model plant for the query. The results were exported and plotted using Cystoscope software.

Gene Ontology (GO) enrichment analysis of watermelon ClHSP70s was conducted using AmiGO 2 (http://amigo.geneontology.org/amigo/landing) with the corrected p-value < 0.05. The GO analysis contains cellular component (CC), biological process (BP), and molecular function (MF).

The sequences of 12 ClHSP70s were downloaded from the watermelon genome, and the bases were set as 2,000 bp before the start codon and 2,000 bp after the stop codon. These sequences were subjected to simple sequence repeat (SSR) mining by MISA (https://webblast.ipk-gatersleben.de/misa/index.php?action=1) (Beier et al., 2017). The minimum standard lengths for dinucleotide, trinucleotide and tetranucleotide repeats were 6, 5 and 5 repeats, respectively (Fu et al., 2023).

Fragments Per Kilobase Per Million (FPKM) values for each ClHSP70s in leaves under cold stress (4°C for 36 h, PRJNA328189) and drought stress conditions (4 days and 8 days, PRJNA454040) were retrieved from CuGenDB database (http://cucurbitgenomics.org/r naseq/home). Watermelon cultivar Y134 was used for cold treatment, while other two varieties Y34 and M20 were used for drought stress treatment. The data for each ClHSP70 were normalized by log2 (FPKM) ( Supplementary Table S2 ) and shown as Heatmap using TBtools.

Total RNA was extracted from samples by Total RNA kit (TIANGEN, Beijing, China), and reverse-transcribed by One-Step gDNA Removal and cDNA Synthesis Super Mix kit (TransGen, Beijing, China) according to the manufacturer’s instructions. The final cDNA was used as the template for subsequent quantitative real-time PCR (qRT-PCR).

Specific primers for detection of ClHSP70 gene expression (Supplementary Table S3) were designed using Primer 5 software. The yellow-leaf-specific proein8 (ClYLS8) and β-actin (ClACT) genes from watermelon were used as the reference genes (Kong et al., 2014). qRT-PCR was conducted using ChamQ universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) with a 20 μL system on the qTOWER³ Series-Real-Time Thermal Cyclers (Analytik Jena, Thüringen, Germany). The experimental procedure of qRT-PCR carried out in accordance with the instruction manual of SYBR qPCR Master Mix, including an initial activation of 95°C for 30 s, then 40 cycles of 95°C for 10 s, 60°C for 30 s. The 2^−△△CT method was used to calculate the relative transcript levels of ClHSP70s (Schmittgen and Livek, 2008).

The microRNA (miRNA) sequences were downloaded from PmiREN2.0 (https://www.pmiren.com/) (Guo et al., 2022). The prediction of binding sites of ClHSP70s targeted by miRNA was performed using psRNATarget (https://www.zhaolab.org/psRNATarget/) (Dai et al., 2018), and the parameters were set to default values. The track predicted ClHSP70 gene binding sites were listed in Supplementary Table S4.

Using HMMER search and BLASTP alignment, 12 HSP70 family members were finally identified (Supplementary Table S5). They were named ClHSP70-1 (62.01 kDa) to ClHSP70-12 (102.80 kDa) based on their molecular weight. The physicochemical properties of ClHSP70s were analyzed by EXPASY online software. The ClHSP70 genes were unevenly distributed on 7 of 11 chromosomes. The minimum length of the encoding sequence was 1,716 bp and the maximum length was 2,766 bp, the number of AAs corresponding to the encoding proteins was 571 and 921, respectively (Table 1). The theoretical isoelectric points were between 4.9 (ClHSP70-6) and 5.78 (ClHSP70-7), indicating that ClHSP70 proteins were acidic. Only the total average hydrophilicity (GRAVY) of ClHSP70-1 was greater than 0, indicating that most ClHSP70s were hydrophilic proteins. Six members (ClHSP70-1, -2, -3, -4, -5, and -11) were mainly located in cytoplasm, one (ClHSP70-10) in nucleus, two (ClHSP70-6 and -9) in chloroplasts, one (ClHSP70-7) in mitochondria, and two (ClHSP70-8 and -12) in endoplasmic reticulum.

In addition, the intron number of all members was analyzed to gain insight into the structural evolution of ClHSP70s (Figure 1). The intron number of ClHSP70s ranged from 0 to 14 (Table 1). The intron number in Group Ⅱ and Ⅲ was higher than other groups, but ClHSP70-1 from Group Ⅲ had no introns. In general, the exon-intron distribution characteristics of ClHSP70s in the same group were similar, while the difference in intron number between different groups was large. Motif analysis showed that the types and amounts of the motifs in Group Ⅰ and Ⅱ were similar, while ClHSP70-10, -11, and -12 in Group Ⅲ did not contain Motif 2, and ClHSP70-11 and -12 did not contain Motif 4. Unlike other ClHSP70s, ClHSP70-12 did not yet contain Motif 6 (Figure 1). The 3D protein structures of ClHSP70s were predicted by SWISS-Model (Supplementary Figure S2) and their GMQE values were higher than 0.8. The number of β folds in the protein structure of each ClHSP70 member was higher than that of α helices, indicating that β folds were the main components of ClHSP70s. The protein structure of ClHSP70-1 was quite different from that of other members, mainly attributing to the fact that the number of α helices contained in ClHSP70-1 (13) was much lower than that of other members (≥18). This structural difference was consistent with the results of the cluster analysis of ClHSP70-1 (Figure 1A). The remaining protein structures were highly similar, such as ClHSP70-2, -3 and -4. They had the same number of α helices, while ClHSP70-4 had only one more β fold than the other two proteins. These results suggest that ClHSP70s are evolutionarily conserved and may have functional redundancy, while structural differences indicate functional specificity among different members.

Based on the annotation file of the watermelon genome, TBtools was used to visualize the distribution of ClHSP70s on chromosomes (Chrs). Twelve watermelon ClHSP70 genes were randomly distributed on 7 Chrs, located on Chr1 (1), Chr4 (1), Chr5 (2), Chr6 (1), Chr9 (4), Chr10 (1) and Chr11 (2), respectively (Figure 2). Tandem duplicated and segmental duplication have certain effects on the formation of gene families during evolution (Cannon et al., 2004). Three pairs of duplication genes in ClHSP70 family were identified, of which ClHSP70-3/-5 was tandem duplicated, ClHSP70-2/-3 and ClHSP70-2/-5 were segmental duplication (Figure 2). To understand whether ClHSP70 genes were affected by natural selection during evolution, the Ka/Ks values of ClHSP70-2/-3, ClHSP70-2/-5 and ClHSP70-3/-5 were calculated. Their Ka/Ks ratios were all less than 1, suggesting that ClHSP70s underwent strong purifying selection during evolution (Li et al., 2022), and duplicated events occurred at 123.4 (ClHSP70-2/-5), 163.6 (ClHSP70-3/-5) and 175.5 (ClHSP70-2/-3) Mya, respectively (Table 2).

A ML phylogenetic tree was constructed using the protein sequences of 157 HSP70s from multiple plants, including A. thaliana (18), S. lycopersicum (19), S. tuberosum (21), C. sativus (13), V. vinifera (15), Z. mays (28), O. sativa (31), and C. lanatus (12) (Figure 3; Supplementary Table S1). Based on their evolutionary relationship, HSP70 proteins were classified into five subfamilies, including cytosol (C), ER, mitochondria (M), chloroplast (Chl) and nucleus (N) subfamily. Nearly half of HSP70s belonged to the cytoplasmic subfamily, which was consistent with the subcellular localization of watermelon ClHSP70s, suggesting that plant HSP70s function mainly in the cytoplasm. HSP70 members from different species of the same subfamilies were more closely related than members from different subfamilies. For example, watermelon ClHSP70-2 and cucumber Cucsa.066080.1 in C subfamily, rice Os06g46600.1 and maize GRMZM2G056766-T01 in N subfamily. Each subfamily had different HSP70 branches formed by monocots and dicots, which also indicated that HSP70 family was formed before the differentiation of monocots and dicots. Twelve pairs of orthologous genes were identified between watermelon and cucumber, suggesting that the ancestral genes of this family were formed before these two species. Additionally, there were several pairs of paralogous genes, indicating that ClHSP70s had undergone multiple replications after watermelon speciation.

To explore the underlying evolution of ClHSP70 family, the collinearity of watermelon with model plants Arabidopsis and monocotyledon rice was analyzed (Figure 4A). Three pairs of homologous gene pairs between watermelon and Arabidopsis (ClHSP70-2/AtHSP70-2, ClHSP70-2/AtHSP70-3, and ClHSP70-8/AtHSP70-12) were identified (Lin et al., 2001), which was much less than 8 pairs between watermelon and rice (ClHSP70-2/OscHSP70-1, ClHSP70-2/OscHSP70-5, ClHSP70-2/OscHSP70-6, ClHSP70-3/OscHSP70-1, ClHSP70-3/OscHSP70-2, ClHSP70-3/OscHSP70-5, ClHSP70-11/OsHSP110-1, and ClHSP70-11/OsHSP110-4) (Jung et al., 2013).

To investigate the expression characteristics and potential functions of ClHSP70 family, analysis of the gene promoter sequences (2,000 bp upstream of the start codon) was performed (Supplementary Table S7). The distributions of cis-acting elements of ClHSP70 promoters were visualized, and 15 functional elements were labeled (Figure 4B). The identified cis-elements were divided into hormone response group, stress response group, and growth and development group. Among them, the hormone response group contained the most elements, up to 66, mainly including ABA (24), methyl jasmonic (15), gibberellin (11), and salicylic (8). The stress response group also had many cis-elements, such as anaerobic induction (19), defense and stress (9), drought stress (8), and low-temperature stress (4). However, only 14 elements in the growth and development group. In summary, ClHSP70s may present different expression patterns during developmental and environmental stress response.

To provide more clues about the function of ClHSP70s, STRING online software was used to predict the protein interaction network of this family members based on the interolog of the Arabidopsis interactome. ClHSP70-8 and -12 both interacted with other proteins, while ClHSP70-2 and -5 did not interact with each other, and neither of them interacted with other three proteins (ClHSP70-1, -3, and -4). Additionally, ClHSP70-7 did not interact with ClHSP70-6 or -9 (Figure 5). These data provide a reference for in-depth study of the biological function of ClHSP70s.

In order to deeply understand the functions of this family, GO annotation and enrichment analysis of ClHSP70s were performed (Supplementary Figure S3). In MF, all ClHSP70s were assigned in binding, such as unfold protein, small molecule, anion and ATP binding. However, only 3 ClHSP70s (25%) were assigned in unfold protein binding. In CC, cellular component was enriched in all 12 genes. Nevertheless, only three members (25%) were assigned in protein folding of BP. These data suggest that ClHSP70s may function as molecular chaperones in cells.

A total of 66 SSRs were discovered from 12 ClHSP70 sequences (Supplementary Table S8). Only the ClHSP70-9 sequence had one SSR, and the remaining members all had multiple SSRs. In addition, 3 SSRs present in compound formation were detected in ClHSP70-6, -11 and -12 sequences. The 2000 bp promoter regions contained 23 SSRs, the gene regions contained 19 SSRs, the post-2000 bp sequences contained 24 SSRs, and ClHSP70-6 contained the most SSR sites (14) (Supplementary Table S8). Among all the SSR types, dinucleotide repeats were the most abundant (53.33%), followed by trinucleotide repeats (40.00%), while there was only one SSR (6.67%) with a tetranucleotide repeat. SSRs with five sequence repetitions (26.67%) were the most abundant (Supplementary Table S9).

To analyze the tissue expression specificity of ClHSP70s, the relative levels of all members in 4 tissues (roots, stems, true leaves and cotyledons) of 30-day-old seedlings were determined by qRT-PCR (Figure 6). The expression levels of ClHSP70-1 and -4 were the highest in stems and the lowest in true leaves. The levels of ClHSP70-2, -3, -11 and -12 in stems were similar to those in roots, but significantly increased in cotyledons and true leaves. The abundance of other six genes (ClHSP70-5, -6, -7, -8, -9, and -10) in stems was the lowest. The expression levels of ClHSP70-5, -8 and -9 in cotyledons and true leaves were significantly higher than those in other tissues, while ClHSP70-6 and -10 were preferentially expressed in cotyledons. In general, ClHSP70s may function in the biological processes of different tissues in watermelon.

HSP70s as molecular chaperones are widely involved in plant heat adaptation, and their roles in other abiotic stresses (such as drought and cold stress) have also been demonstrated in a variety of plants (Wang et al., 2004). However, the response of watermelon ClHSP70s to drought and cold stress remains unknown. To understand the response of ClHSP70 to cold and drought stress in detail, the effects of cold and drought stress on ClHSP70 expression were analyzed based on the RNA-seq data from the CuGenDB database (Figure 7). Under cold conditions, all ClHSP70s except ClHSP70-11 responded to low temperature (4°C for 36 h) to varying degrees. In particular, the levels of ClHSP70-1 and -4 were significantly downregulated (Figure 7A). Three genes (ClHSP70-2, -5, and -9) were not sensitive to drought stress, however, the levels of other nine members were observably regulated by drought treatment (4 or 8 days) (Figure 7B).

Moreover, the dynamic expression patterns of all members in 30-day-old seedlings were analyzed under simulated drought (10% PEG6000) and cold stress (4°C) conditions (treatment for 0, 0.5, 1, 3, 6, 12, and 24 h). As shown in Figure 8A, the expression levels of all ClHSP70s except ClHSP70-1 and -4 were higher after 12 h and/or 24 h of simulated drought treatment. The expression level of ClHSP70-1 was the highest after drought treatment for 6 h, and the level of ClHSP70-4 reached the peak after treatment for 0.5 h. Notably, ClHSP70-1 and -7 were significantly downregulated after 1 h, the levels of ClHSP70-2 and -8 reached the lowest after 6 h, while the lowest level of ClHSP70-6 was found after 3 h. The expression of ClHSP70-1, -3, -4, -5, -7, -8 and -10 were significantly upregulated by short-term (0.5 or 1 h) cold treatment (Figure 8B). Prolonged treatment inhibited the transcription of many genes to varying degrees. For example, both ClHSP70-6 and -8 were significantly downregulated after 3 and 6 h, while the lowest level of ClHSP70-12 appeared after 12 h. The levels of ClHSP70-2 and -3 at 12 and 24 h of treatment were higher than those at other time periods. Different from other ClHSP70s regulated by low temperature, ClHSP70-9 was not sensitive to cold stress. In summary, ClHSP70s may function in the formation of drought and cold stress response networks in watermelon.

ABA signaling pathway is involved in various stresses. The promoters of ClHSP70s contained many ABA-responsive elements, suggesting that the expression of these genes may be regulated by ABA. Therefore, 30-day-old watermelon seedlings were treated with exogenous ABA (100 μM), and the expression patterns of ClHSP70s in leaves at different times after treatment (0, 1, 3, 6, 12, and 24 h) were analyzed. As shown in Figure 9, the levels of five members (ClHSP70-7, -9, -10, -11, and -12) were downregulated by ABA, while the remaining seven genes (ClHSP70-1, -2, -3, -4, -5, -6, and -8) were upregulated at specific time points. The expression patterns of ClHSP70-1 and -6 were similar. Their levels increased gradually with the extension of post-treatment time and reached the highest at 24 h after treatment, which increased 2.7- and 24-fold, respectively. The response of ClHSP70-4 to ABA was relatively slow, and its expression level increased significantly (nearly 30-fold) only at 24 h after treatment. On the contrary, ClHSP70-2, -3 and -5 responded rapidly to ABA, and were strongly induced at 1 h after treatment, which increased by about 300-, 30- and 500-fold, respectively (Figure 9). These data suggest that watermelon ClHSP70s may be widely involved in ABA-mediated regulation of abiotic stresses.

To determine whether ClHSP70s were regulated by microRNAs (miRNAs) through targeted binding sites, the binding sites of miRNAs targeting ClHSP70s were predicted. As shown in Figure 10, 10 ClHSP70s (ClHSP70-1, -2, -3, -4, -6, -8, -9, -10, -11, and -12) might be regulated by 28 miRNAs. Among them, four genes (ClHSP70-2, -3, -10, and -11) contained only one miRNA binding site, while the binding sites of ClHSP70-12 targeted by miRNA was the most. Notably, among the five binding sites regulating ClHSP70-4, two sites each were targeted by miR396a and miR396c, and one was targeted by miRN828. Furthermore, both miR396a and miR396c could regulate ClHSP70-6 and -9, and miR397 regulated ClHSP70-10 and -11 (Figure 10; Supplementary Table S4). These results indicate that one miRNA can regulate multiple ClHSP70s and vice versa.