Temporal dynamics of the tomato rhizosphere microbiome in response to synthetic communities of plant growth-promoting rhizobacteria

The quality and productivity of agricultural crops are strongly influenced by soil health and the efficiency of interactions in the rhizosphere, where plant roots host highly dynamic microbial communities¹. In this context, the root microbiome is not just a passive bystander, but actively contributes to plant metabolism, regulating nutrient uptake, growth, stress response and pathogen tolerance^2,3. However, the growing pressure of intensive agriculture has compromised soil microbial biodiversity, reducing its ecological functionality and increasing dependence on chemical fertilisers and pesticides^4,5.

Plant growth-promoting rhizobacteria (PGPR) have emerged as key biotechnological tools for restoring the microbiological balance of the rhizosphere and improving plant performance under unfavourable conditions^2,6. These microorganisms act through multiple mechanisms: release of phytohormones (auxins, cytokinins), production of siderophores and solubilising enzymes, induction of systemic resistance and competition with pathogens⁶. Their actions can generate systemic effects on plant physiology and the entire rhizosphere ecosystem³.

The application of PGPR in agriculture has shown promising results, both through the use of single strains and microbial consortia^7,8. While single-strain inoculants offer significant agronomic benefits, their effectiveness is often hindered by limited colonization capacity and competition with the resident microbiome^7,9. To address these challenges, research has shifted toward multispecies microbial consortia, or Synthetic Communities (SynComs), which provide greater functional resilience and adaptability^4,10. However, results vary depending on the experimental context and agroecological conditions⁷.

Recent evidence suggests that the functional impact of a consortium does not depend solely on the long-term persistence of the inoculated strains or its genetic detectability in the microbiome^10,11. Even in the absence of clear molecular traceability, consortia can exert indirect and systemic effects on the resident microbiome, contributing to the reorganisation of microbial communities and the modulation of soil ecological functions, mainly mediated by competition, modulation of trophic networks and stimulation of secondary metabolites^8,10,11.

Tomato (Solanum lycopersicum L.) is a leading global vegetable crop, with an annual production exceeding 180 million tons¹². Beyond its relevance as a dietary cornerstone, tomato represents one of the primary model species in plant science and intensive horticulture, particularly in greenhouse systems where yield and fruit quality are maximized under controlled conditions¹³. Accordingly, tomato has become a major target for the development of microbial bioinoculants, and numerous studies have demonstrated the potential of plant growth–promoting microorganisms to enhance plant growth, nutrient acquisition and tolerance to both biotic and abiotic stresses (reviewed in Adedayo et al.; Benaissa et al.)^14,15. In line with trends observed in other crops, tomato research has increasingly moved toward the application of microbial consortia to achieve more robust and reliable plant growth–promoting and stress-mitigating effects¹⁵. Bacterial consortia have been proved to improve growth parameters^16,17,18 and yield¹⁹ of tomato crop. They have also shown the ability to mitigate the adverse effects induced by abiotic stresses such as water^20,21, salt²² , high temperature²³ and chilling stress²⁴. In addition, consortia of bacteria or bacteria and fungi have been shown to reduce the detrimental effect of biotic stressors, reducing the incidence and severity of bacterial wilt²⁵, grey mold¹⁷, Sclerotinia stem rot¹⁸, Fusarium wilt²⁶, nematodes and insect damages²⁷.

Despite increasing interest, only a few studies have investigated the impact of microbial consortia on the resident root-associated microbiome of tomato, and none have focused on consortia composed exclusively of native bacterial isolates^10,28,29. Furthermore, previous research has predominantly characterized the bacterial communities of the rhizosphere, while overlooking fungal communities’ response to bacterial inoculation^10,28,29.

The application of microbial consortia on tomato rhizosphere has been shown to restructure resident bacterial communities by modulating alpha- and beta-diversity through selective effects^10,28,29. Beyond the general trends, it has been observed that microbial responses are highly contextual. Factors such as inoculum composition, experimental settings, and the phenological stage of the plant, ranging from flowering to fruit setting, induce unique ‘microbial signatures’^10,28,29.

Specifically, the temporal behaviour of microbial consortia has not yet been systematically examined. However, among the effects reported in previous studies, initial inoculation can induce shifts in community structure, probably due to long term functional adjustments even after inoculant decline^10,28. Central to this successional response appears to be the “rare biosphere”, low-abundance microbial taxa that act as a vast genetic reservoir. Although these taxa are frequently overlooked, they can be remarkably active in key biogeochemical processes³⁰.

In this study, we aimed to contribute to the understanding of how plant growth–promoting influence the plant microbiome by using three small SynComs (MIX1, MIX2 and MIX3) composed of taxonomically diverse tomato-native bacterial strains isolated from the tomato endosphere. These strains were selected through a top-down, microbiome-informed approach based on the tomato core microbiome, following a microbiome study conducted across the cultivation chain, from seed to nursery to commercial greenhouse^31,32. The selected set comprised strains exhibiting plant growth–promoting and/or biocontrol activity in tomato when tested individually³². These included genera of well-characterized PGPR and BCAs, such as Pseudomonas and Bacillus, which displayed the most pronounced effects, as well as strains from less-studied genera, Leclercia, Chryseobacterium, Glutamicibacter, and Paenarthrobacter, that also contributed measurable benefits³². By combining strains with distinct genomic potentials, we aimed to investigate the complex ‘frenemy’ dynamics and synergistic outcomes of these interactions³³. The three SynComs, comprising four, six, and ten bacterial strains, consistently promoted plant growth in both determinate (Pizzutello) and indeterminate (Proxy) tomato varieties, with the SynComs of six and ten strains—containing pseudomonads—exhibiting the most pronounced effect. Our results showed that temporal factors were the primary drivers of rhizosphere microbial community composition in the ‘Proxy’ variety. Within this overarching temporal framework, SynCom application elicited a distinct, time-dependent response. Initially, bacterial communities underwent treatment-specific shifts before later converging, while fungal communities experienced a progressive, time-modulated restructuring. At the taxonomic and functional level, bacterial responses to the SynComs were both shared and consortium-specific, predominantly involving the ‘rare biosphere’. Notably, predictive functional profiling (PICRUSt) suggested that rhizospheres treated with Pseudomonas-containing consortia converged toward similar functional states, despite the inherent limitations of the inference method.

Analysis of the in vitro compatibility of ten PGPR strains

In this study we used ten genome-sequenced bacterial strains that had previously shown both plant growth promotion and biocontrol activities in tomato plants³². Strain identification, source and GenBank genome sequence ID are described in Table 1. Strains’ reciprocal growth inhibition was evaluated in vitro by cross-streak method, revealing some antagonistic activity on PDA but not on LBA or TSA (Supplementary Table S1). Notably, Pseudomonas strains inhibited the growth of all other bacterial strains, regardless of taxon, and partially inhibited each other as well (Supplementary Table S1). The only exception was recorded for Leclercia sp. S52 that was able to grow together with P. simiae POE78A. B. velezensis PFE11 showed inhibition activity against the other B. velezensis strains (PFE42 and PSE31B). Zones of inhibition were also detected when cross-streaking Leclercia sp. S52 against Paenarthrobacter strains (S54 and S56) and Chryseobacterium sp. POE47 (Supplementary TableS1).

