Root development is maintained by specific bacteria-bacteria interactions within a complex microbiome

Plants grow within a complex web of species interacting with each other and with the plant via a wide repertoire of chemical signals. To model plant-microbe-microbe-environment interactions, we inoculated seedlings with a defined 185-member bacterial synthetic community (SynCom), and manipulated the abiotic environment to enable classification of the SynCom to modules of co-occurring strains. We deconstructed the SynCom based on these modules, identifying a single bacterial genus, Variovorax, which reverts phenotypic effects on root development induced by a wide diversity of bacterial strains and by the entire 185-member community. Variovorax use mechanisms related to auxin and ethylene manipulation to balance this ecologically realistic root community’s effects on root development. We demonstrate metabolic signal interference within a complex model community, defining Variovorax as determinants of bacteria-plant communication networks.

microbiota on agar plates. We inoculated 7-day-old seedlings with a defined 185-member bacterial SynCom (fig. S1 and data S1) composed of genome-sequenced isolates obtained from seedling enrichment trend; modules C and D were composed mainly of Alphaproteobacteria and The mono-association assay confirmed that RGI is a prevalent trait across the plant bacterial microbiota. Thirty-four taxonomically diverse strains, derived from all four modules, induced RGI.
To ascertain the phylogenetic breadth of the Variovorax ability to attenuate RGI, we tested additional Variovorax strains from across the genus' phylogeny ( fig. S9A and data S1; Materials morphological changes in root phenotypes via both auxin-and ethylene-dependent pathways, 165 and both are reverted when Variovorax are present. We conclude that seedlings growing in a 166 realistically diverse bacterial microbiota depend on specific taxa for maintaining a controlled developmental program by tuning the root's response to its microbially-encoded chemical landscape.

Discussion
important are the interactions between these signals, as in the case of quorum quenching (12) or degradation of MAMPs (22). The degradation of bacterially-produced auxin is prevalent in the rhizosphere, including among Variovorax (10,11,23), but its consequences for the plant have  Signaling molecules and other secondary metabolites are products of evolution towards 184 increasing complexity that allows microbes to survive competition for primary metabolites. Our 185 results illuminate the importance of an additional trophic layer of microbes that utilize these 186 secondary metabolites for their own benefit, while providing the unselected exaptation (24) of 187 metabolic signal interference between the bacterial microbiota and the plant host. This potential 188 exaptation, and the consequent homeostasis of plant hormone signals that emerges from it, 189 allows the plant to maintain its root developmental program within a chemically rich matrix.
Moreover, our data suggest that Variovorax accomplish this function while not disrupting the 191 colonization of other taxa, some of which likely provide the host additional ecosystem services.

192
As we proceed to design and develop plant probiotics, we need to consider the maintenance of intra-and inter-taxa microbial homeostasis as potential contributors to satisfactory invasion and persistence of probiotic strains into standing heterogeneous microbial communities.

Fig. 2. Arabidopsis root length is governed by multiple bacteria-bacteria interactions within a community. (A)
Primary root elongation of seedlings grown with no bacteria (NB), with the full 185-member SynCom (Full) or with its subsets: Modules A, B, C and D alone (single modules), as well as all six possible pairwise combination of modules (module pairs). Differences between treatments are denoted using the compact letter display. (B) Primary root elongation of seedlings inoculated with single bacterial isolates. Isolates are colored by taxonomy and grouped by module membership. The strips across the panels correspond to the interquartile range (IQR) as noted at far right. The dotted line represents the cutoff used to classify isolates as root growth inhibiting (cutoff RGI). (C) Binarized image of representative seedlings inoculated with modules A, C and D, and with module combinations AC and AD. (D, E) Heatmaps colored by average primary root elongation of seedlings inoculated with different pairs of bacterial isolates: (D) with four representative RGI-inducing strains from each module (columns) alone (Self) or in combination with isolates from module A (rows), (E) with eighteen RGI-inducing strains (rows) alone (Self) (14). (C) Standardized expression of 12 late-responsive auxin genes across the tripartite and drop-out systems. Each dot represents a gene. Identical genes are connected between bacterial treatments with a black line. Mean expression (95% CI intervals) of the aggregated 12 genes in each treatment is highlighted in red and connected between bacterial treatments with a red line. (D) Primary root elongation of seedlings grown with six hormone or MAMP RGI treatments (panels) individually (Self) or with either Burkholderia CL11 or four Variovorax isolates. Significance between the bacterial treatments is shown using the confidence letter display. (E) GFP intensity of DR5::GFP Arabidopsis seedlings grown with no bacteria, Arthrobacter CL28 and Arthrobacter CL28+Variovorax CL14. Significance within time points is denoted with asterisks. (F) Primary root elongation, standardized to sterile conditions, of wild type (Col-0) auxin unresponsive (axr1-2), ethylene unresponsive (Col-0 + MCP), or auxin/ethylene unresponsive (axr1-2 + MCP) seedlings inoculated with RGI-inducing Arthrobacter CL28 or the Variovorax dropout SynCom (-Variovorax). The blue dotted line marks the relative mean length of uninoculated seedlings. The horizontal shade in each panel corresponds to the interquartile range of seedlings (all genotypes) grown with: Arthrobacter CL28+-Variovorax CL14), or the full 185-member SynCom including 10 Variovorax isolates (Full SynCom). Differences between treatments are denoted using the compact letter display.    20

