Abstract

Background

   Pumpkin (Cucurbita moschata) is an important vegetable crop that often
   suffers from low-temperature stress during growth. However, the
   molecular mechanism involved in its response to chilling stress remains
   unknown. In this study, we comprehensively investigated the effect of
   chilling stress in pumpkin seedlings by conducting physiological,
   transcriptomic, and metabolomic analyses.

Results

   Under chilling stress, there was an overall increase in relative
   electrical conductivity, along with malondialdehyde, soluble sugar, and
   soluble protein contents, but decreased superoxide dismutase and
   peroxidase activities and chlorophyll contents in seedling leaves
   compared with controls. Overall, 5,780 differentially expressed genes
   (DEGs) and 178 differentially expressed metabolites (DEMs) were
   identified under chilling stress. Most DEGs were involved in plant
   hormone signal transduction and the phenylpropanoid biosynthesis
   pathway, and ERF, bHLH, WRKY, MYB, and HSF transcription factors were
   induced. Metabolomic analysis revealed that the contents of salicylic
   acid (SA), phenylalanine, and tyrosine increased in response to
   chilling stress. The findings indicated that the SA signaling and
   phenylpropanoid biosynthesis pathways are key to regulating the
   responses to chilling stress in pumpkins.

Conclusion

   Overall, our study provides valuable insights into the comprehensive
   response of C. moschata to chilling stress, enriching the theoretical
   basis of this mechanism and facilitating the development of molecular
   breeding strategies for pumpkin tolerance to chilling stress.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12864-024-10939-2.

   Keywords: Cucurbita moschata, Chilling stress, Transcriptome,
   Metabolome

Background

   Low temperatures are one of the main environmental factors affecting
   crop plant growth, development, and yields [[34]1]. Cold stress can be
   categorized as chilling (0–15 °C) and freezing (< 0 °C) stress based on
   temperature. Plants from tropical and subtropical regions are often
   sensitive to cold stress, which results in plant damage [[35]2].
   Pumpkin (Cucurbita moschata) is an annual dicotyledonous vegetable,
   with creeping or climbing stems bearing tendrils. Due to its resistance
   to chilling stress, C. moschata has been used successfully as rootstock
   to produce chilling-stress sennsitive watermelon (Citrullus lanatus)
   crops [[36]3]. To overcome chilling stress, plants undergo a series of
   physiological and molecular changes in metabolic processes, ‌signal
   transduction, and stress-sensing mechanisms [[37]4, [38]5]. Therefore,
   unraveling how C. moschata defends itself against chilling
   stress-induced impairment is crucial to understanding the molecular
   mechanisms of resistance to cold stress.

   Recently, valuable insights into cold stress have been reported. For
   example, endogenous plant hormones, including auxin, ethylene (ET),
   gibberellin (GA), abscisic acid (ABA), brassinosteroid (BR), jasmonic
   acid (JA), and salicylic acid (SA) are known to play crucial roles in
   the defense against chilling stress [[39]6–[40]12]. In cucumber,
   RNA-sequencing (RNA-seq) analysis has revealed that genes in auxin, SA,
   JA, and ET signaling pathways are regulated in response to chilling
   stress [[41]13]. In particular, the biosynthesis and regulation of SA
   have attracted extensive attention. SA can be synthesized via the
   isochorismate synthase (ICS) pathway or the phenylalanine ammonia-lyase
   (PAL) pathway. Both require the primary metabolite chorismate [[42]14].
   Cheng et al. revealed that chilling stress induces the accumulation of
   free and conjugated SA in watermelon plants, and that SA production is
   attributed to the PAL pathway, as PAL gene expression and activity in
   watermelon increase significantly under chilling stress [[43]15].
   Although plant hormones are understood to be crucial to signal
   transduction and resistance to chilling stress, the response mechanism
   of pumpkin rootstock under chilling stress is unclear.

   Plant secondary metabolites originating from the phenylpropanoid
   pathway, such as coumarins, lignans, flavonoids, phenolic acids, and
   lignins, are also critical to plant responses to cold stress [[44]16,
   [45]17]. Recent studies have revealed that cold stress induces the
   expression of various key enzymes involved in the phenylpropanoid
   biosynthesis pathway, including PAL, chalcone synthase (CHS),
   peroxidase, 4-coumarate: CoA ligase (4CL), and cinnamoyl CoA reductase
   (CCR) [[46]18–[47]22]. However, whether the phenylpropanoid pathway is
   involved in the response of pumpkins to chilling stress is not unknown.

   In addition, transcription factors (TFs) coordinate the expression of
   downstream target genes. In Arabidopsis and grasses, TFs such as ERF,
   MYB, bHLH, bZIP, AREB/ABF, and DREB1/CBF regulate gene expression under
   abiotic stress [[48]17]. Low-temperature signaling can be mediated by
   the dehydration-responsive element binding factor (DREB) and ERF
   [[49]16, [50]21]. Zhu et al. showed that the expression levels of
   AtCOR413, AtCOR15B, AtCBF1, and AtCBF2 are regulated in transgenic
   Arabidopsis overexpressing AeWRKY32 in Okra (Abelmoschus esculentus L.)
   [[51]22]. Moreover, FvMYB82 from Fragaria vesca is involved in
   regulating the expression of downstream genes associated with cold
   stress [[52]23]. Therefore, TFs play an important role in plant
   responses to cold stress. To date, there have been no reports of TFs
   responding to chilling stress in C. moschata.

   Currently, there is limited information about the genes involved in
   cold tolerance in pumpkins. For example, seven cold stress-related
   glutathione transferase (GST) genes and chromosome condensation 1
   (CmRCC1) gene were characterized when pumpkin seedlings were exposed to
   4–5 °C for 24 h, respectively [[53]24, [54]25]. Furthermore, the
   expression of seven CmoDELLA genes in C. moschata was found to be
   significantly induced under cold stress [[55]26]. These findings
   suggest that these genes play vital roles in the response to cold
   stress. Additionally, multi-omics approaches have emerged as successful
   technologies for plant systems and have enhanced our understanding of
   molecular regulatory networks controlling response to cold stress.

   In this study, we conducted a comprehensive analysis of the
   physiological responses, transcriptomics, and metabolomics of the
   cold-tolerant inbred line C. moschata IL1 at 5 °C for 24 h. The
   objectives of the present study were as follow: (1) to determine the
   gene expression levels and metabolite profiles of IL1 under chilling
   stress; (2) to analyze differentially expressed genes (DEGs) and
   differentially expressed metabolites (DEMs) in response to chilling
   stress; and (3) to identify important genes and pathways that drive
   chilling stress-induced changes. We believe that the results of our
   research provide important insights into the molecular mechanism of
   chilling tolerance in pumpkins, which can be applied in molecular
   breeding programs.

