Multi-Omics Reveals Domestication-Associated Shifts in Phosphorus Adaptation Strategies in Tomato

TL;DR

Wild tomatoes handle low phosphorus better than cultivated ones by using internal strategies, while cultivated tomatoes rely more on external phosphorus sources.

Contribution

This study reveals how domestication has altered phosphorus adaptation strategies in tomatoes through multi-omics analysis.

Findings

Cultivated tomatoes show higher biomass under sufficient phosphorus but are sensitive to deficiency.

Wild tomatoes maintain growth under low phosphorus through efficient internal phosphorus management.

Transcriptomic analysis shows regulatory divergence between wild and cultivated tomatoes in phosphorus response.

Abstract

Phosphorus (P) limitation is a major selective pressure in plant evolution and a persistent constraint on modern crop production. However, how domestication has reshaped P adaptation strategies remains poorly understood. Here, we compared wild (Solanum pimpinellifolium) and cultivated (Solanum lycopersicum) tomatoes under contrasting P conditions using integrated physiological, ionomic, and transcriptomic analyses. Our findings reveal distinct P strategies between the examined genotypes. Cultivated tomatoes achieved higher biomass under sufficient P supply but were highly sensitive to P deficiency, responding through acquisition-driven phenotypic plasticity characterized by extensive root remodeling and enhanced external P mobilization. In contrast, wild accessions maintained growth and higher P use efficiency under low P by relying on an optimized internal P management strategy, including efficient P uptake, preferential allocation to photosynthetically active tissues, and effective remobilization from older leaves. Consistently, ionomic profiling revealed that wild tomatoes preserved coordinated macro- and micronutrient homeostasis under P stress. Tissue-specific transcriptomic analyses further uncovered pronounced divergence in P-responsive regulation, with cultivated tomatoes showing predominantly root-centered responses, whereas wild accessions exhibited strong activation in old source leaves. This tissue-specific specialization was accompanied by a putative regulatory divergence, with HD-ZIP transcription factors enriched in cultivated tomatoes and G2-like and bHLH factors central in wild accessions. Together, our results indicate that modern cultivars exhibit a stronger reliance on external P acquisition and greater growth sensitivity under sustained P limitation compared to wild accessions, which showed relatively more stable internal P allocation patterns, highlighting wild germplasm as a resource for improving crop P efficiency.

Figures

Full text

1. Introduction

Phosphorus (P) is an essential macronutrient for plant growth and development, serving as a core component of nucleic acids, ATP, phospholipids, and signaling molecules [1]. Despite its biological significance, P availability in agricultural soils is severely constrained by fixation into insoluble forms, with approximately 40–60% of global arable land affected by P deficiency [2,3]. Consequently, agricultural systems heavily rely on substantial P fertilizer application to sustain crop productivity, yet this excessive P application diminishes P use efficiency (PUE) and significantly contributes to environmental degradation, particularly aquatic eutrophication driven by P runoff [4,5]. Therefore, understanding the molecular mechanisms of P nutrition and enhancing PUE is crucial for developing nutrient-resilient crops amidst dwindling phosphate rock reserves and environmental concerns.

Plants have evolved sophisticated adaptive strategies to scavenge, remobilize, and efficiently utilize scarce P resources. These strategies involve dynamic adjustments at multiple levels, including morphological changes in root system remodeling, enhanced secretion of organic acids and phosphatases, and strategic reallocation of internal P pools [6,7,8]. These adaptive responses constitute integral components of phosphate starvation responses (PSRs), which coordinate morphological, physiological, and molecular adaptations to phosphate limitation. At the molecular level, transcriptional reprogramming during PSRs is primarily mediated by MYB-CC-type transcription factors of the PHR family, which orchestrate the expression of phosphate starvation-induced genes [9,10,11,12,13]. PHR activity is dynamically regulated by SPX-domain proteins (named after SYG1, Pho81, and XPR1) and inositol polyphosphates (InsPs), ensuring precise control over P homeostasis [14,15,16,17,18,19]. Furthermore, efficient internal P utilization further involves the coordinated activity of vacuolar P transporters and PHT1-family transporters, which facilitate redistribution of phosphate between roots, leaves, and reproductive tissues to sustain growth under limited external supply [1]. Concurrently, a dual-affinity phosphate uptake system, comprising both low- and high-affinity transporters, enables dynamic adaptation to fluctuating external P availability [20]. Together, these mechanisms integrate P acquisition, internal redistribution, and utilization to optimize growth and survival under P-limiting conditions.

Crop domestication has profoundly reshaped plant nutrient acquisition systems. While selective breeding prioritized traits favorable to human needs, adaptive characteristics essential for stress resilience, including nutrient use efficiency, were frequently neglected or lost [21]. In cereals, domestication signatures include reduced seed dormancy, enlarged grain size, and loss of natural seed dispersal mechanisms, reflecting intense selection under resource-rich agricultural regimes [22,23,24,25]. This process coincided with habitat expansion and the adoption of management practices heavily dependent on external inputs such as irrigation, fertilizers, and pesticides. These anthropogenic interventions have substantially altered soil environments, particularly fertility profiles. It is widely recognized that many wild progenitors typically thrive in nutrient-poor soils, whereas modern cultivars tend to be adapted to comparatively fertile, managed conditions. Wild relatives often retain phenotypic plasticity that optimizes fitness under variable environments: under stress, they prioritize reproductive allocation (higher harvest index) at the expense of vegetative biomass, whereas resource abundance favors growth [22,26]. Such adaptive flexibility is often attenuated in modern cultivars, compromising environmental resilience. Root system architecture exemplifies this trade-off; despite its critical role in nutrient foraging, root traits received minimal direct selection during aboveground-focused domestication [27,28]. Wild relatives commonly exhibit superior root adaptations, such as enhanced root-to-shoot biomass partitioning (>60% under drought) and pronounced root apical, conferring stress tolerance yet inadvertently discarded during breeding in fertile, irrigated fields [27,28].

