SlBEL11-RNAi transgenic line was kindly donated by Dr. Daqi Fu, School of Food Science and Nutrition Engineering, China Agricultural University. All tomato plants, wild type (Micro-tom) and SlBEL11-RNAi line used in this experiment were cultivated in a growth incubator under photo-cycle condition of 16-h light (22000 Lux) at 25°C and 8-h dark at 20°C and a maintained humidity at 70%~80%. Fresh fruit samples were collected at breaker stage and used for the subsequent transcriptomic and metabolic analyses. Three biological and technical replicates were implemented for both transcriptome and metabolome profiling.

Fruit firmness was measured using a pointer-type fruit firmness tester (Model GY-3, Aipu Measuring Instruments Co., Ltd., China). The test sample was placed face up on a horizontal experimental bench, and the compression force required to break the fruit was recorded. The value was divided by the surface area of the compressed region, and the pressure required per unit area was taken as the firmness metric of the tomato fruit.

Tomato fruit samples were flash-frozen in liquid nitrogen, freeze-dried using a vacuum freeze-dryer (Scientz-100F), and ground into powder with zirconium oxide beads using a mixer mill at 65 Hz for 1 minute. Total RNA was extracted using a RNA extraction kit (Tiangen Biotech, Beijing, China) according to the manufacturer’s instructions. RNA quantity and purity were measured using a Nano Drop ND-1000 (Thermo Fisher), with acceptable thresholds set as A260/A280 = 1.8-2.1 and A260/A230 ≥ 2.0. RNA integrity was evaluated using an Agilent Bioanalyzer 2100, and only samples with RNA Integrity Number (RIN) ≥ 7.0 were selected for downstream analysis. cDNA libraries were constructed using the Illumina TruSeq Stranded mRNA Library Prep Kit, including mRNA enrichment, fragmentation, double-stranded cDNA synthesis, end repair, adapter ligation, and PCR amplification. After quality validation, libraries were sequenced on an Illumina NovaSeq 6000 system (LC-Bio, Hangzhou, China) in paired-end (PE150) mode, generating ≥6 GB of raw data per sample.

Raw sequencing reads were preprocessed using Fastp to remove low-quality reads (Q < 20), adapter-contaminated sequences, and reads with >5% ambiguous bases (N). Paired-end reads were aligned to the tomato reference genome (SL4.0, downloaded from Sol Genomics Network) using HISAT2 v2.2.1 with parameters: –rna-strandness RF –dta –phred33. Index files were generated using hisat2-build with default settings. Gene expression levels were quantified as Fragments Per Kilobase of transcript per Million mapped reads (FPKM), a widely used metric for estimating transcript abundance. Differential expression analysis was performed using DESeq2 (v1.38.3), with significance thresholds set as |log₂(fold change)| ≥1 and Benjamini-Hochberg adjusted P-value (FDR) < 0.05. Functional enrichment analysis included KEGG pathway analysis via hypergeometric testing (FDR < 0.05) and Gene Ontology (GO) term analysis using Fisher’s exact test, both referenced against the tomato genome annotation database.

The pretreatment process for tomato samples was consistent with transcriptomics protocols. A 50 mg aliquot of the powdered sample was mixed with 1 mL of pre-chilled extraction solvent (methanol/water/formic acid, 15:4:1, v/v/v), vortexed, and sonicated in an ice bath (20 kHz, 5-second intervals, total duration 1 hour). The mixture was centrifuged at 8,000 × g for 5 minutes at 4°C, and the supernatant was collected, vacuum-dried, and reconstituted in 80% methanol. After purification via centrifugation (20,000 × g, 20 minutes, 4°C), the solution was filtered through a 0.22 μm cellulose acetate membrane and stored in HPLC vials at -80°C. Three biological replicates were included per group, with quality control (QC) samples prepared by pooling equal amounts of WT and SlBEL11-RNAi extracts. Three consecutive injections of QC samples were performed prior to formal analysis to stabilize the instrument. Chromatographic separation was carried out on an Agilent SB-C18 column (1.8 μm × 2.1 mm × 100 mm) using a UPLC system (ExionLC™ AD) coupled with a 6500 QTRAP mass spectrometer. The mobile phases consisted of 0.1% formic acid in water (A) and 0.1% formic acid in acetonitrile (B), with a gradient program: 95% A to 95% B over 9 minutes, held for 1 minute, then returned to initial conditions in 70 seconds (flow rate: 0.35 mL/min; column temperature: 40°C). Mass spectrometry parameters included electrospray ionization (ESI) in positive/negative switching mode, ion source temperature of 550°C, and spray voltages of ±5,500/4,500 V.