Results

Physiological changes under chilling stress

   Pumpkin inbred line IL1 seedlings at the two-leaf stage were assessed
   for any physiological changes after 24 h of chilling stress treatment
   (5 °C). The leaves appeared wilt with faded edges and displayed inward
   curling, indicating that the health of seedlings was affected,
   resulting in their slow growth. Before and after 24 h of chilling
   stress, the second-true leaves (about 0.1 g) were sampled to measure
   the relative electric conductivity (REC), leaf greenness index (Soil
   and Plant Analysis Development, SPAD), soluble sugar content, soluble
   protein content, and the activities of superoxide dismutase (SOD),
   peroxidase (POD), and catalase (CAT). The REC and soluble sugar, and
   soluble protein contents in the plants treated with chilling stress
   were significantly higher than those of the control group, while the
   SPAD index, and activities of SOD and POD were significantly lower than
   those of the control group (Fig. [56]1). The MDA content, an indicator
   of lipid peroxidation, was significantly increased in treated plants.
   These results suggest that REC, SPAD index, soluble protein, and sugar
   contents, SOD and POD activities, and MDA content were affected by
   chilling stress. There was no significant difference in CAT activity
   between the control and treatment groups (Fig. [57]1). Our results
   revealed that chilling stress leads to significant physiological
   changes.

Fig. 1.

   [58]Fig. 1
   [59]Open in a new tab

   Physiological changes of IL1 under control and cold treatment
   condition. A, REC; B, SPAD index; C, Soluble sugar content; D, Soluble
   protein content; E, SOD; F, POD, G, CAT; H, MDA. Values are presented
   as mean ± standard deviation (n = 3). Bars with the same small letters
   do not statistically differ by the Duncan’s multiple range test at
   P < 0.05

Transcriptome sequencing data analysis

   A total of six cDNA libraries with three biological replicates were
   constructed for transcriptome sequencing to explore the molecular
   mechanism of the pumpkin response to chilling stress. After removing
   the low-quality reads, each library was represented by over 40 million
   high-quality reads ranging from 47.23 to 48.31 (Table [60]S1).
   Subsequently, we acquired approximately 6.9 Gb of high-quality clean
   reads for each library (Table [61]S1). The mapped reads in six samples
   were relatively high in the range of 97.59–98.28% (Table [62]S1). The
   GC percentage varied between 45.35% and 45.89% for all sequences (Table
   [63]S1). The samples were mapped to the C. moschata genome assembly.
   The Pearson’s correlation coefficients (R^2) between biological
   replicates were all > 0.95, demonstrating a high degree of biological
   reproducibility among the samples (Figure [64]S1). A principal
   component analysis (PCA) of the transcript samples showed that the
   first two principal components explained 98.34% of the total variation
   (Figure [65]S1). These results suggest that the sequencing reads
   generated from the RNA samples were of excellent quality.