The trials were carried out assembling the strains in three synthetic communities (SynComs) of four, six and ten members. In particular, MIX1 comprised four compatible strains (two Bacillus, PSE31B and PFE42, Glutamicibacter halophytocola PFE44 and Leclercia sp. S52); MIX2 was composed of the MIX1 strains supplemented with two Pseudomonas strains (P. salmasensis POE54 and P. simiae POE78A); MIX3 contained all ten strains, adding B. velezensis PFE11, Paenarthrobacter ureafaciens S54, Paenartrhobacter sp. S56 and Chryseobacterium sp. POE47 to the previous consortia (Table 1). Figure S1 shows the co-growth of SynCom strains when cultured together on different substrates (Supplementary Figure S1).

Effect of SynCom treatment on tomato growth

Tomato seedlings ‘Pizzutello’ and ‘Proxy’ were treated by soil drenching (T0) with the SynComs consisting of four (MIX1), six (MIX2), or ten (MIX3) strains (Table 1). The treatments resulted in a significant increase in plant height, compared to a control drenched with water, four weeks after the treatment (Fig. 1). The increase was more pronounced in the indeterminate-growth variety ‘Proxy’, than in the determinate bush-type ‘Pizzutello’. The growth-promoting effects of the SynComs were also impacted the fresh and dry weight of the shoots, which were significantly higher four weeks after treatments (p < 0.0001 and p = 0.001, respectively) (Fig. 2A–B). The effect on root biomass was statistically significant only in ‘Pizzutello’ seedlings treated with MIX1 and MIX2, with a limited effect observed in ‘Proxy’ seedlings treated with MIX2 (Fig. 2A–B). When considering whole-plant data, all the treatments increased the fresh and dry weight of ‘Pizzutello’ plants (p < 0.0001), while in ‘Proxy’ plants MIX2 and MIX3 significantly differed from the control (p < 0.0001 and p = 0.001, for fresh and dry weight respectively) (Fig. 2A–B). In addition, the height of ‘Proxy’ seedlings was monitored weekly during the sampling for the study of the rhizosphere microbial communities. After one week (T1), all SynCom-treated seedlings showed a significant increase in height compared to control plants drenched with tap water (Supplementary Figure S2). Seedlings treated with MIX2 and MIX3 exhibited a significant increase in height at all time points compared to both control plants and those treated with MIX1 (p < 0.0001) (Supplementary Figure S2). Over the four-week period (T4), plants treated with MIX1, MIX2, and MIX3 exhibited a 39%, 94%, and 84% increase in height, respectively, compared to the water-treated control seedlings (Supplementary Figure S2).

Effect on resident rhizosphere bacterial and fungal communities

The impact of SynCom treatments on the rhizosphere microbial communities of tomato ‘Proxy’ plants was assessed compared to water-treated control plants at four distinct time points: after the bacterial suspension was applied to the soil (T0), and subsequently one week (T1), two weeks (T2), and four weeks (T4) post-treatment.

Illumina sequencing of the bacterial 16S rRNA gene and the fungal ITS region produced 39,067,796 and 35,164,584 reads, respectively. After paired-end alignments, quality filtering, and deletion of chimeras and singletons, a total of 14,242,196 bacterial 16S reads and 7,623,912 fungal ITS reads were generated from 47 and 48 samples, respectively (one sample was removed from bacterial analysis due to the low number of reads), and assigned to 32,201 bacterial ASVs and 15,182 fungal ASVs.

Two α-diversity indices, Chao1 and Shannon, were calculated for both bacterial and fungal communities across all treatments and time points. Overall, no significant differences in α-diversity were observed between microbial treatments at any of the sampling points, for either bacterial or fungal rhizosphere communities (Fig. 3A–B; Supplementary Figure S3A-B). However, bacterial communities in control samples exhibited a higher Chao1 index compared to treated samples at all time points, except the final one (T4), where the control samples showed the lowest Chao1 value. At this time point, MIX1 and MIX2 treatments displayed the lowest Shannon index (Fig. 3A–B).

For β-diversity, Principal Coordinate Analyses (PCoA) based on Bray–Curtis dissimilarities indicated that time was a major driver of bacterial community diversity (p value = 0.001: MIX1, R² = 0.44; MIX2, R² = 0.35; MIX3, R² = 0.42; PERMANOVA, 999 permutations), regardless of treatment (Fig. 3D–F). At time points T1 and T2, the bacterial communities in the rhizosphere of plants treated with MIX2, as well as at T2 for those treated with MIX3, exhibited statistically significant differences compared to the communities in control plants (interaction time x treatment, p value = 0.04, R² = 0.12; PERMANOVA, 999 permutations). Regarding fungal communities, only in plants treated with MIX3 a significant influence of time was found (p-value = 0.025, R² = 0.27; PERMANOVA, 999 permutations) while in the other cases, neither time nor treatment appeared to significantly influence fungal community diversity (Supplementary Figure S3D-F).

The relative abundance of bacterial communities at the phylum level reflected the influence of time, as shifts in abundance were observed over the course of the study, despite the overall phylum composition remaining similar. From T0 to T4, a general increase in Verrucomicrobiota and Planctomycetota and a decrease in Proteobacteria, Actinobacteriota, and Bacteroidota was observed across all treatments (Fig. 3C). Linear correlation analysis of phylum-level relative abundances from T0 to T4 confirmed a statistically significant increase in Planctomycetota (Linear Regression test: p < 0.05, R² > 0.80; Supplementary Table S2) in all treatments (Supplementary Figure S4). Likewise, also Firmicutes and Verrucomicrobiota phyla showed an increase across treatments, although only in MIX1 (LR: p < 0.05, R² > 0.53) and MIX2 (LR: p < 0.05, R² > 0.53) were statistically significant, respectively. Meanwhile, a decrease in Proteobacteria, Actinobacteriota, and Bacteroidota was observed across all treatments (LR: p < 0.05, R² > 0.51; Supplementary Table S2). Regarding Patescibacteria phylum, a high trend related to the decrease of its relative abundance in control (LR: R² = 0.66) and MIX3-treated (LR: R² = 0.25) plants and increase in MIX1 and MIX2 treatments (LR: R² > 0.50) was assessed as well, despite statistically significant only in control plants (p < 0.05) (Supplementary Table S2; Supplementary Figure S4).

Looking at the different time points, at T1 the relative abundance of Proteobacteria and Actinobacteriota was higher in the rhizosphere of plants treated with MIX2 than in control plants (29.7% vs 23.4% and 12.4% vs 10.4%, respectively). In addition, the relative abundance of Planctomycetota and Myxococcota was lower in MIX2 compared to the control (12% vs 15% and 1.9% vs 3.7%, respectively) (Fig. 3C). Patescibacteria abundance was higher in the rhizosphere of SynCom-treated plants at T2 (MIX1, 2.4%; MIX2, 3.1%, MIX3, 3.2%) and T4 (MIX1, 5.0%; MIX2, 7.8%, MIX3, 3.6%) compared to the control (1.6%), while Verrucomicrobiota abundance was lower at T2 (MIX1, 19.0%; MIX2, 20.1%; MIX3, 21.7% vs 23%) but higher at T4 (MIX1, 24.3%; MIX2 22.7%; MIX3, 21.0% vs 20.4%) (Fig. 3C). In addition, at T4, the relative abundance of Planctomycetota was higher in MIX3 compared to the control (26% vs 21.9%).

