From microbiome to biostimulants: unlocking the potential of tomato root endophytes

Abstract

Background

Microbe-based biostimulants offer a sustainable and promising alternative to synthetic inputs, potentially reducing or replacing conventional inputs in crop management. Studying the native microbiota, particularly endophytic microbes, helps in selecting those that are naturally adapted to persist and to enhance plant growth under specific environmental conditions. This study aims to define the endophytic microbiota adapted to tomato crops by selecting discriminant amplicon sequence variant (ASVs) that are enriched during key plant growth stages and found in the core microbiota.

Results

This study presents a large-scale analysis of tomato root endophytic prokaryotic microbiota using 16 S sequencing across the most common and widespread conditions used for tomato cultivation, offering comprehensive insight into its structure and dynamics. The results revealed a predominance of the Actinobacteriota and Proteobacteria phyla; less abundant groups included Bacteroidota, Verrucomicrobiota, Patescibacteria, and Firmicutes. Core microbiota analysis and discriminant ASV identification across different plant growth stages enabled the selection of the most abundant and persistent taxa adapted to the tomato endorhizosphere. Streptomyces, Shinella, Devosia, and Pseudoxanthomonas, as well as the lesser known genera Variovorax, Pseudarthrobacter, and Lechevalieria, represented the key genera identified, suggesting long-term host‒microbe associations.

Conclusions

The description of the representative framework of the tomato-associated microbiota and the identification of its most important components provide a basis for developing tailored microbial formulations that can increase crop resilience and reduce dependence on synthetic agricultural inputs, aimed at developing more sustainable environmental management strategies.

Background

In the past few decades, the human population has doubled, thereby requiring a similar increase in food production [1]. This has led to the use of massive amounts of chemical fertilizer in agricultural production, which negatively affects plant, soil, and water ecosystems and translates into risks for human health and environmental preservation [2]. In addition, plant domestication and crop breeding under fertilization have altered the microbiome associated with plants, which may have caused the selection of plant genotypes that have lost the ability to recruit large numbers of microorganisms that support plant growth [3]. The use of plant growth-promoting bacteria (PGPB) offers a sustainable alternative to traditional chemical inputs [4, 5]. The importance of plant–microbe interactions has only recently been studied and better elucidated; thus, reengineering the characteristics of beneficial rhizosphere microbiomes back into agricultural cropping systems is possible [3]. This calls for more research focused on microbial rhizosphere engineering to redesign plant microbiomes to maintain or even improve current productivity.

The identification of PGPB that can be combined to form a specifically tailored community of compatible plant-beneficial strains is of great interest. Compared with synthetic products, their use still has certain drawbacks, mostly due to their greater environmental sensitivity, since the efficacy of PGPB may vary in different crop species and environmental conditions. However, the use of microbial consortia, rather than applications of single microorganism, appears to have many advantages, such as increased stability across a wider range of environmental conditions or the opportunity to create synergistic interactions that can impart multiple desirable characteristics [6]. A critical step in the effective selection of PGPB is the characterization of the indigenous microbiome to identify microbial populations that can persist on plants under several conditions [7]. Among these, endophytic microbial populations live inside plant tissues without harming their host [8]. They are known to increase plant growth by producing growth-stimulating hormones [9, 10] and by exerting typical plant growth-promoting activities (e.g., nitrogen fixation and phosphate adsorption) [11]. Furthermore, endophytic populations can increase plant defenses owing to their ability to generate antimicrobial compounds, such as salicylic acid and ethylene [11,12,13]. Therefore, investigating endophytes is of great interest because of their great potential for sustaining plant growth and thus reducing synthetic inputs in agriculture.

Here, we propose an in silico approach based on the analysis of high-throughput 16 S amplicon sequencing data, which can reveal the composition and behaviors of the bacterial microbiota adapted to the tomato endorhizosphere. We considered a high number of samples and diverse experimental variables, which include the most common and widespread conditions used for tomato cultivation, to obtain information on the endophytic microbial population that is shared between the different conditions. We generated a hypothetical microbial consortium by selecting the amplicon sequence variant (ASVs) from a discriminant analysis of endophytic bacteria enriched in the flowering and harvesting growth stages (exhibiting high abundance and prevalence) and from the core microbiota (i.e., detected in at least 80% of the samples). This approach allowed us to evaluate the presence of specific microbial groups related to certain experimental variables or that are always associated with the host, revealing transient and long-lasting relationships, to lay the groundwork for the development of a tailored microbial formulation for enhancing plant growth.

Results

The endophytic microbiota varied strongly across host plant growth stages

