Co-inoculation with novel nodule-inhabiting bacteria reduces the benefits of legume–rhizobium symbiosis

We chose the M. truncatula genotype DZA 315.6 (hereafter DZA) because we previously characterized its growth when paired with hundreds of rhizobial strains (Batstone et al. 2022a, 2022b), and it has a well-documented ability to form root nodules with diverse soil endophytes (Etienne-Pascal and Pascal 2013). We chose S. meliloti strain 141 (hereafter Sinorhizobium) as our focal rhizobium partner because it was found to be a high-quality symbiont (conferring a high fitness benefit to its host) across several M. truncatula genotypes (A17 and DZA; Batstone et al. 2022b). The four NRE used in this study, along with Sinorhizobium, were previously isolated from the root nodules of natural populations of M. truncatula, as described in detail in Riley et al. (2023). Briefly, strains were isolated from soils surrounding M. truncatula roots from 21 sites spanning the species’ native range: Spain, France, and Corsica. For more detail on how microbial strains were isolated and taxonomically characterized, see Supplementary methods.

We first assigned our four NREs to putative genera by aligning the 16S rRNA gene sequences of each against the NCBI RefSeq 16S rRNA database (O'Leary et al. 2016), and inferring genera based on the top five BLAST hits. Since NRE strain 522 contained 16S sequences of both Paenibacillus and Brevibacillus origin (see Results), Brevibacillus contigs in the genome assembly for strain 522 were removed using VizBin (Laczny et al. 2015) before further analysis. Because species-level identification of bacteria often requires a high-resolution phylogenetic analysis (Lan et al. 2016), we constructed minimum-evolution phylogenomic trees based on whole-genome alignments using the Type (Strain) Genome Server (Meier-Kolthoff and Göker 2019) (Supplementary methods). For the phylogenetic inference, we conducted all pairwise comparisons among the set of genomes using the Genome BLAST Distance Phylogeny approach (Meier-Kolthoff et al. 2013) and we inferred accurate intergenomic distances under the algorithm “trimming” and distance formula d₅. One-hundred distance replicates were calculated each. We calculated digital DNA:DNA hybridization (DDH) values and confidence intervals using the recommended settings of the GGDC 4.0 (Meier-Kolthoff et al. 2013, 2022). We used the resulting intergenomic distances to infer a balanced minimum evolution tree with branch support via FASTME 2.1.6.1 including subtree pruning and regrafting (SPR) postprocessing (Lefort et al. 2015). We inferred branch support from 100 pseudo-bootstrap replicates each, and we rooted trees at the midpoint. We performed type-based species clustering using a 70% dDDH radius around each of the 17 type strains as previously described (Meier-Kolthoff and Göker 2019). The resulting groups are shown in Table S1 and Table S2. We clustered subspecies using a 79% dDDH threshold as previously introduced (Meier-Kolthoff et al. 2014).

We razor scarified and surface-sterilized seeds of M. truncatula genotype DZA, washing them in 95% ethanol for 30 s and then commercial bleach for 7 min, followed by sterile water for 4 min to rinse off any excess bleach. Before planting, we packed the bottoms of SC10R Ray Leech “Cone-tainers” (Stuewe & Sons, Inc., Tangent, OR) with a small handful of autoclaved polyester fiber to cover their drainage holes. We filled each Cone-tainer with ∼200 mL of an autoclave-sterilized mixture of one-part root wash: one-part sand: four-parts turface MVP calcined clay (Profile Products LLC, Buffalo Grove, IL, USA). The resulting potting media composition is suitable for M. truncatula growth because it facilitates root extraction and rinsing upon harvest, which greatly expediates measuring nodule phenotypes (see below). We sowed seeds 0.75 cm deep in the potting media-filled Cone-tainers using sterile forceps, and immediately watered to compact the soil and prevent seed desiccation. We sowed two seeds in each pot in case one failed to germinate, and then covered the potting media surface and seeds with 0.5 cm of autoclaved-sterilized vermiculite to help retain moisture around seeds after sowing. Prior to inoculation, we thinned all pots to one seedling using sterile forceps.

For each of our five strains (four NREs and one Sinorhizobium), we separately streaked out stocks onto sterile Petri dishes with solid Tryptone-Yeast (TY) medium (Vincent 1970). We allowed colonies to grow in a 30 °C dark incubator until ample growth was observed, approximately 48 h. For each strain, we picked a single colony to inoculate a sterile 15 mL Falcon tube filled with liquid TY medium, and then placed the capped tubes in a shaking incubator set to 30 °C and 200 rpm for overnight growth. Between 20 and 24 h later, we combined tubes of the same strain into a single sterile 50 mL falcon tube, gently inverted to mix, and pipetted out 500 µL to estimate cell density (cells per mL) via measuring absorbance at OD₆₀₀ using a NanoDrop 2000c (Thermo Scientific; Waltham, MA, USA). To ensure that each strain started at an equal inoculation density, we added an appropriate amount of sterile liquid medium and culture to reach a final OD₆₀₀ of 0.1, which corresponds to ∼1 × 10⁸ cells/mL.

We tested the direct and indirect effects of NREs on the legume–rhizobium symbiosis by comparing co-inoculations of Sinorhizobium with each of the four NREs to single-inoculations with Sinorhizobium alone or each of the four NREs alone, for a total of 10 treatments: five single-inoculation, four co-inoculation, plus one uninoculated control to monitor contamination levels. The single-inoculation treatments consisted of either a 500 µL dose of Sinorhizobium (hereafter “Sinorhizobium-only”) or a 500 µL dose of one of the four NREs (hereafter “NRE-only”), each dose totaling ∼5 × 10⁷ cells (1 × 10⁸ cells/mL × 0.5 mL). The co-inoculation treatments consisted of one full 500 µL dose of one of the four NRE strains and a full 500 µL dose of S. meliloti 141 (hereafter “co-inoculation”), the combined dose totaling 1 × 10⁸ cells (1 × 10⁸ cells/mL × 1 mL). Given that Sinorhizobium was the plant’s only source of fixed N and was thus expected to be a major limiter on plant growth regardless of NRE presence, we opted to control for the number of Sinorhizobium cells across inoculation treatments rather than the total number of cells. The positions of the racks that held Medicago plants were randomized within the greenhouse (Supplementary methods). We inoculated plants with their respective treatments seven days after seeds had been sown (i.e., when the first true leaf had emerged) to give the plants sufficient time to establish their root systems and begin photosynthesizing. Plants were destructively harvested 4 weeks post inoculation (Supplementary methods).

For NREs that were found to significantly impact plant traits, we wanted to tease apart whether these effects were acting on the plants directly (i.e., observed in the absence of co-inoculation with rhizobia) or indirectly (i.e., via their interactions with rhizobia). However, because we did not supply plants with any external sources of N in our greenhouse experiment (described above), plants grew poorly in all treatments in which Sinorhizobium was absent, precluding our ability to thoroughly test direct versus indirect effects of NREs on plant growth. To address this limitation, we conducted an additional experiment in which plants were supplied with moderate amounts of N and were either inoculated with a single NRE or were left uninoculated (control). If plants grew larger or smaller when inoculated with an NRE compared to uninoculated controls, then we would consider this NRE to incur direct benefits (i.e., mutualistic) or costs (i.e., pathogenic) to the plant, respectively. We prepared 16 replicate plants per inoculation treatment, using similar protocols as described for the greenhouse experiment (above). Briefly, plants were grown in magenta boxes (PlantMedia, Dublin, OH) and were placed in a growth chamber and were treated with N-supplemented media once a week. One week after transplanting, plants were either inoculated with each NRE individually or with sterile media. After 4 weeks, we destructively harvested plants as described above. For more detail on the nitrogen addition experiment, see Supplementary methods.

For data measured in the greenhouse experiment, we constructed a linear mixed-effects model (LMM) using the package “lme4” (Bates et al. 2015) in R (R Core Team 2020) to first test whether measured traits were influenced by inoculation treatments. For each model, we included the fixed effects of inoculation treatment (10 levels) and whether plants were inoculated with rhizobia or not (2 levels: present/absent) to account for effects varying between the two groups. Plant rack was included as a random effect to account for the spatial arrangement of racks in the greenhouse. We performed type II ANOVA on the model generated for each trait using the “lmerTest” package (Kuznetsova et al. 2017). We were most interested in testing how traits responded to co-inoculations with NREs versus those without; thus we used the “emmeans” package (Lenth 2022) to estimate the marginal means of each treatment group for plants inoculated or uninoculated with rhizobia separately. To compare among groups, we used the Dunnet-adjusted “trt.vs.ctrl” option to compare the treatment groups to two different controls: for plants uninoculated with rhizobia, treatment groups were considered to be plants inoculated with each of the four NREs individually while the control was uninoculated plants. For plants inoculated with rhizobia, treatment groups were co-inoculated plants (rhizobia + each NRE) while the control was plants inoculated with rhizobia only. For the nitrogen experiment, we analyzed data using LMMs and a similar post hoc analysis as described above, but with inoculation treatment (control, inoculated) as fixed effect and pot position as a random effect. For each model, we performed residual diagnostics to check for heteroskedasticity and non-independence, and decided that transforming our response variables was not required.

Finally, we wanted to confirm whether NREs that significantly impacted the symbiosis, which were originally isolated from nodules that formed on plants growing in the field, could reinfect plant tissues after inoculation. We conducted an additional co-inoculation experiment of Sinorhizobium paired with each NRE that affected host growth (see Results). Seeds of M. truncatula DZA were prepared, inoculated, and harvested using the same methods as the nitrogen addition experiment (see above), and with slight deviations (Supplementary methods). Following the harvest, we removed nodules from roots as well as root tissue sections without nodules. We then partitioned these tissue samples for DNA sequencing into two categories: endophytes (within tissues) as well as endophytes + epiphytes (on tissue surfaces). In total, we gathered 48 tissue samples for DNA sequencing (3 treatments × 2 tissue sections × 2 sterilization treatments × 4 replicates = 48 samples). To sequence only endophytic bacteria, we randomly chose four samples from each tissue section per inoculum treatment and surface sterilized each individually by adding 1 mL of 30% commercial bleach, thoroughly mixing for 60 s, and then removing bleach thoroughly by washing each sample with 1 mL of sterile DI water. To sequence epiphytic in addition to endophytic bacteria, we randomly chose four additional samples per tissue section per treatment, washing each individually with 1 mL of sterile DI water, and thoroughly mixing for 60 s to remove loosely-associated bacteria on tissue surfaces.

We extracted and purified genomic DNA from tissue samples using a DNeasy Plant Pro Kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions. For each sample, we evaluated the DNA concentration using optical density measurements obtained by a Nanodrop spectrophotometer (NanoDrop Technologies, Rockland, DE, USA) at 260 and 280  nm wavelengths. We submitted DNA samples to the W.M. Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign for 16S rRNA V3-V4 region amplicon sequencing. Briefly, DNA quality was assessed with a Qubit fluorometer (Invitrogen, MA, USA) and the targeted 16S regions were amplified by PCR using primers V3-F357_N (5′-CCTACGGGNGGCWGCAG-3′) and V4-R805 (5′-GACTACHVGGGTATCTAATCC-3′) using a Fluidigm Biomark HD PCR machine (Fluidigm Corporation, South San Francisco, CA, USA). PCR products were quantified with a QuantiT PicoGreen fluorometer (Invitrogen, MA, USA), then cleaned, purified, and size-selected in a 2% agarose E-gel (Invitrogen, MA, USA) followed by gel extraction with a QIAquick Gel extraction kit (Qiagen, Hilden, Germany). Amplicons were sequenced with Illumina MiSeq v2 platform (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions.

For each library (i.e., sample), we inspected raw, demultiplexed reads for quality using DADA2 v1.28.0 (Callahan et al. 2016) in R (R Core Team 2020). Reads were trimmed with DADA2 to remove forward and reverse primers at lengths of 32 and 35, respectively, to account for the length of V3-V4 primers plus the length of Fluidigm-specific CS primers. We filtered the resulting trimmed reads with DADA2 to retain reads with a maximum number of expected errors of two on the forward strand and five on the reverse strand, in addition to removing putative phiX sequences. We truncated reads at the first instance of a quality score less than or equal to two. The resulting filtered and trimmed reads were dereplicated with DADA2 and amplicon sequence variants (ASVs) were inferred from the dereplicated reads with DADA2. We merged ASVs from forward and reverse read pairs and removed putative chimeras with DADA2, using the “consensus” method. The resulting merged, non-chimeric ASVs were used to construct a sequence table, and the taxonomy of the ASVs was inferred with DECIPHER v2.28.0 (Wright 2016) with the SILVA SSU r138 database (Quast et al. 2012). We manually renamed the resulting taxonomic assignments to the genus Ensifer to Sinorhizobium to reflect the most recent modifications to the genera Ensifer and Sinorhizobium (Kuzmanović et al. 2022) Additionally, we removed ASVs without a taxonomic assignment at the domain level and ASVs assigned to the orders Rickettsiales or “Chloroplast” using phyloseq v1.44.0 (McMurdie and Holmes 2013) to exclude ASVs inferred from host DNA from downstream analyses. We constructed a heatmap visualizing the matrix of variance-stabilized ASVs counts per sample using the R package pheatmap v1.0.12 (Kolde 2019) using color palettes generated by RColorBrewer v1.1.3 (Neuwirth 2022) and viridis v0.6.3 (Garnier et al. 2023). The resulting heatmap was manually edited to position color legends and to format the font face of ASV labels. ASV names with polyphyletic or Candidatus genus assignments from SILVA were manually replaced with “NA”.

To determine the relative enrichment of ASVs across inoculum treatment, tissue section, and surface-sterilization treatment, we used the absolute, non-normalized ASVs counts in conjunction with DESeq2 (Love et al. 2014) v1.40.2 for a differential abundance analysis. Briefly, we built negative binomial generalized models using the ASV counts in DESeq2 with inoculum, tissue section, and surface-sterilization treatment included as factors in the models. To determine ASVs whose abundances were significantly affected by each factor individually, we performed likelihood ratio tests using DESeq2 using full models of the three factors together and reduced models lacking one of each factor. For each test, we shrunk the resulting log2 fold changes of every contrast of the factor’s levels using the DESeq2 function “lfcShrink”. Lastly, to determine ASVs differentially abundant across surface-sterilization treatments among the same tissue sections, we performed an additional likelihood ratio test. The full model in this test contained one factor representing the interaction of tissue section and surface-sterilization treatment, while the reduced model included only the intercept. As done with the other tests, we shrunk the log2 fold changes of each level contrast. We generated all plots with ggplot2 (Wickham 2016).

Using NCBI BLAST + v2.14.0, we created a nucleotide database from the 16S rRNA sequences annotated from the genome assemblies of Sinorhizobium, Bacillus sp. 717A, and Pseudomonas sp. 733B. All ASV nucleotide sequences inferred by DADA2 (Table S7) were used as a query in a BLASTn search with BLAST+ (Camacho et al. 2009), retaining alignments with a minimum e-value of 0.01 and at least 80% identity. ASV28 (Table S7) was determined to be a likely representative of Bacillus sp. 717A (see Results) and was aligned against the NCBI rRNA type strain database on the BLAST + online platform (blast.ncbi.nlm.nih.gov) to verify its potential taxonomic assignment. We determined that ASV1 and ASV5 were representatives of Sinorhizobium 141 and Pseudomonas sp. 733B, respectively (see Results). To investigate the presence of strains related to these two ASVs plus ASV28 in other contexts, we aligned ASV1, ASV5, and ASV28 against the 50 most abundant sequences inferred from 16S rRNA amplicon sequencing in another study of symbiotic bacteria of M. truncatula (Brown et al. 2020) using BLASTn with the same parameters.

Our 16S rRNA BLAST results suggested NRE strain 522 belongs to the genus Paenibacillus, 717A belongs to Bacillus, while 702A and 733B belong to Pseudomonas. This process revealed that the whole-genome sequence for Paenibacillus sp. 522 was contaminated with contigs from a Brevibacillus, and so Brevibacillus contigs were removed from the 717A genome assembly before tree construction to avoid errors; however, it is possible that this Brevibacillus was present in the 522 inoculations. We subjected all strains to minimum-evolution species tree construction using whole-genome sequences. Based on the inferred phylogenomic trees (Fig. 1) and the calculated genome distances between the query and type strains (Tables S1–S2), we identified Paenibacillus sp. 522 (Fig. 1A), Bacillus sp. 717A (Fig. 1B), Pseudomonas spp. 702A and 733B (Fig. 1C) each as novel species in their genera.

We found no clear impact on shoot mass between plants inoculated with Sinorhizobium alone compared to plants co-inoculated with Sinorhizobium and either Paenibacillus sp. 522 or Pseudomonas sp. 702A, whereas plants coinoculated with Sinorhizobium plus either Bacillus sp. 717A or Pseudomonas sp. 733B grew significantly more poorly compared to Sinorhizobium-only plants (P values < 0.05), with a 37% and 33% reduction in shoot mass, respectively (Fig. 2A; Table S4). Every group of plants singly inoculated with an NRE without N-fixing Sinorhizobium and uninoculated control plants were in observably poor condition upon harvest (shoot mass ranging from ∼36 to 46 mg on an average compared to ∼166–263 mg for plants inoculated with rhizobia; Table S5), with no clear difference in shoot mass between the NRE-only and uninoculated control groups (Fig. 2A; Table S4).

Uninoculated plants showed negligible signs of contamination; we only found two nodules on a single control plant, while all other controls (n = 29) were nodule-free. Similarly, we only found two plants inoculated with an NRE but not Sinorhizobium that formed nodules (one with 46 and the other with 31 nodules), which were subsequently removed from downstream analyses.

The lack of nodules on most plants (n = 78/80) that had been inoculated only with NREs (Fig. 2B) indicated that these NREs were unable to form nodules on their own. Similarly to shoot mass, we found no clear impact on the average number of nodules formed when plants had been inoculated with Sinorhizobium-only versus those co-inoculated with Paenibacillus sp. 522 or Pseudomonas sp. 702A (Fig. 2B; Table S4). However, both Bacillus sp. 717A and Pseudomonas sp. 733B co-inoculated plants formed significantly fewer nodules on an average compared to Sinorhizobium-only plants, with a 22% and 25% reduction (n = 20), respectively, while none of the co-inoculated groups were significantly different from each other (Fig. 2B; Table S4). This pattern is consistent with the differences observed in shoot mass, above. For a comprehensive list of all pairwise contrasts for each trait measurement and their associated P values, see Table S4.

While co-inoculations of Sinorhizobium with Bacillus sp. 717A or Pseudomonas sp. 733B negatively impacted host–symbiont traits (Fig. 2B; Table S4), the effects on M. truncatula shoot mass could have been either an indirect result of an antagonistic interaction between NRE and Sinorhizobium or a direct result of parasitism by the NREs on the M. truncatula host. Since plants in the greenhouse experiment did not receive supplemental N, we could not distinguish whether the poor performance of M. truncatula plants singly inoculated with NREs (without N-fixing Sinorhizobium) was due to N-deficiency alone or if NRE inoculations also contributed. To address this, we tested whether the NREs were plant pathogens by examining the direct effects of Bacillus sp. 717A and Pseudomonas sp. 733B on host plants that received supplemental N fertilizer (Fig. 3). There was no clear effect on growth when plants were inoculated with either strain compared to the uninoculated controls (P > 0.1, type II ANOVA), suggesting that the NREs 717A and 733B indirectly affected plant fitness in co-inoculations by inhibiting the legume–rhizobium symbiosis. Again, plants grown without supplemental fertilizer (and withoutSinorhizobium) were in poor condition (Table S6), and no nodules were observed on any plant.

Finally, we used 16S rRNA amplicon sequencing of both surface-sterilized and non-surface-sterilized root and nodule tissues from coinoculated plants to locate these NREs in the endosphere and/or rhizosphere. Although we inoculated plants with combinations of only Sinorhizobium, Bacillus sp. 717A, and Pseudomonas sp. 733B in initially axenic conditions, taxonomic assignments of ASVs after 4 weeks in the greenhouse showed that other bacteria were also associated with our tissue samples (Fig. 4). A total of 248 unique ASVs were inferred across all samples (Table S7); however, the vast majority of these ASVs were in very low abundance (Table S7). This low abundance is expected to some degree considering that 16S rRNA amplicon sequencing can report ASVs with abundances < 0.1% (Nikodemova et al. 2023), which either reflect spurious inferences or trace amounts of microbial cells. Nonetheless, ASVs were inferred from all inoculation, tissue, and surface-sterilization groups (Fig. 4).

Our goal was to determine whether we could recover our inoculum strains from nodules in the expected treatments, which would delineate these strains as either nodule endophytes or tissue surface epiphytes. After aligning our inferred ASVs against the 16S rRNA sequences of our isolate genomes, we found that ASV1 and ASV5 represent our Sinorhizobium and Pseudomonas sp. 733B inocula, respectively (Supplementary Results, Table S9). ASV1 (Sinorhizobium based on SILVA references) was abundant in every inoculum treatment, tissue section, and surface-sterilization treatment. The abundance of ASV1 was not different between samples inoculated with Sinorhizobium alone and either 717A + 141 or 733B + 141 samples (P > 0.05, Wald test, false discovery rate [FDR] adjusted) (Fig. 4, Figs. S1A–B, Table S12). As expected, ASV1 was significantly more abundant among nodule samples compared to root samples (P < 0.01, Wald test, FDR adjusted) (Fig. S1C). Thus, we find no evidence that Bacillus sp. 717A or Pseudomonas sp. 733B influence the abundance of Sinorhizobium once they are within the nodules. ASV5 was assigned to Pseudomonas SILVA reference sequences and was only detected in samples from Pseudomonas sp. 733B co-inoculated plants. This ASV was present in non-surface-sterilized nodule samples plus one additional surface-sterilized root sample at trace levels, but was significantly more abundant in nodule samples compared to root samples (P < 0.01, Wald test, FDR adjusted) (Fig. 4, Fig. S1C, Table S12). For Bacillus sp. 717A, we determined that ASV28 is the most likely representative of out of all other ASVs inferred from the amplicon data. This ASV was found only in non-surface-sterilized nodule and non-surface-sterilized root samples from Bacillus sp. 717A co-inoculated plants, with no difference in abundance between the two tissue sections (P > 0.05, Wald test, FDR adjusted) (Fig. 4, Fig. S1A, Table S12).

To determine whether inoculum strains were differentially abundant across the M. truncatula rhizosphere versus the endosphere, we compared surface-sterilized to non-surface-sterilized samples (Figs. S1D–E). Strains that are endosphere-associated should either be enriched in the surface-sterilized groups (indicating preference for the endosphere over the rhizosphere) or not differentially abundant between the two groups (indicating similar levels across the endosphere and rhizosphere). However, strains that are rhizosphere-associated should be enriched in non-surface-sterilized groups. In the nodule samples, ASV5 (Pseudomonas sp. 733B) and ASV28 (Bacillus sp. 717A) were enriched in non-surface-sterilized samples (P < 0.05, P < 0.01, Wald test, FDR adjusted), while ASV1 (Sinorhizobium) was not differentially abundant (P > 0.05, Wald test, FDR adjusted) (Fig. S1E, Table S12). In the root samples, ASV28 was enriched in non-surface-sterilized samples (P < 0.01, Wald test, FDR adjusted). Although ASV1 was more abundant in surface-sterilized root samples compared to non-surface-sterilized root samples and ASV5 was less abundant (P < 0.05, Wald test, FDR adjusted), neither ASV met our minimum fold-change criterion of 1.5 to be considered enriched in either surface-sterilized or non-surface-sterilized root samples (Fig. S1F, Table S12). Collectively, these comparisons show that Sinorhizobium was indeed endosphere-associated in both root and nodule sections. Bacillus sp. 717A, on the other hand, was rhizosphere associated in nodules and root sections, while Pseudomonas sp. 733B was rhizosphere associated in nodules and endosphere-associated in root sections.

To further investigate the ecological relevance of our NREs, we asked whether our three ASVs of interest, representing the three inoculum strains, could also be detected among amplicon data from a separate study of the leaf, root, and nodule endophytes of M. truncatula endophytes grown in native field soil (Brown et al. 2020). Aligning our ASVs of interest against the 50 most abundant sequences generated by Brown et al. (2020) revealed imperfect but numerous alignments to taxa in their study (Table S10). The best alignment of ASV1 was OTU00001, an “Ensifer” bacterium (now reclassified as Sinorhizobium) that was found to be a significant nodule endophyte in the Brown et al. (2020) study, at 99% identity with 92% coverage of the OTU (Table S10). Our ASV28 (Bacillus sp. 717A) aligned to OTU00009 at 100% identity and 92% coverage (Table S10). This OTU was a significant rhizosphere member in Brown et al. (2020) and was labelled as a “Paenisporosarcina” in the order Bacillales, inferred by aligning OTUs against a SILVA 16S reference alignment (v123). Lastly, our ASV5 (Pseudomonas sp. 733B) aligned equally well to OTU00013 and OTU00003, which were labeled by Brown et al. (2020) as “Streptomyces” and “Pseudomonas”, respectively, at 99% identity and 92% coverage (Table S10). Because the OTU sequence lengths in Brown et al. (2020) were over 100 bases shorter than the ASVs generated here, the alignments of ASV5 to both Streptomyces and Pseudomonas OTUs are likely due to the lower sequence resolution in their study. Given that ASV5 (Pseudomonas sp. 733B) aligned perfectly to a region of the 16S rRNA gene in our Pseudomonas sp. 733B genome assembly, we believe that our ASV5 (Pseudomonas sp. 733B) is likely closely related to Pseudomonas OTU00003, a significant member of the leaf phyllosphere in Brown et al. (2020).

Symbionts in a microbiome do not always evolve to benefit their plant host, as microbes have their own fitness interest that do not necessarily align with host fitness (Klein et al. 2022). In contrast to previous studies, the NREs included here were isolated from the same region as our focal Sinorhizobium strain and negatively impacted the legume–rhizobium symbiosis by decreasing the number of nodules formed and shoot mass. The vast majority of NREs described to date have demonstrated positive or neutral effects on plant or rhizobium performance by increasing host shoot mass or the number of nodules formed (Martínez-Hidalgo and Hirsch 2017; Khan 2019). Since two NREs in this study (Bacillus sp. 717A and Pseudomonas sp. 733B) negatively impacted the legume–rhizobium mutualism while the other two (Paenibacillus sp. 522 and Pseudomonas sp. 702A) showed neutral effects, NREs should not generally be considered beneficial. Our experimental approach allowed us to examine the impacts co-inoculation with NREs have on plant growth. But we recognize that NREs might possess other traits from which plants could benefit, including drought resistance, nutrient acquisition, and herbivory resistance. Further studies that examine rhizobia–NRE interactions in the context of a complex soil microbiome under various environmental conditions would be fruitful for understanding the full range of potential effects of NREs on host fitness.

The negative effects of Bacillus sp. 717A and Pseudomonas sp. 733B on plant growth could have been due to several non-mutually exclusive mechanisms, including that these strains: (i) represent plant pathogens, whereby plants perform worse in the presence of these NREs compared to when they are not associated with them regardless of the presence of rhizobia, (ii) reduce Sinorhizobium’s ability to infect and form nodules on plants; and/or (iii) reduce Sinorhizobium’s ability to proliferate and/or fix N within nodules. The results from our nitrogen addition experiment, whereby plants were provided with sufficient levels of N critical for minimal growth (Küster 2013), provided little evidence for (i), that NREs are plant pathogens; plants performed similarly when inoculated by NREs compared to sterile media. Our results also do not support (iii), that NREs impede Sinorhizobium’s ability to inhabit and proliferate within nodules because we found little evidence for changes in Sinorhizobium ASV abundance within nodules due to NRE co-inoculation. Instead, our results are most consistent with (ii), that NREs interfere with Sinorhizobium’s ability to form nodules on host plants, given the significant reduction in nodules formed by plants when co-inoculated with NREs compared to Sinorhizobium alone.

Finding a reduction in nodulation despite plants being N starved in our greenhouse experiment suggests that these NREs inhibited nodulation, negatively impacting N acquisition by the host from rhizobial partners. When the only source of N is that fixed by rhizobia, legumes will continue to form nodules until they acquire sufficient levels of N via a tightly regulated process (Bauer 1981; Caetano-Anollés and Gresshoff 1991). A recent study in which legumes were either co-inoculated with two strains of rhizobia or singly inoculated with each strain individually found a similar reduction in nodulation due to interference competition between strains, whereby the growth of each strain was inhibited by the presence of the other (Rahman et al. 2023). Given their ability to co-inhabit root and nodule tissues, direct interactions between rhizobia and NREs are likely. In other studies, rhizobia-NRE interactions were found to impact nodulation in a variety of ways. For example, S. meliloti was found to cross-utilize siderophores produced by Exiguobacterium NREs in co-inoculations, which stimulated nodule formation and increased the mass of their Trigonella foenum-graecum host (Rajendran et al. 2012). Additionally, production of indole-3-acetic acid by Pseudomonas trivialis 3Re27 and Pseudomonas extremorientalis TSAU20 stimulated nodule formation by Rhizobium galegae on roots of Galega orientalis and increased G. orientalis growth (Egamberdieva et al. 2010). Furthermore, Bacillus subtilis UD1022 was found to antagonistically downregulate quorum sensing and biofilm formation by S. meliloti Rm8530, which is critical to initiate root nodule formation (Rosier et al. 2021). Future studies could determine whether the NREs Pseudomonas sp. 733B and Bacillus sp. 717A act antagonistically against Sinorhizobium via extracellular secretions that inhibit nodule formation, or whether they compete with Sinorhizobium in root tissues. Either or both of these mechanisms could explain the indirect impacts of these NREs on plant growth. Overall, our results highlight that NREs have the potential to negatively impact plants indirectly through their interactions with rhizobia, which in turn, impacts the fitness benefits received by both legumes and rhizobia. Improving symbiotic outcomes in natural or managed fields will therefore require examining the community context in which the symbiosis occurs rather than the focal partners alone.

Because they shift the costs and benefit exchange, biotic players such as NREs may have important roles in shaping mutualism evolution (Strauss 1991; Strauss and Irwin 2004; Keller et al. 2018). Evolutionary theory predicts that mutualisms are susceptible to “cheaters” that reap rewards from their symbionts without paying any costs (Jones et al. 2015). One current objective of coevolutionary research is to understand how genetic variation and abiotic factors stabilize or destabilize mutualisms (Sachs et al. 2004; Heath and Tiffin 2009; Jones et al. 2015; Batstone et al. 2018, 2022b). Some perspectives suggest that maintaining variation in partner quality, even continuing to associate with potential cheaters, offers a selective advantage to the host if the benefits obtained depend on environmental conditions (Batstone et al. 2018). By extension, the benefits exchanged in symbiosis are likely to depend on the biotic context, and more specifically, the microbiome in which the focal symbiosis unfolds. The natural co-habitation of rhizobia and NREs that impact legume–rhizobia symbiosis here could provide yet another mechanism for maintaining variation in partner quality, if quality depends on the identities of the NREs present within a microbiome (Batstone et al. 2018).

Soil bacteria interact with each other and influence the fitness of their host plant with and without infection (de la Fuente Cantó et al. 2020). Whether the microbes behind these interactions inhabit the nodule/root endosphere or rhizosphere has different implications for community composition and the phenotypic outcomes of symbioses (Brown et al. 2020). While the NRE strains used here were originally isolated from cultures of the M. truncatula nodule endosphere, we were not able to confirm reinfection of nodule endosphere by either NRE strain, although we did detect both on surface nodule and root tissues and Pseudomonas sp. 733B in the root endosphere at relatively low levels. Our findings do not rule-out that these NREs occupied nodules alongside Sinorhizobium. Low-abundance bacteria can have significant impacts on the surrounding community (Lynch and Neufeld 2015; Jousset et al. 2017). These community members may be difficult to detect in amplicon sequencing data if their sequence counts in a sample are too low compared to other community members (Huse et al. 2010; Paulson et al. 2013; Escudié et al. 2018). Given Sinorhizobium’s extensive population growth and genome duplication within nodules (Mergaert et al. 2006), sequencing technologies that enrich for particular targets, such as capture sequencing (Hayden et al. 2022), may be required to detect low frequency non-rhizobial nodule occupants.

The Pseudomonas and Bacillus NRE strains studied here were originally isolated from nodules of M. truncatula grown in soils collected from its natural range. By comparing our sequencing results to those of another study that grew M. truncatula in such soils (Brown et al. 2020), we found OTUs representing taxa closely related to Bacillus sp. 717A and Pseudomonas sp. 733B to be among the most abundant community members. This relationship suggests that NREs in this study are likely reflective of NREs in the natural microbiome of M. truncatula. This presents opportunities to understand the mechanisms by which complex root-associated communities establish and influence each other using this as a model multiplayer symbiosis.