Raw data were processed using MS-DIAL for peak alignment, retention time correction, and peak area extraction. Metabolites were identified by matching accurate mass (mass tolerance < 0.01 Da) and MS/MS spectra (mass tolerance < 0.02 Da) against in-house standards, the Human Metabolome Database (HMDB), and MassBank. Features detected in > 50% non-zero measurements within at least one experimental group were retained for downstream analysis. Differential metabolites were identified through a dual-filter approach combining Orthogonal Partial Least Squares-Discriminant Analysis (OPLS-DA) parameters and statistical criteria: (1) variable importance in projection (VIP) scores > 1,(2) absolute fold-change (FC) ≥ 2 with p< 0.05.

Total RNA was extracted from tomato tissues using the RNApure Plant Kit (CWBIO, Beijing, China). For first-strand cDNA synthesis, 2 μg of total RNA was reverse-transcribed using reverse transcriptase and oligo(dT) primers. Quantitative PCR (qPCR) was performed on a qTOWER3/G real-time system (Analytik Jena, Germany). Each reaction (20 μL total volume) contained 25 ng cDNA, 200 nM of each primer, and 4 μL SuperReal PreMix Plus (Tiangen Biotech, Beijing, China; containing DNA polymerase, dNTPs, and optimized buffer components). The thermal cycling program included an initial denaturation at 95°C for 30 s, followed by 40 cycles of 95°C for 5 s (denaturation) and 59°C for 30 s (annealing/extension). Melt curve analysis was performed to verify amplification specificity. Gene expression levels were normalized to the tomato actin gene as an internal control. the 2^-ΔΔCt was rigorously applied for relative quantification of gene expression (Livak and Schmittgen, 2001). The primer sequences used in this study are provided in Supplementary Table S1 .

Data are presented as mean ± standard deviation (SD). Multivariate data analysis and graphical visualization were performed using R (version 4.0.3) and associated R packages.

Observations of developing fruits in wild-type and SlBEL11-RNAi lines revealed that silencing SlBEL11 expression significantly enhanced chlorophyll accumulation in immature fruits (a phenotype previously reported by Meng et al., 2018). No obvious signs of fruit softening were detected during the growth phase ( Figure 1A ). However, fruits began to abscise progressively upon entering the ripening stages (Dong et al., 2024), with noticeable softening observed via tactile evaluation. Subsequent analysis confirmed the silencing efficiency of SlBEL11 in transgenic lines, demonstrating a marked reduction in SlBEL11 transcript levels at the breaker stage fruits ( Figure 1B ). Firmness measurements revealed a 30% reduction in SlBEL11-RNAi fruits at breaker stage ( Figure 1C ).

To elucidate the molecular mechanisms, we conducted comparative transcriptome profiling of wild-type and SlBEL11-RNAi fruits using Illumina NovaSeq 6000 sequencing. As shown in Supplementary Table S2 , a total of 13.43 GB of raw data (267,776,280 paired-ended reads) were generated. Stringent quality control using Fastp v0.23.4 was conducted to remove low-quality reads, adapter sequences and reads containing > 5% ambiguous bases (N), yielding 39.72 GB of high-quality clean data (263,556,950 valid reads) with Q30 > 95.97%, and GC content of 42%-45%.

The biological repeatability of the samples was evaluated using Pearson correlation coefficient ( Supplementary Figure S1A ). Intra-group sample correlations exceed R² > 0.9, revealing the reliability and reproducibility of the experimental design. Gene expression levels were normalized using the FPKM method and visualized via violin plots ( Supplementary Figure S1B ) and density distribution map ( Supplementary Figure S1C ). These analyses revealed similar gene expression patterns between groups, with log₁₀(FPKM) values concentrated in the range of -2 to 2, indicating that SlBEL11 silencing did not induce global transcriptional alterations.

Differentially expressed genes (DEGs) were further detected using DESeq2 v1.38.3 with a threshold of |log₂ Fold Change| > 1 and FDR-corrected P < 0.05 ( Figure 2A , Supplementary Table S3 ). Only 665 DEGs were identified, including 417 up-regulated and 248 down-regulated genes. Hierarchical clustering heatmap ( Figure 2B ) revealed distinct intergroup segregation and tight intragroup clustering of DEGs. To verify the transcriptomic results, 14 DEGs were selected for RT-qPCR analysis ( Supplementary Figure S2 ). The expression patterns of the tested DEGs were consistent with that in the transcriptome.

