Transcriptomic Analysis of Early Fruit Development in Micro-Tom Tomato Reveals Conserved and Cultivar-Specific Mechanisms

1. Introduction

Tomato (Solanum lycopersicum) is an important model species for studying fleshy fruit development due to its extensive phenotypic diversity, particularly in fruit weight and tissue morphology [1]. This diversity provides valuable insights into the fundamental biological processes of fruit growth, ripening, and hormonal regulation. Fruit development in tomato follows a biphasic pattern comprising an initial phase of active cell division post-anthesis, followed by a phase of cell expansion without clear boundaries between these stages [2,3]. Hormones, notably gibberellins (GAs) and auxins, play critical roles in controlling fruit set and development. Gibberellins promote cell division, while auxins primarily promote cell expansion, influencing overall fruit size and morphogenesis [4].

The early developmental stages are critical for fruit formation and final fruit size. Shortly after anthesis, fertilization determines whether fruit development will continue or abort, a process known as fruit set [5,6]. The size and weight of the fruit largely depend on the number of cells in the pericarp, established during the cell division stage [7]. Following fertilization, the embryo within the ovule undergoes rapid cell division, reaching a multicellular embryonic stage around 3 to 6 days post-anthesis (DPA) [8]. Simultaneously, cells in the ovary wall begin dividing around 2 DPA, with cell division in the ovary wall nearly completed by 5 DPA, marking the onset of the cell expansion phase [8,9,10].

Several genes have been identified that affect fruit size by regulating cell division during these early developmental stages. Remarkably, despite the wide diversity in fruit size and shape observed across cultivated tomatoes, a limited set of genes seems to drive the major variations in fruit morphology among domesticated varieties [11,12]. Alleles of the major genes SUN, OVATE, LOCULE NUMBER, and FASCIATED individually explain up to 71% of the fruit shape variations in cultivated tomato (Rodríguez et al. 2011). Quantitative trait loci-associated genes Cell Number Regulator and SlKLUH (from QTLs fw2.2 and fw3.2, respectively) also influence fruit size through modulation of cell division rates [13].

The dwarf tomato cultivar Micro-Tom has emerged as a popular model system for genetic and developmental studies, due to its compact size (~15 cm in height), short life cycle (~70 days from germination to fruiting), and adaptability to high-density cultivation [14,15]. Originally bred as an ornamental variety, Micro-Tom carries distinctive mutations, including the dwarf (d) and self-pruning (sp) alleles, that underlie its reduced size and determinate growth habit [15]. These traits, along with the uniform ripening (u) gene, make it an ideal controlled experimental model, particularly for studying plant hormone interactions and responses to environmental stimuli [16]. The uniform ripening phenotype is especially advantageous for fruit development research, as it enables synchronized sampling of large quantities of fruit tissue for downstream analyses. Furthermore, Micro-Tom possesses a relatively compact diploid genome of approximately 900 Mb and is supported by extensive genetic and genomic resources, which have collectively facilitated major advances in elucidating the molecular and genetic regulation of fruit formation [17,18].

Transcriptomic studies in tomato, particularly those utilizing the Micro-Tom model, have provided valuable insights into gene expression and regulatory networks throughout different developmental stages. Tissue-specific transcriptome profiles have been mapped across essential reproductive organs, shedding light on the dynamics of gene regulation influencing fruit formation and maturation [19,20]. These studies have highlighted the complexity of gene regulatory networks driving the transition from ovary development to fruit maturation, including the roles of specific loci such as SUN and OVATE in determining fruit shape and the processes leading cellular expansion and endoreduplication [1].

More recently, cutting-edge approaches such as spatial transcriptomics [5,13] have begun to uncover the cellular heterogeneity and spatial organization of gene expression during fruit development. However, most studies focus on fruit ripening [2][23] and stress responses [24], with fewer transcriptomic analyses of early fruit development. An RNA-seq-based investigation of early tomato fruit development was performed in the 'Moneymaker' cultivar [13], providing insights into gene expression dynamics during the initial stages of fruit formation. This study identified key regulators involved in pericarp- and ovule-specific development post-fertilization, highlighting the roles of auxin biosynthesis genes and tissue-specific transcription factors. Nonetheless, comprehensive transcriptomic analyses focusing on the early stages of fruit development in the Micro-Tom cultivar remain lacking.

