In this study, the watermelon cultivar black jade (cultivated by the Institute of Agro-Biological Resources of Fujian Academy of Agricultural Sciences, Chinese Variety Registration Number: GPD Watermelon (2021) 350012) was chosen for research. The varietal characteristics were as follows: hybrid variety, medium-sized watermelon, fruit development period of 26 days, and full growth period of 80 days. The first female flower position is 8, with an interval of 4 between the appearance of female flowers. The fruit is oval-shaped, with waxy skin, a single fruit weight of 4.5 kilograms, a red flesh color, a fruit skin thickness of 1 centimeter, a green base color of the skin, and hidden skin stripes. The soluble solids content in the center is 11.2%, while that in the edge is 9.1%. The fruit skin is hard, has fine, crisp and juicy flesh, has good toughness, and is suitable for transportation. The plant is susceptible to wilt disease and shows strong drought resistance. The yield per mu in the first growth cycle was 2150.7 kilograms, an increase of 11.03% compared to that of the control, black beauty. The yield per mu in the second growth cycle was 2200.9 kilograms, an 11.20% increase compared to that of the control, black beauty. Seeds with good maturity and plump and uniform grains were selected, soaked in warm water at 55°C for 4 h, and subsequently placed in the dark at 30°C for accelerated germination. Seedlings were grown in a greenhouse and watered with 1/2 Hoagland nutrient solution every day. The day and night temperatures of the greenhouse were 28/20°C, the air humidity was 65%, the light conditions were 16 hours of light and 8 hours of darkness, and the light source was fluorescent lamps. After the second true leaf of the plants unfolded, plants with consistent growth were selected, the substrate adhering to the roots was carefully removed, and the plants were subsequently placed in a hydroponic box for hydroponic cultivation. The hydroponic method was performed as follows: a 30 L plastic box was used as a hydroponic box, the culture nutrient solution was changed once a week (Hoagland nutrient solution), foam boards were cut to match the length and width of the incubator, the incubators were covered with foam boards, and spaces were evenly made on the foam boards. When watermelon plants had 5-6 true leaves, the nutrient solution in the hydroponic tanks was replaced with 20% PEG6000. The net photosynthetic rate (Pn), transpiration rate (Tr), intercellular carbon dioxide (Ci) and stomatal conductance (Gs) were measured with an LI-6400 portable photosynthetic instrument after 20% PEG6000 treatment for 0 h, 1 h, 6 h, 12 h and 24 h (Lawrence et al., 2019). Moreover, 11 replicates of each sample (3 for RNA-seq, 5 for metabolome sequencing, and 3 for qRT−PCR) were collected, wrapped in aluminum foil, quickly placed in liquid nitrogen and subsequently stored at -80°C.