Tomato (Solanum lycopersicum L.) serves as a pivotal model for dissecting P nutrition due to its well-annotated genome, rich wild germplasm diversity, and global agricultural relevance. Wild relatives (such as S. pimpinellifolium and S. habrochaites), native to the P-impoverished Andean soils, frequently demonstrate superior P-foraging capacity and stress resilience compared to cultivated varieties [29]. Tomato underwent a two-stage domestication trajectory: from berry-sized S. pimpinellifolium (SP) to cherry-sized S. lycopersicum var. cerasiforme (SLC) in South America, and subsequently to big-fruited S. lycopersicum var. lycopersicum (SLL) in Mesoamerica [29]. This entire domestication process was characterized by severe genetic bottlenecks during various migratory events [29,30], leading to a substantial reduction in genetic diversity and extensive gene loss throughout its trajectory [31]. This progression entailed intense selection for fruit size loci, while traits like nutrient use efficiency were largely overlooked. Contemporary breeding priorities, such as yield, shelf life, and disease resistance, have further marginalized nutrient efficiency. Accordingly, cultivated tomatoes often exhibit attenuated root hair development, reduced organic acid exudation, and lower phosphatase activity under P deficiency relative to wild accessions [31,32]. Critically, genotypic variation in PUE directly influences productivity, and P-efficient tomato genotypes achieve up to 77% greater biomass than inefficient counterparts under low-P conditions [31,32].

Despite these observed phenotypic differences, the systemic mechanisms driving the divergence in P adaptation strategies between wild and cultivated tomatoes remain poorly explored. Specifically, how tissue-specific P allocation, internal homeostasis, and molecular regulation are coordinated to balance growth and stress resilience is largely unknown. Here, we address this knowledge gap through an integrated analysis combining phenotypic, ionomic, and tissue-specific transcriptomic profiling of wild and cultivated tomato accessions under low- and sufficient-P conditions. Based on our results, we hypothesize that wild tomatoes may employ a strategy centered on high P-use efficiency by maintaining stable internal P concentrations, enhancing tissue-specific P remobilization, and coordinating a broader suite of micro- and macro-nutrients, supported by conserved phosphate starvation response transcriptional networks. In contrast, cultivated accessions appear to prioritize high-acquisition strategies and root-centric growth plasticity, often at the expense of internal P stability. By dissecting these physiological and molecular divergences, this study reveals evolutionary trade-offs in P adaptation and highlights candidate genes and regulatory modules that could be harnessed to breed P-efficient tomato cultivars. Our findings provide a foundation for leveraging wild germplasm to improve crop resilience and sustainability under P-limited agroecosystems.

2. Results

2.1. Divergent Responses of Tomato Accessions to P Gradients

To investigate adaptive strategies of cultivated and wild tomatoes under varying P availability, we evaluated the growth and P-nutritional status of six selected tomato accessions, including four cultivated (S. lycopersicum; ZYF, ZF, MM, DY) and two wild relatives (S. pimpinellifolium; LA, PI), under low P (0.05 g kg^−1^ P, designated as LP) and moderate P (1.0 g kg^−1^ P, designated as MP) supplies (details in Materials and Methods). As expected, soil P availability significantly affected plant growth, with low P supply significantly inhibiting the plant height and shoot fresh weight (FW) of all accessions (Figure 1A–C). Notably, a clear divergence in growth pattern was observed between genotypes. Under MP, the four cultivated accessions (ZYF, ZF, MM, and DY) exhibited robust growth, with shoot FW ranging from 18.8 to 22.8 g. In contrast, the two wild accessions (LA and PI) showed significantly lower biomass accumulation, reaching only 11.5 to 16.8 g. Under LP conditions, shoot growth was severely inhibited across all genotypes, reducing the absolute biomass differences between cultivated and wild accessions (Figure 1C).

Two-way ANOVA revealed highly significant genotype by P level (G × P) interactions for plant growth (plant height and shoot fresh weight) as well as Pi and TP concentrations in both leaf and root tissues (Table S1, p < 0.001), indicating that these genotypes exhibit distinctly differential responses to varying P availability. Consistent with these statistical interactions, despite the overall reduction in biomass under LP, a genotype-specific response was evident in internal P concentrations. Inorganic phosphate (Pi) content in leaf and root, while generally decreasing sharply with low P supply, was maintained at higher levels in wild accessions (leaf: 0.033 and 0.025 mg g^−1^ FW; root: 0.013 and 0.01 mg g^−1^ FW) compared to most cultivated lines under LP (Figure 1D,E). This trend was further confirmed by total phosphorus (TP) concentrations in both leaf and root tissues (Figure 1F,G). Under LP conditions, wild accessions consistently exhibited significantly higher leaf and root TP concentrations than the cultivars. Leaf TP concentrations in wild accession LA were approximately 2.3 mg g^−1^ DW under LP, which was notably higher than the range observed in cultivated genotypes (1.161.88 mg g^−1^ DW) (Figure 1F). Similarly, root TP concentrations for wild accessions were 1.39 and 1.55 mg g^−1^ DW, substantially exceeding cultivated genotypes (1.191.26 mg g^−1^ DW) under LP conditions (Figure 1G). Conversely, under MP supply, cultivated tomatoes achieved higher TP concentrations than their wild relatives in both roots and leaves, aligning with their greater biomass (Figure 1F,G).

To quantify the genotypic differences in maintaining P status under stress, the ratios of P-related traits under LP versus MP conditions were calculated (Figure 2). These ratios consistently demonstrated that wild accessions LA and PI maintained significantly higher relative shoot and root Pi content, as well as shoot and root TP concentrations, compared to cultivated lines. For example, the LP/MP ratio for root TP concentration in wild accessions was around 0.5~0.6, while for cultivated accessions it was lower than 0.4 (Figure 2D). These results suggest that wild accessions exhibit a more stable internal P concentration profile when faced with P deprivation.