Analysis of differentially expressed genes in the transcriptome

   DEGs were identified according to the parameters by the conditions of
   |log[2][fold-change (FC)] | ≥ 1 and P-value ≤ 0.05. Overall, 5,780 DEGs
   were identified in the treatment versus control samples, with 2,713
   upregulated and 3,067 down-regulated (Table [66]S2, Figure [67]S2).
   Among the DEGs, 325 were identified as TFs under chilling stress (Table
   [68]S3). Twenty-nine (27 upregulated and two down-regulated) were
   identified and assigned to various TFs (|log[2] FC| ≥ 4), indicating
   that most TFs were significantly induced by chilling stress
   (Fig. [69]2, Table [70]S4). Among 27 TFs, the ERF family was
   predominant, with nine ERF genes upregulated under chilling stress
   (Fig. [71]2), followed by the WRKY and MYB families. Interestingly, the
   MYBAS2 and WRKY49 genes were significantly downregulated, but other TFs
   genes were significantly upregulated compared to the control (Table
   [72]S4). These results suggest that ERF, WRKY, and MYB TFs are involved
   in regulating gene expression under chilling stress.

Fig. 2.

   [73]Fig. 2
   [74]Open in a new tab

   The heat map of expression patterns of high expression of transcription
   factors

   To investigate the functions of the DEGs, we performed a GO functional
   enrichment analysis (count ≧ 10, padj < 0.01). Of the three core GO
   annotation categories, biological processes (BP) comprised 72.2% of the
   total assigned annotations, whereas molecular functions (MF) and
   cellular components (CC) comprised 20.4% and 7.4%, respectively (Table
   [75]S5), leading us to focus on terms associated with biological
   processes. The GO terms in the biological processes category included
   transcription regulatory region DNA binding, response to ABA, positive
   regulation of transcription, DNA-templated, response to cold, and
   ethylene-activated signaling pathway (Fig. [76]3). GO terms, such as
   response to GA, JA, SA, and PAL activity were also significantly
   enriched (Table [77]S5).

Fig. 3.

   [78]Fig. 3
   [79]Open in a new tab

   Functional classification of DEGs by using Gene Ontology (GO)

   To further investigate the functions of DEGs in the response to
   chilling stress in pumpkins, a KEGG pathway enrichment analysis was
   performed. A padj ≤ 0.05 was used as the threshold, resulting in
   significant enrichment of 5,780 DEGs in 18 pathways. Notably, pathways
   related to plant hormone signal transduction and phenylpropanoid
   biosynthesis were closely linked to chilling stress (Fig. [80]4, Table
   [81]S6).

Fig. 4.

   [82]Fig. 4
   [83]Open in a new tab

   Analysis of DEGs by using KEGG terms

   To verify the accuracy of the transcriptome sequencing results, we
   selected 14 DEGs associated with stress, ABC1K1, RVE8, PLDZETA1, ERF1B,
   dnaJ protein homolog, SAMDC, XTH22, ATX1, At1g62810, GH3.6, AED3, TPST,
   PAL, and DREB3. The expression levels (|log[2]FC|) of these genes,
   detected by quantitative reverse transcription (qRT) PCR ranged from 2
   to 6 compared to the control group (Figure [84]S4; melting curves are
   shown in Figure [85]S3). The trends in the expression of the 14 genes
   were consistent with the transcriptome sequencing results, indicating
   that the sequencing results were highly reliable and could be used to
   examine the response mechanism to chilling stress in pumpkins (Figure
   [86]S4). The qRT-PCR results revealed ABC1K1, RVE8, PLDZETA1, ERF1B,
   dnaJ protein homolog, SAMDC, and PAL to be upregulated, whereas AED3,
   XTH22, ATX1, Atlg62810, GH3.6, TPST, and DREB3 were downregulated in
   the treatment samples compared to the control. The qRT-PCR data showed
   a strong positive correlation (R^2 = 0.961) with the transcriptome
   sequencing data.

Analysis of differentially expressed metabolites in the metabolome

   In this study, based on high-resolution mass spectrometry (HRMS)
   technology, a non-targeted quantitative analysis of metabolites in six
   pumpkin samples was performed. To evaluate the differentially
   accumulated metabolites between the control and treatment samples, the
   OPLS-DA model was applied. The established OPLS-DA model showed
   goodness-of fit (Figure [87]S5). A PCA of the metabolomic profiles was
   conducted to provide information about distinctly separate groups
   (Figure [88]S6). Overall, 178 metabolites were detected and are shown
   in volcano plots (Figure [89]S7). Compared with the control samples, 63
   metabolites were upregulated and 115 downregulated in the chilling
   stress-treated samples (Table [90]S7). Of these metabolites, 64.5% were
   found to be annotated based on the HMDB database. These metabolites
   were mainly carboxylic acids and derivatives, followed by organooxygen
   compounds and fatty acyls (Figure [91]S8). Some of the metabolites, the
   carboxylic acids, and derivatives, including mytilin A, L-dopa,
   pyroglutamic acid, L-tyrosine, argininosuccinic acid, L-phenylalanine,
   L-leucine, and L-valine were significantly upregulated, which indicated
   that amino acids are important metabolites involved in the chilling
   stress response (Table [92]S7). Furthermore, the only plant hormone to
   accumulate in pumpkin leaves under chilling stress was SA (Table
   [93]S7).

   To identify the significant changes in metabolite expression under
   chilling stress, a KEGG pathway analysis of DEMs was performed to
   identify significantly enriched metabolic pathways (Fig. [94]5, Table
   [95]S8). Our metabolomic analysis revealed significant changes in
   several metabolic pathways including isoquinoline alkaloid
   biosynthesis, aminoacyl-tRNA biosynthesis, nitrogen metabolism,
   tyrosine metabolism, alanine, aspartate and glutamate metabolism, and
   phenylalanine metabolism. These outcomes indicate that chilling stress
   mainly regulates the metabolism of some amino acids. Furthermore,
   several key metabolites in the phenylpropanoid pathway, including
   phenylalanine and tyrosine, accumulate in response to chilling stress.