Fungal communities of tomato rhizosphere, at phylum level, were dominated by Ascomycota, followed by Basidiomycota, unidentified Fungi and Mucoromycota (Supplementary Figure S3C). The relative abundance of Basidiomicota increased while the abundance of Ascomycota and Mucoromycota decreased over time (Supplementary Figure S3C). Linear correlation analysis of phylum-level relative abundances from T0 to T4 confirmed that the shift of the most abundant fungal phyla was influenced by time (Supplementary Table S3; Supplementary Figure S5). The relative abundance of Basidiomycota in the rhizosphere of SynCom-treated plants increased over time in all treatments, although not significantly, whereas in the control plants it showed the opposite trend (Supplementary Table S3; Supplementary Figure S5). Conversely, the relative abundances of Ascomycota declined over time, although significantly only in MIX1 (LR: p < 0.05, R² > 0.85). The same trend was observed for Mucoromycota (not significant) and Mortierellomycota, with a significant decline observed in MIX3 (LR: p < 0.05, R² > 0.88) (Supplementary Table S3; Supplementary Figure S5). In MIX3, a significant increase of Rozellomycota was also observed (LR: p < 0.05, R² > 0.58) (Supplementary Table S3; Supplementary Figure S5).

Comparing the treatments at the different time points, MIX2 and MIX3 at T2 showed a lower abundance of Ascomycota compared to the control (77.7% and 80.6% vs 84.4%), while in turn Basidiomycota (14% and 10.6% vs 6.3%) and Mucoromycota (4.0% and 5.1% vs 1.8%) were more abundant (Supplementary Figure S3C). At T4 all the plants treated with the consortia had a higher abundance of Basidiomycota (MIX1, 11.6%; MIX2, 10.3%; MIX3, 11.3%) than the control (6.4%). In addition, in MIX2 a higher relative abundance of Mucoromycota (2.6%) was observed compared to the rhizosphere of control (1.7%) (Supplementary Figure S3C). At the genus level, all samples were characterized by the high relative abundances (from 46 to 65%) of Pseudogymnoascus (Ascomycota) (data not shown).

Treatments with SynComs affected bacterial communities more than fungal communities

Differential abundance analysis was performed to compare the bacterial and fungal genera of the rhizosphere of SynCom-treated plants with the rhizosphere microbial communities of control plants at each time point.

The analysis revealed significant alterations in the composition of tomato rhizosphere bacterial communities following the treatment with synthetic consortia. The number of differentially abundant taxa compared to the untreated control varied both over time and across treatments (Fig. 4).

At T1 (one-week post-treatment), the highest overall impact was observed, with 136 significant (p value < 0.05) differentially abundant genera identified (Fig. 4A). The MIX2 consortium induced the broadest response with 53 unique taxa, followed by MIX1 (37 taxa) and MIX3 (21 taxa) (Fig. 4C). Seven taxa were shared across all three treatments (Fig. 4C). Fold-change analysis revealed that among these core taxa, Glutamicibacter, Bacillus, Azospirillum, and MM2 were consistently enriched, whereas Demequina, Phenylobacterium, and Nitrobacter were suppressed across all treatments (Supplementary Figure S6A). Pairwise analyses further identified specific shared signatures: MIX1 and MIX2 shared 12 taxa, of which 9 enriched, including several involved in the sulphur cycle such as Desulfosporosinus, Sulfurovum, and Woeseia, while Niveispirillium was depleted (Supplementary Figure S6A). MIX1 and MIX3 shared five taxa, including the nitrogen-fixer Paramesorhizobium, which was significantly enriched. MIX2 and MIX3 shared only the taxon Ahniella, which was suppressed in both (Supplementary Figure S6A). Additionally, MIX2 showed strong specific enrichments of Robiginitalea, Tyzzerella, Syntrophus, Candidatus Competibacter and Solibacillus. MIX1 specifically recruited anaerobic taxa such as Clostridium sensu stricto 10 and Anaerocolumna, along with Apilactobacillus, while MIX3 exhibited a unique enrichment of Piscinibacter and a depletion of Snodgrassella (Supplementary Figure S6A).

At T2 (two weeks post-treatment), the total number of significant differentially abundant taxa decreased to 83 (Fig. 4A). Three taxa were common to all treatments: Solimonas (enriched), and Cellvibrio and Tistrella (both depleted) (Fig. 4C; Supplementary Figure S6B). MIX1 showed the highest number of differentially abundant genera (27), with strong enrichments of Lactobacillus and Methanosaeta (Supplementary Figure S6B). MIX2 and MIX3 differential abundant taxa accounted for 22 and 19 genera, respectively (Fig. 4C). Shared signatures at T2 demonstrated high consistency: all 4 taxa shared by MIX1/MIX2 and all 3 taxa shared by MIX2/MIX3 showed congruent patterns; among the 5 taxa shared by MIX1/MIX3, only one showed an opposing trend (Supplementary Figure S6B).

By T4 (four weeks post-treatment), 108 significant differentially abundant genera were identified (Fig. 4A). An increase in shared responses was observed, with 18 DA ASVs common to all treatments (Fig. 4C). Fold-change analysis revealed that 14 of these 18 shared taxa were depleted, with particularly marked suppression of Thermoanaerobaculum and Methylobrevis (Fig. 5). Treatment-specific taxa accounted for 8 genera for MIX1 and 22 genera for both MIX2 and MIX3 (Fig. 4C). MIX1 and MIX2 exhibited the highest similarity at T4, sharing 26 taxa with congruent patterns. Among them, eight genera showed enrichment, including Glutamicibacter and 18 genera were suppressed, including Desulfobacter, Sulfurovum and Woeseia (Fig. 5). MIX1 and MIX3 shared 5 taxa with consistent patterns. MIX2 and MIX3 shared 7 taxa with congruent trends, including a strong suppression of Snodgrassella (Fig. 5).

Overall, most of the differentially abundant genera in the rhizosphere of SynCom-treated plants, compared to the control, were rare taxa, with very low relative abundance in the rhizosphere microbiome (Fig. 5). At T4, the only exception was the genus Gemmata, which stood out as one of the most abundant bacterial taxa across all treatments and was significantly enriched in the rhizosphere of plants treated with MIX2 and MIX3 (Fig. 5).

Statistically significant differentially abundant genera were detected only at two and four weeks after SynCom treatments and their number was lower than that observed for bacteria (Fig. 4B). At T2, the number of depleted genera was higher than enriched ones, while the opposite trend was observed at T4, with the rhizosphere of plants treated with MIX2 showing a higher number of differentially abundant genera (19) than the other treatments, compared to the rhizosphere of control plants (Fig. 4D). At T2, all SynCom treatments showed a significant enrichment of Mucor (Mucoromycota) and depleted Cladosporium (Ascomycota) in the rhizosphere (Supplementary Figure S7A). A few genera showed treatment-specific depletion or enrichment (Supplementary Figure S7A).