The prokaryotic endophytic microbiota associated with tomato roots varied strongly with the host plant growth stage, showing three different groups, especially for the Sicily samples (Fig. 1A).

Fig. 1

Influence of the experimental variables on endophytic microbiota diversity. (A), Beta diversity (multidimensional scaling MDS using Bray‒Curtis distance) for the n = 192 samples highlights variations with respect to tomato plant growth stage (post-transplant, flowering, and harvesting) and sampling region (Campania and Sicily). (B), Percentage of the total microbiota variation influenced by the experimental variables: sampling time, sampling region, tomato cultivar, sampling area, and treatment

Indeed, sampling time, that is, post-transplant (Pt), flowering (F) and harvesting (H), influenced 25.1% of the total bacterial microbiota variability, and the combination of other variables (i.e., sampling region, tomato type, sampling area, and treatment) contributed to 7.3% (Fig. 1B). Of note, no differences were detected in the different varieties of elongated tomatoes (i.e., Pietrarossa F1 and Maraskino), which we therefore grouped in the same tomato type. The sample biodiversity varied significantly with respect to host plant development, with a statistically significant reduction from the Pt stage to the F stage in both the Campania (Fig. 2A-B) and Sicily (Fig. 2C-D) samples.

Fig. 2

Endophytic microbiota biodiversity and structural variations across host plant growth stages. (A-D), Differences in terms of alpha diversity for the Campania (A and C) and Sicily (B and D) samples according to the Shannon index (C and D) and observed richness (A and B). (E-F), Average relative abundance of the main phyla detected in the Campania (E) and Sicily (F) samples and stratified across tomato growth stages (post-transplant, flowering, and harvesting)

The main phyla followed concordant variations in both regions. The most abundant phylum overall, Proteobacteria, decreased significantly in abundance from the Pt stage to the F stage (average relative abundance from 56.3 to 41.5% in the Campania (Fig. 2E) and from 60.2 to 43.8% in the Sicily (Fig. 2F) samples). Similarly, the abundance of Verrucomicrobiota and Planctomycetota decreased significantly from Pt to F in both regions, whereas that of Bacteroidota and Chloroflexi decreased in both regions but significantly only in the Campania samples. These findings indicate temporal shifts in microbial composition, although causal relationships would require additional experimental validation or mechanistic studies. The second most abundant phylum, Actinobacteriota, exhibited a significant increase in abundance from Pt to F (from 23.1 to 37.2% in Campania and from 23.8 to 40.5% in Sicily). Another phylum of interest was Patescibacteria, whose abundance increased significantly from Pt to F (from 3.0 to 8.7% in the Sicily and from 1.7 to 6.8% in the Campania samples), and that of Firmicutes also increased significantly in the F stage in both regions. When comparing the later timepoint H with F, most of the main phyla showed significant variations and tended to return to the Pt condition. Indeed, the abundance of Actinobacteriota and Patescibacteria decreased from the F to H stages (but only significantly so in Sicily), whereas that of Proteobacteria and Verrucomicrobia significantly increased in the H stage in both sampling regions.

Several hundred ASVs were differentially abundant during plant development

We performed a more in-depth analysis and identified specific ASVs discriminating between the growth stages of tomato plants. After removing low-abundance ASVs (< 0.005% average relative abundance across all samples), we identified 1,191 ASVs that were differentially abundant (adjusted p value < 0.1) between Pt and F in Campania and 1,157 in Sicily. When comparing F and H, only 396 ASVs were differentially abundant (adjusted p value < 0.1) in Campania and 859 in Sicily. Finally, 1,142 (in Campania) and 1,031 (in Sicily) ASVs were discriminant between Pt and H. This resulted in 1,715 ASVs that were significant in at least one of the comparisons, and 87 of them were detected as abundant in our set of samples (i.e., average relative abundance ≥ 0.7% in at least one of the 12 groups represented in Fig. 3). More specifically, 49 abundant ASVs increased in abundance from the Pt stage to the F stage, whereas the other 38 decreased in abundance from the Pt stage to the F stage (Fig. 3, Supplementary Data 2, and Supplementary Data 3). Among the 49 ASVs enriched in F, 21 belonged to the Actinobacteriota phylum, including Streptomyces (n = 5), Amycolatopsis (n = 4), Micromonosporaceae (n = 3), Lechevalieria (n = 2), and Kribbella (n = 2). Another highly represented phylum was Proteobacteria, with 18 ASVs, including Gammaproteobacteria (n = 3, Pseudomonas 87, Microbulbifer okinawensis, and Stenotrophomonas 44), Betaproteobacteria (n = 3, Variovorax 15, and Variovorax paradoxus 1 and 2) and Alphaproteobacteria (n = 12, including Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium 50, 51 and 35, Aminobacter 20, Bosea 34, Ensifer 11, Devosia riboflavina, Shinella 17, and Sphingobium herbicidovorans). The remaining 10 ASVs belonged to the Patescibacteria (n = 3; LWQ8 family and TM7a genus), Firmicutes (n = 3, all assigned to Bacillus), Verrucomicrobia (n = 2, all assigned to Luteolibacter) and Bacteroidota (n = 2, all assigned to Flavobacterium) phyla (Fig. 3). These 49 ASVs tended to increase throughout plant development, as the abundance of most of them was significantly higher in the H stage compared with the Pt stage in both geographical regions (e.g., Streptomyces 130, Bacillus 46 and 54, and Lechevalieria 2).