Together, these results reveal a clear divergence of cultivated and wild tomatoes in P adaptation: cultivated tomatoes exhibited maximized growth under optimal P but were highly sensitive to deficiency, whereas wild accessions maintained more stable internal P status under low-P conditions.

2.2. Differential Root Physiological Strategies for P Acquisition

To elucidate the physiological mechanisms driving the observed genotypic divergence in P adaptation, we selected four representative accessions (ZYF, DY, LA, and PI) for in-depth analyses under LP and MP conditions (Figure 3A). Two-way ANOVA indicated significant genotype × P level (G × P) interactions for root Pi concentration (p = 0.004), root TP concentration (p < 0.001), phosphorus uptake efficiency (PUpE; p < 0.001), and the activities of both acid phosphatase (ACP; p = 0.02) and alkaline phosphatase (ALP; p < 0.001) (Table S2). These interactions suggest distinct adaptive strategies of these genotypes to varying P availability, which we further explored.

We first investigated the root physiological responses to clarify the basis of P uptake. Root Pi content of some cultivated accessions like ZYF maintained comparable or slightly higher Pi levels than wild types under LP, contrasting with the MP conditions where cultivated lines accumulated significantly more root Pi (Figure 3B). Root TP concentrations in wild accessions, especially PI, were significantly higher than those in cultivated lines under LP conditions, whereas under MP conditions, root TP concentrations were higher in cultivated lines (Figure 3C). This divergence extended to P uptake efficiency (PUpE) (Figure 3D). Under LP conditions, wild accessions demonstrated significantly higher PUpE values (2.57 and 2.13 mg P g^−1^ for LA and PI, respectively) compared to cultivated lines (ZYF and DY, with 1.57 and 1.55 mg P g^−1^). However, under MP conditions, this trend was reversed or minimized, with cultivated lines showing comparable or higher PUpE values. These results indicate that wild tomatoes possess a superior capacity to efficiently acquire phosphorus per unit of root biomass, specifically under P-deficient environments.

To understand the P acquisition capability, we further analyzed root phosphatase activities. Interestingly, cultivated accessions exhibited a more aggressive acquisition strategy. Under P deficiency, the activities of both acid phosphatase (ACP) and alkaline phosphatase (ALP) were significantly induced in cultivated accessions, reaching levels markedly higher than those in wild accessions (Figure 3E,F). These results indicate that cultivated tomatoes may rely on a high-acquisition strategy involving enhanced rhizospheric activation, whereas wild tomatoes appear to adopt a different strategy that does not rely on high enzymatic exudation.

2.3. Genotypic Variation in P Allocation and Ionomic Profiles

Beyond root acquisition, the ability to optimize internal phosphorus distribution is critical for P adaptation. We therefore examined the P allocation patterns across different organs, which revealed that wild accessions allocated a greater proportion of approximately 69% (LA) and 57% (PI) accumulated P to the leaves under LP conditions, which was notably higher than the 51% (ZYF) and 50% (DY) observed in cultivated lines (Figure 4A). Moreover, wild accessions exhibited a lower proportion of P retained in old leaves compared to cultivated lines (Figure 4B), indicating a more efficient remobilization of P from senescing tissues to active growing organs.

To synthesize the physiological findings and systematically quantify the adaptive divergence, we calculated the Low Phosphorus Tolerance Coefficients (LPTC) for key traits (Figure 4C). This comprehensive analysis illuminated two contrasting strategies adopted by the different genotypes. Cultivated accessions exhibited superior resilience in structural metrics, maintaining higher LPTCs for morphological traits such as stem diameter (0.493 and 0.5) and root DW (~0.17), compared to wild relatives (stem diameter ~0.40 and root DW ~0.1). In contrast, wild accessions prioritized the maintenance of physiological P status. They displayed remarkably higher LPTCs for TP concentrations in both root (up to 0.807 in PI) and shoot (up to 0.492 in LA) tissues. Crucially, the tolerance coefficient for root PUpE in wild accessions was 0.306 and 0.364, which was approximately double that of cultivated accessions (0.154 and 0.163). Collectively, these quantitative indices suggested differential tolerance patterns, namely, that cultivated lines relied heavily on morphological adjustments to cope with P deficiency, while wild tomatoes prioritized the maintenance of internal P homeostasis and uptake efficiency.

We further explored whether this divergence in P efficiency is linked to broader nutrient homeostasis. Heatmap analysis of the ionomic profiles (Figure 5A) visualized a broad shift in nutrient accumulation patterns. Under MP conditions, cultivated and wild accessions exhibited distinct but relatively balanced element distributions. While under LP stress, wild accessions (LA and PI) maintained relatively higher Z-scores for a broad spectrum of elements in shoot tissues, including macronutrients P, N, K, and Na, as well as micronutrients such as B, Mn, Zn, and Fe, compared to cultivated lines, which showed generally lower relative abundances. This suggests that wild tomatoes possess a superior capability to maintain ionomic homeostasis and potentially synergistic nutrient uptake even under P-limiting conditions.

PCA results further revealed a fundamental and systematic restructuring of the ionomic profile driven by genotypes (Figure 5B). In shoot tissues under LP conditions, the first two principal components clearly separated the wild accessions from the cultivated ones, indicating distinct elemental signatures. Notably, the vectors for P, N, K, Mg and Mn were closely associated with the wild accession, whereas cultivated accessions clustered in the opposite direction, showing weaker associations with these key nutrients. This systematic ionomic divergence was not limited to low-P shoots. Root tissues and MP conditions confirmed that the ionomic signatures of wild and cultivated tomatoes are fundamentally distinct regardless of tissue type or P availability. Specifically, in P-deficient roots, wild accessions, particularly PI, clustered in a direction associated with P, K, N, Mn, B, and Fe vectors, distinguishing them from the cultivated group, which showed closer associations with Na, Mg, Zn, and Ca.