At T4, a consistent enrichment of several Basidiomycota genera was observed in the rhizosphere of SynCom-treated plants, including, Clitopilus, across all treatments, and Apiotrichum and Myriococcum in MIX2 and MIX3 compared with the control (Fig. 4D; Supplementary Figure S7B). Only a few additional fungal genera showed significant differential abundance at this time point. Pseudochaetosphaeronema was depleted in the rhizosphere of plants treated with MIX1 and MIX2, while Pseudophialocephala was enriched in plants treated with MIX1 and MIX3 (Supplementary Figure S7B). The rhizosphere of plants treated with MIX1 also displayed higher Saccharomyces and lower Gomphillus (consistently with T2) (Supplementary Figure S7A-B). MIX2 showed the largest number of unique differentially abundant taxa, including Smittium, which exhibited the highest log₂FoldChange (Supplementary Figure S7B). In MIX3 Cladosporium remained depleted, confirming the trend already detected at T2 (Supplementary Figure S7A-B).

Predicted biological functions of rhizosphere bacterial communities

To investigate the potential functional characteristics of the rhizosphere bacterial communities of control and SynCom-treated plants, we used the PICRUSt2 tool. It was predicted that the bacterial communities participate in various functions, categorized into four major groups encompassing 23 functional subcategories, most of them associated with metabolic pathways (approximately 71% in relation to the total genes) (Supplementary Figure S8). Predicted genes related to carbohydrate, amino acid and energy metabolism were the most abundant in all treatments, representing approximately 16%, 14% and 9% of the total genes, respectively (Supplementary Figure S8). A comparative analysis of the predicted KEGG pathways between the rhizosphere bacterial communities of control and SynCom-treated plants was performed, highlighting significant differences (p value < 0.05, FDR) (Fig. 6). The predicted genes of the bacterial communities in the rhizosphere of plants treated with MIX1 were overall enriched compared to control plants, while those ones of plants treated with MIX2 were mostly depleted (Supplementary Figure S9). In particular, the depleted predicted gene families in the metabolic pathways were involved in terpenoid and polyketide metabolism as well as metabolism of carbohydrates and amino acids (Fig. 6). In turn, enrichment of predicted genes associated with xenobiotic biodegradation and metabolism was found (Fig. 6). The opposite trend was observed in the rhizosphere bacterial communities of plants treated with MIX1 (Fig. 6).

Targeting SynCom strains in the microbiome

The 16S rRNA gene sequences of the inoculated strains in the SynComs were aligned against the 16S rRNA gene amplicon-based metagenomic data of rhizosphere samples. Sequences with ≥ 99.5% similarity were designated targeted ASV sequences according to Hu et al.¹⁰. Eight ASVs with the highest match with the inoculated strains were found and their absolute abundance was employed as a metric to determine the over-time strains’ abundance in the rhizosphere microbiome (Supplementary Table S4). The three Bacillus strains, all belonging to the species B. velezensis ³², showed the highest sequence similarity with ASV_407 (Supplementary Table S4).

The absolute abundance of the ASVs belonging to the same genera as the inoculated strains decreased over time in the rhizosphere of both control and SynCom-treated plants (Fig. 7). On the whole, we were able to detect most of the ASVs of our inoculated strains up until the end of the trial, although to a different extent depending on the strain (Fig. 7; Supplementary Table S4). In particular, ASV_2430, matching the inoculated strain Leclercia sp. S52, and included in all SynComs, was not or was detected only at very low levels at the end (Fig. 7; Supplementary Table S4). The same occurred for ASV_2685 (matching the strain Chryseobacterium sp. POE47) and ASV_6620 (matching the strain P. simiae POE78A), inoculated in the MIX3 and in both MIX2 and MIX3, respectively (Fig. 7; Supplementary Table S4).

Overall, the ASVs of the inoculated strains were detected consistently within the treatments although, in some cases, some of the ASVs putatively matching the inoculated strains (i.e. Paenarthrobacter ASV_1004 and Bacillus ASV_407) were also detected in the rhizosphere of the untreated control plants, as well as in some treatments in which they were not inoculated (Fig. 7; Supplementary Table S4).

Three synthetic bacterial communities (SynComs), comprising four, six, and ten native bacterial strains isolated from tomato endosphere (seed and root), denominated MIX1, MIX2 and MIX3 respectively, consistently promoted plant growth in both determinate (‘Pizzutello’) and indeterminate (‘Proxy’) tomato varieties. The SynComs were assembled using strains identified as part of the core microbiome of tomato through a longitudinal study tracking microbiome assembly from seeds to root compartments throughout nursery cultivation and post-transplantation^31,32. Among the three SynComs tested, all of which included Bacillus strains, only MIX2 and MIX3—which also harboured Pseudomonas species—significantly promoted seedling growth. MIX1, which lacked pseudomonads, did not produce comparable plant-growth benefits.

Rhizosphere microbiome dynamics were analyzed in the ‘Proxy’ variety, where temporal factors were the primary drivers of community composition and SynCom treatments induced detectable shifts in resident communities, with the strongest differential abundance changes primarily affecting low-abundance bacterial taxa. The strongest microbiome response was detected one week after treatment and was characterized by the enrichment of multiple bacterial taxa. This initial response was followed by intermediate community changes at the second sampling time. At the final time point, four weeks after SynCom application, community responses converged across treatments, with a higher number of shared responses among the three SynComs. These responses were predominantly characterized by taxon depletion, impacting both taxa shared across all treatments or between pairs of SynComs and those uniquely responsive to individual consortia, with the depletion largely confined to rare taxa. Differential abundance analyses revealed a greater overlap of responsive taxa between SynCom MIX1 and MIX2. In contrast, functional predictions generated with PICRUSt, although based on prediction, indicated a stronger similarity between MIX2 and MIX3, both of which included Pseudomonas spp., the SynComs that exhibited the highest efficacy in promoting tomato seedling growth. While SynCom inoculation did not lead to a pronounced restructuring of the rhizosphere fungal community in treated plants, we identified trends indicative of a subtle, yet coherent, modulation. Collectively, these shifts indicate that the primary effect was a targeted modulation of the fungal succession dynamics over time.

The bacterial strains used to assemble our SynComs were selected through a top-down approach, by matching the data from the tomato core microbiome³¹ with culturable isolates obtained from the same samples³². Although core taxa are not necessarily the most abundant members of the community, the core microbiome represents the set of taxa consistently associated with a host across contexts and enriched through evolutionary processes to support holobiont fitness^34,35. Native core taxa are expected to colonize soil rapidly, persist over time, and contribute to the assembly of stable microbial communities³⁵. Building on this principle, multispecies microbial consortia can enhance inoculum survival and establishment by mimicking key features of native core communities, while also providing a broader and more functionally diverse set of benefits than single-strain inoculants^10,11.

This set of strains included not only well studied biocontrol agents such as Pseudomonas and Bacillus but also less conventional bacterial genera including Leclercia, Chryseobacterium, Glutamicibacter, and Paenarthorbacter³². Because the isolates were selected through a top-down approach rather than on the basis of predefined functional traits, their phenotypic and genomic characterisation proved essential for uncovering their actual properties. These analyses showed that, although they displayed in planta PGPR and biocontrol properties, the in vitro antimicrobial activity of Chryseobacterium sp. POE47 and Micrococcaceae strains (G. halophytocola PFE44, P. ureafaciens S54 and Paenarthrobacter sp. S56) was actually limited ³². Consistent with these results Plant Bacterial Interaction Factor (PIFAR) analysis showed that these strains clustered separately from Pseudomonas, Bacillus and Leclercia, which exhibited the highest percentages of toxin-related factors³².