In this study, we used the most recent, high-quality assembly of the Micro-Tom tomato genome [17,18] to investigate transcriptomic changes associated with early fruit development. We selected three critical time points: 3 DPA, representing the cell division stage; 5 DPA, corresponding to the peak of cell proliferation; and 8 DPA, marking the transition to cell expansion. Focusing on predicted protein-coding genes, we identified nine distinct expression clusters and functionally correlated these transcriptional patterns with the biological processes occurring during each developmental stage. These findings were compared to known patterns from other tomato cultivars to generate a similarity map, thereby enhancing our understanding of conserved and cultivar-specific regulatory mechanisms governing early fruit development in Micro-Tom. Additionally, integration of untargeted metabolomic data derived from previously published GC-MS and LC-MS analyses provided complementary insights, revealing a strong concordance between transcript abundance and metabolite profiles during this critical developmental window.

2. Results

2.1. RNA-Seq Quality and Alignment Statistics

RNA-seq statistics indicate high-quality data suitable for downstream analysis. On average, 80.7 ± 4.2% of raw reads passed quality filtering and trimming, resulting in read libraries with a mean Q30 of 94.86 ± 0.17% (Table 1). These libraries were successfully aligned to the Micro-Tom reference genome [18] with an average of 91.8 ± 0.9% of reads mapping uniquely (Table 1). Gene body coverage analysis shows a characteristic bell-shaped distribution, typical of high-quality poly(A)-selected libraries [25] and does not indicate degradation or technical artifacts (Supplementary Figure S1). In addition, Principal Component Analysis (PCA) before and after Trimmed Mean of M-value (TMM) normalization show consistent clustering by developmental stage, with normalization having a minimal impact on overall sample relationships (Supplementary Figure S2). This indicates that the dataset was well balanced and free from major biases related to library size or composition.

2.2. Differential Expression Shows Stage-Specific Transcriptional Changes

To characterize gene expression dynamics during early fruit development, we performed pairwise differential expression analysis between 3 DPA vs. 5 DPA, 5 DPA vs. 8 DPA, and 3 DPA vs. 8 DPA (Supplementary File S1). Applying a threshold of |log2FC| > 1.5 and FDR < 0.01, we identified 368 upregulated and 120 downregulated genes between 3 and 5 DPA, 915 upregulated and 1,611 downregulated genes between 5 and 8 DPA, and 2,158 upregulated and 2,056 downregulated genes across the 3 to 8 DPA interval (Figure 1A). These differentially expressed genes (DEGs) represent 2.25%, 11.63%, and 19.4% of all genes expressed in the dataset, respectively, based on a total of 21,725 genes retained after filtering for low expression.

Venn diagram analyses of significant DEGs show high overlap between the 5 vs. 8 DPA comparison and the 3 vs. 8 DPA set: 702 upregulated and 1,050 downregulated genes are shared (Figure 1B,C). This suggests that the major transcriptional reprogramming occurs between 5 and 8 DPA, and the differences observed in the 3 vs. 8 DPA comparison largely reflect changes during the latter stages.

Gene Ontology (GO) enrichment analysis of Biological Process terms supports this interpretation (Supplementary File S1). Between 3 and 5 DPA (Supplementary Figure S3), upregulated genes are enriched in processes such as response to wounding, fatty acid biosynthesis, and cellulose production, suggesting early activation of membrane lipid metabolism and cell wall synthesis. No strong enriched downregulated terms were detected at this stage, consistent with a modest overall transcriptional shift.

In the 5 to 8 DPA transition (Supplementary Figure S4), upregulated genes were strongly associated with lipid storage, methylation, and pectin catabolism, reflecting increased energy allocation and active cell wall remodeling. Concurrently, downregulated genes were enriched in defense responses, microtubule-based movement, and carbohydrate metabolism, suggesting a reduction in stress-related and structural reorganization pathways as the tissue progresses toward functional specialization.