Fig. 3

Variation in the most abundant discriminant ASVs across host plant growth stages. Identification of the ASVs with an average relative abundance ≥ 0.7% within each tomato cultivar sample group and discriminant (according to the Wilcoxon‒Mann‒Whitney test, FDR-corrected p value < 0.1) in at least one of the three comparisons among the tomato growth stages (post-transplant, flowering, harvesting), in at least one sampling region (Campania and Sicily). The taxa are sorted in descending order according to the average fold change (flowering/post-transplant) between the two sampling regions. The Wilcoxon‒Mann‒Whitney test results are shown (left heatmap): red denotes ASVs significantly enriched in the post-transplant stage, green denotes those significantly enriched in the flowering stage, blue denotes those significantly enriched in the harvesting stage, and gray denotes nonsignificant (Ns) comparisons. Prevalence (central heatmap) and average relative abundance (right heatmap) values were calculated by grouping samples with respect to tomato cultivar (Ciliegino and Pietrarossa F1/Maraskino), sampling time, and sampling region

Among the 38 ASVs whose abundance decreased from Pt to F (Fig. 3), most (n = 26) were assigned to the Proteobacteria phylum: 18 to class Alphaproteobacteria (particularly to genus Devosia: 9 ASVs, Shinella: 4 ASVs, Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium: 2 ASVs), 6 to class Gammaproteobacteria (Pseudoxanthomonas Mexicana 2 and 4, Pseudomonas 115, Cellvibrio fibrivorans, Cellvibrio 86, Lysobacter 230) and 2 to Betaproteobacteria (Acidovorax facilis and Acidovorax 47). The least represented phyla were Actinobacteriota (n = 6, including Lechevalieria 3, Promicromonospora 4, Streptomyces 43, Glycomyces algeriensis, Microbacteriaceae 54 and Salinibacterium 1), Planctomycetota, Patescibacteria, Chloroflexi, Bacteroidota and Verrucomicrobiota. These ASVs tended to decrease in abundance throughout plant development, as their abundance was significantly lower in H than in Pt in at least one sampling region, with the exception of Pseudomonas 115.

The core microbiota included 25 ASVs prevalent across several conditions

We identified 25 ASVs that were detected in more than 80% of the samples as core microbiota (Fig. 4, Supplementary Data 4, and Supplementary Data 5). The most represented phylum was Proteobacteria (n = 16), followed by Actinobacteriota (n = 8) and Planctomycetota (n = 1). Focusing on Proteobacteria, only one ASV belonged to class Gammaproteobacteria (i.e., Pseudoxanthomonas 41, average relative abundance of 0.80%), whereas the others were assigned to the class Alphaproteobacteria: three belonged to genus Devosia (including Devosia 78, with an average relative abundance of 0.53%), two to genus Shinella (i.e., Shinella 10: 1.26; Shinella 17: 0.48%), and the remaining ones were scattered across multiple genera. These genera also included several highly abundant ASVs, such as Bosea 34 (1.46%) and Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium 50 (0.82%). For Actinobacteriota, the most abundant ASVs were Streptomyces 43 (3.53%) and Streptomyces 77 (2.39%), followed by Pseudarthrobacter 5 (0.79%). Finally, the only prevalent ASV belonging to the phylum Planctomycetota was assigned to Pir4 lineage 181 (0.36%).

Fig. 4

Identification of the core microbiota of tomato root endophytic prokaryotic communities. The core microbiota was determined by considering the ASVs detected in at least 80% of the samples. Prevalence (left heatmap) and average relative abundance (right heatmap) values were calculated by grouping samples with respect to tomato cultivar (Ciliegino and Pietrarossa F1/Maraskino), sampling time (post-transplant, flowering, harvesting) and sampling region (Campania, Sicily)