Collectively, these data demonstrated that wild tomatoes optimized P allocation and remobilization and also maintained a coordinated ionomic homeostasis under P stress, whereas cultivated lines exhibited morphological plasticity, reflecting divergent strategies in nutrient management.

2.4. Transcriptomic Reprogramming Under P Starvation

To unravel the molecular networks underlying the distinct physiological strategies, we conducted RNA sequencing analysis on 72 leaf samples (covering young and old leaves) and 24 root samples from the four selected tomato accessions (ZYF, DY, LA, and PI) under low and moderate P conditions. After removing low-quality sequences (quality scores < 25), a total of 655 Gb of high-quality data were obtained, with each sample yielding over 5 Gb of data. Filtered sequencing reads were aligned to the tomato reference genome (SL4.0, Solanum lycopersicum), achieving alignment rates ranging from 74.36% to 98.36% (Table S3). Principal component analysis (PCA) of the transcriptomic data confirmed high intra-group reproducibility for both leaf and root samples, indicating the reliability of our RNA-seq dataset in capturing specific transcriptional responses to varying P levels (Figure S1).

Gene expression levels were quantified as Transcripts Per Million (TPM) and normalized for downstream analysis. Differentially Expressed Genes (DEGs) were identified using the DESeq2 package in R, with thresholds of False Discovery Rate (FDR) < 0.01 and absolute fold-change (FC) > 2 (FC > 2 for up-regulation, FC < 0.5 for down-regulation). In total, 33,158 expressed genes were detected, comprising 27,622 known and 5536 novel genes. Similarly, 69,365 expressed transcripts were identified, including 27,067 known and 42,298 novel transcripts, which were not previously annotated in the reference genome.

Analysis of DEGs revealed a striking genotype-dependent, tissue-specific bias in transcriptional reprogramming. In leaf tissues, the transcriptional response exhibited strong positional specificity (Figure 6A). For up-regulated genes, the proportion of activated genes in old leaves relative to young leaves was notably higher in wild accessions compared to cultivated ones (Figure 6A). Specifically, in wild accession PI, approximately 79% of its total up-regulated leaf DEGs were found in L3 (3450 DEGs), whereas cultivated lines like DY typically showed a substantially lower proportion of around 28% in L3 (505 DEGs). In contrast, for down-regulated genes, cultivated accessions ZYF (85%) and DY (68%) and wild accession LA (91%) showed a larger proportion of DEGs in young leaves (L1 and L2) (Figure 6A). Notably, wild accession PI still presented a significant number of down-regulated DEGs in old leaves (L3: 1644). The number of common DEGs across different leaf positions within each variety was relatively low for both up- and down-regulated sets. In root tissues, the cultivated accession DY stood out with the most extensive transcriptional changes, showing 1328 up-regulated and 915 down-regulated exclusive DEGs, significantly exceeding those in other accessions. This data indicates that while cultivated tomatoes like DY undergo predominantly root-localized transcriptional reprogramming, wild tomatoes like PI exhibit a more pronounced response in senescing source leaves.

Further clustering of leaf DEGs from the four tomato varieties under different P supplies grouped all differentially expressed genes into six distinct sets (G1–G6; Figure S2), with DY and PI showing a relatively larger number of DEGs. Gene set G3 in ZYF and PI and G1 in DY and LA showed an up-regulated trend from young to old leaves, with higher expression in older leaves under LP conditions. In contrast, G1 in ZYF and PI and G3 in LA were more highly expressed in older leaves under MP conditions. Conversely, G4 in ZYF, G2 in DY, and G2 and G6 in PI were up-regulated from old to young leaves and showed higher expression in young leaves under LP, while G6 in ZYF and G3 in DY exhibited elevated expression in young leaves under MP. G5 in ZYF, DY, and PI, together with G4 in LA, displayed the highest expression in leaf L2. Clustering of root DEGs (Figure S2) also identified six gene sets. Sets G1, G2, and G3 were up-regulated under MP but suppressed under LP and were predominantly enriched in varieties DY, LA, and PI, respectively. In contrast, sets G4, G5, and G6 were induced by low P stress and showed lower expression under moderate P, with these responses being prominent in varieties ZYF, LA and PI, and PI, respectively.

To further characterize the low P response across tomato varieties, leaf and root DEGs were grouped into 15 major clusters representing dominant tissue-specific expression patterns (Figure S3). In leaves, clusters such as 1 (ZYF/DY), 10 (LA), and 12 (PI) showed low expression in young leaves but high expression in old leaves, especially under LP. Clusters 2 and 4 (ZYF), 8 (DY), 7 (LA), and 5 (PI) decreased from young to old leaves. Young leaves under MP displayed high expression in clusters 3 (ZYF/DY), 1 (LA), and 9 (PI), whereas young leaves under LP were marked by clusters 7 (ZYF), 6 (DY), 14 (LA), and 10 (PI). Clusters 10 (ZYF), 11 (DY/LA), and 2 and 3 (PI) were highest in the second young leaf. Stress-specific induction was observed in old leaves (L3) of clusters 6 (ZYF), 15 (LA), and 12 (PI), and in young leaves (L1) of cluster 15 (DY) under LP only. In roots, clusters 2 and 8 (DY) and 15 (LA) were highly expressed specifically under MP. Low P stress induced root-specific expression in cluster 7 (ZYF) and cluster 9 (PI), highlighting distinct responses between cultivated and wild tomatoes. Cluster 6 showed consistently higher expression under LP in DY, LA, and PI, suggesting a core set of low-P responsive genes.

Beyond the overall transcriptional changes, we sought to identify the core transcriptional regulators orchestrating these divergent P responses. To analyze the regulatory patterns of tomato leaves under LP stress, we constructed gene co-expression networks for each accession. Using the GENIE3 R package, a random forest approach with 1000 decision trees was employed to infer transcription factor (TF) prioritization networks based on co-expression profiles from our tissue-specific transcriptomic data. For enhanced visualization and analysis, networks were filtered to display and annotate only the top 1000 co-expressed interactions by weight (Figure 6B, Table S5).