Across the 3 to 8 DPA interval (Supplementary Figure S5), upregulated genes were predominantly enriched for photosynthesis, chloroplast biogenesis, and fatty acid biosynthesis, highlighting the metabolic shift toward organellar development and energy capture. In contrast, downregulated genes were associated with protein phosphorylation, defense response to fungi, and cell wall modification, indicating a broad transcriptional reprogramming from stress adaptation and structural turnover to sustained biosynthesis and tissue maturation.

2.3. Distinct Expression Clusters Mark Developmental Transitions

To better understand the molecular processes and key factors involved early fruit development, we performed k-means clustering analysis for the three developmental stages (Figure 2, Supplementary File S2). This analysis identified nine distinct expression patterns (clusters), which highlight the dynamic and coordinated regulation of genes as the fruit transitions from cell proliferation to cell expansion. Additionally, GO enrichment analysis provided further insight into the biological processes underlying these expression patterns, which represent diverse trends of gene upregulation, downregulation, or stage-specific peaks in the 3-8 DPA interval (Figure 3, Supplementary File S2).

Cluster 1, consisting of 955 genes, shows high expression levels at 3 and 5 DPA, followed by a notable decrease at 8 DPA. This expression pattern suggests that these genes are strongly involved in early developmental processes, particularly those associated with DNA metabolic activities. GO enrichment analysis revealed that a substantial portion of these genes, including at least 55 proteins involved in DNA replication, play a crucial role in cell division.

Similarly, clusters 6 and 7, containing 1,252 and 2,038 genes, respectively, also exhibit a declining expression trend towards 8 DPA, although the decline is more gradual compared to the sharp drop observed in cluster 1. These clusters are enriched in genes related to RNA metabolic processes regulation. The observed decline in both DNA- and RNA-related processes at 8 DPA reflects a reduced need for high proliferative activity as the fruit shifts into the cell expansion phase. Taken together, the coordinated downregulation in clusters 1, 6, and 7 indicates a developmental transition from proliferation-related activities to processes more aligned with cell expansion and maturation.

A distinct set of genes shows a peak at 5 DPA (cluster 8, 605 genes), with the cell cycle being the top-enriched category, including 50 genes (Figure 4A). This peak coincides with the culmination of cell proliferation, suggesting that these genes play a key role in orchestrating the final rounds of cell division prior to the fruit's transition into the expansion phase. The predicted protein-protein interaction network of genes in this cluster, based on conserved interactions in the STRING database (Figure 4B) shows their potential roles in DNA replication, particularly within the S phase of the cell cycle, as well as mitosis, including members of the Anaphase Promoting Complex/Cyclosome (APC/C).

Clusters 4, 5, and 9 represent the largest expression groups, all exhibiting an upward trend in gene expression toward the 8 DPA timepoint (Figure 2). Cluster 5 (1,349 genes) display a more pronounced increase between 5 and 8 DPA and is enriched for protein transport functions, with 105 genes associated with this GO term. In contrast, cluster 4 (1,344 genes) shows a more gradual transition and is enriched in photosynthesis-related processes (76 genes) as well as broader cellular metabolic functions (667 genes). Cluster 9, the largest cluster with 2,332 genes, is enriched for small molecule metabolic processes. The late activation of genes in these clusters likely reflects a ramp-up in metabolic activity required to support fruit growth, highlighting the importance of energy production and metabolic support in ensuring proper maturation, fruit quality, and ripening.