GO enrichment analysis of DEGs are shown in Figure 2C and Supplementary Table S4 . In the category of Biological Process, DEGs are significantly enriched in the pathways of single-organism process like single-organism metabolic process, single-organism localization, single-organism transport, suggesting that SlBEL11 regulates basal physiological functions. The enrichment of DEGs in other processes, such as oxidative-reduction process, localization and transport related processes, are also detected. In the category of molecular function, the significant enrichment of oxidoreductase activity, cofactor binding, and coenzyme binding, further supported the alteration of oxidative-reduction process in SlBEL11-RNAi tomatoes. The detection of binding and transport activities, such as tetrapyrrole binding, heme binding, fructose 1,6-bisphosphate 1-phosphatase activity, benzoate and xenobiotic transporters, hinted at potential changes in secondary metabolism. In the category of cellular component, however, only the “photosystem II oxygen evolving complex” was significantly enriched, indicating a potential impact on chloroplast function. This finding aligns with the result of KEGG enrichment analysis ( Figure 2D , Supplementary Table S5 ) where DEGs clustered in photosynthesis-antenna protein pathways. Additionally, enrichment in linoleic acid metabolism, brassinosteroid biosynthesis and ABC transporter were also detected.

As the chromatography system/mass spectrometer is in direct contact with the samples, the accumulation of residues in the chromatographic column and the mass spectrometry ion source may cause signal drift or system errors with increasing sample load (Hao et al., 2023). To ensure data reliability and repeatability, three quality control (QC) samples were used for continuous monitoring of the instrument in this study. The superimposed analysis of total ion chromatograms in both positive and negative ion modes showed that the peak intensities and time reproducibility of the QC samples were highly consistent ( Supplementary Figure S3 ), demonstrating excellent signal stability of the instrument. Further pearson correlation analysis of the QC samples showed that the correlation coefficients were greater than 0.9 ( Supplementary Figure S4 ), confirming the stability of the experimental procedure and the optimal performance of the instrument.

Metabolites were structurally identified by matching retention time, molecular mass (mass error <10 ppm), MS/MS fragmentation patterns, and collision energy against both in-house and public databases. All identifications were subjected to rigorous manual verification. Metabolites with a coefficient of variation (CV) <30% in QC samples were retained for subsequent analysis. A total of 714 metabolites were identified in wild-type (WT) and SlBEL11-RNAi tomato samples, spanning 22 metabolic categories, including alcohols(16), alkaloids(41), amino acid and derivatives(92), anthocyanins(12), carbohydrates(20), flavanone(21), flavone(51), flavonoid(18), flavonol(29), indole derivatives(6), isoflavone(5), lipids(75), nucleotide and derivates(59), organic acids and derivatives(106), phenolamides(27), phenylpropanoids(62), polyphenol(7), proanthocyanidins(1), quinones(2), sterides(5), Terpene(13), Vitamins and derivatives(16) and unclassified compounds(30) ( Supplementary Table S6 ).

Multivariate analyses of 714 metabolites revealed distinct metabolic profiles between WT and SlBEL11-RNAi tomato lines. Principal component analysis (PCA) separated the two groups along the primary axis (PC1, 66.93% variance), with WT and SlBEL11-RNAi samples clustering negatively and positively, respectively ( Figure 3A ). While PCA validated experimental stability and intergroup variability, its unsupervised nature limited sensitivity to subtle biological differences. To address this, supervised orthogonal partial least squares-discriminant analysis (OPLS-DA) was employed, yielding an enhanced group discrimination ( Figure 3B ). The model exhibited high reliability (permutation test: R²Y > 0.5, Q² >0.5) with no overfitting ( Figure 3C ), confirming robust metabolic distinctions between genotypes.

A three-tiered screening strategy (absolute FC > 2, P < 0.05, OPLS-DA-derived VIP > 1) was implemented to identify metabolically significant features. A total of 189 differential metabolites were identified in the WT and SlBEL11-RNAi tomato samples. As shown in Figure 4A , compared with WT, 115 metabolites were up-regulated and 74 were down-regulated in SlBEL11-RNAi tomatoes compared to WT. These differential metabolites include 26 lipids, 25 organic acids and derivatives, 20 phenylpropanoids, 18 amino acids and derivatives, 16 phenolic amines, 16 flavonoids, 11 nucleotides and derivatives, 10 flavonols, 10 alkaloids, 7 flavones, 7 flavanones, 4 terpenoids, 3 alcohols, 3 vitamins and derivatives, 2 polyphenols, 2 anthocyanins, 2 isoflavones, 2 indoles and derivatives, 1 carbohydrate, 1 proanthocyanidin and 3 unclassified compounds ( Figure 4B ).