Fig. 5.

   [96]Fig. 5
   [97]Open in a new tab

   The KEGG pathway enrichment of DEMs in the chilling stress treatment
   groups compared to the control group. The P-value is the significance
   of enrichment of this metabolic pathway. The ordinate is the name of
   the metabolic pathway; the abscissa is the impact (impact = the number
   of significantly different metabolites/the total number of metabolites
   in this pathway). The color from red to orange indicates that –ln
   P-value decreases in turn

Integrative analysis of genes and metabolites

   For the integrative analysis of genes and metabolites, nine common KEGG
   enrichment pathways of the DEMs and DEGs were used to show the overlap
   of these two omics (Fig. [98]6, Table [99]S9). The DEMs and DEGs were
   mainly enriched in plant hormone signal transduction, phenylpropanoid
   biosynthesis, and phenylalanine metabolism (Fig. [100]7).

Fig. 6.

   [101]Fig. 6
   [102]Open in a new tab

   Venn diagram of DEGs and DEMs involved in the KEGG pathway

Fig. 7.

   [103]Fig. 7
   [104]Open in a new tab

   The top 10 KEGG pathways with the highest number of DEGs and DEMs
   identified by trancriptome and metabolome

Responses of plant hormone signal transduction to chilling stress

   GO and KEGG pathway enrichment analyses revealed that plant hormone
   signal transduction was significantly enriched when comparing the
   treatment and control groups. In the GA signaling pathway, three GID1
   genes and one DELLA gene were upregulated (Fig. [105]8). In the ABA
   signaling pathway, the orthologs of three ABA receptors (PYL/PYR), five
   protein phosphatases (2 C/PP2C), and two ABF genes were upregulated
   (Fig. [106]8). In the ET signaling pathway, three ETRs, two SIMKKs, two
   EIN3-binding F-box proteins/EBF1/2, and three ethylene-responsive TF
   1ERF1B genes were upregulated in response to chilling stress
   (Fig. [107]8). In the JA signaling pathway, one JAR1, eight JAZs, and
   TF MYC2 genes were significantly upregulated (Fig. [108]8). In
   particular, ERF1B and MYC2 showed very high expression levels (Table
   [109]S4). Moreover, the accumulation of SA induced the upregulated
   expression of downstream regulatory protein NPR3, NPR3-like, and TF
   TGA21 genes under chilling stress (Fig. [110]8). Although the contents
   of ABA, ET, JA, and GA were nearly unchanged in pumpkin in response to
   chilling stress (Table [111]S7), the SA content was significantly
   increased after exposure to 5 °C for 24 h (Table [112]S7). Although the
   contents of many metabolites varied, only the transcripts and the
   content of SA varied simultaneously. Furthermore, PAL genes were
   significantly induced, and enrichment of PAL activity in the GO
   analysis revealed its contribution to SA synthesis (Table [113]S5).
   These findings strongly support the involvement of the SA signaling
   pathway in the response to chilling stress in C. moschata.

Fig. 8.

   [114]Fig. 8
   [115]Open in a new tab

   The heat map of expression patterns of DEGs involved in plant hormone
   signaling transduction pathway. Red boxes indicate upregulated DEGs
   under chilling stress in C. moschata. Blue boxes indicate
   down-regulated DEGs. Arrows with solid lines indicate direct
   activation. Lines ending with a bar indicate negative regulation. The
   block represents the log[2]FKPM of genes, respectively