Various microbial combinations have been reported to enhance tomato growth, but these consortia typically originated from highly heterogeneous sources, including commercial formulations (e.g., Micomix, Maxi Soil, composed of bacteria and fungi)^27,28,29 and/or strains taken from culture collections^10,16,17,36. However, to our knowledge, no studies have evaluated synthetic consortia assembled exclusively from native bacterial isolates.

Bacillus species, among the most widely utilized PGPR as single-strain inoculants, have been successfully integrated into consortia with other beneficial microbes, such as Rhizobium³⁶, Enterobacter^16,18, Achromobacter^16,37, Acinetobacter¹⁷, and cyanobacteria¹⁹, successfully improving tomato growth. Consistent with this evidence, our SynComs included Bacillus strains as core component. MIX1 (4 strains) combined two Bacillus velezensis representatives with Glutamicibacter and Leclercia. MIX2 (6 strains) expanded this consortium by adding two Pseudomonas strains. MIX3 (10 strains) further extended MIX2 by incorporating two Paenarthrobacter isolates, one Chryseobacterium, and an additional Bacillus strain. The strains of the species Bacillus velezensis (PFE11, PFE42 and PSE31B) and Pseudomonas simiae POE78A and P. salmasensis POE54 used in this study were previously found to exhibit biocontrol and PGP activities in vivo in tomato³². In particular, Bacillus representatives showed the highest growth-promotion effects with strain PSE31B being the most effective in controlling Fusarium crown and root rot³². P. simiae POE78A was the most effective at reducing disease symptoms of bacterial leaf spot, supporting the hypothesis that its effect involves the induction of systemic resistance³². Genome analysis highlighted the presence of various beneficial traits related to nutrient acquisition and stress relief. Metabolites predicted by AntiSMASH included surfactin and fengycin lipopeptides in B. velezensis among others and obafluorin in P. salmasensis POE54. PIFAR analysis clustered Pseudomonas and Bacillus strains together, showing the highest percentages of toxin-related factors compared to the other strains³².

In planta studies have reported that synergistic combinations of Bacillus and Pseudomonas outperform individual strains in the biocontrol of Cephalosporium in maize³⁸, Fusarium in beans³⁹, and Xanthomonas axonopodis pv. malvacearum in cucumber, radish, and cotton^40,41. In tomato, consortia containing Bacillus and Pseudomonas species have been shown to induce systemic resistance against Sclerotium rolfsii⁴², and combinations of Bacillus and Pseudomonas with Trichoderma effectively controlled Fusarium oxysporum and induced systemic resistance to Botrytis cinerea²⁶.

However, interactions between Bacillus and Pseudomonas are often described as ‘frenemy’-type relationships, characterized by a complex balance between competitive and cooperative dynamics^33,43. For instance, negative interactions and biofilm inhibition against B. subtilis were associated in vitro with species such as Pseudomonas capeferrum, P. entomophila and P. protegens often driven by 2,4-diacetylphloroglucinol (DAPG), but not exclusively, since no significant difference in the abundance of encoded BGCs between inhibiting and non-inhibiting isolates was found⁴⁴. On the opposite in dual-species biofilms, B. velezensis and P. stutzeri have been shown to coexist synergistically by occupying distinct spatial and nutritional niches, with Pseudomonas favouring nutrient-rich, oxygen-limited environments and Bacillus thriving in oxygen-rich and nutrient-limited conditions⁴⁵. Andrić and coworkers⁴⁶ showed that Pseudomonas employs toxic cyclic lipopeptides (CLPs) such as sessilins to impose inhibitory pressure, to which Bacillus responds by overproducing surfactin, a lipopeptide that mitigates Pseudomonas toxicity and promotes motility. Building on this work, Andrić and colleagues⁴⁷ further demonstrated that Bacillus can actively perceive Pseudomonas through a siderophore pyochelin-mediated signalling, triggering a robust secondary metabolic response that includes the upregulation of polyketides and bacteriocins with broad activity against both Gram-negative and Gram-positive bacteria⁴⁷. From a colonization perspective, as shown by Andrić et al.⁴⁷, Bacillus is excluded from regions of the root elongation zone already occupied by the sessilin-producing Pseudomonas strain CMR12a, and its cells are absent or undetectable in close proximity to Pseudomonas colonies. Nevertheless, Bacillus is still able to persist: the production of bioactive secondary metabolites and surfactin, enhancing both chemical defense and surface motility, enables the bacterium to relocate to less competitive micro-niches and maintain its fitness along the root surface^46,47. The fact that we observed benefits in planta despite antagonism in vitro does not contradict the literature on Pseudomonas and Bacillus interactions. Rather, it suggests that in the plant environment, their interaction may shift towards a dynamic where effective niche colonization and high metabolite production by the consortium facilitate the exclusion of competitors, ultimately benefiting the host.

Limited information exists on the PGPR activity of the other genera included in our SynComs; however, some have been reported to exert plant-beneficial effects in earlier studies^48,49,50,51. Chryseobacterium balustinum was used in combination with different Pseudomonas species and cyanobacteria, resulting in enhanced wheat growth in a hydroponic system⁵². Moreover, mixtures of different Chryseobacterium species, or combination of Chryseobacterium with other plant-beneficial bacteria (e.g. Pseudomonas sp., Streptomyces sp., Sphingomonas sp.), protected wheat from Rhizoctonia solani and promoted root growth in Arabidopsis⁵³. Paenarthrobacter nitroguajacolicus, used together with Pseudomonas and Bacillus, accelerated post-fire soil recovery by improving aggregation, nutrient availability (N, P, K), and Bituminaria bituminosa germination and development⁵⁴. Much less is known about the performance of Glutamicibacter in multispecies mixtures. Nevertheless, previous studies showed that this genus can enhance tomato growth, salinity tolerance^51,55, and mitigate biotic stresses³². Although Leclercia sp. S52 did not show plant growth promotion activity in tomato³², other representatives of this genus have been reported to promote the growth of rice and soybean^48,56. The combined evidence from single-strain and consortium assays, along with genome-encoded functional traits, underscores that strains belonging to lesser-known genera merit deeper investigation to uncover their molecular mechanisms.

Rhizosphere bacterial community dynamics in the ‘Proxy’ tomato variety were primarily driven by temporal factors. Linear correlation analysis at the phylum level revealed a significant increase over time in Planctomycetota, Verrucomicrobiota, and Firmicutes, accompanied by a reduction in Proteobacteria, Actinobacteriota, and Bacteroidota, regardless of the treatment. These shifts likely reflect the natural selection process exerted by the host plant, which acts as a selective filter for its associated microbiome⁵⁷. As the community transitions from bulk soil to the rhizosphere, it undergoes refinement driven by root-imposed gradients, including carbon source availability, oxygen levels, pH fluctuations, and nutrient depletion⁵⁷. The increase of Patescibacteria, also known as Candidate Phyla Radiation (CPR), only observed in the tomato rhizosphere of MIX2 treated plants and decreasing in the control plants, is difficult to interpret functionally, as little is known about their ecology in agricultural soils. Nevertheless, their detection reflects the presence of non-cultivable members of the microbial “dark matter” revealed by metagenomic surveys⁵⁸. Recent studies show that the only cultivated representatives, the Saccharibacteria, grow as obligate epibionts on Actinobacteria and depend on specialized secretion systems and pili-mediated adhesion for host-associated growth, suggesting that shifts in CPR abundance may mirror changes in host populations or microbial interactions rather than direct plant-related functions^58,59.