Visual inspection of these networks immediately suggested a distinct organizational principle between cultivated and wild accessions (Figure 6B). Quantifying the top TFs by their degree (connectivity) further confirmed this divergence (Table S5). In cultivated accessions, the core transcriptional regulators primarily belonged to the HD-ZIP, MYB_related, and ARF families. Among these, HD-ZIP TFs consistently ranked high with degree values up to 131 in DY and 124 in ZYF. Conversely, wild accessions predominantly utilize G2-like and bHLH TF families as core regulators. In LA, G2-like TFs exhibited the highest degree (145), followed by WRKY and bHLH. Similarly, in PI, bHLH and G2-like TFs were among the top regulators, with degrees of 161 and 80, respectively.

Together, core transcriptional regulator analysis revealed genotype-specific differences in their composition patterns under P limitation. Specifically, cultivated tomatoes exhibited root-centric transcriptional networks with HD-ZIP dominance, whereas wild tomatoes showed activated old-leaf-centric networks coordinated by G2-like and bHLH TFs.

2.5. Expression Dynamics of Key Phosphate Starvation Response Genes

To further dissect the molecular mechanisms underlying the diverse P responses observed between cultivated and wild tomato accessions, we identified 19 core gene families involved in P homeostasis and signaling, comprising 63 genes, based on homology to known phosphate starvation response (PSR) genes in Arabidopsis thaliana (Table S4). Expression profile analysis revealed that the induction of these genes was generally more robust in roots than in leaves, and within leaves, more pronounced in old leaves (L3) compared to young leaves (L1/L2) (Figure 7). This tissue-specific pattern suggests a coordinated strategy where old leaves act as primary sensors and sources for P remobilization under deficiency, while roots actively respond to external P scarcity.

Specifically, key P signaling and homeostasis regulators exhibited distinct tissue- and accession-specific expression patterns. SlPHR2 was predominantly expressed in old leaves and roots, serving as a central hub for P-responsive transcription. SlPHR4 was strongly induced only in roots, with minimal expression in leaves. SlPHR1 and SlPHR3 were specifically induced in old leaves but not in roots and showed stronger induction in wild accessions. Concurrently, SPX family genes, negative regulators of PHR, were significantly up-regulated under P deficiency, particularly in old leaves and roots, while SlPHO2 was mainly induced under MP conditions in old leaves and roots. SlPHF1 showed higher expression in old leaves relative to young leaves and was more strongly induced in roots under P deficiency, with notably higher expression in wild accessions. Genes regulating iron content, such as SlHRZ1, SlHRZ3, and SlHRZ4, were predominantly induced in roots, with higher expression observed in wild tomatoes, indicating coordinated regulation of P and Fe homeostasis.

Further insights into P redistribution and efficient utilization were gained through analyzing coordinated regulation of vacuolar P transport systems and membrane transporters. Vacuolar P efflux genes (SlGlpTs) were significantly induced in old leaves under P deficiency, with SlGlpT3 notably showing strong induction in DY roots (Figure 7), providing transcriptional evidence for efficient P remobilization. Additionally, transporter genes involved in xylem loading of phosphate (SPX-EXS) exhibited stronger induction in old leaves of cultivated accessions under low P conditions compared to wild types. Conversely, vacuolar P influx transporter genes (SlSPX-MFSs) showed higher expression under MP conditions in leaves. Two E3 ubiquitin ligases (SlSPX-RING), which can degrade phosphate transporters, were mainly expressed in roots and exhibited stronger expression under MP conditions compared to LP conditions. Interestingly, the phosphate transporter SlPHT1;5 was primarily induced in roots under P deficiency, with remarkably stronger expression in the wild accession PI. Consistently, the acid phosphatase gene family (SlPAP) also showed P-deficiency induction in roots, and specific members like SlPAP6, SlPAP9, and SlPAP17 exhibited notably stronger transcriptional abundance in the roots of wild accessions (LA and PI) compared to cultivated lines (Figure 7). Additionally, SlPAP genes also showed relatively higher expression in the old leaves of wild accessions.

Collectively, these expression dynamics revealed that wild tomatoes exhibited high gene expression related to root P acquisition and internal P remobilization, particularly from old leaves, whereas cultivated accessions showed higher expression of specific transporters in leaves under deficiency. This coordinated regulation of PSR genes, vacuolar transporters, and phosphatases underpins the contrasting P adaptation strategies between wild and cultivated tomatoes, linking transcriptional responses to the observed physiological and ionomic divergences.

3. Discussion

P is a critical macronutrient that frequently limits plant growth and agricultural productivity globally, especially in naturally P-deficient soils [8]. Low P availability is known to induce a cascade of adverse physiological and morphological symptoms in plants, including suppressed growth, leaf chlorosis, slender stems, and impaired reproductive development, ultimately leading to significant yield losses [33,34]. Recognizing that modern agricultural reliance on intensive P fertilization for maximal yields is unsustainable due to environmental degradation and finite phosphate reserves, improving crop P use efficiency stands as a critical goal for sustainable agriculture. Our study aimed to unravel how the domestication process has influenced P adaptation strategies in tomato by comparatively analyzing phylogenetically distinct wild (S. pimpinellifolium) and cultivated accessions (S. lycopersicum) across varying P availabilities. Our findings reveal that wild and cultivated tomatoes have evolved divergent P adaptation strategies: cultivated lines primarily exhibit robust morphological plasticity and P acquisition capacity, whereas wild lines prioritize efficient physiological P utilization and internal P homeostasis under P limitation. These insights are crucial for informing future breeding efforts to develop P-efficient tomato varieties.

3.1. Domestication-Associated Divergence in Growth and Morphological Adaptation