Proposal of microbial consortia

We finally defined a microbial consortium (Fig. 5A) of 12 ASVs by selecting the most prevalent and abundant ASVs from the discriminant analysis as well as from the core microbiota. More specifically, 8 ASVs were selected from the core microbiota (prevalence greater than 90% and average relative abundance greater than 0.4% across all samples): Streptomyces 43, Shinella 10, Pseudoxanthomonas 41, Pseudarthrobacter 5, Devosia 78, Sphingobium herbicidovorans 1, Streptomyces 77 and Bosea 34. From the discriminant analysis results, we selected 4 additional ASVs (present in at least 60% of the samples overall and with an average relative abundance greater than 3% in at least one group of samples represented in Fig. 3), and only ASVs that significantly increased in abundance from Pt to F and from Pt to H in at least one sampling region were selected. The resulting ASVs were Streptomyces 77, Bosea 34, Lechevalieria 2, Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium 50, Variovorax paradoxus 1, and Mycobacterium 34. Interestingly, the first two ASVs were also selected from the core microbiota. This defined microbial consortium was thus selected from the ASVs from a discriminant analysis of endophytic bacteria enriched in the flowering and harvesting growth stages (exhibiting high abundance and prevalence) and from the core microbiota. This comprised 7 ASVs belonging to the phylum Proteobacteria and 5 ASVs belonging to the phylum Actinobacteriota. Correlation analysis (Fig. 5B) revealed that the ASVs selected from the discriminant analysis presented strong positive correlations. On the other hand, most of the ASVs selected from the core samples presented slightly positive or neutral correlations (i.e., Streptomyces 43, Shinella 10, Pseudoxanthomonas 41, Pseudarthrobacter 5, and Devosia 78), although they presented some significant variations among the sampling times, particularly an enrichment in Pt.

Fig. 5

Proposal of a microbial consortium for tomato crops. (A), The selected ASVs were identified by choosing the most prevalent (> 90% overall) and abundant (> 0.4% of average relative abundance overall) ASVs from the core microbiota and from the discriminant analysis results (> 60% of prevalence and > 3% of average relative abundance in each group of samples defined by tomato type), excluding from the latter those that decreased in abundance from the post-transplant stage to the flowering and harvesting stages. The membership of the ASVs is shown in the first bar on the left. The Wilcoxon‒Mann‒Whitney test results are shown (left heatmap): red denotes ASVs significantly enriched the post-transplant stage, green denotes those significantly enriched in the flowering stage, blue denotes those significantly enriched in the harvesting stage, and gray denotes nonsignificant (Ns) comparisons. Prevalence (central heatmap) and average relative abundance (right heatmap) values were calculated by grouping the samples with respect to tomato cultivar (Ciliegino and Pietrarossa F1/Maraskino), sampling time (post-transplant, flowering, and harvesting) and sampling region (Campania and Sicily). (B), Correlation matrix highlights positive (in blue) and negative (in red) correlations among the ASVs of the consortium. The * denotes a statistically significant correlation (p value < 0.05). Red and green squares indicate the two main subconsortia identified

Discussion

Assessing the dynamics of plant-associated microbiota across different locations and plant growth stages may lead to powerful insights into microbial community behavior, making it possible to develop precise biotechnologies for crop management. Here, we propose an in silico approach based on the analysis of 16 S amplicon high-throughput sequencing data, which can reveal the composition and behaviors of the bacterial microbiota adapted to the endorhizosphere of tomato. We considered a high number of samples and diverse experimental variables, which include the most common and widespread conditions used for tomato cultivation, to obtain information on the endophytic microbial population that is shared between the different conditions (Fig. 6A). Endorhizospheric bacteria are a subset of the rhizospheric microbiota that manage to pass through the biological filter of the plant roots [14]. It is known that plant root exudates play an important role in microorganism recruitment, and plants can allocate a nonnegligible portion (11%, according to some studies [15, 16]) of photosynthesized fixed carbon to rhizodeposition activity. Microbial soil communities are characterized by high variability, which can depend on soil tillage, host plant genotype, and environmental and edaphic conditions [17]. This makes establishing direct relationships between biotic/abiotic factors and microbiota dynamics challenging.

Fig. 6