Beta diversity analysis further confirmed the temporal shift in bacterial community composition, with samples forming distinct clusters based on sampling time. These findings align with a previous meta-analysis that pointed out that archaeal and bacterial communities exhibit rapid temporal variability, particularly within the first month⁶⁰. In tomato, the richness and evenness of rhizosphere bacterial communities have been reported to be influenced by time^28,29, though those studies focused on distant phenological stages, such as flowering and fruit setting. Flowering, in particular, showed the highest biodiversity in the rhizosphere communities, that was suggested to be associated to increased root exudation, especially in carbon compounds that may favour bacterial proliferation²⁹.

Our study focused on the early stages of tomato growth up to one-month post consortia treatment. By analysing this specific window, we were able to capture the ‘priming’ effect of the SynComs on the resident community before the major phenological shifts associated with reproductive growth. Consistent with Nasuelli and colleagues²⁸, we observed a reduction in alpha diversity indices (Observed and Shannon) in the rhizosphere of consortia-treated plants, which may be attributed to the selective pressure exerted by the specific species present in the inoculum⁶¹.

The most pronounced response in the rhizosphere bacterial communities of treated plants was observed one week after consortium application (T1). A core set of seven differentially abundant bacterial genera was identified across all treated plants; among these, Bacillus and Glutamicibacter, components of all three consortia, were significantly enriched. Furthermore, each consortium induced strong, specific responses. Notably, MIX2 elicited the greatest effect, characterized by high enrichment of rare taxa such as Robiginitalea, Tyzzerella, and Syntrophus. MIX2 was also the only treatment to significantly alter beta diversity compared to the control at T1, alongside a drastic reduction in alpha diversity. A shared signature between MIX1 and MIX2 was the enrichment of Desulfosporosinus, Sulfurovum, and Woeseia. The simultaneous enrichment of these three genera in the rhizosphere may indicate a robust activation of the sulphur cycle, where they can act enhancing the turnover between sulphur reductive and oxidative states⁶². Sulphur is an essential macronutrient that plays a crucial role in synthesizing proteins, chlorophyll, enzymes, and vitamins, which are critical for the plant growth promotion⁶³. Interestingly, at T4 the shared signature across treatments was a depletion of sulphur-cycling bacteria, such as Sulfurimonas, Sulfurovum, Desulfobacter and Desulfatiglans, which may indicate that the treatment-induced perturbation had largely subsided.

A common signature across the three SynCom treatments was the transient enrichment of Azospirillum, together with a clear depletion of the nitrite-oxidizing bacterium Nitrobacter, most pronounced in MIX2, a pattern absent from the untreated control. Nitrospira showed a similar enrichment trend at T1 in all SynComs, although it reached statistical significance only in MIX1. By T4, a broader decline was evident, with the depletion of Nitrobacter, Nitrospira, Nitrosomonas and Nitrosospira, involving different species across the individual SynComs. Azospirillum is a well-characterized rhizosphere-associated diazotroph with documented roles in nitrogen fixation, phytohormone production and root development⁶⁴. In contrast, Nitrobacter and Nitrospira are nitrite-oxidizing bacteria (NOB) that catalyse the aerobic conversion of nitrite to nitrate in the final step of nitrification, and their abundance in agricultural soils is known to respond to changes in nutrient and organic matter availability⁶⁵. Likewise, ammonia-oxidizing bacteria (AOB) such as Nitrosomonas and Nitrosospira contribute to the first step of nitrification and represent key taxa in soil nitrifying communities⁶⁵. The concurrent decline of AOB and NOB at T4 in the SynCom treatments likely reflects a temporary shift in rhizosphere community structure rather than a direct effect on individual functional groups.

By two weeks post-treatment (T2), the impact on the rhizosphere had generally attenuated. The number of differentially abundant bacterial taxa decreased compared to T1, with only three genera common to all treatments, while alpha diversity showed a partial recovery in MIX2 and MIX3. At four weeks (T4), the bacterial community responses converged among the treatments, sharing 18 differentially abundant genera. Most of these shared taxa were depleted. Although beta diversity did not show significant differences between treatments and control at this stage, T4 was characterized by numerous shared responses between MIX1 and MIX2, alongside individual responses in MIX2 and MIX3, primarily defined by the depletion of rare taxa. Despite their low relative abundance, research suggests that rare microbial taxa may be more relevant to ecosystem functioning than previously introduced³⁰. Rare taxa can grow abundant under changing conditions⁶⁶ and respond to perturbations by degrading toxic compounds^67,68. In addition, the high diversity of low abundance taxa represents a large reservoir of genetic traits underpinning a wide spectrum of both known and potentially novel microbial functions⁶⁹.

Hu and colleagues¹⁰ reported that the beneficial effects of a Pseudomonas consortium on tomato growth were more closely linked to changes in resident community diversity and composition, and particularly with an increase in the abundance of initially rare taxa, rather than the direct introduction of plant-beneficial traits in the consortia¹⁰. Similarly, the co-inoculation of Rhodopseudomonas palustris and Bacillus subtilis in rice significantly altered soil bacterial structure, with the rare biosphere exhibiting more pronounced shifts than the abundant community⁷⁰. Structural equation modelling identified the rare bacterial community as the predominant driver of increased crop yields⁷⁰. In our study, this ‘rare biosphere’ may be particularly sensitive to the initial effect of inoculation, where SynCom members, even if they do not persist at high densities, alter the trophic network or competitive landscape^10,11. The specific enrichment or depletion of these genera suggests that the SynCom members can act perhaps through the secretion of secondary metabolites or the modulation of signalling molecules⁷¹.

The taxonomic shifts of bacterial communities observed at T4 were reflected in the predicted functional profiles of the rhizosphere microbiomes. Although MIX1 and MIX2 shared the largest number of differentially abundant bacterial genera, PICRUSt2 indicated that MIX2 and MIX3 displayed more similar functional signatures clearly distinct from MIX1. In these two treatments, pathways related to xenobiotic degradation were proportionally more represented compared to the control, whereas pathways associated with terpenoid and polyketide metabolism, and the turnover of carbohydrates and amino acids showed relatively lower contributions. These contrasting profiles suggest that the presence of Pseudomonas in MIX2 and MIX3 may have contributed to redirecting the community’s predicted functional potential toward more specialized degradation-oriented pathways. Several of the taxa showing the strongest changes at T4 belonged to low-abundance groups, a pattern consistent with the broader notion that rare taxa can exert disproportionate effects on microbiome functioning and contribute to processes such as degradation of complex compounds^30,67, although our data do not allow conclusions about their functional contribution. Overall, these patterns indicate a redistribution of predicted functional potential rather than demonstrated functional changes and should therefore be interpreted as relative shifts in potential functions within the SynCom-influenced communities.