A clustering heatmap was generated to visualize sample relationships and the differences of metabolite intensity, based on the normalized expression values of differential metabolites. As shown in Supplementary Figure S5 , a distinct hierarchical clustering of metabolite among groups was observed. The top 10 up-regulated and down-regulated differential metabolites were selected using fold change as a criterion. As shown in Figure 4C , the top 10 up-regulated differential metabolites included fumaric acid, N-caffeoyl spermidine, geniposide, syringic acid, tricin O-hexosyl-O-syringin alcohol, N-sinapoyl cadaverine, O-p-coumaroyl quinic acid O-rutinoside derivative, 3-O-p-coumaroyl shikimic acid, cinnamoyl tyramine, phosphatidylcholine acyl 19:2/16:0. The top 10 down-regulated differential metabolites were O-feruloyl coumarin, D-erythro-sphinganine, coumarin O-rutinoside, tricin 5-O-hexoside, 3-(4-hydroxyphenyl)propionic acid, eriodictiol C-hexosyl-O-hexoside N-acetyl-L-tyrosine, sakuranetin, hesperetin O-hexosyl-O-hexoside, N-p-coumaroyl hydroxyagmatine.

KEGG pathway enrichment analysis of the differentially expressed metabolites was performed using Metaboanalyst 4.0. The top 20 significantly enriched metabolic pathways are presented in Figure 4D , including flavonoid biosynthesis, phenylpropanoid biosynthesis, biosynthesis of phenylpropanoids, ubiguinone and other terpenoid-guinone biosynthesis, longevity regulating pathway, toluene degradation, dopaminergic synapse, stilbenoid, diarylheptanoid and gingerol biosynthesis, asthma, betalain biosynthesis, biosynthesis of enediyne antibiotics, biosynthesis of vancomycin group antibiotics, bisphenol degradation, fc epsilon RI signaling pathway, folate biosynthesis, histamine H2/H3 receptor agonists/antagonists, monoterpenoid biosynthesis, phosphatidylinositol signaling system, PI3K-akt signaling pathway, puromycin biosynthesis.

KEGG enrichment analysis of differential genes and metabolites identified 25 co-enriched KEGG-enriched pathways ( Figure 5A ). To further explore the relationship between DEMs and DEGs and determine the pathways affected by SlBEL11, we overlaid p-values thresholds on KEGG histograms, prioritizing pathways enriched by both DEMs (p < 0.05) and DEGs (p < 0.01) ( Figure 5B ). This approach identified six key pathways, including ABC transporters, ascorbate and aldarate metabolism, glycine/serine/threonine metabolism, glycolysis/gluconeogenesis, phenylpropanoid biosynthesis, and pyruvate metabolism.

Expression and regulatory patterns of differential metabolites and genes associated with glycolysis/gluconeogenesis, ascorbate/aldarate metabolism, and phenylpropanoid biosynthesis are summarized in Figure 6 . In glycolysis/gluconeogenesis, salicin decreased twofold, accompanied by downregulation of ADH1 (3.3-fold) and PK (2.4-fold) ( Figure 6A , Supplementary Table S7 ). For ascorbate/aldarate metabolism, inositol declined 2.7-fold, while APX (24-fold), ALDH (5.7-fold), and GME (2.1-fold) were upregulated, contrasting with the marked suppression of AO (5.7-fold) ( Figure 6B , Supplementary Table S7 ). In phenylpropanoid biosynthesis, seven metabolites, including coniferyl alcohol (5.1-fold), sinapyl alcohol (4.7-fold), L-tyrosine (4.6-fold), Scopoletin (2.9-fold), caffeate (2.6-fold), coniferyl aldehyde (2.6-fold) and cinnamic acid (2.2-fold), showed elevated abundance, whereas syringin declined 5.9-fold. Concurrently, UGT72E (69.2-fold), CCR (11.3-fold) and E1.11.1.7 (2.3-fold) were upregulated, opposing the 1.9-fold downregulation of PAL ( Figure 6C , Supplementary Table S7 ). Figure 7 illustrates coordinated metabolic and transcriptional interactions across ABC transporters, pyruvate metabolism, and glycine/serine/threonine metabolism. In the category of ABC transporters, ornithine (5.3-fold) and biotin (4.3-fold) accumulated, while inositol decreased 2.7-fold alongside the upregulation of ABCB14 and ABCB3 (8-fold and 7.6-fold, respectively) ( Figure 7A , Supplementary Table S7 ). Pyruvate metabolism exhibited fumaric acid accumulation with four upregulated genes, including ALDH (5.7-fold), maeB (5.4-fold), DLD (4.8-fold) and chMDH (2.2-fold), contrasting with the suppression of PK (2.4-fold) ( Figure 7B , Supplementary Table S7 ). Glycine/serine/threonine metabolism featured elevated L-tryptophan (4.0-fold) and phosphoserine (2.2-fold), concurrent with upregulation of gcvH(5.4-fold), DLD (4.8-fold) and AGXT (2.9-fold), opposing the downregulation of glyA (2.2-fold). ( Figure 7C , Supplementary Table S7 ).