Two luffa varieties were used in this study: browning-sensitive LQY and browning-tolerant JD. All luffa plants were cultivated in the experimental field of Chengdu Agricultural College, China. Luffa fruits were harvested at 15, 20 and 45 days after flowering (April to May 2021), designated as S1, S2 and S3, respectively. These stages were selected based on the life cycle of luffa to capture the onset, progression, and culmination of physiological and biochemical changes potentially associated with fruit browning. At the same time, flowers, leaves and stems were collected at the corresponding time points. The samples were frozen immediately in liquid nitrogen, followed by separation of the seeds and receptacles under frozen conditions.

Fresh luffa fruits were collected 45 days after flowering, and then manually sliced through the middle. The slices were left at room temperature for 15 minutes to observe browning in the cross-section.

PPO activity was measured using spectrophotometry, following the method of González et al. (2020). Specifically, 1 g of frozen sample was homogenized in 5 mL of cold PPO extraction buffer. The homogenate was transferred to a 2-mL centrifuge tube and centrifuged at 8,000 rpm at 4°C for 10 min. Next, 150 µL of the supernatant was mixed with 750 µL of PPO assay buffer (50 mM sodium phosphate buffer, pH 6.0; 0.1% SDS; 15 mM 4-methylcatechol) in a new tube. The mixture was incubated at 37°C for 10 min, and immediately transferred to a 97°C boiling water bath for 5 min. After centrifugation at 10,000 rpm for 10 min at room temperature, the supernatant was transferred to a quartz cuvette, and absorbance was measured at 420 nm (A420 nm). PPO enzymatic activity was defined as one unit (1 U), representing the amount of enzyme required to increase absorbance at 420 nm by 0.01 per min in 1 g of tissue within a 1-mL reaction system.

ROS level measurement was performed as previously described with slight modifications (AbdElgawad et al., 2016). Frozen samples were ground in liquid nitrogen, and 1 g of the resulting powder was added to 9 mL of sodium phosphate buffer (100 mM, pH 7.0). The mixture was vortexed at 3500 rpm, and centrifuged for 10 min to obtain the supernatant as a 10% homogenate. Subsequently, 20 µL of the supernatant and 20 µL of the enzyme solution were added to 200 µL of substrate application solution. The mixture was incubated at 37°C for 20 min. Absorbance was then measured at 450 nm using a microplate reader (MK3; Thermo Fisher Scientific, Waltham, MA, USA). ROS concentration was calculated as: c=A/(ε×b) where ‘c’ represents the sample concentration, A is the absorbance value, ϵ is the molar extinction coefficient, and b refers to the path length of the sample. PPO and ROS measurements were conducted in triplicate, including three biological and technical replicates.

Total RNA was extracted from each frozen fruit sample using an RNA prep Pure Plant Plus Kit (Tiangen, Beijing, China), and mRNA was obtained using oligo (dT) beads (Beyotime, Shanghai, China). The mRNA was reverse transcribed to cDNA, and the library was constructed using an NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, USA) according to manufacturer’s instructions. Three biological replicates were included to ensure the reliability of transcriptome data. Qualified libraries were paired-end sequenced on the Illumina HiSeq 6000 sequencing platform at Yunnan Pulis Biotechnology Co. Ltd. (Kunmin, China), yielding over 6 GB of raw data per sample.