The fungal communities of the tomato rhizosphere followed patterns consistent with those typically observed in soil-root systems, with Ascomycota dominating and Basidiomycota, Mucoromycota and Mortierellomycota contributing smaller fractions^27,31,72,73. Over time, all samples showed a general decline of Ascomycota and Mucoromycota in all treatments, including the control, suggesting that their reduction reflects an underlying temporal dynamic of the rhizosphere community rather than a SynCom-specific effect. In contrast, Basidiomycota displayed a treatment-dependent pattern: while they decreased slightly in the control, they increased in all SynCom-treated soils, indicating that the treatments modulated the relative balance among dominant fungal phyla by sustaining or enhancing groups that would otherwise decline. Although longitudinal studies on fungal dynamics in the tomato rhizosphere are scarce, the temporal patterns observed here align with the notion that fungal communities undergo predictable successional shifts driven by plant development and resource turnover³¹. The β-diversity analysis showed an arch-shaped distribution in the PCoA plot, a pattern that may reflect ordination artefacts and could partly obscure treatment-related differences. This phenomenon is frequently observed in microbial environments with high turnover rates along environmental or temporal gradients, where a single dominant taxon influences the distribution of samples⁷⁴. In our study, this role was likely played by the genus Pseudogymnoascus. Its high relative abundance observed in our dataset is in line with several studies reporting this genus as a common and sometimes dominant component of tomato-associated fungal communities^75,76. Its known ecological traits, like psychrotolerance, saprotrophic ability, and strong colonization of organic substrates, fit well with the growth-chamber environment and with the substrate used in our experiment, which likely favored its proliferation^77,78.

The effects of SynComs on fungal communities were more evident at the genus taxonomic level. At T2, all treatments showed higher abundance of Mucor compared to the control, despite the overall decline of Mucoromycota over time. This indicates that SynComs slowed the natural decrease of this group. In contrast, several Ascomycota genera, including Cladosporium, were consistently depleted across treatments, suggesting that SynComs accelerated the decline of taxa already trending downward. Other Ascomycota responded in a treatment-specific manner, such as Gomphillus, Trichoderma, Paecilomyces, and Aspergillus, which were depleted in individual SynComs. By T4, clearer patterns emerged in Basidiomycota. While this phylum showed a slight decline in control plants, its abundance was sustained or enhanced in all SynCom-treated soils. This suggests that bacterial inoculation may mitigate the natural decline of certain fungal groups by altering niche availability or root exudate profiles⁷¹. Clitopilus was enriched across all SynComs, while Apiotrichum and Myriococcum were enriched particularly in MIX2 and MIX3, points toward a progressive ecological replacement. As dominant Ascomycota declines, bacterial consortia may facilitate the release of niche space, allowing for the expansion of Basidiomycota. Although overall fungal richness and diversity were not significantly altered by the treatment, coherent trends were observed. Together, these patterns suggest that bacterial inoculants influence the fungal assembly not by restructuring overall diversity, but by altering competitive interactions and niche availability in the rhizosphere, an impact likely mediated by shifts in bacterial community composition. Indeed, when introduced microorganisms interact with the resident microflora, they not only directly contribute to plant growth but can also induce indirect effects on the microbial community⁸. These indirect effects are mediated by changes in the diversity, composition, and function of the existing microbiome, that not necessarily requires the persistence of the strains within the community^10,11,28.

In our study, while most inoculated strains remained traceable four weeks post-treatment, a notable decline in their absolute abundance was observed over time. The low detectability or disappearance of specific taxa, such as Leclercia sp. S52 (ASV_2430) and P. simiae POE78A (ASV_6620) by T4, may be attributed to several factors. Biologically, this decline likely reflects a limited rhizosphere colonization fitness of certain strains when faced with competition from the complex resident microbiome, as native soil communities often exhibit strong resistance to invasion⁷⁹. Furthermore, plant-mediated microbial selection, driven by specific root exudate profiles, may have favoured the persistence of certain taxa while filtering out others^57,71. The persistence of significant growth-promoting effects despite the reduction in SynCom abundance suggests that these consortia may act through early-stage indirect microbiome modulation rather than requiring long-term, high-density persistence^10,28.

The comparison of the ASV/OTU sequences with those of the 16S rRNA gene of the bacterial inoculant has already been used by other authors, although with different aims and results^10,22,80. In general, as expected, the matches between a bacterial strain and a community member were not always unambiguous; therefore, a similar sequence could be detected in the control samples, as also observed in other studies^22,80. Moreover, in our study, multiple Bacillus strains were found to be highly homologous to the same ASV. However, it is noteworthy that, even just a few hours after soil drenching commercial plantlets in a commercial horticultural substrate with billions of bacterial cells, the introduced bacteria were not strongly represented in the samples, as compared to the substrate microbial communities. Higher relative abundances of SynCom bacterial strains have, however, been observed in other studies in the tomato root environment²². Innovative approaches, such as metagenomics applying long-read sequencing to the 16S rRNA gene, other housekeeping genes, or metagenomic whole-genome shotgun sequencing, could provide more accurate insights for precisely monitoring the fate and impact of introduced microbial strains or specific rhizosphere resident bacterial members within a complex community context^2,81.

Here, three tomato-native SynComs (MIX1, MIX2 and MIX3) enabled us to characterize plant growth responses alongside time-resolved shifts in the rhizosphere microbiome. All consortia enhanced growth in both determinate and indeterminate tomato varieties, with greater efficacy when Pseudomonas was included (MIX2 and MIX3). By selecting native strains identified as part of the tomato core microbiome, spanning the cultivation chain from seed to commercial greenhouse, we moved beyond traditional culture collections and commercial products to harness host-evolved microbial players. A detailed analysis of the rhizosphere reveals that SynCom application triggers dynamic, time-dependent responses characterized by early treatment-specific shifts and a later convergence of community structure. Notably, the most significant microbial changes occurred in the bacterial communities and were localized within the “rare biosphere”, emphasizing the role of low-abundance taxa as key mediators of plant–microbe interactions. The enrichment of rare genera involved in biogeochemical cycles (e.g., sulphur and nitrogen) and the predicted functional redirection toward specialized metabolic pathways, such as xenobiotic degradation, underscore the complexity of these successional dynamics. Predicted functional profiles also suggested a redistribution of metabolic potential, and fungal communities showed coherent, treatment dependent modulation without wholesale restructuring.

Despite these insights, certain limitations must be noted, including the use of a commercial potting substrate and controlled growth-chamber conditions, which may influence the persistence and behaviour of the inoculants. Furthermore, while PICRUSt2 provided valuable functional predictions, these relative shifts in metabolic potential require future validation through omics-based approaches. This work opens several directions that could be explored in future studies, particularly those aiming to connect early SynCom effects with plant responses and microbiome dynamics in different soils and conditions. Further investigations using complementary functional approaches and measurements of plant physiology could clarify the processes underlying the early shifts we observed. Over time, evidence from multiple studies may help define a more general model of how small native consortia influence microbiome succession and support plant performance.

Bacterial strains used in this study

The ten bacteria used in this study were selected as bioinoculants amongst the core microbiome genera of tomato plants in the cultivation chain³². The strains used belonged to six genera, namely Bacillus, Pseudomonas, Glutamicibacter, Paenarthrobacter, Chryseobacterium and Leclercia. Further details about the strains used and their isolation source are presented in Table 1. The genomes of the bacterial strains are available under BioProject ID: PRJNA1096641.

In vitro evaluation of bacterial strain cross-compatibility