Standardized primary root elongation
Primary root elongation (cm) NS NS *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** **  Network of statistically significant gene ontology terms contained in the 18 genes cohesively overexpressed in RGI treatments in the tripartite and dropout systems. See Figure  4A and 4B. The network was computed using the emapplot function from the package clusterProfiler in R. A p-value for terms across the gene ontology was computed using a hypergeometric test, additionally the size of each point (Gene ontology term) denotes the number of genes mapped in that particular term.   Differences between treatments were indicated using the confidence letter display (CLD) derived from the Tukey's post hoc test implemented in the package emmeans (37).

Co-innoculation with
Cultures from each strain in the SynCom were grown in KB medium and washed separately, and

372
OD 600 was adjusted to 0.01 before spreading. We performed the experiment in two independent All strains were grown in separate tubes, then washed, and OD 600 was adjusted to 0.02 before 394 mixing and spreading. In vitro growth conditions were the same as in Materials and Methods 2b.

395
Upon harvest, root morphology was measured (Materials and Methods 2c) and plant RNA was harvested and processed from uninoculated samples, and from samples with Variovorax CL14, c. Primary root elongation analysis.
We fitted ANOVA models for each RGI-inducing strain tested.

530
We directly compared differences between the full SynCom and Variovorax drop-out treatment using a t-test and adjusting the p-values for multiple testing using false discovery rate. (Fig 2A, fig. S1A, S4, S7 and S9A-B)

Phylogenetic inference of the SynCom and Variovorax Isolates
To build the phylogenetic tree of the SynCom isolates, we used the super matrix approach isolate genomes from the SynCom utilizing the hmmsearch tool from the hmmer v3.1b2 (44).
Then, we selected 47 markers that were present as single copy genes in 100% of our isolates.

539
Next, we aligned each individual marker using MAFFT (45) and filtered low quality columns in the 540 alignment using trimAl (46). Then, we concatenated all filtered alignments into a super alignment.

541
Finally, FastTree v2.1 (47) was used to infer the phylogeny utilizing the WAG model of evolution.

542
For the Variovorax relative's tree, we chose 56 markers present as single copy across 124

628
We applied a variance stabilizing transformation to the raw count gene matrix. We then 629 standardized each gene expression (z-score) along the samples. From this matrix, we subset 12 genes that in a previous study (15) exhibited the highest fold change between auxin treated and 631 untreated samples. Finally, we calculated the mean z-score expression value of each of these 12 genes across the SynCom treatments. We estimated the statistical significance of the trend of these 12 genes between a pair of SynCom treatments (Full SynCom vs Variovorax drop-out, we estimated a p-value by randomly selecting 12 genes 10000 times from the expression matrix 13. Measuring the ability of Variovorax to attenuate RGI induced by small molecules 641 ( Figure 4D and  treatments v and vi were prepared as described above (Materials and Methods 6a).

682
GFP fluorescence in the root elongation zone of 8-10 plants per treatment were visualized using a Nikon Eclipse 80i fluorescence microscope at days 1, 3, 6, 9 and 13 post inoculation. The 684 experiment was performed in two independent replicates.
From each root imaged, 10 random 30 X 30 pixel squares were sampled and average GFP

693
We transferred four 7-day old wild type seedling and four axr1-2 seedlings to each plate in this

721
We used local alignment (BLASTp) to search for the presence of the10 genes (iacA, iacB, iacC, 722 iacD, iacE, iacF, iacG, iacH, iacI, and iacY) of a previously characterized auxin degradation 723 operon (50) across the 10 Variovorax isolates in our bacterial synthetic community and searched 724 for hotspots of clustered hits within genomes (<10kb between genes). Across the 10 Variovorax 725 isolates scanned, we could not reconstruct any hotspot that recapitulated the previously described 726 operon. Additionally to the iac operon, we scanned the auxin degradation operon described in

727
(51) and could not identify it across the Variovorax isolates. Another piece of evidence that 728 supports the fact that Variovorax lack these degradation operons was the weak (< 30%) identity 729 between the spurious hits we did find to some of the components of these operons across the 730 Variovorax genomes.

863
Data and materials availability:

864
The 16S amplicon sequencing data associated with this study is deposited in the NCBI SRA