P deficiency imposed a universal constraint on vegetative growth across all studied tomato genotypes, but the magnitude of this growth penalty revealed a distinct domestication-associated divergence in adaptive strategies between cultivated and wild tomatoes (Figure 1). Although cultivated accessions possessed superior biomass potential under optimal P, they suffered a more precipitous decline in growth under P scarcity compared to their wild relatives (Figure 4). This divergence aligns with the concept of the domestication syndrome, where selection for high yield in nutrient-rich agricultural environments often inadvertently compromises traits essential for stress tolerance in resource-poor natural habitats [24]. Consistent with this, recent studies have documented similar trade-offs between growth potential and stress resilience in other crops, suggesting that modification of nutrient acquisition strategies may be a common consequence of domestication [35].

The physiological basis of this divergence likely involves distinct resource allocation strategies. The rapid growth typical of cultivated lines demands substantial P resources, rendering them highly vulnerable when P is scarce. In contrast, wild types, often adapted to nutrient-poor soils, appear to employ a more conservative growth strategy, prioritizing the maintenance of internal P homeostasis over maximizing biomass, which is a classic survival mechanism in oligotrophic environments [8,36]. This conservative strategy, characterized by efficient P remobilization and reduced growth under limitation, has been well documented in native species from phosphorus-impoverished habitats [37]. From an applied perspective, the stress-resilient traits preserved in wild germplasm represent valuable genetic resources for breeding programs aimed at improving phosphorus efficiency in future cultivars [38,39].

3.2. Contrasting Root-Level P Acquisition and P Status Management

Root physiological responses further distinguished the adaptive strategies. Under P limitation, wild accessions maintained significantly higher root TP concentrations, indicating a greater capacity for root P retention or storage, whereas cultivated lines exhibited higher activities of both ACP and ALP in roots, reflecting an aggressive investment in external P acquisition (Figure 3). However, this high acquisition capacity in cultivated lines did not translate into superior tolerance under sustained P deficiency, as their root P uptake efficiency was notably lower than that in wild accessions. This trade-off between short-term acquisition capacity and long-term sustainability has been observed in other crop species, suggesting that selection for rapid growth under high fertility may inadvertently uncouple investment in acquisition from maintenance of internal P status [40]. The observed contrasting strategies reflected distinct evolutionary solutions to this shared adaptive challenge.

Beyond root-level acquisition, our data highlight the critical role of internal P management. Wild accessions maintain higher shoot TP concentrations through superior P remobilization from old to young leaves, as evidenced by their higher young/old leaf P ratio (Figure 4). This implied efficient translocation of P from senescing tissues to active growing organs, prioritizing P supply for new leaves and maintaining their basic photosynthetic capacity, whereas cultivated tomatoes displayed P retention in old leaves, potentially linked to photosynthetic utilization in those tissues [41,42]. Efficient P use within the plant, particularly through superior P remobilization, is a hallmark of P-efficient plants [43]. The ability to effectively transport P from roots to shoots can also enhance plant immunity and bolster resistance to various abiotic and biotic stresses [7,44]. These findings collectively underscore that the artificial and natural selection pressures have led to distinct survival strategies in tomato [29,31].

These physiological differences were reflected in distinct ionomic signatures (Figure 5). The robust genetic separation indicates that domestication has fundamentally reshaped the entire nutrient homeostasis network. The ionome, serving as a comprehensive fingerprint of a plant’s elemental composition, reflects its nutrient acquisition, transport, and utilization strategies [45,46]. Under P deficiency, wild root tissues clustered along vectors associated with P, K, N, Mn, B, and Fe, whereas cultivated roots were associated with Na, Mg, Zn, and Ca (Figure 5). This coordinated regulation of multiple nutrients in wild tomatoes suggested that P efficiency is an integrated component of a system-wide ionomic strategy shaped by long-term adaptation to low-P soils [47]. Recent comparative studies across multiple crop species have similarly documented that wild relatives engage in broader nutrient network coordination, while domesticated lines exhibit more specialized but potentially less resilient ionomic profiles [48,49]. In lablab, wild accessions maintained higher leaf concentrations of zinc and iron than cultivated lines under drought stress, with differences linked to phosphorus metabolic responses [48]. Similarly, analysis of a genetically diverse set of wild soybean accessions revealed that average concentrations of multiple elements were higher in wild than domesticated soybean, suggesting erosion of ionomic breadth during domestication [49].

3.3. Molecular Mechanisms Underpinning Divergent P Adaptation Strategies

Integrated transcriptomic analysis revealed distinct tissue-specific patterns of transcriptional reprogramming underpinning the varying P adaptation strategies between cultivated and wild tomatoes. The cultivated accession DY exhibited a root-centric transcriptional response, aligning with its physiological strategy of enhanced external P acquisition. In contrast, the wild accession PI mounted a pronounced transcriptional response specifically in older leaves, providing molecular evidence for its prioritized internal P remobilization and conservation (Figure 6). This tissue-specific divergence in transcriptional investment was observed under P limitation. Underlying these patterns was a differential engagement of core transcriptional regulators. Co-expression network analysis indicated that cultivated accessions showed an enrichment of HD-ZIP TFs, whereas wild accessions were enriched for G2-like and bHLH TFs as central regulators (Figure 6). This suggests a putative divergence of the P-stress signaling network, where these distinct TFs may be associated with root-centric developmental adjustments or nutrient recycling and stress tolerance pathways [50,51].

A detailed expression analysis of key PSR genes further demonstrated the specific molecular components of each strategy. Tomato exhibited a functionally diversified PHR-SPX regulatory module. While SlPHR1/2/3 were active in older leaves, SlPHR4 was specifically and strongly induced in roots (Figure 7), contrasting with the single dominant PHR orthologs in Arabidopsis and rice and suggesting evolved tissue-specific regulatory partitions [9,10,11,12,13]. For P acquisition, genes associated with high-affinity uptake (SlPHT1;5) and organic P mineralization (e.g., SlPAP6, SlPAP9) were more robustly induced in roots of wild accessions, aligning with their efficient root foraging phenotype [1]. Conversely, the cultivated strategy appeared less reliant on transcriptional upregulation of these specific uptake genes. Regarding internal P management, wild accessions showed stronger induction of vacuolar phosphate efflux transporters (SlGlpTs) in older leaves, facilitating remobilization of stored P [52,53,54,55], while vacuolar influx transporters (SlSPX-MFSs) were more associated with P-sufficient conditions. Furthermore, the enhanced expression of P translocation genes (e.g., SlPHO1 homologs) in wild types likely supports their observed elevated leaf P content [56].