• Figure 4. Heatmap and predicted interaction network of cell cycle–related genes from cluster 8. A- Heatmap showing the expression profiles of 50 cell cycle related genes from cluster 8 across the three developmental stages. Gene symbols are provided in parentheses for annotated genes The heatmap was generated using the Morpheus tool (https://software.broadinstitute.org/morpheus). B- Predicted protein–protein interaction network based on conserved interactions from homologs in other organisms, as retrieved from the STRING database. Line thickness represents the strength of evidence supporting each interaction.

2.4. Different Transcription Factors May Regulate Early Fruit Development

To investigate the transcription factors (TFs) potentially regulating early fruit development, we analyzed the putative promoter regions (2 kb upstream of coding sequences) of co-expressed genes within each cluster for enrichment of TF binding motifs (Supplementary File S3). Interestingly, only one TF with significantly enriched motifs (Bonferroni-Hochberg adjusted p-value < 0.05) was the MYB family member SLM2ch04g16973, whose binding motifs were found in the upstream regions of 13 genes in Cluster 6 (Figure 5). This TF corresponds to the gene Solyc04g077260 in the ITAG4.0 annotation of the Heinz tomato genome, and its function remains uncharacterized. The predicted targets of SLM2ch04g16973 include genes with diverse functions, such as protein kinases, RNA-binding proteins, proteasome-related proteins, and a histone acetyltransferase. This TF was filtered out in the DE analysis due to its higher adjusted p-value (0.023) in the 3 vs. 5 DPA comparison, but it is slightly more expressed at 3 DPA than at both 5 and 8 DPA (Supplementary Figure S6). The predominant expression of its putative target genes at 3 DPA further supports a potential role for SLM2ch04g16973 in early fruit development, deserving further investigation.

To gain insights into the transcriptional regulatory landscape during early fruit development in Micro-Tom tomato, we focused on the expression patterns of TFs across developmental stages. Of the 1,623 predicted TFs in the Micro-Tom genome, 586 (36.1%) were detected in our filtered dataset, and these were distributed among the nine expression clusters (Supplementary Figure S7, Supplementary File S3). Among the detected TFs, several families such as bHLH, MYB, NAC, WRKY, ERF, C2H2, and HD-ZIP were represented by many members (Supplementary File S3), reflecting their broad regulatory functions. For instance, the Auxin Response Factor (ARF) and C3H zinc finger families showed strong representation in Cluster 1, which is characterized by higher expression at 3 DPA (Supplementary Figure S7). This expression pattern may indicate a regulatory role early in fruit development, potentially contributing to the onset of fruit-specific transcriptional programs.

Although many TF families were broadly represented across clusters, a hypergeometric enrichment analysis revealed that only members of the MIKC-MADS family were significantly overrepresented in a specific cluster, Cluster 9, which comprises genes with increasing expression over time (Supplementary Figure S7). This result suggests that MIKC-MADS TFs may play a particularly important role during the later stages of early fruit development, possibly regulating key transcriptional switches involved in floral organ identity and fruit maturation.

2.5. Integration of Metabolite Clusters with Gene Expression

To better understand the coordination between metabolic and transcriptional programs during early fruit development, we integrated our gene expression clusters with the metabolite profiles reported by [26]. In that study, 72 metabolites identified via LC-MS and 46 via GC-MS were quantified at 3, 5, and 8 DPA in the same biological system. By leveraging this dataset, we aimed to identify temporal relationships between metabolite accumulation and gene expression, which could point to co-regulated pathways or feedback mechanisms linking metabolism to developmental gene regulation. Using hierarchical clustering, we organized these metabolites into seven distinct clusters based on their temporal patterns (Figure 6, Supplementary File S4). These metabolite profiles closely mirrored gene expression dynamics during tomato fruit development. By correlating the metabolite clusters with their corresponding gene expression clusters, we uncovered distinct temporal patterns that further highlight key developmental transitions and suggest potential feedback regulation between metabolism and gene expression.

Cluster M1 comprises compounds such as serine, asparagine, several flavonoid derivatives, tryptophan, and succinic acid. These metabolites correspond to gene cluster 5, which follows a “mid/low–low–high” expression pattern, peaking at 8 DPA. GO enrichment analysis of this gene cluster highlights functions related to protein localization and vesicle-mediated transport, suggesting a temporal association between metabolite accumulation and gene expression.

In cluster M2, primary metabolites, including quinic acid, GABA, several amino acids, and energy-related compounds such as pyruvate, 2-oxoglutarate, fumarate, and glyoxylate, are found. These metabolites align with gene cluster 9, which exhibits a “low–mid/low–high” expression pattern. Genes in this cluster are enriched for biosynthetic and energy metabolism processes, including translation and small molecule synthesis, which coincide with the later accumulation of these metabolites.

Cluster M3 includes hydroxybenzoic acid-hexose, caffeoyl-hexoside 5, malic acid, and various glycosides such as feruloyl-glucoside 3. These metabolites correspond to gene cluster 2, which follows a “high–low–high” expression pattern. GO enrichment of genes in this cluster identifies functions related to redox activity and secondary metabolism, consistent with the observed fluctuations in metabolite levels.

Metabolite cluster M4 is composed of caffeoyl derivatives, quinic acid isomers, flavonoids, and coumaric acid derivatives. These metabolites are associated with gene cluster 6, which exhibits a “high–low–low” expression pattern, with strong early expression that declines at later developmental stages. Genes in this cluster are enriched in regulatory processes, indicating a potential relationship between early gene expression and metabolite accumulation.

Cluster M5, containing beta-alanine, caffeoyl-glucarate, fructose, and feruloyl-glucoside 1, corresponds to gene cluster 4, which follows a “low–mid–high” expression pattern. GO enrichment analysis highlights functions related to photosynthesis and chloroplast organization, suggesting that gene expression and metabolite accumulation are temporally aligned during fruit maturation.

Metabolite cluster M6 includes Que-hexose, quinic acid 3-caffeoyl trans, and quercitrin, which are associated with gene cluster 1. This cluster follows a “high–high–low” expression profile, with early gene expression linked to DNA metabolism and cell proliferation. The corresponding metabolites show a similar early peak, followed by a decline, suggesting synchronization between metabolite levels and cell division activity.

Finally, in cluster M7, which includes salicylic acid and sucrose, corresponds to gene cluster 8. This gene cluster follows a “low–high–low” pattern with a transient peak at 5 DPA. GO enrichment identifies genes related to cell cycle regulation, coinciding with the temporary accumulation of these metabolites.

Across all clusters, hierarchical clustering of metabolite profiles and gene expression profiles reveals strong temporal coordination, were early-stage metabolites align with clusters enriched in regulatory and proliferative functions, and late-stage metabolites associated with genes driving biosynthetic and metabolic processes. This coordinated behavior reinforces the functional integration between transcriptional programs and metabolic transitions during early tomato fruit development.

3. Discussion

Our comprehensive transcriptomic analysis of Micro-Tom tomato revealed a dynamic and tightly regulated gene expression landscape during early fruit development. The nine expression clusters reflect a clear developmental progression from cell division to cellular expansion, with early stages enriched for DNA and RNA metabolic processes, consistent with active cell proliferation. In contrast, later stages showed upregulation of genes involved in biosynthesis, energy metabolism, and cell wall modification, marking the fruit’s transition toward maturation.

Comparisons with transcriptomic data from the “Moneymaker” cultivar [13] reveal both conserved developmental pathways and cultivar-specific regulatory features. In both cultivars, hormone signaling genes, particularly those involved in auxin and gibberellin responses, are key regulators during early fruit development. Similarly, genes associated with cell cycle progression and DNA replication are upregulated in early stages in both datasets, highlighting a shared transcriptional program supporting fruit initiation. However, carbohydrate metabolism genes peaked earlier in Moneymaker (2 to 5 DPA), while their expression was delayed in Micro-Tom (8 DPA), suggesting that the timing of metabolic transitions may differ between cultivars and potentially influence fruit growth dynamics and physiological traits.

A comparison of transcription factors expression supports this pattern of both overlap and divergence. Both studies identified MADS-box and AP2 family genes as central players in early development. In our dataset, MIKC-MADS TFs showed later upregulation, consistent with the findings of [13], who reported M-type MADS expression in ovules at 5 DPA. In contrast, TF families such as TALE, which were prominent in their pericarp-specific analysis, were less represented in our whole-fruit dataset, likely reflecting differences in tissue resolution or genetic background. These findings suggest a conserved regulatory core modulated by cultivar-specific transcriptional programs.

Several insights also emerged from the integration of transcriptomic and metabolomic profiles. Gene expression Cluster 8, which peaks at 5 DPA, aligns with the accumulation of sucrose and salicylic acid (SA). Sucrose is known to act as a metabolic checkpoint, promoting the G1/S transition during cell cycle, and in Arabidopsis, can rapidly induce D-type cyclins in sugar-starved cells [27]. This supports the idea that sucrose availability at this stage may act as a developmental cue for proliferation in tomato fruit.

Salicylic acid may also play a regulatory role at this stage. While high SA levels are associated with stress responses and growth inhibition, moderate concentrations can promote cell division, particularly through crosstalk with auxin and ROS [28]. Given its peak at 5 DPA, SA may act alongside sugars to fine-tune the balance between proliferation and differentiation.

Most amino acids accumulated at 8 DPA, coinciding with the onset of cellular expansion and increased metabolic demand. Flavor-associated amino acids such as glutamate and aspartate, though abundant in later ripening stages [29], also showed early fluctuations, suggesting that amino acid metabolism is already active during the initial phases of fruit development - even if its sensory impact becomes more pronounced later.

Metabolite cluster M2, which includes key carbon intermediates like 2-oxoglutarate, fumarate, and pyruvate, correlated with gene cluster 9, enriched in genes linked to carbon metabolism. This association points to increased mitochondrial activity and metabolic flux, supporting both biosynthesis and energy production during cell wall development and storage compound accumulation [30].

Although our analysis identified several candidate regulators (e.g., MYB family member SLM2ch04g16973), experimental validation is necessary to confirm their roles in development. Approaches such as ChIP-seq or CRISPR-Cas9-mediated knockouts could help determine their downstream targets and regulatory impact during fruit development.

Finally, future studies incorporating spatial transcriptomics and metabolomics could resolve tissue-specific gene-metabolite interactions that are obscured in whole-organ profiling. Such data may clarify whether the observed differences in hormone signaling, sugar and amino acid metabolism are adaptive responses, tissue dilution effects, or inherent to the cultivar. A deeper understanding of these regulatory mechanisms could inform targeted breeding strategies to improve fruit quality, size, and stress resilience in tomato [23].

4. Materials and Methods

4.1. RNA Sequencing

Sequencing libraries were obtained from our previous work [26]. Briefly, stranded total RNA libraries were prepared from samples of developing Micro-Tom fruits collected at 3, 5, and 8 DPA, using three biological replicates per stage, each consisting of a pool of 12 fruits. Sequencing was performed by Fasteris Co., Ltd. (Switzerland) using the Illumina NovaSeq 6000 platform, generating 150 base pairs (bp) paired-end reads. Raw reads were processed with TrimGalore (version 0.6.10) (github.com/FelixKrueger/TrimGalore), using default parameters for adapter trimming and quality filtering. Quality assessment of the processed reads was conducted using FastQC (version 0.12.1) (bioinformatics.babraham.ac.uk/projects/fastqc) and RSeQC (version 5.0.1) [31] before proceeding to downstream analyses.

4.2. Read Mapping and Transcript Quantification

Processed reads were mapped to the most recent assembly of the Micro-Tom genome [18] (NCBI accession GCF_036512215.1) using the STAR aligner (version 2.7.11b) [32], with default parameters optimized for paired-end reads. Gene quantification was achieved by enabling the “--quantMode GeneCounts” option during mapping, which generates read counts per gene directly from the aligner. The count table was then filtered to eliminate low-expressed genes by eliminating genes with less than 10 total reads across at least 3 samples.

4.3. Differential Expression Analysis

Differential expression analysis was conducted using the RNFuzzyApp package (gitlab.com/habermann_lab/rna-seq-analysis-app, cloned in June 5th, 2024) in R (version 4.3.3). Raw gene counts were normalized using the Trimmed Mean of M-values (TMM) method [33]. The resulting expression matrix was analyzed using the TCC method [34], with statistical testing performed by DESeq2 [35]. A False Discovery Rate (FDR) cutoff of 0.01 was applied to identify significantly differentially expressed genes, with |log2FC| > 1.5.

4.4. Expression Clustering

Clustering of expression patterns was performed using a separately filtered set optimized for clustering. Starting from the TMM-normalized counts, genes were retained if they had a minimum of 5 read counts in every sample, resulting in 11,035 consistently expressed genes. The expression values were transformed using a log2(x + 1) function to stabilize variance and then standardized using Z-score normalization. The average expression value for each gene at each time point was calculated for clustering purposes. K-means clustering was executed in R using the kmeans function, specifying nine clusters (k = 9) and setting the maximum number of iterations to 50. A random seed of 123 was used to ensure reproducibility of the clustering results. The k = 9 value was selected based on elbow criterion and cluster stability across 100 random starts (testing k = 2–15).

4.5. Functional Enrichment Analysis

Functional enrichment analysis was conducted using the STRING database (version 12.0) [36] and the clusterProfiler R package [37]. The proteome encoded by the Micro-Tom genome was annotated within STRING to facilitate the analysis and is accessible at https://version-12-0.string-db.org/organism/STRG0A81OUP. Gene IDs from each cluster were uploaded to the STRING workspace to create gene sets. For the DE genes, the enricher function of clusterProfiler was used with the STRING annotation. The total number of expressed genes from the normalized, filtered expression matrices served as the background for statistical enrichment analysis. Overrepresented Gene Ontology (GO) terms and pathways were identified for each cluster using built-in enrichment functions, applying statistical significance cutoffs of p < 0.05 and Benjamini–Hochberg FDR < 0.05 in all cases.

4.6. Transcription Factor Prediction and Analysis

Transcription factors (TFs) were identified through BLASTP similarity searches against known TFs in the PlantRegMap database [38] and via de novo prediction using the planttfhunter R package (version 1.6.0) [39]. Expression values of the identified TFs were transformed using a robust Z-score (subtract median and divide by median absolute deviation). Hierarchical clustering was performed using one minus Pearson correlation as the distance metric, enabling the assessment of expression pattern similarities among TFs during fruit development stages.

4.7. Transcription Factor Motif Mapping and Enrichment Analysis

This analysis was performed similarly to what is described by [40]. Known TF binding motifs for S. lycopersicum were obtained from CisBP v2.00 [41] and JASPAR 2022 [42]. These motifs were mapped to the 2-kb upstream regions of all annotated genes in the genome using FIMO v5.5.7 [43] with default parameters. We retained the top 7,000 scoring motif matches, discarding lower-scoring hits. For each TF motif, hypergeometric tests were used to assess whether the number of motif instances found in genes of a given cluster (Clusters 1–9) was significantly higher than expected based on its total number of occurrences across all expressed genes. Benjamini-Hochberg correction was applied to control the false discovery rate (FDR < 0.05). Representative sequence logos showing nucleotide frequencies at each position were generated using WebLogo (http://weblogo.berkeley.edu/logo.cgi).

4.8. Integration of Metabolomics Data

Normalized metabolite intensities at 3, 5, and 8 DPA [26] were scaled to Z-scores. We performed k-means clustering on the Z-scored matrix, testing k = 2–10 and selecting k = 7 based on the elbow criterion and cluster stability across 100 random starts. To link metabolite clusters (1–7) to gene clusters (1–9), we computed the Spearman correlation between the centroid profiles (mean Z-score at each timepoint) of each metabolite cluster and each gene cluster. Clusters were associated by high Spearman correlation coefficients (>0.8) supported by visual pattern confirmation.

5. Conclusions

In summary, our integrative analysis of transcriptomic and metabolomic data in Micro-Tom tomato revealed distinct gene expression clusters associated with key developmental transitions during early fruit development. These clusters highlight a tightly regulated progression from cell division to expansion, driven by coordinated changes in hormone signaling, transcription factor activity, and metabolic pathways. Our findings also uncover both conserved regulatory mechanisms shared with other tomato cultivars and cultivar-specific differences in transcriptional and metabolic programs.

By linking metabolite accumulation to gene expression dynamics, this study provides new insights into the molecular mechanisms regulating fruit formation. The observed associations between sugars, amino acids, and cell cycle regulators suggest potential metabolic checkpoints that influence early fruit growth. These findings lay the groundwork for future functional studies, such as mutant validation (CRISPR), chromatin accessibility assays (ChIP-seq), and spatial transcriptomics, to further dissect critical regulatory networks.

These insights can directly inform breeding programs aimed at enhancing fruit size, nutritional quality, and resilience to environmental stress, as well as guide biotechnological strategies using precision breeding techniques such as CRISPR/Cas.