Experimental design and processing steps. (A) The sampling plan included a total of 192 samples from two Italian regions, Sicily (n = 96) and Campania (n = 96), with the following experimental variables: three sampling times for each sampling region; post-transplant (n = 32), flowering (n = 32), and harvesting (n = 32); two tomato cultivars per sampling time, Ciliegino (n = 16) and Pietrarossa F1/Maraskino (n = 16); two sampling areas per tomato cultivar, open field (n = 8) and greenhouse (n = 8); and two crop management types per sampling area, organic (n = 4) and conventional (n = 4). (B) Main steps comprising the acquisition of tomato root endophytic microbiota samples, DNA extraction, 16 S rRNA sequencing, data analysis, and the proposal of microbial consortia for crop management (dashed line)

We showed that the host plant growth stage strongly influenced the bacterial endophytic microbiota associated with tomato roots (Fig. 1), which may reflect an association between different mixtures of root exudates produced by the host plant during its growth and the recruitment of the microbiota. Indeed, this hypothesis is well supported by the literature, as during plant development, exudates can vary, resulting in different mixtures [18, 19], and these variations may lead to differential recruitment of root-associated microbiota [20, 21]. The variation in alpha diversity (Fig. 2A-D) clearly revealed that the bacterial endophytic microbiota underwent strong selective pressure during the plant growth stages. Indeed, the observed richness as well as the Shannon index decreased significantly from Pt onward in both sampling regions, suggesting that some microbial taxa tended to become more dominant in the F and H stages. The observed biodiversity reduction was in line with previous findings; for example, some studies [22] reported that the biodiversity of the rhizospheric microbiota associated with Arabidopsis thaliana decreased during plant development, and others [20] showed that different mixtures of root exudates collected from two growth stages of A. thaliana selected distinct microbial communities when applied to soil, suggesting active and differential recruitment of the microbiota by the host plant during its development. In both sampling regions, the Actinobacteriota, Firmicutes and Patescibacteria phyla significantly increased in abundance from the Pt stage to the F stage, whereas other main phyla, such as Proteobacteria, Planctomycetota and Verrucomicrobiota, significantly decreased in abundance (Fig. 2E-F). In the H stage, the general bacterial microbiota structure tended to return to the Pt condition: the abundance of the Proteobacteria, Verrucomicrobiota and Planctomycetota phyla increased significantly again in both regions, whereas that of the Actinobacteriota and Patescibacteria phyla decreased. Therefore, we observed that the bacterial microbiota underwent strong selection in the F stage of tomato plants, while this selective pressure was likely weaker in the later growth stages, probably following the decrease in the plant’s metabolic activities.

We also revealed the most abundant ASVs responsible for the variations in the microbial communities across host plant growth stages. In general, these discriminant ASVs exhibited quite common behavior in the two tomato cultivars and in the two geographical regions (Fig. 3); therefore, some of these taxa may be of interest for further investigations at the species level to develop microbial formulations that are effective across multiple conditions. In general, the abundance values of the ASVs were positively correlated with their prevalence. This is expected; however, it is still not clear if this relationship is universally applicable across all environments and microbial taxa [23]. Many of the ASVs enriched in the F stage belong to the Actinobacteriota phylum (Fig. 3 and consistent with the results presented in Fig. 2E-F), suggesting a key role of these microorganisms in tomato plant development. The Streptomyces genus had five distinct ASVs significantly enriched from the Pt stage to the F stage in at least one sampling region. This genus is well known for having many plant growth-promoting (PGP) activities, such as the production of indole-3-acetic acid (IAA), 1-aminocyclopropane-1-carboxylate deaminase (ACCD) and siderophores, and phosphate-solubilizing activity [24,25,26]. It is also widely recognized that Streptomyces species are efficient biocontrol agents, probably due to the high rate of lytic enzyme and antimicrobial compound production [27], and this genus has already been identified as an endophyte of tomato roots [28]. We also identified the Amycolatopsis genus (n = 4) as enriched in the F stage in at least one sampling region. Amycolatopsis species are known for having PGP activities (siderophores and ACCD production [29]), and this genus has been reported as a tomato root endophyte [30]. Many other ASVs belonging to the Actinobacteriota phylum with significantly increased abundance in the F stage are known for exerting PGP activities. This is the case for the Kribbella genus, with some species described as high-IAA producers and nitrogen-fixing bacteria [31], as well as potential biocontrol agents since they can produce antifungal compounds known as kribellosides [32]. Additionally, the Nocardioides genus has been reported as an IAA producer [33]. Among the other ASVs enriched in F, many were assigned to the Proteobacteria phylum, with several taxa related to PGP traits. More specifically, Variovorax paradoxus has been reported to have interesting characteristics, such as the ability to metabolize many different substrates [34] and the ability to interfere with the quorum sensing signals of several bacterial species [35, 36], in addition to classic PGP traits such as ACCD and siderophore production [37]. The Stenotrophomonas genus is also known for having PGP activities, such as IAA production [38] and saline stress protection [39]. In agreement with the findings shown in Fig. 2E-F, most of the ASVs that were significantly enriched in Pt belonged to the Proteobacteria phylum, with a high number of ASVs associated with the Devosia genus (n = 9) and genera related to the organic matter cycle, such as Cellvibrio and Lysobacter [40, 41]. The genus Lysobacter is also known for its important role as a biocontrol agent [42, 43]. Therefore, our results suggested that the endophytic microbiota probably shifted toward augmented PGP characteristics in the F stage of tomato plants, with the host plant recruiting specific microorganisms during their development.