The bacterial strains were tested for compatibility by cross-streak method. Each strain was preliminarily grown on Nutrient Agar (NA, Oxoid, Italy) supplemented with 1% (w/v) dextrose (NDA) at 27 ± 1 °C for 48 h and then streaked in three rows on different media: Luria–Bertani Agar (LB, Laboratorios Conda, S.A., Spain) supplemented with 1% agar (LBA), Tryptic Soy Agar (TSA, Oxoid, Italy) and Potato Dextrose Agar (PDA, Oxoid, Italy), with three replicates per medium. Then, the other strains were streaked perpendicularly to the previous one forming a “cross”. Ten strains per plate were tested and plates incubated at 27 ± 1 °C for 48–96 h. The presence of an inhibition zone at the intersection of the paired strains was recorded as antagonistic activity.

Bacterial growth conditions and inoculum preparation

Strains were routinely maintained on NDA at 27 ± 1 °C and long-term stored in LB supplemented with 20% (vol/vol) glycerol at − 80°C. For the inoculum preparation, all bacterial strains were individually grown in LB broth for 24 h at 27 ± 1 °C in a rotary shaker (180 rpm). Bacterial cultures were centrifuged at 5000 rpm for 15 min, and after discarding the supernatant, the pellets containing the bacterial cells were resuspended in sterile distilled water and the density was normalized to an OD₆₀₀ of 0.1, containing approximately 10⁸ colony forming units (cfu)·mL^− 1. SynComs were assembled mixing bacterial suspensions in equal proportions.

Plant growth promotion assays

One-month-old seedlings of Solanum lycopersicum cv. ‘Proxy’ (indeterminate-growth variety) and cv. ‘Pizzutello’ (determinate-growth variety) were produced under standard conditions in a commercial nursery in Ragusa, Italy. The seedlings were individually transplanted into pots (8 cm Ø) filled with a commercial potting substrate (Professional Mix, Vigorplant: pH: 5.5–6.5; porosity: 90%; manuring: 1000 g/m³ NPK + 1500 g/m³ of organic nitrogen + microelements and slow-release nitrogen). After transplanting, 20 mL of the SynCom suspensions (water for the control) were applied to each pot by soil drenching close to the plant crown. Nine replicates were used for each treatment. The trial was carried out in a growth chamber under controlled conditions (25 °C, RH 60–70%, 16h-light/8h-dark). Plants were monitored regularly and watered on a daily basis. The pots were rearranged randomly every two–three days. Four weeks after the treatment, the following growth attributes were recorded for both cultivars: plant height, shoot and root fresh and dry weights. Plant height of ‘Proxy’ seedlings was also measured weekly.

Analysis of rhizosphere microbial communities

Experimental design and DNA extraction

Rhizosphere soil from the ‘Proxy’ plants was collected at four time points: T0, beginning of the trial (a few hours after SynCom treatment) and at one (T1), two (T2) and four weeks (T4) post-treatment. At each sampling time, nine plants per treatment were removed from their pots and the roots were gently shaken to remove excess soil particles. Each sample comprised five grams of roots with adhering soil from three randomly selected plants, thus forming three replicates per treatment. Samples were suspended in 20 ml of sterile saline buffer, vortexed and then centrifuged (13,500 rpm, 20 min at 4 °C) ³¹. Pellets were stored at − 80°C for subsequent experiments. DNA extraction was performed with the DNeasy PowerSoil Pro Kit (Qiagen, Hilden, Germany) following the manufacturer’s protocol. DNA concentration and purity were determined with a NanoDrop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA) prior to downstream analyses.

16S rRNA and ITS amplicon sequencing

Library preparation and amplicon sequencing were conducted at IGA Technology Services (I-33100 Udine, Italy). For bacterial community profiling, the V3–V4 hypervariable region of the bacterial 16S rRNA gene was amplified with primers 16S-341F and 16S-805R⁸². Peptide nucleic acid (PNA)-clamping was applied during the first 16S rRNA gene amplification step to block amplification of host chloroplast and mitochondrial 16S rRNA gene sequences. For fungal community profiling, primers ITS1 and ITS2⁸³ were used to amplify part of the ITS1 region of the fungal rRNA operon. The 16S and ITS libraries were sequenced on an Illumina NovaSeq6000 instrument (Illumina, San Diego, CA, USA) using 250-bp paired-end mode. Read ends were overlapped to generate high-quality full-length sequences and ensure accurate taxonomic classification.

Bioinformatic analysis

Merging forward and reverse reads, quality filtering and trimming, and Amplicon Sequence Variant (ASV) generation were performed using DADA2 (v. 1.26.0)⁸⁴ in R (v. 4.0.2)⁸⁵. The taxonomic assignment of ASVs was performed using the 16S SILVA 138 ⁸⁶ database for 16S reads, whereas the UNITE database was considered for ITS reads (version 9.0, all eukaryotic dynamic)⁸⁷. Plant-related (e.g., chloroplast and mitochondria) and unassigned ASVs were filtered out from the 16S ASV table. α- and β-diversity analysis was performed using the phyloseq package (version 3.17) in R (v. 4.0.2). In order to evaluate the diversity within each sample, α diversity was estimated on the basis of the Chao1 richness and Shannon diversity indices. Statistical significance of α diversity was analyzed using the Kruskal–Wallis test. β diversity was assessed by means of the Bray–Curtis dissimilarity values and depicted with a Principal Coordinate Analysis (PCoA) to evaluate diversity within each group and among the groups of samples. The PERMANOVA test was conducted to assess statistical significance between each group of samples through vegan⁸⁸ in R. Differential abundance analysis was performed with DESeq2 (v. 1.40.2)⁸⁹ to detect the enriched and depleted bacterial and fungal genera in the SynCom-treated plants compared to the control plants.

In addition, the PICRUSt2 tool was used to predict the putative metagenome functions based on the ASVs. The genes of the predicted metagenomes were functionally annotated on the basis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, using ggpicrust2⁹⁰. The resulting data set was filtered at level 1 to exclude categories not relevant in plant samples (i.e. organismal systems and human diseases). Differential abundance analysis was performed with DESeq2 to assess the enriched and depleted putative functions in the SynCom-treated plants compared to the control plants.

Identification of the inoculated strains in the microbiome

The 16S rRNA gene sequences of the inoculated strains were extracted from the whole genome sequences from Nicotra et al.³². To cross-reference whether the 16S rRNA gene sequences could be found in microbiome sequences of the rhizosphere of the SynCom-treated plants, these were compared with the 16S rRNA gene amplicon-based metagenomic data of rhizosphere samples using the Basic Local Alignment Search Tool BLASTN (http://www.ncbi.nlm.nih.gov). Sequences with ≥ 99.5% similarity were assigned to the same ASV¹⁰. This cut off point (≥ 99.5% sequence similarity) was chosen in order to avoid overestimating the percentage of matching sequences due to the shorter length of 16S rRNA amplicon sequences as compared to the whole 16S rRNA gene sequence of the inoculated strains. In case of multiple matches, the ASV with the highest identity percentage was selected from those exceeding the set threshold.

Statistical analysis

Data from the PGP experiments were examined by analysis of variance (ANOVA) using Minitab 20 statistical software (Minitab, Inc., State College, PA). Means were separated using Tukey’s post-hoc HSD test with a significance threshold of α = 0.05. For each of the most abundant bacterial and fungal phyla, linear regression analyses were conducted in R to evaluate the relationship between relative abundance and sampling time points across treatments.