A polysaccharide and polyphenol plant total RNA extraction kit was used to extract RNA (each sample was approximately 100 mg). The extraction process was carried out according to the manufacturer’s instructions. After RNA extraction, a Nanodrop was used to determine whether the purity (OD260/280), concentration, and nucleic acid absorbance peak of the RNA were normal. An Agilent 2100 was used to accurately determine the integrity of the RNA. The indicators included the RIN, 28S/18S ratio, whether the spectral baseline was increased and the presence of the 5S peak. The concentration of RNA in all the samples was ≥600 ng/μl, and the total amount of RNA was ≥2 μg. The samples that met the criteria were shipped to Beijing Biomics Biotech Co., Ltd. (Beijing, China) for sequencing. Oligo (dT) magnetic beads were used to enrich eukaryotic mRNA through the binding of AT complementary pairs to the poly(A) tail of the mRNA. Then, fragmentation buffer was added to fragment the mRNA into short fragments. The mRNA was used as a template to synthesize first-strand cDNA with random hexamers, followed by the addition of buffer, dNTPs, and DNA polymerase I to synthesize second-strand cDNA, which was subsequently synthesized by AMPure. Double-stranded cDNA was purified using XP beads. The purified double-stranded cDNA was end-repaired and A-tailed, and sequencing adapters were ligated. The fragments were subsequently selected using AMPure XP beads, after which PCR enrichment was performed to obtain the final cDNA library. An Agilent 2100 was used to determine the size of the insert in the library. The Q-PCR method was used to accurately quantify the effective concentration of the library (effective concentration of the library>2 nM) to complete the library check. The constructed libraries were sequenced on the Illumina HiSeq 2500 sequencing platform. After obtaining the raw sequences, fastp software (version 0.23.4) was used to remove adapter sequences and filter out low-quality N sequences with a proportion greater than 5%, etc. to obtain clean reads for use in subsequent analyses (Chen et al., 2018). The clean reads were aligned to the watermelon reference genome (http://cucurbitgenomics.org/v2/organism/16) using HISAT2 software (Kim et al., 2019). The fragments per kilobase of exon model per million mapped fragments (FPKM) were used to determine the expression level. The P value and fold change were calculated with DESeq2 software, with a P value ≤ 0.05 and |log2fold change|>1 serving as the screening criteria for identifying DEGs according to differential expression (Liu et al., 2021). The DEGs were analyzed based on the KEGG database (http://www.genome.jp/kegg/) for the annotation of gene function.

After the tissue samples were stored at -80°C, placed at -20°C for 30 min and subsequently placed in a refrigerator at 4°C for thawing, 25 mg of sample was weighed and added to precooled methanol/acetonitrile/water solution (2:2:1). After vortexing, the mixture was sonicated at low temperature for 30 min, held at -20°C for 10 min, and centrifuged at 14,000 × g for 20 min at 4°C. The supernatant was collected and dried under vacuum. During mass spectrometry analysis, the solution was reconstituted by adding 100 μL of acetonitrile aqueous solution (acetonitrile:water = 1:1) followed by vortexing. The mixture was centrifuged at 14,000 × g for 15 min at 4°C, after which the supernatant was collected for sample analysis. The samples were separated on an Agilent 1290 Infinity LC ultrahigh-performance liquid chromatography (UHPLC) C-18 column (column temperature, 40°C; flow rate, 0.4 mL/min; injection volume, 2 μL; mobile phase composition, A: water + 25 mM acetic acid; ammonium + 0.5% formic acid; and B: methanol). The gradient elution program was as follows: 0.5 min, 5% B; 0.5–10 min, B was changed linearly from 5% to 100%; 10.0–12.0 min, B was maintained at 100%; B was changed linearly from 100% to 5% at 12.0–12.1 min; B was maintained at 5% at 12.1–16 min; and the samples were placed in an autosampler at 4°C during the entire analysis. The primary and secondary spectra of the samples were collected using an AB Triple TOF 6600 mass spectrometer. The raw data in the Wiff format were converted to the.mzXML format by ProteoWizard, and MS-DAIL software was subsequently used for peak alignment, retention time correction and peak area extraction. Using the data extracted by MS-DAIL, the first metabolite structure was identified, the data were preprocessed, the experimental data quality was evaluated, and data analysis steps were performed (Tsugawa et al., 2015).

In this study, the metabolomics database of the New Life Plants of the Chinese Academy of Sciences was used [in-house database (Shanghai Applied Protein Technology)]. The structures of the metabolites in the biological samples were identified by matching the molecular mass (molecular mass error within <10 ppm), second-order fragmentation spectrum, and retention time of the metabolites in the database, and the identification results were analyzed, checked and confirmed. Based on the metabolite content data matrix, principal component analysis (PCA) was performed on each sample using the R 4.2.3 language PCAtools software package (R package version 2.14). The Human Metabolome Database (HMDB, Version 5.0) and the KEGG database (https://www.genome.jp/kegg/) were used to classify the metabolites and functionally annotate the pathways to determine the main biochemical metabolic pathways and signal transduction pathways associated with the metabolites (Kanehisa and Goto, 2000). Partial least squares regression was used to construct a model of the relationship between the expression levels of metabolites and sample categories to predict sample categories (Pi et al., 2023). The fold change in the expression of each metabolite in each comparison group was calculated, and a significance test was performed on the expression level of each metabolite in each comparison group using Student’s t test to obtain the q value (corrected p value for the t test). A |log2fold change|>1 and q value < 0.05 were used as the standards for identifying differentially abundant metabolites.

The WGCNA package of the R package was used for coexpression analysis of DEG expression profiles via the dynamic branch and cut method (Langfelder and Horvath, 2008). The weighting factor β should satisfy a correlation coefficient close to 0.8. In this study, β=4 was chosen as the weighting factor. The network was constructed using the automatic network construction function “Blockwise Modules” to obtain gene coexpression modules. Each module contained different numbers of genes. The modules with minModuleSize = 30 and merge cut height = 0.25 were used as standards to merge modules with a similarity of 0.75. The correlation coefficients between the characteristic vector ME (module eigen) of the module and Pn, Tr, Ci, Gs, ABA, SA and JA were calculated. Cytoscape (version 3.10.1) software was used for visualization of the coexpression network (Shannon et al., 2003).

Total RNA was extracted using an RNAprep Pure Polysaccharide and Polyphenol Plant Total RNA Extraction Kit (Tiangen, China). The concentration of each RNA sample was determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). cDNA was obtained by reverse transcription of RNA using an M-MLV RTase cDNA Synthesis Kit (TaKaRa, Japan). qRT−PCR analysis was performed using a Roche LC480 instrument (Roche Diagnostics GmbH, Mannheim, Germany) and SYBR Green (Takara Bio, Inc.). The volume of the reaction mixture was 20 μL. The reaction program was set as follows: predenaturation at 95°C for 30 s; denaturation at 95°C for 5 s; annealing at 60°C for 5 s; and extension at 72°C for 30 s (35 cycles). The results were quantitatively analyzed by the 2^–ΔΔCt method. The Ct value of the internal reference gene was subtracted from the Ct value of the target gene to obtain the ΔCt. The mean ΔCt value of the control group (0 h) was subtracted from each ΔCt of the treatment group to obtain the ΔΔCt. The final expression level was determined based on the equation 2^–ΔΔCt (Livak and Schmittgen, 2001). The internal reference gene was ACTIN (Cla97C02G026960), and three biological replicates were performed for each procedure. All primers used in this study are listed in Supplementary Table S1 .

In this study, Excel 2010 was used for Pn, Tr, Ci, Gs and qRT−PCR data statistics, and the R ggplot2 software version was used for light and data visualization (Pn, Tr, Ci and Gs). The data before 0 h drought stress (CK) and after stress were tested using the ggpubr package to calculate the significance level. Transcriptome and metabolome heatmaps were drawn using the R pheatmap software package, and correlation analysis was performed using the corrplot software package. The correlation and p values for genes and metabolites were calculated using the R Hmisc software package, and the network diagram was visualized using Cytoscape (version 3.10.1) software (Shannon et al., 2003).

After drought stress, plants undergo a series of morphological and physiological changes to adapt to external conditions. The openings of the leaf stomata decrease in size or even close, thus hindering the entry of CO₂ into the leaves, affecting the utilization of CO₂ in the carboxylation center, and reducing the rate of photosynthesis. Therefore, measuring photosynthetic indicators and analyzing the response of plants to drought stress can help elucidate the process and mechanism underlying the impact of drought stress on plant photosynthesis; this information can be used to scientifically formulate strategies for breeding drought -resistant plants. Therefore, the Pn, Tr, Ci and Gs of watermelon leaves were measured after exposure to drought stress for 0 h, 1 h, 6 h, 12 h, and 24 h. Compared with the Pn, Tr, Ci, and Gs in watermelon leaves before treatment, the Pn, Tr, Ci, and Gs in leaves exposed to drought stress decreased significantly with increasing drought stress time, reaching a minimum at 24 h ( Figure 1 ). Pn and Tr decreased the most, with Pn at 29.64% of the pretreatment level and Tr at 29.07%. Ci and Gs decreased by 48.39% and 54.04%, respectively, of the pretreatment levels. This study showed that plants adapt to drought by changing their stomatal density, which affects photosynthesis through changes in intercellular CO₂ concentrations. Watermelon mainly reduces Pn levels by regulating Gs and Tr and closing stomata to reduce CO₂ intake, thereby mitigating the damage to its photosynthetic system caused by drought stress. To further identify key genes and metabolites involved in watermelon drought stress, RNA-seq and metabolomics sequencing were performed on leaves at five time points (0 h, 1 h, 6 h, 12 h and 24 h) before and after drought stress.

A total of 102.53 Gb of raw RNA-seq data were obtained from 15 watermelon samples collected at five time points (0 h, 1 h, 6 h, 12 h and 24 h) before and after drought stress. After filtering, a total of 99.99 Gb of clean data were obtained. The amount of clean data from each sample exceeded 5.81 Gb, the Q30 base percentage exceeded 90.67%, the GC content exceeded 41.71%, and the alignment rate with the reference genome exceeded 91.96% ( Supplementary Table S2 ). The correlation coefficients between the same sample and different biological replicates were all greater than 0.95. The PCA results showed that the biological replicate s were clustered together, indicating that the transcriptome data were reliable and repeatable ( Figure 2A ; Supplementary Figure S1 ). To analyze the transcriptome dynamics of watermelon drought stress, cluster analysis and principal component analysis (PCA) were performed on RNA-seq data from 15 samples collected at 5 time points (0 h, 1 h, 6 h, 12 h and 24 h). The RNA-seq data was divided into 3 groups. The samples at the earliest time points (0 h and 1 h) formed the first group, which represented the early stage of stress. The samples at 6 h and 12 h formed the second group, representing the middle stage of stress, during which the most DEGs were produced. The 24 h sample belongs to the third group, corresponding to the later stage of stress. The clustering of the earliest stress time points (0 h and 1 h) and 24 h was closer, indicating that their transcriptional regulation patterns were highly similar.

To study the regulatory mechanisms of transcription of watermelon plants at different stages of drought stress, DEGs at each time point during drought stress were identified ( Figure 2B ). There were 1610 DEGs at 1 h compared to that at 0 h. There were 4552 DEGs at 6 h, 5047 DEGs at 12 h, 4614 DEGs at 24 h, and a total of 701 shared DEGs at the four time points ( Figure 2B ). To elucidate the functions of all 7829 DEGs, DEG enrichment was determined by GO term and KEGG pathway analyses ( Supplementary Table S3 ). The results showed that the most significantly enriched terms in the cellular component subclass were related to the membrane and cell. The most prominent terms in the molecular function subcategory were involved in NAD and carbohydrate binding, while those in the biological process category were involved in cellular glucose homeostasis and carbohydrate metabolism ( Figure 2C ). KEGG enrichment analysis of all DEGs revealed that the DEGs were involved mainly in starch and sucrose metabolism, glutathione metabolism, plant hormone signal transduction, and photosynthesis ( Figure 2D ).

To further investigate the functional transformation among genes associated with watermelon drought stress, the 7829 genes that were differentially expressed in all periods were clustered into five clusters using the k-means clustering algorithm, and then KEGG analysis was performed to annotate the functional category of each cluster ( Figure 3 ). Cluster 1 was specifically and highly expressed at 0 h and was enriched mainly in the starch and sucrose metabolism, galactose metabolism, glutathione metabolism and flavonoid biosynthesis pathways. Cluster 2 was highly expressed at 0 h and 1 h and was mainly enriched in pyruvate metabolism, photosynthesis, galactose metabolism and carbon fixation in photosynthetic organism pathways. Cluster 3 was highly expressed at 6 h and was mainly enriched in terms related to phenylalanine metabolism, carbon fixation in photosynthetic organisms, glycolysis and gluconeogenesis, and the mitogen-activated protein kinase (MAPK) signaling pathway. Cluster 4 was highly expressed at 12 h and was mainly enriched in the starch and sucrose metabolism, plant hormone signal transduction, alpha-linolenic acid metabolism and photosynthesis pathways. Cluster 4 was highly expressed at 24 h, and the main annotations involved phenylalanine metabolism, starch and sucrose metabolism, the pentose phosphate pathway and photosynthesis pathways.

Based on the 949 metabolites obtained by UPLC−MS, PCA revealed the segregation of 25 watermelon samples at five time points (0 h, 1 h, 6 h, 12 h and 24 h) before and after drought stress ( Figure 4A ). To understand the classification and functional properties of different metabolites, the identified metabolites were annotated. The annotated metabolites were divided into 10 main categories. There were 303 metabolites related to lipids and lipid-like molecules; 103 metabolites related to organoheterocyclic compounds; 97 metabolites related to phenylpropanoids and polyketides; and 14 metabolites related to lignans, neolignans and related compounds ( Supplementary Figure S2 ). To study the changes in metabolites in watermelon at different stages of drought stress, a total of 479 differentially regulated metabolites (DRMs) were identified; compared to those at 0 h, there were 154 DRMs at 1 h, 176 DRMs at 6 h, 239 DRMs at 12 h, 428 DRMs at 24 h, and 22 DRMs at all four time points ( Figure 4B ; Supplementary Table S4 ). Lipids, lipid-like molecules and organoheterocyclic compounds exhibited the greatest differences (164 and 64, respectively) ( Figure 4C ). The changes in DRM content showed that these metabolites were affected by drought stress, and they were grouped into four clusters using the k-means clustering algorithm ( Figure 4D ). Cluster 1 had the highest levels at 6 h and 24 h and included carbohydrates, fatty acyl glycolosides and fatty amides. At 6 h, Cluster 2 exhibited the greatest increase in carbohydrate, fatty acid and terpene glycoside contents. Cluster 3 had the highest content at 0 h, 6 h and 12 h and included glycosylglycerols, quinoline carboxylic acids and fatty acids. Cluster 4 had the highest content of diterpenoids, fatty acids and carbohydrates at 24 h.

Through joint analysis of the transcriptome and metabolome, the relationships between 7829 DEGs and 479 DRMs under drought stress in watermelon were systematically explored. KEGG pathway enrichment analysis revealed that DEGs and DRMs were mainly enriched in 11 KEGG pathways, such as starch and sucrose metabolism, plant hormone signal transduction, and photosynthesis pathways ( Figure 5A ). First, changes in the levels of metabolites in the starch and sucrose metabolism pathways were analyzed. The levels of most metabolites increased with the prolongation of drought stress, while the D-sucrose content reached a minimum at 1 h of drought stress and then gradually increased, reaching a maximum at 12 h ( Figure 5B ). The expression levels of the glycogen biosynthesis pathway genes 1,4-alpha-glucan branching enzyme (GEB) and UTP-glucose-1-phosphate uridylyltransferase (GalU) both showed an upward trend, reaching a maximum level at 24 h, while the expression levels of starch synthase (ss), granule-bound starch synthase (GBSS) and glucose-1-phosphate adenylyltransferase (GlgC) decreased ( Figure 5C ). The expression of all genes in the glycogen degradation pathway (including glycogen phosphorylase (GP), 4-alpha-glucanotransferase (4GT) and phosphoglucomutase (PGM)) tended to increase. This result indicates that the increase in sugar content under drought conditions may be due to a decrease in glycogen degradation. A network of DEGs and DRMs in the starch and sucrose metabolism pathway was constructed using the screening criteria of PCC ≥ 0.75 and P < 0.05. The 15 genes that were most strongly regulated by uridine 5-diphospho-D-glucose and the other metabolites were regulated by 2-6 genes ( Figure 5D ). GBE-1 was a key gene associated with eight starch and sucrose compounds.

The levels of two metabolites (plastoquinol-1 and adenosine 5-diphosphate) in the watermelon photosynthesis pathway began to decrease 12 h after drought stress, and the level of adenosine 5-triphosphate began to decrease after 1 h ( Figure 6A ). The expression levels of the photosynthesis pathway genes photosystem I subunit II (PSI-II), photosystem I subunit III (PSI-III), photosystem I subunit IV (PSI-IV), cytochrome b6 (Cyt b6) and F-type H+-transporting ATPase subunit a (F-type H+-ATPase) decreased after 6 h of drought stress ( Figure 6B ). A network of DEGs and DRMs in the photosynthesis pathway was constructed using the screening criteria of PCC ≥0.75 and P < 0.05. PSI-II regulated only adenosine 5-diphosphate and adenosine 5-triphosphate, and the remaining seven genes were shown to interact with three metabolites. These parameters were significantly correlated ( Figure 6C ).

UPLC−MS revealed significant changes in the levels of four hormones (ABA, IAA, JA and SA) in watermelon plants under drought stress ( Figure 7A ). The levels of all hormone-related metabolites increased except for those of 2-azaniumyl-3-(1H-indol-3-yl) propanoate. First, the expression levels of genes related to the ABA signaling pathway were analyzed. The expression levels of four pyrabactin-resistant/pyrabactin-resistant 1-L like (PYR/PYL) genes tended to decrease, while those of two genes tended to increase ( Figure 7B ). The expression levels of protein phosphatase 2C (PP2C) and ABA-responsive element (ABRE)-binding factor (ABF) both showed increasing trends. The expression of the three SNF1-related protein kinase 2 (SnRK2) genes decreased, while that of SnRK2-4 tended to increase, reaching a maximum at 6 h. These findings indicate that different genes in the ABA pathway respond to drought stress at different times. These genes may promote watermelon tolerance to drought stress. In contrast, the expression of most genes in the IAA pathway tended to decrease, which was contrary to the changes in IAA content ( Figure 7C ). Most genes in the JA pathway (except jasmonate ZIM-domain (JAZ)-5) had upregulated expression, and genes in the SA pathway, including nonexpressor of pathogenesis-related genes 1 (NPR1) and pathogenesis-related protein 1 (PR1) showed a downward trend, except for TGACG-binding (TGA)-3 and -4 which had upregulated expression ( Figure 7D ). Using the screening criteria of a PCC ≥ 0.75 and P < 0.05, a network of DEGs and DRMs in the plant hormone signal transduction pathway was constructed. Specifically, 21 genes regulated ABA, 18 genes regulated IAA, and eight genes regulated JA. SA was regulated by seven genes ( Figure 7E ).

Based on the expression levels of the 7829 DEGs, β=4 was chosen for network construction, and the dynamic pruning tree method was used to merge modules with similar expression levels. A total of 7 coexpression modules were obtained, and different colors were used to represent different modules ( Figure 8A ). The calculated modules correlated with Pn, Tr, Ci, Gs, ABA, SA and JA ( Figure 8B ). The turquoise module was significantly correlated with Pn, Tr, Ci and Gs (r > 0.8, p < 0.05); the brown module was significantly correlated with ABA (r = 0.9, p < 0.05); and the green module was significantly correlated with ABA, SA and JA (r > 0.86, p < 0.05). For each module, the five genes with the highest connectivity were determined to be hub genes, and 15 hub genes were ultimately identified ( Figure 8C ). The 15 genes included one encoding a bHLH transcription factor (Cla97C03G068160), one encoding a MYB transcription factor (Cla97C01G002440), one encoding an HSP transcription factor (Cla97C02G033390), one encoding a GRF transcription factor (Cla97C02G042620), one encoding the ABA pathway-related gene SnRK2-4 (Cla97C10G186750) and one encoding GP-2 (Cla97C05G105810) in the starch and sucrose metabolism pathway ( Table 1 ).

To confirm the accuracy of the transcriptome expression profile, 10 genes were randomly selected for qRT−PCR analysis. Three independent replicates were performed, and the correlation between the qRT−PCR results and the fold difference in RNA-seq data was calculated. The results showed that the transcriptome data and qRT−PCR data were significantly correlated (R=0.94, p < 0.01), and these comprehensive results indicated that the transcriptome sequencing data were reliable ( Figure 9 ).