While the identification of the ASVs discriminant among plant growth stages has suggested transient relationships between the host plant and endophytic microbiota, analysis of the core microbiota allowed us to characterize more long-term associations. We identified a set of ASVs (Fig. 4) that were characterized by high persistence across growth stages and geographical distribution in both sampling regions. Interestingly, Shinella, Streptomyces, and Sphingobium have already been described as core taxa of the endophytic microbiota associated with tomato roots [44], and Devosia, Bosea, and Bradyrhizobium were identified as some of the most abundant taxa in the tomato rhizosphere and endorhizosphere in previous studies [45]. Notably, we identified Pseudarthrobacter, not previously reported as a tomato root endophyte, as one of the most prevalent taxa (detected in 99.5% of the samples) along with its high relative abundance (0.79%).

We ultimately selected a total of 12 ASVs, the most abundant and prevalent of which were from both the core microbiota and the discriminant analysis results. Among the latter, we excluded ASVs with decreased abundance from Pt to F or from Pt to H, since they were probably not selected by the metabolic activities of the host plant during its development (Fig. 5A). Correlation analysis revealed a separation between five ASVs (i.e., Streptomyces 43, Shinella 10, Pseudoxanthomonas 41, Pseudarthrobacter 5, and Devosia 78) and the others, suggesting the possible selection of two subconsortia. These five ASVs (Fig. 5B, red squares) were selected from the core microbiota: they presented slightly positive or neutral correlations among themselves and negative correlations with the other selected ASVs (except for Pseudarthrobacter 5). The other ASVs of the consortium (i.e., Sphingobium herbicidovorans 1, Streptomyces 77, Bosea 34, Lechevalieria 2, Allorhizobium-Neorhizobium-Pararhizobium-Rhizobium 50, Variovorax paradoxus 1, and Mycobacterium 34) showed strong positive correlations among themselves (Fig. 5B, green squares), as they were selected from the discriminant analysis results (except for S. herbicidovorans), meaning that they significantly increased in abundance from Pt to the F and H stages in at least one geographical region. Therefore, we suggest a first consortium with the role of establishing long-term relationships with host plants and thus ensuring phytosanitary stability and a second consortium that instead may be better-suited to promote tomato growth at later stages, as these ASVs probably exerted some PGP activities specifically recruited from the host plant during its development. Interestingly, most of the selected ASVs reported in the literature have PGP activities, such as those of the Streptomyces genus and Variovorax paradoxus, as discussed above [24, 37], in addition to Shinella (IAA and hydrogen cyanide production [46]), Pseudoxanthomonas, Pseudarthrobacter, Mycobacterium (IAA and ACCD production, phosphate solubilization, and nitrogen fixation [47,48,49,50]), and Devosia (nitrogen fixation [51]). Other selected ASVs (i.e., Sphingobium herbicidovorans, Bosea or Lechevalieria) are less known as PGP candidates, although Lechevalieria is known for antibiotic production [52]. Finally, the selected ASVs belonging to the Allorhizobium-Pararhizobium-Neorhizobium-Rhizobium group could not be assigned to a specific genus according to the molecular data used in this study. Furthermore, several ASVs selected in the consortium have already been reported in the literature as tomato root endophytes, such as Streptomyces, Shinella, Devosia, and Bosea, as well as Pseudoxanthomonas, Sphingobium herbicidovorans, and Mycobacterium [44, 45, 53,54,55]. Conversely, other selected ASVs, such as Variovorax paradoxus, Pseudarthrobacter, and Lechevalieria, are likely reported here for the first time as important components of the endophytic microbiota of tomato roots. Since these ASVs showed interesting behavior in the microbial community in terms of persistence at high abundance during plant development and across different tomato cultivars and geographical regions or enrichment in the F stage, suggesting recruitment by the host plant, further investigations into the role of these ASVs in the holobiont and their possible PGP activities are needed.