Ultimately, this understanding is vital for sustainable agriculture. The rich genetic diversity present in wild tomato relatives offers invaluable genetic resources for breeding P-efficient cultivated varieties. The identification of putative P-adaptive genes and mechanisms from wild types provides a basis for the future verification and development of new tomato cultivars that could be less dependent on intensive P fertilization, thereby potentially mitigating environmental impacts and contributing to enhanced agricultural sustainability. This approach can help overcome the domestication syndrome, where selection for yield may have inadvertently compromised stress tolerance.

3.4. Study Limitations and Future Directions

While this study provides an integrative multi-omics comparison of P adaptation between cultivated and wild tomato accessions, several limitations should be acknowledged.

First, our conclusions are derived from comparative analyses of four cultivated and two wild accessions and therefore demonstrate group-level divergence rather than direct genetic causation. While the observed physiological and regulatory differentiation is consistent with divergence accompanying tomato domestication, we did not perform phylogenetic reconstruction, selection analyses, or genetic mapping to identify specific domestication loci. Accordingly, the patterns reported here should be interpreted as domestication-associated rather than strictly domestication-driven.

Second, the wild accessions analyzed belong to Solanum pimpinellifolium, the closest progenitor of cultivated tomato, but do not represent the full diversity of wild Solanum species. Although the two wild accessions exhibited consistent responses under P limitation, broader sampling across additional taxa will be necessary to assess the generality of the identified regulatory and physiological patterns.

Third, transcriptional regulatory divergence was inferred from co-expression network modeling and differential expression analysis without functional validation. Thus, the identified hub transcription factors represent testable regulatory hypotheses rather than experimentally confirmed causal regulators. Future experimental approaches such as CRISPR-mediated perturbation, promoter-reporter assays, chromatin immunoprecipitation, and downstream target validation will be necessary to delineate the complete mechanistic pathway linking transcriptional regulation to specific physiological traits. In addition, responses were assessed at a single developmental stage under sustained stress, precluding resolution of temporal dynamics. Future genetic perturbation and time-resolved multi-omics analyses will be required to establish mechanistic hierarchies and response kinetics.

Collectively, these limitations do not diminish the robustness of the comparative patterns reported here but delineate the scope of inference. The present study provides a systems-level framework and testable regulatory hypotheses that establish a foundation for future genetic, functional, and ecological investigations into phosphorus adaptation in tomato.

4. Materials and Methods

4.1. Plant Materials and Culture System

Six tomato genotypes were used in this study: four cultivated varieties (S. lycopersicum var. cerasiforme ZheYingFen No. 1, ZheFen No. 702, Money-Maker, and DongYangJuXing) and two wild relatives (S. pimpinellifolium LA1589 and PI438898) (hereafter termed ZYF, ZF, MM, DY, LA, and PI, respectively). Seeds were surface-sterilized with 5% NaClO for 10 min, rinsed thoroughly with distilled water, and pre-germinated on moist filter paper in Petri dishes at 25 °C in darkness for 72 h. Uniform seedlings were transplanted into soil in a greenhouse with a 16 h light/8 h dark photoperiod, with photosynthetically active radiation maintained at approximately 150 μmol m^−2^ s^−1^, a day/night temperature regime of 28 °C/24 °C, and 65% relative humidity. To minimize positional effects and ensure uniform light exposure, pots were randomly arranged and repositioned weekly throughout the experimental period. For each physiological parameter measured, three biological replicates were established per treatment combination. Each biological replicate comprised 3–6 independently grown plants.

4.2. Experimental Design and P Treatments

The experiment was conducted using a completely randomized design. The experiment comprised two P treatment levels: moderate phosphorus (MP, 1 g P kg^−1^ soil) and low phosphorus (LP, 0.05 g P kg^−1^ soil) supply. The test soil was prepared by mixing low-P soil, vermiculite, and sand in a 1:1:1 volume ratio, then portioned into containers. Based on soil weight, corresponding amounts of KH_2_PO_4_ were applied to achieve the MP and LP levels (1 and 0.05 g P kg^−1^ soil, respectively). For the LP treatment, K_2_SO_4_ was added to supplement an equimolar amount of potassium. Additionally, all treatments were supplemented with a P-free Hoagland nutrient solution applied in sufficient quantity to ensure an adequate supply of other essential mineral elements.

4.3. Phenotypic Documentation and Physiological Indicators Determination

Plant phenotypes were documented at the end of the experiment using a digital camera (Canon EOS 5D). For biomass determination, plants were separated into roots, stems, and leaves. Fresh weights were recorded immediately after harvesting, followed by drying at 70 °C for 72 h to obtain dry weights. The following traits were calculated to assess P acquisition, utilization, and tolerance:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The Low Phosphorus Tolerance Coefficient (LPTC) for each trait was calculated as the ratio of its value under low P conditions to its value under moderate P conditions [57,58,59].

4.4. Elemental Analysis and Soil Phosphatase Activity Quantification

For inorganic phosphate (Pi) determination, frozen tissue samples (100 mg) were homogenized in 1 mL of 0.5 M H_2_SO_4_ and incubated for 30 min. After centrifugation (12,000× g, 10 min), the supernatant was used for Pi quantification using the molybdenum-antimony colorimetric method. For elemental analysis, dried samples were digested in concentrated H_2_SO_4_-H_2_O_2_, and the elemental concentration was determined by inductively coupled plasma optical emission spectrometry (ICP-OES; Thermo Scientific, Waltham, MA, USA).