Integrated analysis of genes and metabolites related to the phenylpropanoid
biosynthesis pathway in response to chilling stress

   Based on the integrated analysis of gene expression patterns and
   metabolite profiles, among the 66 DEGs in the phenylpropanoid
   biosynthesis pathway, 12 PALs, one cytochrome P450 CYP73A100-like, and
   five glucosyltransferase-likes genes were significantly upregulated
   compared to the control (Fig. [116]9, Table [117]S6). Seven upregulated
   and 14 downregulated genes were associated with peroxidase activity.
   The accumulation of phenylalanine, tyrosine, and
   trans-2-hydroxycinnamate (Table [118]S8) led to the upregulated
   expression of PAL, 4CL, CCR, and POD, key genes in the phenylpropanoid
   biosynthesis pathway. These expression changes due to chilling stress,
   also facilitated the synthesis of key metabolites (Fig. [119]9). Our
   results indicate that chilling stress activates the phenylpropanoid
   biosynthesis pathway in pumpkin leaves.

Fig. 9.

   [120]Fig. 9
   [121]Open in a new tab

   Diagram of phenylpropanoid biosynthesis pathway with their related
   differentially expressed genes and differentially expressed
   metabolites. The below heatmap shows the expression changes of genes
   involved in the phenylpropanoid biosynthesis pathway. PAL,
   phenylalanine ammonia-lyase; BGLU, beta-glucosidase; 4CL, 4-coumarate:
   CoA ligase; CCR, cinnamoyl CoA reductase. Metabolites labeled in red
   indicate the content of metabolites that significantly increased in C.
   moschata. The block represents the log[2]FKPM of genes, respectively

Discussion

   Under low-temperature conditions, plants enhance their cold tolerance
   by inducing or repressing gene expression [[122]27, [123]28]. The
   release of the C. moschata reference genome has facilitated a
   comprehensive transcriptomic and metabolomic analysis, providing
   insights into chilling stress tolerance [[124]29]. In this study, 5,780
   DEGs were activated in pumpkin leaves following exposure to chilling
   stress. Furthermore, TFs such as MYB, WRKY, ERF, and HSF were highly
   expressed, and metabolites such as phenylalanine, tyrosine, and SA
   accumulated, leading to the change in physiological changes in the
   pumpkin seedling leaves. The integrated transcriptomics and
   metabolomics analysis revealed SA signaling pathway and phenylpropanoid
   biosynthesis contribute to the chilling stress response in pumpkins.

   Regulatory genes, particularly TFs, play an essential role in plant
   stress response [[125]30]. For example, MYB-like17, a cold-inducible
   gene in Brassica napus L., and AhMYB30 from peanuts increase the
   resistance of these plants to freezing stress [[126]31, [127]32], while
   PmWRKY57 in Prunus mume and VvWRKY28 in Vitis vinifera improve cold
   tolerance in Arabidopsis thaliana [[128]33, [129]34]. In addition,
   NtbHLH123 in Nicotiana tabacum and IbbHLH79 in sweet potato demonstrate
   a positive regulation effect during cold stress [[130]35, [131]36].
   Furthermore, CdEFR1 in Cynodon dactylon positively regulates the
   plant’s cold response probably by activating stress-related genes,
   PODs, CBF2, and lipid transfer proteins [[132]37]. HSFs such as Oryza
   sativa OsHSFA2a, OsHSFA3, and Cucumis sativus HSFA1d are induced after
   exposure to low temperatures [[133]38–[134]40]. In the present study,
   numerous TFs were identified (Table [135]S3), and those from the MYB,
   WRKY, ERF, and HSF families were highly expressed (Fig. [136]2).
   Therefore, they may be involved in the response to chilling stress in
   C. moschata. In the future, these TFs may be used as candidate genes in
   plant breeding programs to improve the cold tolerance of pumpkins.

   GO and KEGG pathway enrichment analyses revealed that ABA, ET, GA, JA,
   and SA signaling pathways were activated under chilling stress (Table
   [137]S5, Table [138]S6). Among the plant hormones, SA plays diverse
   regulatory roles in plant growth and stress resistance. It has been
   demonstrated that SA can be synthesized via PAL in watermelon, with
   conversion efficiency from benzoic acid to SA [[139]15, [140]41]. In
   our study, the accumulation of SA and the GO enrichment of PAL activity
   (Table [141]S5, Table [142]S7) suggest that SA synthesis occurs via the
   PAL pathway in C. moschata. Furthermore, the NPR3 and TGA21 genes
   associated with the SA signaling pathway in C. moschata were induced by
   chilling stress (Fig. [143]8). The functions of the SA receptor NPR3
   and TGAs/WRKYs TFs in SA signaling have previously been established
   [[144]42]. NPR3 functions with TGAs as transcriptional corepressors of
   SA response genes and is inhibited by SA [[145]43]. SA acts through its
   signaling pathway ultimately resulting in the buildup of secondary
   metabolites [[146]44]. Therefore, the application of SA enhances plant
   resilience to chilling stress.

   Phenylpropanoid biosynthesis is a crucial pathway in the synthesis of
   secondary metabolites, playing a significant role in plant development
   and stress responses [[147]45]. Previous studies have found that cold
   stress influences the expression of genes involved in the
   phenylpropanoid biosynthesis pathway, including PAL, cinnamate
   4-hydroxylase (C4H), and 4CL genes and the accumulation of coumarins,
   flavonoids, and lignin facilitates the adaptation to low-temperature
   environments in Arabidopsis, cucumber, and winter rapeseed
   [[148]46–[149]48]. PAL is the first enzyme that catalyzes the
   deamination of l-phenylalanine to trans-cinnamate to initiate the
   synthesis of all phenylpropanoids. Moreover, POD catalyzes the
   polymerization of coumarin, terpineol, and mustard alcohol to produce
   macromolecular lignin [[150]49]. In this study, we found that the
   expression levels of PAL, CCR, 4CL, and POD were upregulated under
   chilling stress, with increased phenylalanine, tyrosine, and
   trans-2-hydroxycinnamate accumulation seen in C. moschata seedlings
   (Fig. [151]9). Based on these results, phenylpropanoid biosynthesis may
   contribute to chilling stress resistance in pumpkins.