Conclusions

This research provides an in-depth characterization of the tomato root endophytic bacterial microbiota, revealing significant temporal shifts during plant growth stages (post-transplant, flowering, and harvesting). These shifts are associated with changes in bacterial community structure, with Actinobacteriota and Proteobacteria being the dominant phyla. Therefore, it appears that host‒microorganism associations are subject to precise mechanisms of coevolution that drive the assembly of the plant-associated microbiota rather than other environmental factors. The highly representative framework of the tomato root endophytic microbiota determined by discriminating taxa between sampling times and the core microbiota provides insight into Streptomyces, Shinella, Devosia, and Pseudoxanthomonas, as well as the lesser known genera Variovorax, Pseudarthrobacter, and Lechevalieria. These persistent microbial taxa, the most important components of the tomato root endophytic microbiota, suggest long-term host‒microbe associations. These findings provide a basis for developing tailored microbial formulations that may increase crop resilience, reduce dependence on synthetic agricultural inputs, and improve overall agricultural sustainability. Further studies are needed to validate the proposed microbial consortia under real agricultural conditions and assess their effectiveness in promoting plant growth and soil health. However, given the huge biodiversity of bacterial populations present in the soil rhizosphere, and the time and resources required for microbial isolation and selection using culture-dependent methods, understanding specific key persistent bacterial genera could streamline the process of developing a commercial microbial-based biostimulant.

Methods

Experimental design and soil measurements

The experimental design is summarized in Fig. 6A, and sample information is reported in Supplementary Data 1. In detail, the sampling plan included a total of 192 samples from two Italian regions: Sicily (Vittoria, RG; 36°57′ N 14°32′ W; n = 96) and Campania (Parete, CE; 40°57′35.57″ N 14°09′53.87″ W; n = 96). For each sampling region, three sampling times were considered based on known changes in the level of phytohormones that occur during the transition to the different growth stages: post-transplant (n = 32), flowering (n = 32), and harvesting (n = 32). Two main tomato types were considered per sampling time: the small-rounded fruited type Solanum lycopersicum var. cerasiforme, commonly called Ciliegino, a putative feral hybrid between wild currant-type and domesticated garden tomatoes (n = 16) and the elongated type S. lycopersicum cv. Pietrarossa F1 (Campania region) and Maraskino (Sicily region), domesticated garden tomatoes (n = 16). For each tomato type, two crop treatments were carried out: organic (cattle manure at 29 t ha^− 1; n = 8) and conventional (mineral fertilizer as ammonium nitrate at 26% N at a dose of 440 kg ha^− 1; n = 8). Finally, for each treatment, the plants were grown in two agricultural areas: a greenhouse (tunnel polyethylene greenhouse with a proper ventilation and dehumidification systems to regulate moisture levels and with a drip irrigation system; n = 4) and an open field (n = 4). The tomato seedlings used in this study were sourced from local markets in the Sicily and Campania regions. Before tomato transplantation, soil sampling was carried out at 0–20 cm for physical and chemical characterization. The soil of the Campania region was a clay sandy loam soil (USDA classification) with a pH of 7.4, an EC of 4.7 dS m^− 1, total nitrogen content of 1.6 g kg^− 1, organic matter content of 22.9 g kg^− 1, phosphorus content of 127.3 ppm, potassium content of 382.3 ppm, calcium content of 2196 ppm, magnesium content of 102.3 ppm, sodium content of 13.3 ppm, organic carbon content of 13.3 and a C/N ratio of 8.3. The soil of the Sicily region was a clay sandy loam soil (USDA classification) with a pH of 7.3, an EC of 4.9 dS m^− 1, total nitrogen content of 1.5 g kg^− 1, organic matter content of 21.1 g kg^− 1, phosphorus content of 133.3 ppm, potassium content of 322.3 ppm, calcium content of 2068 ppm, magnesium content of 116.7 ppm, sodium content of 147.7 ppm, organic carbon content of 12.2 and a C/N ratio of 8.1.

Sampling and DNA extraction