For the soil phosphatase activity, the rhizospheric soil, approximately 1 mm of surface soil adhering to the root, was used. The samples (0.5 g) were incubated in 2 mL of assay buffer (pH 8.8 for ALP and pH 6.5 for ACP containing 5 mM p-nitrophenyl phosphate) for 1 h at 37 °C. The reaction was terminated by adding 1 mL of 0.5 M NaOH, and the absorbance of p-nitrophenol released was measured at 410 nm using a TECAN spectrophotometer (Infinite^®^ F50 microplate reader), (Tecan Austria GmbH, Grödig, Austria). Specific activity was expressed as nmol p-NP h^−1^ g^−1^ soil.

4.5. RNA Extraction, Sample Library Preparation and Sequencing

For tissue-specific transcriptomic analysis, roots and leaves were washed with clean water, and tissues were excised using a surgical blade. Tissue samples (young leaves, old leaves, and roots) were collected from three biological replicates per treatment and immediately frozen in liquid nitrogen. A total of 96 samples of leaves and roots of cultivated (S. lycopersicum; ZYF and DY) and two wild relatives (S. pimpinellifolium; LA and PI) tomatoes were collected. Total RNA was extracted using TRIzol reagent (Invitrogen, Waltham, MA, USA) according to the manufacturer’s protocol. RNA quantity and quality were assessed using a NanoDrop spectrophotometer and Agilent Bioanalyzer 2100 system (Thermo Fisher Scientific, Waltham, MA, USA). RNA of each sample was sequenced by Illumina NovaSeq 6000 (Illumina, Inc., 5200 Illumina Way, San Diego, CA, USA).

4.6. Transcriptome Analysis

The high-quality sequencing reads were analyzed using pipelines as described previously [60]. The genome version SL4.0 of S. lycopersicum was used as a reference for transcriptome mapping. In brief, fastp was used to perform adapter trimming and quality filtering, and high-quality sequencing reads were mapped to the tomato reference genome using HISAT2. Transcript quantification was performed by featurecount. Transcripts Per Million (TPM) values were used to normalize gene expression levels by StringTie (v2.2.1), and a gene was defined as an effectively expressed gene if its average TPM values of biological replicates were more than 3.

Differentially expressed genes (DEGs) analysis was performed by the R package DESeq2 with raw counts of each gene. For each gene, the p-value and fold change in expression between the samples were obtained. The p-value was adjusted by the BH method to control the FDR and correct multiple tests (Benjamini and Hochberg, 1995). Genes with FDR < 0.01 were considered a significant threshold for detecting DEGs, and fold-change > 2 and fold-change < 0.5 for up-regulation and down-regulation, respectively. Gene clustering analysis of time-course RNA-seq data was carried out by Mfuzz (Kumar and E. Futschik 2007) with DEGs identified above.

4.7. Transcription Factor (TF) Identification and Co-ExpWression Network Analysis

The full list of transcription factors in tomato was obtained from PlantTFDB v5.0 (http://planttfdb.gao-lab.org/) [61]. A TF-prioritized network was constructed using the GENIE3 R package [62]. The random forest method with 1000 decision trees was used for network construction. The top 2000 items filtered by weight were extracted for network visualization. The resulting network was visualized using the Cytoscape v3.9 software [63].

4.8. Identification of Phosphate Starvation Response Genes in Tomato

Identification of core genes in phosphate starvation response in tomato was performed by a similar method as our previous study [54]. In brief, nine core gene families were identified by BLASTP with an E-value cutoff at 1 × 10^−10^, using typical genes of each family from Arabidopsis. Candidate subject genes validated by InterProScan 5 [64] must contain the completed core domain of each family. Finally, a total of 63 core phosphate starvation response genes were identified in the tomato genome, covering phosphate signaling, phosphate transport, and phosphorus metabolism.

4.9. Statistical Analysis

All data were subjected to analysis of variance (ANOVA) using R software (v4.1.0). Significant differences between treatments were determined by Tukey’s honest significant difference (HSD) test at p < 0.05. To evaluate the main effects of genotype, P treatment, and their interaction on plant growth and physiological parameters, a two-way ANOVA was performed. When significant interactions were detected, differences among group means were further analyzed using Tukey’s HSD post hoc test at a significance level of p < 0.05. All statistical analyses were conducted using GraphPad Prism 9.0. Principal component analysis (PCA) and Pearson correlation analysis were performed to evaluate relationships between physiological parameters and gene expression patterns. Heatmaps of differentially expressed genes were generated using the pheatmap package in R. All figures were prepared using GraphPad Prism 9.0.

5. Conclusions

By integrating physiological, ionomic, and tissue-specific transcriptomic analyses, this study reveals contrasting strategies of P management in cultivated and wild tomatoes. Cultivated accessions rely primarily on an acquisition-driven strategy, with enhanced root remodeling and rhizospheric phosphatase activity, resulting in higher biomass under sufficient P but greater growth sensitivity under P limitation. Wild accessions, in contrast, maintain more stable growth under low P through efficient internal P allocation and remobilization, preserving coordinated nutrient homeostasis. Consistent with these physiological patterns, ionomic analyses revealed that wild accessions preserved coordinated macro- and micronutrient homeostasis under P stress, indicating that cultivated and wild accessions differed in the extent of system-level buffering responses under P limitation. These physiological differences were accompanied by a putative regulatory divergence, with HD-ZIP transcription factors enriched in cultivated tomatoes and G2-like and bHLH factors central in wild accessions. Together, our results indicate that modern cultivars exhibit a stronger reliance on external P acquisition and greater growth sensitivity under sustained P limitation compared to wild accessions, which showed relatively more stable internal P allocation patterns. We note that the limited number of genotypes and the inference-based nature of regulatory networks preclude definitive attribution of these patterns to domestication or functional confirmation of the predicted regulatory mechanisms. Future studies with broader sampling and functional validation will be required to substantiate these domestication-associated differences.