Conclusions

   The physiological, molecular, and metabolic response of C. moschata to
   chilling stress were investigated. TFs such as ERF, bHLH, WRKY, and MYB
   may contribute to molecular breeding for chilling stress tolerance of
   pumpkin in the future. It was found that the SA signaling pathway and
   phenylpropanoid biosynthesis play key roles in regulating the responses
   of pumpkins to chilling stress. In addition, POD activity decreased
   significantly. The findings of this study offer valuable insights into
   the comprehensive response of C. moschata to chilling stress. However,
   further research is needed to elucidate the complex regulatory
   networks, and facilitate the development of breeding programs to
   develop resilient cultivars.

Materials and methods

Plant material and treatments

   The C. moschata inbred line IL1, purified through more than eight
   consecutive generations of self-pollination, was obtained from the
   Qingdao Academy of Agricultural Sciences of Shandong Province, China.
   The IL1 seedlings were grown in a growth chamber under the following
   conditions, temperature: 25/18°C (day/night); light/dark: 16 h/8 h;
   relative humidity: 80%; and light intensity: 24,000 lx [[152]26]. After
   three weeks of growth, the “second-true-leaf stage” seedlings, were
   subdivided into chilling stress treatment (5 °C) and control (25 °C)
   groups [[153]26]. Each group contained three biological replicates.
   After 24 h of treatment, the second expanding leaves were collected for
   physiological, qRT-PCR, and metabolite analyses and mRNA-sequencing.

Physiological index measurements and statistical analysis

   The second-true leaves of pumpkin seedlings were sampled before and
   after 24 h of chilling stress for physiological analysis. The REC was
   measured according to a published method [[154]50], with minor
   modifications, using a PHS-25 conductivity meter (Yueping, Shanghai,
   China). The pumpkin leaves were sliced into thin segments and shaken at
   150 rpm for 2 h in 10 mL of ddH[2]O, at room temperature. The REC (C0)
   value of ddH[2]O was measured as the control using a DDS-IIA detector
   (Shanghai, China). The initial REC (C1) was taken before boiling, and
   the final REC (C2) was measured after boiling for 30 min. The REC was
   calculated according to the following equation: REC (%) = (C1–C0)/(C2
   –C0) × 100.

   The index of total chlorophyll, indicated by the Soil Plant Analysis
   Development (SPAD) values, was measured by averaging three readings
   taken from the midpoints of pumpkin leaves using a portable chlorophyll
   meter (SPAD-502 Plus, Minolta Camera Co. Ltd., Japan) [[155]51].

   The MDA content was measured using an MDA assay kit (Nanjing Jiancheng
   Bioengineering Institute, China). A 100 mg sample of pumpkin leaves was
   weighed and 1 mL of extraction solution was added before centrifugation
   at 8000 × g for 10 min at 4 °C. The supernatant was removed, and the
   reagents were added sequentially according to the manufacturer’s
   protocol. The solution was kept in a 100 °C water bath for 60 min and
   then centrifuged at 10,000 × g. After centrifugation at room
   temperature for 10 min, the absorbance at 532 nm was measured, and the
   MDA content was calculated according to the manufacturer’s
   instructions.

   SOD activity was determined by spectrophotometry (Jinghua, Shanghai,
   China) using a SOD activity assay kit (Nanjing Jiancheng Bioengineering
   Institute, China). A 100 mg sample of pumpkin leaves was weighed, and 1
   mL of extraction solution was added for ice bath homogenization before
   centrifugation at 8000 × g for 10 min at 4 °C. The supernatant was
   removed, and the reagents were added sequentially according to the
   manufacturer’s protocol. The solution was placed in a 37 °C water bath
   for 30 min, and the absorbance was then measured at 560 nm. The SOD
   activity was then calculated according to the manufacturer’s
   instructions.

   CAT activity was determined by spectrophotometry (Jinghua, Shanghai,
   China) using a CAT assay kit (Nanjing Jiancheng Bioengineering
   Institute, China). A 100 mg sample of pumpkin leaves was weighed and
   0.9 mL of 100 mM phosphate buffer, pH 7.0, was added for ice bath
   homogenization. The sample was centrifuged at 4000 × g for 10 min at
   4 °C before removing the supernatant and adding the reagents in
   sequence according to the manufacturer’s instructions. The solution was
   mixed immediately, then incubated at 37 °C for 1 min. Further reagents
   were added to the mixture according to the subsequent steps in the
   protocol. Finally, the absorbance was measured at 405 nm and the CAT
   activity was calculated according to the manufacturer’s instructions.

   POD activity was measured at 420 nm using a POD assay kit (Nanjing
   Jiancheng Bioengineering Institute, China). A 100 mg sample of pumpkin
   leaves was weighed and 0.9 mL of 100 mM phosphate buffer, pH 7.0, was
   added for ice bath homogenization. The sample was centrifuged at 3500 ×
   g for 10 min at 4 °C before removing the supernatant and adding the
   reagents in sequence according to the manufacturer’s instructions. The
   solution was incubated at 37 °C for 30 min. The absorbance was measured
   at 420 nm, and the POD activity was calculated according to the
   manufacturer’s instructions.

   The soluble sugar content was measured using a Plant Soluble Sugar
   Content test kit (Nanjing Jiancheng Bioengineering Institute, Nanjing,
   China). First, 100 mg of pumpkin leaves were weighed and combined with
   1 mL of dH[2]O for homogenization. The sample was centrifuged at 4000 ×
   g for 10 min, the supernatant was removed, and the reagents were added
   according to the manufacturer’s instructions. Finally, the mixture was
   heated at 100 °C for 10 min and left to cool naturally. The absorbance
   was measured at 620 nm, and the content of soluble sugars was
   calculated according to the manufacturer’s instructions.

   Soluble proteins were measured according to the Bradford [[156]52]
   method using albumin as the standard. The total soluble protein
   absorbance was recorded at 595 nm using a spectrophotometer (Jinghua,
   Shanghai, China), and the concentration is presented as mg g^− 1 fresh
   weight.

   The obtained data were subjected to a two-tailed t-test (5%
   probability) in SPSS 22.0 (version 20.0; IBM Inc., USA). P < 0.05 was
   considered statistically significant.