Rhizosphere sampling of the tomato plants was performed as previously described [56]. To obtain the endophytic bacterial microbiota, representative subsections of root tissue were placed directly into 50 mL of sterile chilled water (4 °C) to facilitate deep washing. The samples were shaken vigorously for one minute to eliminate the ectorhizosphere and centrifuged at 6500 rpm for 10 min at room temperature. The washing step was repeated until visible soil on the root surface was eliminated. The roots were subsequently sterilized with 70% (v/v) ethanol for 10 min, washed with sterile distilled water, frozen in liquid nitrogen and crushed in a mortar to obtain a homogeneous powder. The absence of bacterial and fungal epiphytes was assessed by placing the surface-sterilized roots on plate count agar (PCA; Oxoid, Milan, Italy) and potato dextrose agar (PDA; Oxoid, Milan, Italy) plates at 28 °C for 3‒4 days. Total genomic DNA was extracted using a FastDNA SPIN Kit for Soil (MP Biomedicals, Illkirch Cedex, France) according to the manufacturer’s instructions and stored at -80 °C until analysis. The main steps comprising DNA extraction, 16 S rRNA sequencing and data analysis are summarized in Fig. 6B.

Illumina miseq sequencing and Raw data processing

The primers S-D-Bact-0341-b-S-17 (5’-CCTACGGGNGGCWGCAG-3’) and S-D-Bact-0785-a-A-21 (5’-GACTACHVGGGTATCTAATCC-3’) were used to amplify the V3-V4 hypervariable region of the 16 S rRNA gene. The amplicon libraries were sequenced on a MiSeq platform (Illumina, USA), and 250-bp paired-end reads were obtained. The primers were removed from the raw reads using Cutadapt [57] with the command-line option -g to detect nonanchored and partial primer sequences. Preliminary quality control was performed with FastQC [58]. Reads were further processed through the DADA2 pipeline [59], which included the filtering of low-quality reads through the filterAndTrim function with the following parameters: the maxN = 0, maxEE = c(2,2) and truncQ = 2 options were used, and the recommended settings were selected, although the truncLen option was omitted since the reads were of good quality at all lengths. After the paired-end reads were merged, a total of 7,695,390 sequences were generated, with an average of 40,080 merged sequences per sample. Chimeras were also removed via the removeBimeraDenovo option, and a total of 6,128,328 sequences were retained (with an average of 31,918 per sample). Taxonomic profiling and assignment were performed using the SILVA database [60] (v138.1).

Statistical analyses

Downstream statistical analyses were performed in the R environment [61] (v4.3.3). The sample biodiversity was covered by sequencing depth (Supplementary Fig. 1). The original dataset comprised a total of 49,924 distinct ASVs. We removed ASVs that were not assigned at the kingdom level or that were assigned to Eukaryota, Mitochondria, or Chloroplast. This resulted in a total of 49,085 ASVs, which were further filtered to 19,733 ASVs by removing the ASVs with fewer than 5 occurrences overall. Beta diversity in terms of MDS was computed through the vegan package [62] and the Bray‒Curtis distance. This MDS plot and the bar plots were represented using the ggplot package [63]. Alpha diversity in terms of the Shannon index and observed richness was calculated using the microbiome package [64] and represented with the ggpubr package [65]. Statistical tests between groups of samples in terms of alpha diversity were performed using the Kruskal‒Wallis test, which uses the function available in the ggpubr package. The ASVs statistically discriminant among host plant growth stages were calculated by removing ASVs with an average relative abundance < 0.005% overall. This resulted in 1812 ASVs. Discriminant ASVs were identified through the Wilcoxon‒Mann‒Whitney test; multiple hypothesis correction was performed using the false discovery rate; only the ASVs with an adjusted p value < 0.1 and an average relative abundance > 0.7% in at least one group of samples are reported in the heatmap in Fig. 3. The discriminant ASVs were ordered in the heatmap on the basis of the average fold change between the two sampling regions, which was calculated as follows:

Fold change = log2(Average abundance of the ASV in F + c / Average abundance of the ASV in Pt + c), where the average abundance of the ASV is the average computed between the two tomato cultivars at a given sampling time and in a sampling region and c is a constant equal to 2.2e-16. The fold change was calculated separately for the Campania and Sicily regions, and the average was then calculated as the final fold change value. The core microbiota was calculated by considering the ASVs detected in at least 80% of the samples without defining a threshold in terms of relative abundance. Heatmaps were generated using the ComplexHeatmap package [66]. The microbial consortium was selected by defining two criteria, resulting in 12 ASVs: (i) 6 ASVs were selected from the discriminant analysis results by choosing those that increased in abundance from Pt to F and H and exhibited high abundance (> 3% average relative abundance in at least one group of samples) and prevalence (> 60% overall); and (ii) 6 additional ASVs were selected from the core microbiota by choosing the most prevalent (> 90% overall) and abundant (> 0.4% average relative abundance overall) ones. The correlation matrix was computed using the function cor (Spearman method), and p values were calculated using the function rcorr.