Clean reads obtained after filtering the raw data were aligned to the luffa reference genome (https://www.ncbi.nlm.nih.gov/, PRJNA596077) using HISAT2 (v2.1.0) (Kim et al., 2019). The number of reads that mapped to each gene was quantified using HTSeq (v0.6.1) (Anders et al., 2015). The fragments per kilobase of transcript sequence per million base pairs sequenced (FPKM) of each gene were then calculated based on gene length and the number of mapped reads. Principal component analysis (PCA) was performed using the prcomp function (Maćkiewicz and Ratajczak, 1993) in R to compute the variance ratio explained by each principal component (PC). DEGs were identified using the DESeq2 package (v4.0.4) (Love et al., 2014) in R. P-values were adjusted using the Benjamini-Hochberg procedure to control the false discovery rate (FDR). DEGs were considered significant with P-adjusted value < 0.05 and |log2(fold change)| > 1. Functional enrichment analysis of the DEGs was performed using the clusterProfiler package (Yu et al., 2012) in R, with the KEGG database as reference.

Freeze-dried fruit material was crushed for 1.5 min at 30 Hz using a mixer mill (MM 400, Retsch, Haan, Germany) with zirconia beads. A total of 100 mg of the powder was extracted overnight with 1.2 mL of 70% aqueous methanol at 4°C. The extracts were then filtered (SCAA-104, 0.22 μm pore size; ANPEL, Shanghai, China) after centrifugation at 12000 rpm for 10 min. Three biological replicates were included.

The samples were examined using a UPLC–ESI–MS/MS system (UPLC, Nexera X2, SHIMADZU, Japan; MS, 4500 Q TRAP, Applied Biosystems, USA). The UPLC system used an Agilent SB-C18 column (1.8 µm, 2.1 mm*100 mm) (Agilent, Santa Clara, CA, USA). The mobile phase consisted of solvent A (0.1% formic acid in water), and solvent B (0.1% formic acid in acetonitrile). A gradient program was adopted, starting with 95% A and 5% B. A linear gradient was applied over 9 min, gradually reaching 5% A and 95% B, which was maintained for 1 min. The composition was adjusted back to 95% A and 5% B within 1.1 min and held for 2.9 min. The column oven was maintained at 40°C, and the injection volume was 4 µL. The effluent was then introduced into an ESI-triple quadrupole-linear ion trap (QTRAP)-MS. PCA was performed using the prcomp function in R, with the data subjected to unit variance scaling prior to analysis. Metabolites with VIP values ≥ 1 and absolute Log2 fold changes (|Log2FC|) ≥ 1 were considered significantly regulated across groups. The annotated metabolites were linked to the KEGG database. Pathways with substantially upregulated metabolites were then loaded into metabolite set enrichment analysis (MSEA), and their significances were established using the p-values of the hypergeometric test.

Pearson correlation coefficients between DEGs and DAMs were calculated using the cor function in R. DEGs and DAMs with correlation coefficients greater than 0.80 and p-values less than 0.05 were selected for further analysis. Enrichment analysis of the significantly correlated DEGs and DAMs was subsequently conducted using the KEGG database to identify associated biological pathways.

Total RNA was extracted from samples of fruits, flowers, leaves and stems of the two luffa varieties using RNA simple total RNA kit (Tiangen Biotech Beijing Co.,

Ltd., Beijing, China). The extracted RNA was reverse transcribed into cDNA and passed through EasyScript^R superMix (Transgen Biotech, Beijing, China) according to the manufacturer’s instructions. q-PCR was performed on an SYBR Green system (TaKaRa, Dalian, China). The cycling conditions were 3 min of predenaturation at 95°C, then 95°C for 10 s, 60°C for 30 s, and 72°C for 90 s (40 cycles). To ensure accurate normalization in our q-PCR analysis, we selected Maker00006543 (40S ribosomal protein SA-like) as the reference gene based on its stable and relatively high expression across all samples in our transcriptome dataset. While commonly used reference genes (Chen et al., 2021) failed to exhibit stable and reliable amplification in our experimental system, Maker00006543 demonstrated minimal variation in expression and consistent amplification. The FPKM values across samples are provided in Supplementary Table 1 . Primers were designed using Primer Premier 6.0 (Premier Biosoft, International, Palo Alto, California, USA) and are listed in Supplementary Table 2 . Three biological and technical repeats were performed. q-PCR results were analyzed using the 2^−ΔΔCt method (Kim et al., 2021) and statistical significance was determined using Student’s t-test.

To investigate differences in browning susceptibility between LQY and JD, fruit samples were cross-sectioned. After 15 min in temperature, LQY slices exhibited varying degrees of browning, while browning-tolerant JD showed no visible signs of browning ( Supplementary Figure S1A ). Moreover, LQY demonstrated notably higher PPO enzyme activity ( Supplementary Figure S1B ) and ROS levels ( Supplementary Figure S1C ) than JD across S1, S2, and S3. These results suggest that LQY may contain higher concentrations of polyphenols and flavonoids, which are more susceptible to oxidative changes.

PPO activity and ROS levels from earlier experiments indicated that LQY and JD may exhibit distinct polyphenol metabolism and antioxidant-related gene profiles. Therefore, a comparative transcriptome analysis was conducted at three developmental stages for both luffa varieties, with most samples generating more than 6 GB of clean data ( Supplementary Table 3 ). PCA showed clear clustering within each variety as well as significant differentiation between the two varieties ( Figure 1A ). To investigate transcriptome changes across developmental stages, we performed differential expression analysis. The results revealed a higher number of DEGs in the JD variety at both S2 and S3 compared with S1, with more than double the number detected than in LQY. This may suggest more pronounced transcriptome changes in the JD variety. In addition, more DEGs were observed between S2 and S1 than between S3 and S2, implying that the most significant transcriptome changes occur during the transition from S1 to S2 ( Figure 1B ). A global expression analysis of all DEGs revealed variety-specific differences between LQY and JD ( Figure 1C ). Functional enrichment analysis of the DEGs indicated that the S1-specific DEGs were primarily associated with pathways such as phenylpropanoid biosynthesis, glutathione metabolism, oxidative phosphorylation, and plant hormone signaling pathways. At S2, DEGs were predominantly involved in beta-alanine metabolism, propanoate metabolism, and the degradation of valine, leucine, and isoleucine. At S3, the enriched pathways included phenylalanine, tyrosine, and tryptophan biosynthesis, alongside phenylpropanoid biosynthesis ( Figure 1D ). The results highlighted that both S1 and S3 exhibited active phenylalanine metabolic processes, suggesting a potential link between the phenylalanine metabolic pathway and the regulation of polyphenol and flavonoid metabolism.

Due to the dynamic variations in the expression profiles between LQY and JD, the k-means clustering method was applied to perform cluster analysis of all DEGs ( Figure 2 ). A total of 12 modules were identified, with 5 modules (Module 5, 6, 9, 10, and 11) showing opposing or divergent expression patterns between LQY and JD. Functional enrichment analysis revealed distinct characteristics among the five modules. Notably, Module 5 was upregulated in JD but not in LQY during S2, with genes in this module associated with phosphotransferase activity, response to oxidative stress, response to chitin, and response to oomycetes. This suggested that at S2, a considerable number of genes linked to stress tolerance were upregulated in JD, which may contribute to its browning tolerance properties. Protein-protein interaction (PPI) analysis of the DEGs in Module 5 revealed that SYP121 and CLF may be crucial genes during the process. Furthermore, the DEGs in Module 6 were upregulated at both S2 and S3 in JD but not in LQY. Analysis of their roles suggested that these genes may be associated with cellular anatomical entities, intracellular organelles, and related processes. HEN2, AT1G26370, and NRPB7 may play key roles in these functions. Module 9 and Module 11 revealed comparable change tendencies, with DEGs being more abundant in LQY. Functional analysis displayed that the DEGs may be associated with pathways such as oxidoreductase activity, galactose metabolism, and folate biosynthesis. UBP26 and AT3G11960 were identified as key genes in Module 9, while AT3G62870 and AT5G02960 were recognized as key genes in Module 11. In Module 10, the JD variety exhibited higher expression at all stages compared with LQY, peaking at S2 before decreasing slightly. This module may be associated with pyrophosphatase activity, ATPase activity, and ncRNA processing pathways, with AT2G40360 and AT2G20635 as pivotal genes.

Transcriptome analysis suggested that genes involved in polyphenolic biosynthesis may differ significantly between LQY and JD, potentially leading to variations in polyphenolic compound accumulation. Subsequent metabolomic analysis identified 13, 5, and 8 variety-specific polyphenol metabolites at S1, S2, and S3, respectively. Similarly, 18, 14, and 9 variety-specific flavonoid metabolites were detected during the same stages. In LQY, 12 out of 13 polyphenolic metabolites were upregulated at S1, while 7 out of 8 were upregulated at S3. In contrast, 13 of 18 flavonoid metabolites were upregulated in JD at S1 ( Supplementary Figure S2 ). These periodic fluctuations in polyphenol and flavonoid metabolite levels indicate that LQY accumulated large amounts of polyphenolic metabolites at S1 and S3, while JD accumulated substantial amounts of flavonoid metabolites at S1.

To investigate the molecular mechanisms underlying secondary metabolite production, we explored expression changes in the structural genes of the polyphenol and flavonoid metabolism pathways. A total of 18 genes were identified to be involved in tyrosine and hydroxytyrosol production. PPO expression peaked at the S3 stage in both LQY and JD, exhibiting stage-specific regulation. Moreover, copper amine oxidase (CuAO) and tryptophan decarboxylase (TDC) showed higher expression in JD, with a distinct increase at S3. The aldehyde dehydrogenase (ALDH) genes, ALDH-1, ALDH-3, ALDH-4, ALDH-5, ALDH-6, and ALDH-7 showed significant upregulation at S3 in LQY, potentially due to the higher polyphenol metabolic activity in this variety ( Figure 3 ).

In the flavonoid metabolism pathway, we identified 57 genes from eight families, including 5 HCT (hydroxycinnamoyltransferase) genes, 17 C3H (cinnamate 3-hydroxylase) genes, 4 CSE (caffeoyl shikimate esterase) genes, 17 PAL (phenylalanine ammonia-lyase) genes, 4 4CL (4-coumarate-CoA ligase) genes, 4 CHS (chalcone synthase) genes, 4 F3H (flavanone 3-hydroxylase) genes, and 4 UGT (uridine diphosphate glycosyltransferase) genes. Further examination of structural gene expression changes in the flavonoid metabolic pathway revealed significant differences in the key genes between LQY and JD. Notably, HCT, CHS, CSE, and UGT were highly expressed in LQY at S1, while F3H was highly expressed in JD at S3. Other genes, such as PAL, C3H, and 4CL, exhibited varying expression patterns between species, which may contribute to the significant accumulation of flavonoid metabolites in JD ( Supplementary Figure S3 ).

To identify TFs linked to polyphenol and flavonoid metabolic pathways, as well as potential pathways involved in luffa browning, we performed a transcriptome-based analysis to assess the differential expression of TFs that may regulate these processes. The results revealed differential expression of TFs between LQY and JD, including members of the MYB, WRKY, bHLH, and WD40 families. Specifically, 34 MYB, 15 bHLH, 19 WD40, and 14 WRKY TFs were identified. Notably, bHLH2–bHLH6 were highly expressed at S1 in LQY, whereas bHLH8–bHLH12 were predominantly expressed at S1 in JD. WRKY1–WRKY9 showed peak expression at S3 in JD, while WRKY6–WRKY14 were primarily expressed at S1 in the LQY variety. Furthermore, MYB1–MYB18 were strongly expressed at S2 and S3 in JD, whereas the remaining MYB TFs exhibited stage-specific expression in LQY. The WD40 TFs displayed more variable expression patterns, reflecting both variety- and stage-specific regulation ( Figure 4 ).

To investigate the mechanisms and patterns of metabolite changes associated with luffa browning, we performed a comparative metabolomic analysis of LQY and JD fruit samples. PCA was able to distinguish the metabolic profiles of the two luffa varieties across the three developmental stages. PC1, which accounted for 72.67% of the variation, primarily separated the samples by variety. The close clustering of replicates within each group indicated strong reproducibility. Specific metabolites, such as CHCHON233 and CHCHON343, were major contributors to the observed variation, highlighting metabolic differences between the two varieties ( Supplementary Figure S4A ). Differential analysis detected 202, 84, and 150 DAMs between LYQ and JD at S1, S2, and S3, respectively. Metabolomic analysis revealed a clear distinction between LQY and JD at S1 ( Supplementary Figure S4B ), which was consistent with transcriptome profile. KEGG functional enrichment analysis of the DAMs indicated that metabolic changes at S1 were primarily associated with phenylalanine metabolism, monoterpenoid biosynthesis, flavonoid biosynthesis, tryptophan metabolism, and tyrosine metabolism ( Supplementary Figure S4C ). At S2, glycine, serine, and threonine metabolism, along with sphingolipid and flavonoid metabolism were implicated in the DAMs between varieties ( Supplementary Figure S4D ). At S3, metabolic pathways including phenylalanine metabolism and tyrosine metabolism, were highly conserved between varieties ( Supplementary Figure S4E ). Phenylalanine metabolism was highly active during at S1 and S3, aligning with DEG-enriched metabolic pathways. Metabolomic and transcriptomic evidence indicated that S1 and S3 were crucial for polyphenol production.

We applied k-means clustering to categorize the metabolites into 12 subgroups, aiming to identify similar variation patterns within the metabolomic data. Module 4, Module 5, and Module 6 were identified as distinct metabolic modules, each exhibiting significant differences between the two varieties. In Module 4, metabolites decreased sharply at S2 in LQY, while remaining stable at all stages in JD. These metabolites were found to be associated with the biosynthesis of flavonoids, anthocyanins, and unsaturated fatty acids. Gene-metabolite interaction analysis identified PRL, COMT, and PPARG as key genes involved in these metabolic pathways. In Module 5, metabolites in LQY were upregulated slightly at S2 and S3, whereas in the JD variety, they exhibited notable increases. This module was primarily involved in the metabolism of glycine, serine, and threonine, as well as nicotinate, niacinamide, and thiamine. Genes such as PTGS2, PYGM and NT5E were identified as key genes. In contrast, Module 6 remained stable in LQY but was significantly downregulated in the JD variety at S2. This module was primarily involved in the biosynthesis pathways of flavonoids, flavones, and flavonols. COMT, FOS, and PRL were identified as key genes in Module 6 ( Figure 5 ).

The expression patterns of POD and MYB genes in the fruits ( Figure 6 ), flowers ( Supplementary Figure S5A ), stems ( Supplementary Figure S5B ), and leaves ( Supplementary Figure S5C ) of LQY and JD were analyzed by q-PCR. The results indicated that the expression levels of these two gene families were highest in the fruits and lower in all other tissues. Furthermore, expression of both gene families was generally higher in LQY than in JD, with the exception of Maker00036505 in fruit. Gene expression increased with developmental stages, typically peaking at S2 in LQY, while in JD, the highest expression was observed at S3. These findings were consistent with RNA-seq results and suggest that MYB TFs and POD enzymes may regulate the expression of browning-related genes in LQY and JD, potentially contributing to browning in luffa.