Transcriptome sequencing and data analyses

   Three biological replicates of control and treatment samples stored
   frozen in liquid nitrogen were used to isolate total RNA using a
   Tiangen RNAprep Pure Plant Kit (Tiangen, China). A TruSeq Stranded mRNA
   LTSample Prep Kit (Illumina, San Diego, CA, USA) was used to construct
   the cDNA libraries. The libraries were sequenced using an Illumina
   HiSeqTM 2500 sequencing platform. The RNA-seq samples were mapped to
   the C. moschata genome
   ([157]https://www.ncbi.nlm.nih.gov/assembly/GCA_002738365.1) as the
   reference genome using Tophat [[158]53].

   Gene expression calculations used fragments per kilobase per million
   (FPKM) values, which can be compared directly between samples to assess
   gene expression differences. The DESeq R package (1.10.1) was then used
   to identify genes exhibiting differential expression levels across
   sample groups. P < 0.05 and |log[2]FC| ≥ 1 served as the criteria for
   screening DEGs. The function of DEGs was annotated using the public
   KEGG and GO databases [[159]54]. All sequencing data were deposited in
   the NCBI Short Read Archive (SRA) database under the BioProject ID:
   PRJNA1044531.

qRT‑PCR validation

   To verify the reliability of the transcriptome results, 14 DEGs
   involved in cold stress, ABC1K1, RVE8, PLDZETA1, ERF1B, dnaJ protein
   homolog, SAMDC, XTH22, ATX1, At1g62810, GH3.6, AED3, TPST, PAL, and
   DREB3, were selected to validate the RNA-seq results. These genes were
   analyzed by qRT-PCR. The specific primers for the selected genes were
   designed using Primer Premier 5.0. All primers are listed in Table
   [160]S10. The qRT-PCR reactions were normalized using C. moschata
   Actin7 as a reference gene for all comparisons [[161]26, [162]55]. Each
   reaction contained 10 µL of 2 × SYBR Premix ExTaq™, 0.4 µL of 10 µM of
   each gene-specific primer, 8.2 µL of ddH[2]O, and 1 µL of cDNA template
   in a final volume of 20 µL. The PCR conditions were as follows: initial
   denaturation at 95 °C for 30 s, 40 cycles of 95 °C for 5 s, 60 °C for
   30 s. Three biological replicates and three technical replicates were
   performed for each of the analyzed genes. The cycle threshold (Cq)
   value was calculated using the 2 − ^ΔΔCt method to determine the
   relative expression levels of genes. A melting curve was created for
   each PCR product to avoid non-specific amplification. The expression
   levels of the genes in different samples were calculated and compared.

Metabolite extraction and detection

   The second-true leaves of pumpkin seedlings were sampled after
   treatment for 0 h (control) and 24 h and quickly placed in liquid
   nitrogen. Five independent biological replicates were set up for each
   treatment. The freeze-dried samples were crushed using a mixer mill at
   a frequency of 45 Hz for 4 min. Then, 50 mg of the sample was added to
   a solution containing extracted material [[163]56]. Final samples were
   stored at -80 °C until use in an ultra-high performance liquid
   chromatography (UHPLC)-tandem mass spectrometry (MS/MS) analysis.
   Multiple quality control (QC) sample aliquots were prepared from a
   single QC sample.

Metabolomic analysis

   The PCA and OPLS-DA were performed. Further identification and
   quantification of metabolites were performed. The significant DEMs with
   variable importance in the projection (VIP) > 1 and P < 0.05 were
   determined based on the VIP obtained by the OPLS-DA model and the
   P-value generated by Student’s t-test. Liquid chromatography-mass
   spectrometry (LC-MS) was used to detect the metabolites [[164]56]. The
   Lipidmaps database ([165]https://www.lipidmaps.org/), KEGG database
   ([166]http://www.genome.jp/kegg/), and HMDB database
   ([167]http://www.hmdb.ca/) were used to annotate the different
   metabolites. DEMs from the treatment and control groups were submitted
   and mapped onto their biochemical pathways through metabolic enrichment
   and pathway analyses based on the KEGG database
   ([168]http://www.genome.jp/kegg/).

Statistical analysis

   All experiments in this study included three biological replicates. All
   data were analyzed using SPSS 26.0. Data are expressed as the
   mean ± standard error (SE) of the three replicates. Data were analyzed
   using Student’s t-tests. The significance of differences between groups
   was analyzed by variance analysis using GraphPad Prism software
   (version 8.0). Origin 8.0 (Microcal Software Inc., Northampton, MA,
   USA) was used to construct the figures. Transcriptome and metabolome
   heatmaps were drawn using the R pheatmap software package.

Electronic supplementary material

   Below is the link to the electronic supplementary material.
   [169]12864_2024_10939_MOESM1_ESM.docx^ (3MB, docx)

   Supplementary Material 1: Figure S1. Correlation and principal
   component analysis (PCA) of all samples. (A) Pearson correlation among
   six samples. IL1C1, IL1C2, and IL1C3 represent three biological
   replicates of the control group; IL1T1, IL1T2, and IL1T3 represent
   three biological replicates of the treatment group. (B) PCA of all
   transcripts of treatment and control groups. Figure S2. A Volcano map
   of DEGs in the IL1. The x-axis is the log[2] scale of the fold change
   of gene expression (log[2](fold change)). The y-axis is the minus
   log[10] scale of the P-values, which indicates the significant level of
   expression difference. Genes significantly differentially expressed are
   indicated by red dots (upregulated) and green dots (downregulated).
   Genes with no significant differential expression are indicated by gray
   dots. Figure S3. The melting curve of actin, AED3, ABC1K1, RVE8, XTH22,
   ATX1, At1g62810, GH3.6, PLDZETA1, TPST, ERF1B, dnaJ protein homolog,
   SAMDC, PAL, and DREB3. Figure S4. qRT-PCR validation of gene expression
   in treatment vs. control, linear relationship between qRT-PCR and
   RNA-seq. Relative quantification was obtained through 2^−ΔΔCT method
   using Actin reference gene. Data represent the average from three
   biological replicates and the error bars indicate the standard
   deviation (± SD). Figure S5. The scores plots of the OPLS-DA model
   showing an obvious separation between IL1T (green squares) and IL1C
   (bule circles). Figure S6. Principal component analysis (PCA) of
   metabolites identified by chilling stress-treatment (IL1T, green
   squares) and control (IL1C, bule circles). Figure S7. Volcano plots for
   differentially expressed metabolites (DEMs) between chilling stress
   steatment and control.
   [170]12864_2024_10939_MOESM2_ESM.xlsx^ (864KB, xlsx)

   Supplementary Material 2: Table S1. Summary of sequencing data of C.
   moschata leaves under chilling stress. Table S2. Complete list of
   differentially expressed genes (DEGs) regulated by chilling stress in
   C. moschata leaves. Table S3. Transcription factors of differential
   expression under chilling stress. Table S4. High expression of
   transcription factors. Table S5. Gene Ontology (GO) analysis of all
   annotated DEGs. Table S6. Kyoto Encyclopaedia of Genes and Genomes
   (KEGG) pathway enrichment analysis of all annotated DEGs. Table S7.
   Differential metabolites under chilling stress. Table S8. KEGG pathway
   enrichment analysis of DEMs. Table S9. The overlapped pathways from
   DEGs and DEMs involved in the KEGG pathway. Table S10. Information of
   gene specific primers for qRT-PCR.

Acknowledgements