Design of synthetic bacterial communities for predictable plant phenotypes

1The composition of plant-associated microbial communities influences plant health and development. This has raised interest in the use of microbes for biotechnology and agriculture. However, it is challenging to measure the contribution of individual microbes from a complex microbiota to host health. Thus, a number of in vitro screening strategies are commonly applied to identify candidate plant-interacting microbes; however, none of the traits typically screened are correlated with a plant-beneficial outcome. Another common prescreening strategy involves performing plant–bacterium binary-association assays, but only a few have been successfully translated into agricultural settings, suggesting that these assays also fail to capture critical aspects of nature’s complexity. Moreover, it is well established that microbial consortia can produce strong and unexpected effects on host health, and such emergent properties are hard to predict, hindering the rational design of microbial consortia with desired host outputs. Previous strategies to address this conceptual problem included the exhaustive study of possible communities assembled from a small number of microbiota constituents in zebrafish and the analysis of randomized combinations of bacteria in mice. Other approaches often begin with an exhaustive evaluation of all combinations of, for example, the nutrients nitrate and ammonium and the hormones auxin, cytokinin, and abscisic acid on plant root growth and development. Although exhaustive approaches can provide a complete picture of interactions within complex systems, they are unfeasible for systems with more than a handful of variables, given the astronomical number of possible factorial combinations. Even in the rare cases in which functional microbial consortia have been assembled, most studies focus on a single community that is considered a treatment, and rarely is an effort made to dissect the contribution of its constituents. This makes it impossible to establish predictable generalizations beyond the tested communities or conditions used. An instance that dissected the components of a consortium consisted of only 2 bacterial strains. These findings reinforce the necessity for reduced complexity and modular model systems to associate microbial community composition with host phenotypes.

2Our approach is summarized in Fig 1. In short, we first characterized the relationship between in vitro bacterial assays and plant–bacterium binary-association assays, and we used the latter to define functional bacterial blocks. These blocks are groups of bacteria that, by themselves, have a similar influence on a host phenotype. Then, we defined a subset of all the possible communities by constructing partially overlapping synthetic communities (SynComs) of 2 blocks each, tested the effect of these consortia on multiple plant phenotypes, and characterized the plant transcriptional response to these consortia. We evaluated the predictive performance of different statistical models on communities that the models had not seen before. We selected a neural network (NN) because it maximized predictive performance and used this model to design novel synthetic communities that maximized the change in 1 plant phenotype. Finally, we tested the model designs by constructing the novel communities it suggested and validated nearly all of the predicted host phenotypic outputs.

3We focus on bacterial manipulation of the plant response to phosphate (Pi) starvation, a commonly limiting nutrient for plant growth. Pi is an essential macronutrient for plants and also for microbes and is limited in soil. Microbial communities living in the proximity of the plant take up Pi from the environment using a highly efficient Pi transport system. Therefore, the available Pi in the close vicinity of plants is subject to direct and intense competition for uptake between microbes and plants. Although the response of Arabidopsis thaliana seedlings in axenic conditions to phosphate starvation is well characterized, the elucidation of the regulatory mechanisms of this response in the presence of the plant microbiome is only recently emerging.

4We systematically evaluated the performance of a large collection of root bacterial isolates using in vitro screening and binary plant–bacterium association assays as predictors for the effect of derived bacterial consortia on plant phenotypes in response to phosphate starvation. We confirm that bacterial in vitro assays have no correlation with bacterial effects on plant phenotypes. However, we found that plant–bacterium binary-association assays are informative for designing small synthetic communities. Surprisingly, the effects of bacterial consortia on host physiology were mostly additive and independent of bacterial abundances, suggesting that functional stacking within a microbial consortium can determine its effect on host phenotypic response. Finally, we successfully validated novel synthetic communities designed by an NN that led to predictable changes in plant shoot Pi content. Our results provide a useful road map from binary host–microbe assays to the design and testing of useful small consortia to predictably alter host phenotypes.

5In response to Pi deficiency, plants change root exudate metabolite profiles and root architecture to explore Pi-rich soil patches. This may lead to bacterial soil community shifts. In order to learn how root exudate profiles change in response to Pi, we harvested root exudates from A. thaliana plants in response to 2 short and complementary nutritional transitions that mimic the dynamics of Pi stress (S1A Fig; Materials and methods 1a, 1b, 1d, 1e). We demonstrated that our Pi transitions were sufficient to induce a reconfiguration of plant exudate primary metabolic profiles (S1C and S1D Fig, and S1 Table).

6We next tested whether these exudates modified the in vitro growth capacity of a collection of 440 bacterial strains isolated from the roots of Brassicaceae grown in soil that is not overtly Pi deficient (nearly all from A. thaliana) (, S1B Fig, and Materials and methods 1c). We identified a range of bacterial growth behaviors (Materials and methods 1f, 1g) and found that the bacterial growth differences between phosphate conditions are much weaker than the differences between strains (S2 Fig). As expected, phylogeny explained most of the growth differences between strains (S2A Fig and Materials and methods 1h). Most of the bacterial growth parameters provided the same information, so we selected the area under the growth curve (optical density [OD] versus time) (AUC) as a growth marker for subsequent analyses.

7Hierarchical clustering of AUC differences between in vitro conditions identified 10 groups of bacteria that represented different response patterns to exudates derived from roots grown in different Pi concentrations and media supplemented or not with Pi (S2B Fig and S2 Table). We found that root exudates could enhance or inhibit bacterial growth and that this effect could be either general or specific to one type of exudate (S2B Fig). Thus, consistent with previous findings, plant-derived root exudates modulated the growth of bacterial root isolates depending on the plant’s Pi starvation status.

8We selected a subset (n = 183) of the strains from the in vitro assays for determining whether they exerted a functional role on the plant under different phosphate conditions. We selected bacterial isolates that belonged to all of the different response patterns (S2B Fig) and that were most responsive to both Pi levels and the presence of exudates (Materials and methods 1g, S2 Table). We measured the change in plant shoot Pi content, a direct marker of phosphate starvation responsiveness, in response to the presence of each of 183 individual strains, when compared to axenically grown plants. We evaluated shoot Pi content under 4 Pi conditions that represented a 2 × 2 design matrix of 2 Pi levels used for plant germination (full; 1 mM; and depleted, about 5 μM Pi) and 2 Pi concentrations (30 μM Pi and 100 μM Pi) to which seedlings were switched, concomitant with the application of each bacterial strain (Fig 2A and Materials and methods 2). The use of 2 germination conditions in the experimental design allowed us to evaluate the effect of the activation of the phosphate starvation response and the shoot Pi content on the plant–bacterium interaction under different Pi concentrations.

9On average, bacteria had a slightly negative effect on plant shoot Pi content, visualized as a small tail in the bacterial treatment graphs (pink) in Fig 2B. This effect was stronger when the environmental Pi concentration was lower (Fig 2B and 2C and S3 Table). These findings are consistent with our previous results that a bacterial synthetic community drives a context-dependent competition with the plant for Pi. Overall, we found that more strains had a negative than a positive effect on shoot Pi content (Fig 2B, S3 Table, and Materials and methods 2e). Specifically, there were significantly more strains that had a stronger negative effect on plant shoot Pi content in the most limiting Pi conditions (germination in Pi depleted, followed by transfer to 30 μM Pi) (Fig 2B and S3 Table), in which the phosphate starvation response should be active. Conversely, the least Pi-deprived condition (germination in full Pi, followed by transfer to 100 μM Pi) exhibited a significant enrichment of strains that positively affected shoot Pi content (Fig 2B and S3 Table). These results are consistent with bacterial effects on Pi content in the shoot being modulated by the nutritional status of the plant. Importantly, germination conditions did not alter bacterial colonization (Fig 2D), and the effect of individual strains on plant shoot Pi content was independent of the ability of root-inoculated bacteria to colonize the shoot and independent of bacterial titers in different plant organs (Fig 2D; S3A Fig and Materials and methods 2c).

10We detected a weak phylogenetic signal in the ability of bacterial strains to modulate plant Pi content that was significant in only 2 of the 4 conditions (Fig 2C and Materials and methods 2f). Accordingly, we found no correlation between the effect of individual bacterial isolates on shoot Pi content and their in vitro growth phenotype in response to switched Pi levels and root exudates (S4 Fig). Overall, our survey of plant–bacteria binary associations and the resulting distribution of bacterial effects on shoot Pi content argue that the majority of plant–bacteria interactions are competitive, at least in the context of phosphate starvation response.

11We recently demonstrated that the A. thaliana phosphate starvation response is largely antagonistic to immune system function. We therefore asked whether activation of the plant phosphate starvation response could modulate the outcome of binary bacteria–plant interactions. We analyzed shoot Pi content in plants pretreated with phosphite (Phi) (KH2PO3) and then transferred to either 30 μM Pi or 100 μM Pi in the presence of each of 30 selected bacterial strains that either reduced, increased, or had no effect on the shoot Pi content (10 strains per class; S3E Fig). Phi is a nonmetabolizable analog of Pi and its accumulation delays the phosphate starvation response, resulting in low accumulation of Pi in the shoot(S3B and S3C Fig). We found that germinating plants on Phi (low shoot Pi, phosphate starvation response off) dramatically reduced the number of bacterial isolates that diminished shoot Pi content, compared to germination on low Pi (low shoot Pi, phosphate starvation response on) (S3D and S3E Fig). Additionally, we observed that under Phi pretreatment, none of the strains significantly increased shoot Pi content compared to germination on high Pi (high shoot Pi, phosphate starvation response off) (S3D and S3E Fig). Importantly, Phi treatment did not alter bacterial colonization (Fig 2D, S3A Fig and Materials and methods 2c). These findings indicate that activation of the plant phosphate starvation response results in different modes of bacterial interactions with the plant that are independent of shoot phosphate content. These results indicate that an active phosphate starvation response can modulate the outcome of both positive and negative interactions with bacteria, likely mediated via coregulation of the plant immune system. An analogous mechanism has been described for the interaction between A. thaliana and a beneficial fungus.

12We sought to establish whether the results from binary associations are indicative of bacterial effects when a more complex bacterial community is present. We used a microcosm reconstitution approach, in which we inoculated plants with defined complex bacterial synthetic communities (Materials and methods 3a). A subset of 78 strains analyzed in the binary-association experiments was grouped into 3 functional groups consisting of positive (P1-P3), indifferent (I1-I3), and negative (N1-N3) bacteria, depending on their effect on shoot Pi accumulation. For the positive and negative groups, we focused on strains that had a statistically significant effect on plant Pi accumulation, after correcting for multiple testing (Materials and methods 3i). Each functional group was further divided into another 3 blocks of 8–9 bacterial strains, according to the magnitude of their individual effects (Fig 3A, S3 and S4 Tables, and Materials and methods 3i). We then combined pairs of these blocks to define 14 partially overlapping bacterial synthetic communities (Fig 3B). This scheme was designed to maximize the probability of observing extreme plant phenotypes by stacking functionally similar blocks and to gain information from combining the most extreme phenotypic blocks defined in the binary-association assays.

13We evaluated shoot Pi content, primary root elongation, shoot size, and total root network in A. thaliana plants grown in association with the 14 bacterial synthetic communities in the same growth conditions used for the binary-association analysis (Fig 2A and Materials and methods 3b). We found that synthetic communities, like individual bacterial strains, were more likely to reduce plant shoot Pi content, and that synthetic communities made of negative blocks led to lower shoot Pi accumulation than those composed of positive blocks (Fig 3C and 3D and S5 Fig). For example, in general, the estimated negative effect on shoot Pi accumulation for negative blocks is significantly larger than for positive or indifferent blocks (Fig 3C and 3D, S5a Table, and Materials and methods 3j). At the synthetic community level, the effect of the negative strains was clearly dominant; only 2 communities containing negative blocks (I3N1 and N2N3) showed a nonsignificant reduction in shoot Pi content and this in only 1 of the tested conditions. Importantly, the only significantly positive effect with respect to no bacteria involved 2 positive blocks and was weak (P1P3) (Fig 3C, S5b Table, and Materials and methods 3j). We also observed that the majority of the cases in which a synthetic community did not significantly reduce the shoot Pi accumulation occurred under the less Pi-restricted condition (100 μM) (Fig 3C, S5b Table, and Materials and methods 3j), consistent with the results from the individual strains (Fig 2B). This trend was generally consistent for the other plant phenotypes analyzed (Fig 3C and 3D and S5 Fig). Overall, the reduction in shoot Pi content associated with negative blocks correlated with less shoot area, shorter primary roots, and bigger root networks (top and bottom rows in Fig 3C and 3D, and S5 Fig), morphological changes that match the canonical phosphate starvation response in axenic conditions. In contrast, positive bacterial blocks caused less intense plant phosphate starvation response phenotypes. These effects were more obvious in plants grown at low environmental Pi concentration (Fig 3C and 3D and S5 Fig). Thus, the binary-association assays were generally informative with regard to the behavior of bacteria in more complex biotic backgrounds.

14Interestingly, we observed that a number of synthetic communities, for example P2P3 and P1P2, led to increased shoot area compared to axenically grown plants, despite exhibiting reduced shoot Pi content (Fig 3C and S5 Fig). In contrast, plants treated with P1P3 in +Pi_100 μM Pi condition, had shoot Pi content similar to Pi-sufficient plants but unexpectedly exhibited a reduced shoot area (Fig 3C and S5 Fig). Thus, bacterial consortia can decouple shoot Pi-content accumulation from the growth inhibition responses typically associated with the canonical phosphate starvation response.

15We estimated the common (additive) effects of each block of strains across different bacterial backgrounds (e.g., in different synthetic communities) (Materials and methods 3j). Surprisingly, we found that additive contributions of the bacterial functional blocks are sufficient to explain most of the plant phenotypic variation observed (Fig 4). We found that synthetic community membership (i.e., ignoring bacterial relative abundances) typically explained more than 50% of the plant phenotypic variance (Fig 4). This indicates that intrablock bacterial interactions contribute at least as much as interblock interactions to the plant phenotypes tested. Furthermore, the effect of bacterial blocks on the phenotypes analyzed is generally consistent across different synthetic communities, despite each strain’s relative abundance being dependent on the microbial context (S6 and S7 Figs).

16We found that the bacterial abundances in either agar substrate or in the root endophytic compartment were poorly correlated with plant phenotypes (Materials and methods 3e, 3k). Despite the consistent taxonomic profiles of the inoculum, we observed that bacterial communities of agar and root samples were dominated by variable bacterial taxa, depending on the specific combination of bacterial blocks present (S6A and S6B Fig). This suggests that bacteria–bacteria interactions are important in shaping the final community. Furthermore, we found clear taxonomic differences between root and agar samples. Most notably, Streptomyces strains (order Actinomycetales) were particularly good root colonizers despite their limited success on agar, while Pseudomonadales strains were relatively more successful in agar than in root samples (S6A and S6B Fig). These results recapitulate previous findings in natural soils, indicating that Actinobacteria are enriched in A. thaliana roots. Phosphate concentrations in the media had only a minimal effect on the final community composition (S6B Fig).

17We then quantified the information gained by incorporating relative abundance data (S6 Fig) into our additive model (Materials and methods 3k). Surprisingly, in all cases (16/16) the plant phenotypic variance explained by microbiota composition decreased when we incorporated relative abundance (S7 Fig). While in some cases, the differences might not be statistically significant, together, this result demonstrates significantly better performance by the model that ignores relative abundance (p-value = 0.000481; 2-tailed Wilcoxon signed-rank test). Our results indicate that bacterial blocks disproportionately modulate shoot Pi content with respect to their strain abundances, an observation analogous to that seen in bacteria modulating zebrafish immune responses.

18The synthetic communities differentially modulated plant phenotypes related to phosphate starvation response. Therefore, we examined the transcriptomes of plants growing with different synthetic communities. We first explored the expression of a literature-based core set of 193 phosphate starvation response transcriptional markers. Plants did not exhibit induction of phosphate starvation response markers in axenic conditions, even when Pi was low (Fig 5A). However, some synthetic communities induced the canonical transcriptional response to Pi starvation in plants grown on 30 μM Pi (Fig 5A). Plants that showed transcriptional activation of the phosphate starvation response displayed lower shoot Pi accumulation. However, we also observed that some synthetic community treatments lead to low shoot Pi content and no activation of the transcriptional phosphate starvation response (S8 Fig). The effect of synthetic communities was in general dependent on the presence of negative bacterial strain blocks (Fig 5A). In contrast, synthetic communities consisting of only positive blocks of bacteria did not induce the phosphate starvation response transcriptional signature in any condition analyzed (Fig 5A). No induction of the phosphate starvation response genes was observed when the Pi stress was released (following transfer to 100 μM Pi) except for the bacterial combination P3N3, which exhibited induction on 100 μM Pi (Fig 5A). In accordance with the shoot Pi content data (Figs 3C, 3D and 4), we found that additive effects of bacterial blocks could explain the level of transcriptional induction (Fig 5B). The specificity in the bacterial modulation of plant phenotypes suggests that the changes observed in the plant in response to the synthetic communities are linked to bacterial block activities.

19We next explored the overall plant genome-wide transcriptional response to bacteria consortia, Pi conditions, or both. Our design allowed us to test both the response to synthetic communities and to individual bacterial blocks between and within conditions (Materials and methods 3h). As anticipated, plants growing with bacterial synthetic communities on low Pi generally induced phosphate starvation responsive genes and modified the expression of immune system–related genes (S9 Fig, S6 and S7 Tables). Overall, there was not a common response to bacterial presence, with only 45 and 35 genes being significantly up- or down-regulated by more than half of the bacterial blocks, respectively (S6 Table). The number of genes differentially expressed in response to different bacterial blocks did not correspond with the strain composition of the blocks; blocks P2, N3, and I1 altered the expression of the most genes, and blocks I3, N1, and P3 influenced the least (S6 Table). In particular, block I3 only altered the expression of 17 genes, despite being detected in plant roots and surrounding agar (S6 Fig). At the functional level, most of the bacterial blocks induced the expression of the plant defense response, specifically up-regulating genes for salicylic acid biosynthesis (S10 Fig), consistent with overall Bacteria versus No Bacteria comparisons (S9C Fig).

20We also investigated differences between the genes induced by different bacterial blocks. Comparison of genes differentially expressed between positive and negative blocks across all conditions showed that positive blocks had higher expression of genes involved in energy production, while negative blocks specifically induced abiotic stress–responsive gene sets, specifically abscisic acid–related genes (Fig 6A, S6 and S7 Tables). Negative blocks of bacteria also increased the expression of a specific sector of the jasmonic acid response involved in glucosinolate biosynthesis (Fig 6A–6C, S6 and S7 Tables). The glucosinolate pathway modulates the interaction of A. thaliana with a beneficial fungus at low Pi, and its expression is regulated by the master regulator of phosphate stress response, PHR1 (PHOSPHATE STARVATION RESPONSE1). When the environmental Pi was low (30 μM Pi), we observed many more differentially expressed genes between positive and negative blocks (Fig 6A), with negative blocks driving higher expression of genes of both the phosphate starvation and defense response.

21We then focused our analysis on the Pi-limiting conditions (30 μM). In this condition, synthetic communities containing negative blocks showed a strong induction of the phosphate starvation response (Fig 5A). We asked whether the different negative blocks (N1, N2, and N3) differed in their effects. There were almost no expression differences between the 2 most negative blocks (N2 and N3), but we identified 103 genes differentially regulated by bacterial blocks N1 and N3. These genes were mostly stress-related genes, including general abiotic stress and defense response, the expression of which was comparatively reduced in the phenotypically more negative block N3 (Fig 6D, S6 and S7 Tables). This result indicates that under phosphate starvation, all negative blocks activate a similar set of phosphate starvation response genes but differentially suppress other stress responses.

22We found that some genes were induced in response to specific block combinations. For example, we found that PHOSPHATE2 (PHO2), a ubiquitin-conjugating E2 enzyme in A. thaliana required for the degradation of Pi transporters at high Pi, is highly expressed only in plants exposed to the synthetic community P3N3 in all Pi conditions analyzed (Fig 6E). This finding may explain the strong transcriptional response to Pi starvation caused by this synthetic community (Fig 5A). The auxin-regulated gene AUXIN-REGULATED GENE INVOLVED IN ORGAN SIZE (ARGOS) showed a weak positive correlation with the induction of the phosphate starvation response, and it was induced in plants grown with the synthetic communities P3N3 and N2N3 (Fig 6F). ARGOS controls organ size in A. thaliana and its transgenic expression results in enlarged aerial organs. This could serve to counterbalance the negative effect on shoot size that low Pi typically causes.

23Our design of synthetic communities emphasized placing every bacterial functional block into at least 2 microbial backgrounds; therefore, we should be able to estimate bacterial effects that are independent of background. In principle, this estimation could be used to design novel synthetic communities with predictable outputs. We found that additivity of bacterial effects could explain most, but not all, of the host phenotypic variation. Therefore, we built 3 different quantitative predictive models capable of capturing different levels of complexity and evaluated their performance. We constructed a simple linear model (LM), a linear model that included pairwise interactions between bacterial functional blocks (INT), and a Neural Network model (Fig 7A and Materials and methods 4). We focused on shoot Pi content, which had the strongest signal-to-noise ratio (SNR) of all plant phenotypes tested (S11 Fig, Materials and methods 4b). To evaluate the predictive performance of each model, we used a form of cross validation in which the data from each synthetic community were held out one at a time, and the remaining synthetic communities were used to train each of the 3 models; those trained models were then used to predict the plant phenotypic output of the held-out consortium. We found that the NN had the lowest cross-validated prediction error of the 3 models and that the difference was statistically significant (p-value = 0.0073) (Fig 7B). Neural Networks are popular predictive models because they can capture more complex and nuanced relationships that simpler (linear) models cannot; however, this can come at a cost of reduced interpretability. We performed a sensitivity analysis (Materials and methods 4f) by calculating the effect that changing each variable would have on shoot Pi content according to the NN and the 2 linear models (LM and INT). We found a general agreement between the 3 models; for example, all models showed that Pi level in the media and the presence of negative bacterial blocks had the strongest effect on shoot Pi content, but the NN produces much more fine-grained results, because it is able to predict the change differentially across each condition (Fig 7C).

24In order to validate the prediction accuracy of the NN, we chose the 25 bacterial block replacements that were predicted to result in the largest increase in shoot Pi content and experimentally tested whether an increase was produced by these synthetic communities (Fig 7D, S8 Table, and Materials and methods 4g). There was a significant correlation (ρ = 0.42, p-value = 0.0375) between predicted and observed shoot Pi content changes caused by the bacterial block replacements (Fig 7E). Strikingly, we found that 23 out of 25 bacterial block replacements increased shoot Pi content on average (p-value = 0.004; 1,000 permutation tests with synthetic community labels randomly permuted) (S9 Table). Moreover, the improvement in shoot Pi content was statistically significant in 16 out of 25 bacterial block replacements (p-value = 0.032; 1,000 permutation tests with synthetic community labels randomly permuted) (S9 Table). Only 1 out of 25 bacterial block replacements significantly decreased Pi content (S9 Table). Again, we noted little correlation between bacterial abundances and their effect on Pi content (S12 Fig and S9 Table). Compared to linear models (LM and INT), the NN had significantly lower prediction errors (p-value ≤ 4.65 × 10−7) (Fig 7F). In summary, we were able to rationally design novel synthetic communities that lead to predictable plant phenotypic outputs.

25While plant responses to stress have been shown to be influenced by associated microbial communities, causal relationships in plant–microbe interactions in a community context and measured phenotypes have proven difficult to establish. This limitation is, in fact, generally true across complex host–microbial interaction systems. Here, we demonstrate that binary-association assays can inform the design of synthetic bacterial communities that lead to predictable plant phenotypes, an observation seen only once previously, in one animal system. The host phenotypic output of the bacterial synthetic communities was consistent with the output expected from binary interactions. Validation of predictions from an NN confirmed that we could predictably alter certain plant phenotypes by changing the plant’s microbiota membership.

26Other phenotypes and host–microbiota systems can likely be studied with this approach. The only requirements are that a microcosm reconstitution system is available and that functional bacterial blocks can be defined, so that synthetic communities that maximize the expected range of phenotypic variance can be constructed (Fig 3A and 3B). In practice, other aspects that are likely to influence the tractability of a system are the functional bacterial diversity and the SNR of the phenotypes being measured. In the case of plant phosphate starvation, we found that bacterial abundances provided no information, and while it is too early to say if this is a general feature, the only other work that directly manipulated a well-defined microbiota to establish its effect on a host phenotype reached a similar conclusion. If this trend continues across other host–microbiota systems, then our approach has the added advantage that strains need not be distinguishable by a specific marker gene. While a simple additive model typically explained more than 50% of the host phenotypic variation, we found value in utilizing an NN approach that was able to capture more complex relationships but remained interpretable and significantly increased our prediction accuracy for novel communities. Our framework is based on empirical validation and thus remains flexible enough to allow for simpler or more complex models, depending on the case.

27We achieved high prediction accuracy across an untested set of synthetic communities, thus demonstrating that selecting a subset of the possible communities by partial overlap of bacterial functional block pairs is sufficient to characterize this system. This method requires no design of specific heuristics. Thus, this methodology provides an opportunity to expand the capacity for mechanistic understanding not only of biological networks that control plant phenotypes but of other complex ecological systems.

28Furthermore, by focusing on block replacements as testable hypotheses, we provide a simple outcome that can be extracted from both linear and nonlinear predictive models. This can guide the next set of experimental designs, thus providing nonlinear methods like deep learning a stronger empirical grounding, rendering them less of a “black-box.”

29We demonstrate the utility of our approach by defining mechanistic aspects of the plant phosphate stress response in the presence of combinations of bacterial blocks. We observed that bacteria range in their effect on phosphate content in the plant between severely decreasing and moderately enhancing it. These data are consistent with our previous findings that bacterial interactions with the plant are controlled by negative regulation exerted by the phosphate starvation response on the plant immune system. A similar mechanism was described for the interaction under phosphate-limiting conditions of A. thaliana with the beneficial fungus Colletotrichum tofieldiae. Thus, our results provide additional evidence for mechanisms by which plants and bacteria compete in times of nutritional stress.

30The use of multiple bacterial synthetic communities led us to define interesting particular aspects of the phosphate stress response. We observed that certain synthetic communities, such as P2P3 and P1P2, drive an increase in the shoot area compared to axenically grown plants, despite the low shoot Pi content that they engender. These data recapitulate a previous observation on the effect of altering the activity of PHOSPHATE1 (PHO1), a gene required for Pi loading into the xylem. These authors found that shoot Pi content could be uncoupled from the developmental responses typically linked to Pi scarcity. We corroborated that reduced shoot growth is not necessarily a direct consequence of Pi limitation. The observations that both bacterial activity and the modulation of PHO1 expression can uncouple plant phenotypes during the response to low Pi leads us to hypothesize that microbes could interdict PHO1 transport activity, thus modifying Pi translocation from roots to shoots.

31Additionally, we found that synthetic community P3N3 uniquely induced a strong transcriptional response to phosphate starvation in the majority of the conditions tested. Plants exposed to this bacterial combination showed a high-level induction of PHO2, a ubiquitin-conjugating E2 enzyme required for the degradation of Pi transporters at high Pi. This discovery may explain the intense transcriptional response to Pi starvation caused by this particular synthetic community.

32We observed much more variability in the bacterial colonization patterns than in their effects on plant phenotypes. Synthetic communities tended to be dominated by 1 block, but the identity of that block did not correlate with plant phenotype. On the other hand, synthetic communities had remarkably consistent effects on plant phenotypes, and synthetic community membership was sufficient to predict host phenotype. These observations suggest that bacteria–bacteria interactions are critical for microbial community assembly, which is probably a highly dynamic process in which the microbial background determines which bacteria perform well. On the other hand, the effect of bacteria on plant phenotypes is probably due to functional stacking, in which many phenotypically redundant strains with potentially different niches maximize the chance of attaining the desired host phenotypic output. This “lottery model” has been proposed as a major driving mechanism of host colonization by its microbiota at the taxonomic level, and it would be interesting to test whether a similar process governs functional assembly.

33In conclusion, we provide a general method for the study of various biological host–microbiome systems through rational selection of a tractable subset of the possible combinations of bacteria from a defined collection. We demonstrate that complex relationships among host phenotypes, the microbiota, and the abiotic environment can be captured using deep learning techniques. By testing each block of bacterial strains in multiple synthetic communities, and successfully validating predictions derived from an NN, we demonstrated that it is possible to both infer causality and attain generality when it comes to predicting host phenotypes in this complex system. This approach contributes to the rational design and deployment of microbes to improve responses of hosts to biotic and nutritional stresses.

34All seeds were surface sterilized with 70% bleach and 0.2% Tween-20 for 8 minutes, followed by 3 rinses with sterile distilled water. This treatment eliminates any seed-borne microbes on the seed surface. Seeds were stratified at 4°C in the dark for 2 days.

35Col-0 seeds were germinated on Johnson medium (KNO3 [0.6 g/L], Ca(NO3)2.4H2O [0.9 g/L], MgSO4.7H2O [0.2 g/L], KCl [3.8 mg/L], H3BO3 [1.5 mg/L], MnSO4.H2O [0.8 mg/L], ZnSO4.7H2O [0.6 mg/L], CuSO4. 5H2O [0.1 mg/L], H2MoO4 [16.1 μg/L], FeSO4.7H2O [1.1 mg/L], Myo-Inositol [0.1 g/L], MES [0.5 g/L], pH 5.6–5.7, 1% bacto-agar [BD, Difco]), 0.5% sucrose, solidified with 0.6% agar and supplemented or not with 1 mM Pi, in a horizontal position (approximately 160 seedlings per plate). After 7 days of growth, seedlings were transferred to a 12-well plate. Each well was filled with 3 mL of liquid Johnson medium and between 50 and 60 seedlings. Seedlings were transferred to the opposite Pi concentration from the germination conditions (i.e., plants that were initially grown in 1 mM Pi were transferred to liquid medium with no supplementation of Pi and vice versa).

36Plants were grown in liquid media with agitation (30 rpm) for 24 hours in a growth chamber, in a 16-hour light/8-hour dark regime (24°C/21°C). In these conditions, we incubated additional 12-well plates containing Johnson medium alone, supplemented with or without 1 mM Pi. We used these samples as controls for the next experiments. Liquid supernatants containing root exudate and control samples were collected, filtered (0.22 μm), and used for the next experiments.

37All bacterial strains used in this work were isolated from roots of Brassicaceae grown in 2 well-studied wild soils from North Carolina, US. For the screening of the bacterial strain collection in different plant root exudates, bacteria from −80°C glycerol stocks were grown on LB plates at 28°C. A single colony was then inoculated in 200 μL of 2xYT medium (16 g/L Tryptone, 10 g/L yeast extract, 5 g/L NaCL, about 5.5 mM Pi) in a 96-well polystyrene plate (Costar) and covered with a breathable sealing film (Excel) to prevent contamination. Bacterial cultures were grown with agitation at 28°C. After 24 hours, all cultures were diluted 1/10 into the different plant exudates and control media and grown at 28°C with agitation. Exudate and media samples without bacteria were included as controls of contamination (uninoculated controls). To analyze the growth curves of the different isolates, the OD at 600 nm (OD600nm) was measured every 3 hours during the day and every 14 hours during the night for 3 days using a microplate reader (Tecan, GENios).

38Primary metabolites profiling was performed using an ALEX-CIS GCTOF MS in the NIH West Coast Metabolomics Center (University of California, Davis, CA). Plant root exudates and control samples were extracted following Fiehn et al.. 30-μL aliquots of each sample were extracted by 1 mL of degassed acetonitrile:isopropanol:water (3:3:2, v/v/v) at –20°C, centrifuged, and decanted, with subsequent evaporation of the solvent to complete dryness. A cleanup step with acetonitrile/water (1:1) removed membrane lipids and triglycerides. The cleaned extracts were aliquoted into 2 equal portions and the supernatants were dried down again. Internal standards C08-C30 FAMEs were added, and the samples were derivatized by methoxyamine hydrochloride in pyridine and subsequently by N-methyl-N-trimethylsilyltrifluoroacetamide for trimethylsilylation of acidic protons. Data were acquired using the chromatographic parameters published in Fiehn et al.. We used a column from Restek corporation: rtx5Sil-MS (30 m length × 0.25 mm internal diameter, with 0.25-μm film made of 95% dimethyl/5% diphenylpolysiloxane. Helium was used as the mobile phase with a column temperature of 50–330°C and a flow rate of 1 mL min−1. 0.5 μL of sample was injected with 25 splitless time into a multi-baffled glass liner at 50°C ramped to 250°C by 12°C s−1.

39Mass spectrometry parameters were used as follows: a Leco Pegasus IV mass spectrometer was used with unit mass resolution at 17 spectra s−1 from 80–500 Da at −70 eV ionization energy and 1,800 V detector voltage with a 230°C transfer line and a 250°C ion source.

40Raw data were normalized according to NIH West Coast Metabolomics Center (University of California, Davis, CA) quality standards.

41For the metabolite analysis, we standardized the metabolite abundances by dividing the abundance of each metabolite by the total abundance in each sample. To visualize the general patterns in our dataset, we independently applied hierarchical clustering on metabolites and samples. We clustered samples according to their Bray-Curtis dissimilarities, which were calculated with function vegdist in the R vegan package. Metabolites were clustered according to a Pearson correlation dissimilarity matrix, computed using the formula 1−cor(x) via the function cor in R. In both cases, the hclust function in R was used with the method “ward.D” to create the corresponding dendrograms. The resulting dendrograms were visualized together with the metabolite relative abundances using the function heatmap.2 from the package gplots, with the parameter scaling by row, which corresponded to metabolites.

42We calculated the log(fold-change) for each metabolite in each exudate sample, in relation to its respective media control sample (Table 1).

43Then, for each sample, we calculated the log2(fold-change) of each metabolite with respect to its control by subtracting the log2-transformed mean relative abundance of that metabolite in the control samples to the log2-transformed relative abundance value of the metabolite in a given sample.

44To visualize the matrix derived from the log(fold-change) calculation, we independently applied hierarchical clustering on the metabolites and on the samples. We clustered exudate samples based on the euclidean distance calculated with the dist function in R. Metabolites were clustered according to a Pearson correlation dissimilarity matrix, computed using the formula 1−cor(x) via the function cor in R. In both cases, the hclust function in R was used with the method “ward.D” to create the corresponding dendrograms. The resulting dendrograms were visualized together with the metabolite log2(fold-changes) using the function heatmap.2 from the package gplots, with the parameter scaling by row, which corresponded to metabolites.

45All the scripts developed for this analysis are available at: https://github.com/isaisg/primary_metabolites_phosphate.

46Bacterial growth curves were quality filtered by removing strains that had growth profiles that were highly similar to uninoculated samples. All the following operations were done with functions available via the PGCA R package (https://github.com/surh/PGCA). For each strain and condition, the median growth curve was obtained by calculating the median OD600nm per time point. The 4 resulting growth curves (for 4 conditions) were concatenated and grouped by hierarchical clustering based on their correlation distance according to the formula dxy = 1 − ρxy, where ρxy is the Pearson correlation coefficient between strains x and y. The resulting clustering dendrogram was cut at a height of 0.5, based on visual inspection. Clusters that had more than 40% uninoculated controls were discarded together with any strains that fell within them. The remaining uninoculated controls were also discarded. We then calculated the AUC for each strain and condition by adding the median OD600nm per time point using PGCA (https://github.com/surh/PGCA).

47Four other growth curve features were extracted for each strain and condition. First, a moving average filter smoothed all data. For this, we applied a regular moving average filter with window size of 5 and step size of 1 per individual replicate. We defined the following metrics for each strain under each condition using the smoothed growth curves:

48Maximum OD (MAX). MAX is the maximum OD reached among all 4 replicates and all time-series points.

49Mean density over the last 3 measurements (L3M). L3M is used to give an approximate measurement of the OD at saturation by averaging the last 3 measurements in the time series over all replicates.

50Mean time to reach half maximum (HMT). The time for a strain to reach half maximum is the time stamp of the earliest sample in the time series data that has an OD greater than or equal to half of the MAX. We take the mean over all replicates to calculate HMT. We use HMT to show how fast the strain is growing.

51Maximum growth rate (MGS). The growth rate at a particular time point is calculated by (log2 xtk+1 − log2 xtk)/(tk+1 − tk), where xtk represents the OD of the strain at time tk. MGS is the maximum growth rate among all replicates and all time points.

52Each of these features was calculated for the 4 conditions: Johnson media supplemented or not with Pi (+Pi or −Pi, respectively); Johnson media supplemented with root exudates from 2 short and complementary nutritional Pi transitions (+Pi_−Pi or −Pi_+Pi). Furthermore, we also normalized each growth curve by their control by calculating log2(fold-change) on each specific feature between −Pi_+Pi and +Pi (MLF), and +Pi_−P and −Pi (PLF).

53In general, several features analyzed contain redundant information (S2A Fig). This redundancy allowed us to select 1 parameter, the AUC, as an aggregate measure of bacterial performance for further analysis.

54For grouping strains according to their growth performance, each strain’s AUC was log transformed and then standardized per condition by subtracting its mean and dividing by the standard deviation. Hierarchical clustering was used with the euclidean distance and the complete linkage method in R. In order to select strains to test in binary-association assays, an ANOVA model was fit on each strain, using the AUC values as dependent variables and condition (media) as the only independent variable. We calculated the R2, which indicated which proportion of the variation in the bacterial growth performance (AUC) is attributable to media—in other words, how much of the bacterial growth is affected by plant exudates and phosphate levels. We prioritized testing multiple strains per cluster that had the highest R2 values. The code and data to perform these analyses are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

55For all strains with an available Sanger generated 16S rRNA gene sequence (395 strains out of 440), we used MUSCLE to perform a multiple sequence alignment with default parameters. We then filtered out positions that had more than 99% gaps as well as the top 10% most entropic sequences using QIIME. The resulting filtered alignment was used to build a maximum likelihood tree with FastTree, using midpoint rooting.

56We standardized all the phenotypes to allow for simultaneous visualization and easier comparison. We used the phylosig function from the phytools R package to test Pagel’s λ for phylogenetic signal. Results were visualized with the ggtree R package. The code and data to perform these analyses are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

57Plants were germinated in axenic condition on Johnson medium plus 0.5% sucrose with 1 mM Pi, about 5 μM Pi (traces of Pi from the agar [Difco]), or supplemented with 1 mM phosphite in a vertical position for 6 days. Seedlings (10 seedlings per plate) were then transferred to 30 μM Pi and 100 μM Pi media (without sucrose), alone or with the bacterial culture, at 105 c.f.u/mL of medium for another 7 days. We harvested shoot samples for Pi quantification. Plants were grown in a growth chamber in a 16-hour light/8-hour dark regime (24°C/21°C). We used 3 replicas per condition and bacterial strain; the experiment was repeated once.

58To confirm that plants germinated in different Pi regimens or with phosphite differentially activated the phosphate starvation response, plants carrying the phosphate stress reporter construct IPS1:GUS were grown in Johnson medium containing 1 mM Pi, 1 mM phosphite or traces of Pi (about 5 μM Pi). After 7 days, the expression of the reporter constructs IPS1:GUS was followed by GUS staining (S3C Fig).

59To study the colonization of the plant shoot by root-inoculated bacteria (Fig 2D and S3A Fig), we germinated Col-0 on Johnson medium 1 mM Pi, about 5 μM Pi (traces of Pi from the agar [Difco]), or 1 mM phosphite for 7 days. Seedlings (5 seedlings per plate) were then transferred to 2-compartment plates for a week. In this system, root and shoot were placed in different compartments separated by a plastic barrier to prevent microbe diffusion through the medium. The root compartment was previously filled with Johnson medium 30 μM Pi solidified with 1% agar containing bacteria, and the shoot compartment was filled with a solution of agar (water + 1% agar). Plants were grown in a growth chamber in a 16-hour light/8-hour dark regime (24°C/21°C). We harvested shoot and root separately, and bacteria accumulation in shoots, roots, and root-surrounding agar was analyzed in mock-inoculated plants and plants colonized by bacteria. We used 5 replicas per condition and bacterium; the experiment was repeated once with similar results.

60For binary association and synthetic community experiments, a single colony was inoculated in 4 mL of 2xYT medium in a test tube. Bacterial cultures were grown at 28°C with agitation overnight. Cultures were then rinsed with a sterile solution of 10 mM MgCl2 followed by a centrifugation step at 2,600 g for 8 minutes. This process was repeated twice to eliminate any additional nutrient supplementation in the media. The OD600nm was measured and, assuming that 1 OD600nm unit is equal to about 109 c.f.u/mL, we equalized individual strain concentration to a final value of 105 c.f.u/mL of medium. The medium was cooled (to 40–44°C) near the solidification point and then the bacterial strain or the bacteria strain mix was added to the medium with agitation.

61For bacterial colonization analysis in binary-association assays, roots, shoots, and agar samples were harvested and weighed. Roots and shoots were rinsed 3 times with sterile distilled water to remove agar particles and then placed in a sterile tube with 200 μL of 10-mM MgCl2 and glass beads. Plant material and agar samples were then crushed using 2 cycles of 30 seconds, frequency 30 s−1, in a TissueLyser II (Qiagen). These samples were serially diluted, plated on LB, and c.f.u/mL of harvested samples were determined.

62The method of Ames was used to determine the free phosphate concentration in the shoots of seedlings grown on different Pi regimens and treatments.

63We identified which strains had a statistically significant effect on shoot phosphate accumulation. We compared the phosphate content between plants treated with a given strain and axenically grown plants from the same experiment. To determine significance, we log transformed the phosphate content and used an ANOVA model in R with terms for bacterial treatment and biological replicate. Each of the 4 phosphate conditions was analyzed independently. By testing the bacterial treatment term in the model, we determined whether the effect was significant. We corrected all the obtained p-values using the Benjamini-Hochberg method, and considered significant all effects that had a q-value < 0.1. The code and data to perform this statistical analysis are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

64To facilitate the simultaneous visualization and to control for batch effects, we scaled all the bacterial effects with respect to the no bacteria controls in the same experiments. We did this by using the results from the ANOVA analysis (Materials and methods 2e) as input for the phylogenetic analysis. These ANOVA results represent the log(fold-change in Pi content, with respect to no bacteria) caused by each strain. We used the phylosig function from the phytools R package to test Pagel’s λ for phylogenetic signal. We utilized the same tree that we constructed for the phylogenetic analysis of in vitro growth (Materials and methods 1h). 177 out of 183 strains were present in this tree. Results were visualized with the ggtree R package. The code and data to perform these analyses are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

65For synthetic community experiments, plants were germinated in axenic condition on Johnson medium plus 0.5% sucrose with 1 mM Pi or about 5 μM Pi (traces of Pi from the agar [Difco]) in a vertical position for 6 days, then transferred (10 seedlings per plate) to 30 μM Pi and 100 μM Pi media (without sucrose), alone or with the synthetic community, at 105 c.f.u/mL of medium (see Materials and methods 2b) for another 7 days. Plants were grown in a growth chamber in a 16-hour light/8-hour dark regime (24°C/21°C). Plant material (root) was collected for 16S profiling, for transcriptional analysis (seedlings), and shoots were used for Pi quantification. This experiment was repeated once with similar results, and each repetition included 3 biological replicates (10 plants each) per condition and per synthetic community analyzed.

66Primary root length was measured using ImageJ, and shoot area and total root network were measured with WinRhizo.

67For bacterial colonization analysis using 16S, roots were surface sterilized with freshly made 10% household bleach with 0.1% Triton-X100 for 12 minutes. Following the bleaching, roots were rinsed once in sterile distilled water, then placed in 2.5% sodium thiosulfate to neutralize the bleach for 2 minutes, and rinsed once more with sterile distilled water. Roots were then freeze-dried and powdered in a 2-mL tube with glass beads using the MPBio FastPrep for 20 seconds at 4.0 m/s. These samples were used for DNA extraction using 96-well format MoBio PowerSoil kit (SDS/mechanical lysis).

68To quantify bacteria from agar samples, a freeze and squeeze protocol was used. Syringes with a square of sterilized miracloth on the bottom were completely packed with agar samples and kept at −20°C for a week. Samples were thawed at room temperature and syringes were squeezed gently into 50-mL tubes. Samples were centrifuged at max speed for 20 minutes, and most of the supernatant was discarded. The remaining 1–2 mL of supernatant containing the pellets was moved into clean microfuge tubes. Samples were centrifuged again, supernatants were removed, and pellets were used for DNA extraction with the 96-well format MoBio PowerSoil kit (SDS/mechanical lysis). DNA was extracted simultaneously for both agar and root samples. We performed randomization of the sample order using a mechanical method.

69We amplified the V3–V4 regions of the bacterial 16S rRNA gene using primers 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). Libraries were created using a modified version of the Lundberg et al.. The molecule-tagging step was changed to an exponential amplification to account for low DNA yields, with the following reaction and temperature cycling conditions (Table 2).

70Following PCR cleanup to remove primer dimers, the PCR product was indexed using the same reaction and 9 cycles of the cycling conditions described in Lundberg et al.. Sequencing was performed at The University of North Carolina at Chapel Hill, Chapel Hill, NC, on an Illumina MiSeq instrument, using a 600-cycle V3 kit. The raw data from these sequencing experiments are available at the ENA Sequence Read Archive (accession number PRJEB22060).

71Synthetic community sequencing data were first processed with MT-Toolbox. Categorizable reads from MT-Toolbox (i.e., reads with correct primer and primer sequences that successfully merged with their pair) were quality filtered with Sickle by not allowing any window with a Q-score under 20, and trimmed from the 5′ end to a final length of 350 bp. The resulting sequences were matched to a reference set of the strains in the synthetic community generated from Sanger sequences and A. thaliana organellar sequences. Sequence mapping was done with USEARCH7.1090 with the option “-usearch global” at a 98.5% identity threshold (which translates to 4 mismatches for our sequence length). Sequences from each sample were only mapped to the expected sequences, given the synthetic community block composition. Matching sequences were used to produce an abundance table that was input into AMOR. Host nuclear and organellar-derived sequences were filtered out, and samples that had less than 400 reads were removed. The resulting table was collapsed by strain order or strain block in AMOR. Aggregation was performed by arithmetic addition of all the counts belonging to a given group. Taxonomic and block-level distributions were visualized with the phylogram function in AMOR, and constrained ordination was performed with vegan’s capscale function and visualized with ggplot2. The code and data to perform this analysis are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

72Total RNA was extracted from A. thaliana seedlings according to Logemann et al.. Frozen seedlings were pulverized in liquid nitrogen. Samples were homogenized in 400 μl of Z6-buffer; 8 M guanidinium-HCl, 20 mM MES, 20 mM EDTA pH 7.0. Following the addition of 400 μL phenol:chloroform:isoamylalcohol, 25:24:1, samples were vortexed and centrifuged (20,000 g, 10 minutes) for phase separation. The aqueous phase was transferred to a new 1.5-mL tube and 0.05 volumes of 1 N acetic acid and 0.7 volumes 96% ethanol were added. The RNA was precipitated at −20°C overnight. Following centrifugation (20,000 g, 10 minutes, 4°C), the pellet was washed with 200 μl sodium acetate (pH 5.2) and 70% ethanol. The RNA was dried and dissolved in 30 μL of ultrapure water and stored at −80°C until use.

73We performed an RNA-Seq experiment to evaluate the effect of the different synthetic communities on the transcriptional profile of A. thaliana seedlings. Wild-type Col-0 was germinated under Pi-sufficient and Pi-depleted conditions for 6 days and then transferred to 30 μM Pi and 100 μM of Pi for another 7 days, with or without bacteria. At this point, seedlings were harvested for RNA-Seq. This experiment was repeated once, and each repetition included 3 biological replicates per genotype per condition per synthetic community analyzed. Raw reads and read counts are available at the NCBI Gene Expression Omnibus under accession number GSE102248.

74Illumina-based mRNA-Seq libraries were prepared from 1 μg RNA. Briefly, mRNA was purified from total RNA using Sera-mag oligo(dT) magnetic beads (GE Healthcare Life Sciences) and then fragmented in the presence of divalent cations (Mg2+) at 94°C for 6 minutes. The resulting fragmented mRNA was used for first-strand cDNA synthesis using random hexamers and reverse transcriptase, followed by second-strand cDNA synthesis using DNA Polymerase I and RNAseH. Double-stranded cDNA was end repaired using T4 DNA polymerase, T4 polynucleotide kinase, and Klenow polymerase. The DNA fragments were then adenylated using Klenow exo-polymerase to allow the ligation of Illumina Truseq HT adapters (D501–D508 and D701–D712). All enzymes were purchased from Enzymatics. Following library preparation, quality control and quantification were performed using a 2100 Bioanalyzer instrument (Agilent) and the Quant-iT PicoGreen dsDNA Reagent (Invitrogen), respectively. Libraries were sequenced using Illumina HiSeq2500 sequencers to generate 50-bp single-end reads.

75Initial quality assessment of the Illumina RNA-Seq reads was performed using the FASTXToolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html). Cutadapt was used to identify and discard reads containing the Illumina adapter sequence. The resulting high-quality reads were then mapped against the TAIR10 A. thaliana reference genome using Tophat, with parameters set to allow only one mismatch and discard any read that mapped to multiple positions in the reference. The Python package HTSeq was used to count reads that mapped to each one of the 27,206 nuclear protein-coding genes. Raw sequencing data and read counts are available at the NCBI Gene Expression Omnibus accession number GSE102248.

76For expression of the phosphate starvation response core transcriptional markers and the clustering analysis, we converted the count table into a table of reads per kilobase per million (RPKM) and standardized these values per gene by subtracting the mean gene expression and dividing by the standard deviation of each gene’s expression. Hierarchical clustering was performed with the R function hclust using the complete linkage method. Gene ontology enrichment analysis was performed on the PlantGSEA online platform.

77For specific hypothesis tests, we used edgeR to fit a quasi-likelihood negative binomial model with the function glmQLFit, after estimating tagwise dispersion parameters. We then applied a quasi-likelihood ratio test with the function glmQLFTest in edgeR using different sets of contrast. Our model specification included indicator terms for each bacterial block, as well as terms for phosphate condition, biological replicate, and interaction between phosphate condition and each block. This was defined according to the following formula:

78The first 9 terms are indicator variables for the bacterial blocks of the synthetic communities that take the value of 1 when that block is present. Phosphate has 4 levels that correspond to the 2 × 2 phosphate condition experimental design. The experiment has 5 levels that correspond to independent biological replicates of the synthetic community experiments (each community was in 2 biological replicates). A total of 41 coefficients, including intercept, are generated from this design. The definitions of the contrasts used for the different hypotheses are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

79The code and data to perform this analysis are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

80Due to the fact that the majority of the strains analyzed in binary interaction had a negative effect on shoot Pi accumulation, we first identified strains that had a significant positive effect (q-value < 0.1 from ANOVA) in the 2 conditions that ended on 100 μM phosphate. We found 26 strains and we labeled them as positive strains. Many strains are negative in 2 conditions or more, so we first identified strains that had a statistically significant (q-value < 0.1 from ANOVA) negative effect on shoot phosphate accumulation in at least 3 of the 4 conditions. We then identified strains that had a statistically significant negative effect in at least 2 conditions, but with higher statistical confidence (q-value < 0.05). We removed 2 strains that did not come from our Brassicaceae cultivation efforts in 2 natural soils (strains Pseudomonas fluorescens WCS417r and Rhizobium sp. R219). We combined the 2 sets of negative strains and obtained 26 strains that we termed negative strains. We finally identified strains that had no statistically significant effect (q-value > 0.1 from ANOVA) on phosphate accumulation across all conditions, and from this set we randomly and computationally selected 26 strains that we termed the Indifferent strains. We finally calculated the mean effect of each strain on shoot phosphate accumulation by averaging the coefficients from the ANOVA in all conditions. We sorted the strains within each group according to that mean, and we divided each group into 3 blocks of bacteria by taking groups of 9, 9, and 8 strains (9 + 9 + 8 = 26). The bacterial blocks defined this way are the basis for the synthetic community design. The code and data to perform this statistical analysis are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

81We tested a set of 14 partially overlapping synthetic communities. Each of the synthetic communities was made of a combination of 2 bacterial blocks. We made 9 synthetic communities by combining adjacent blocks (i.e., blocks that are next to each other when sorted by their mean effect). Each of these 9 communities is represented as an outer arc in Fig 3B, and they are composed mostly of strains that have similar effects on shoot phosphate accumulation when tested in binary association; however, they represent the widest possible range of mean effects, and so we expect they might produce different plant phenotypes. We constructed another 5 synthetic communities, which are represented as inner arcs in Fig 3B that represent extra combinations of the most extreme blocks (i.e., P1 and N3), in order to test how strong their effects are in a variety of backgrounds.

82To determine the level of consistency of bacterial block effects on different plant phenotypes, we first compared plants inoculated with each synthetic community versus axenically grown controls. We then estimated the main effects of each block using multiple regression and compared the coefficients obtained from both methods. Phosphate content and shoot area measurements were log transformed to reduce heteroscedasticity, and thus the coefficients from this analysis should be interpreted as log(fold-change) between inoculated plants and axenic controls. The LMs adequately represented measurements from primary root elongation and total root network, and so coefficients for these 2 phenotypes should be interpreted as the difference between inoculated plants and axenic controls. Appropriateness of the LMs and the need and correctness of transformation was determined via inspection of residual plots. Only the 14 original synthetic communities were included in the analysis.

83To estimate synthetic community effects, we fit a linear model per phenotype and condition, with only the samples of one synthetic community at a time, plus the axenic controls performed in the same experiment. We had only bacterial treatment and experiment variables, according to the following formula:

84The resulting SynCom coefficient was denoted as the “measured” synthetic community effect: To find the expected phenotypic effect of each synthetic community, we first estimated each block’s additive (main) effect. We did this by fitting all the data from each media condition into 1 linear model containing only terms for each block and experiment as independent variables, using the following formula:

85Each of the block variables is encoded as an indicator variable, and each has the value of “1” when the corresponding block is present and “0” when the corresponding block is absent. We finally obtained the “predicted” community effect by arithmetically adding the coefficients for the 2 blocks that make each community. To compare the effects of both synthetic communities and blocks together, we rescaled them by dividing them by the standard deviation of all coefficients from the same phenotype in the same condition (i.e., by column). The code and data to perform this statistical analysis are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

86In order to determine if bacterial relative abundances are correlated with host phenotypes, we repeated the analysis above (Materials and methods 3j) but substituted the indicator variables of each block for the relative abundance of that block.

87Because the relative abundance data are obtained from a pool of plants and require destructive sampling, we averaged all the measurements from each independent biological replicate and used those averages as input for the analysis. This effectively reduces the number of observations and increases the standard errors. In order to achieve a fair comparison, we also recalculated the purely additive contributions of each block (without considering abundance), and we equalized the norm of the coefficients between the 2 versions of the analysis. Results are presented in S7 Fig, and the code and data to perform this statistical analysis are bundled in the R package (wheelP) [https://github.com/surh/wheelP].

88Plants were germinated in axenic condition on Johnson medium plus 0.5% sucrose without supplementation of Pi in a vertical position for 6 days, then transferred to 30 μM Pi medium (without sucrose), alone or with the synthetic communities at 105 c.f.u/mL of medium for another 7 days. Plants were grown in a growth chamber in a 16-hour light/8-hour dark regime (24°C/21°C). Plant material (shoot) was collected for Pi quantification and roots were used for 16S profiling. This experiment was repeated once.

89To determine the feasibility of predicting the phenotypic effect of novel communities, we calculated the SNR for each phenotype. There are 120 different input conditions per phenotype (experiment, Pi condition, and synthetic community input). For each condition, there are 3 replicate measurements for plant shoot area, total root network length, and Pi content and around 20–30 replicates for primary root elongation. Let xc = [x1c, x2c, …, xrc] be the measurement of all replicates for a particular phenotype under condition c. We define the signal of a specific input condition as the mean measurement over all replicates in that condition.

90The signal variance is the variance of signals over all input conditions.

91We define the noise for a specific input condition as the residual between the mean and the actual measurement of the replicates:

92The noise variance is the variance of the noise over all input conditions and replicates (Table 3).

93We focused on shoot phosphate accumulation, which had the largest SNR of all phenotypes and serves as a bedrock test for our approach. Input data consist of a biological replicate ID b ∈ {1,2}, a technical replicate ID r ∈ {1,2,3}, pretreatment indicator p ∈ {−Pi, +Pi} and phosphate-level q ∈ {30 μM, 100 μM}, and a SynCom setting S ⊆ {P1,P2,P3,I1,I2,I3,N1,N2,N3}. We call the combination of p and q the phosphate condition. A design is a combination of phosphate condition and synthetic community. An input condition is a combination of design and biological replicate ID. Let zb,p,q,S,r be the standardized Pi-content measurement under biological replicate b, technical replicate r, pretreatment p, phosphate-level q, and SynCom setting S. The mean of Pi content across 3 technical replicates is and its variance is

94An input xb,p,q,S is constructed as a binary vector of length 12 (Table 4). A model learns a function for the prediction of mean Pi content.

95An LM has the following form: where β is the bias term. w is the linear coefficients for each features in x. If interaction features are considered, we have the following form: where ϴ is an upper triangular matrix for which diagonal entries are 0, and each entry ϴi,j indicates the linear coefficient for the interaction term xi xj. Compared to LM, INT is able to capture condition-specific influences. For example, a synthetic community can have different impact on Pi content under different phosphate conditions. Elastic net regularization is used in learning the parameters for both linear model and linear model with interaction features. An elastic net regularization of LM and INT yields the following optimization objective: where λ1, λ2 are the elastic net regularization penalty weights.

96We used a multilayer feed-forward NN, a typical framework in the deep NN structure family, in which input data are combined and transformed nonlinearly through multiple layers of artificial neurons. An NN contains an input layer, an output layer, and L hidden layers. The depth of the network is the number of hidden layers, and the width is the number of nodes in each hidden layer. For convenience, the input layer is defined as h0(x) = x, and the output of lth hidden layer is defined as hl(x). The number of nodes in layer l is ml. The input into the lth layer of the network is defined as where Wl is a real value weight matrix of ml-1 by ml, and βl is a vector of length ml-1. The output of lth hidden layer is where leaky Rectified Linear Unit activation function is:

97Finally, a linear output layer is on top of the last hidden layer:

98The prediction error is evaluated using sum of squares Error(W, β) = ∑t(yt − fNN(xt))2 regularized by weight decay.

99To train the NN, we minimized the regularized prediction error using RMSprop with momentum 0.9. Only structures with equal width in each layer were considered in this paper. The final network was selected by cross validation among the networks trained with the following hyper-parameter settings: depth of 1–4; width of 100, 200, 300, 400, 500; weight-decay parameter α ranging over 0.001, 0.005, 0.01, 0.05, 0.1, 0.5; and the number of training epochs of 100, 200, 300, 400, 500. The final network had depth of 3, width of 200, weight-decay parameter of 0.05, and was trained over 100 epochs.

100Next, we estimated the prediction error associated with our NN with a leave-SynCom-out cross-validation experiment; in this scheme, the model was trained on all but one synthetic community and evaluated on that held-out synthetic community. We calculated the mean Pi content prediction error on each held-out synthetic community as the cross-validation error for a single fold. Based on the cross-validation experiment, we chose the best NN model architecture, which had 3 hidden layers and 200 neurons in each layer, as shown in (Fig 7A).

101We note that there are more parameters than observations in the NN model. However, the number of parameters in a model is not a proper estimate of the complexity in the model. By choosing proper regularization parameters (early stopping, weight decay), the model can be trained without overfitting, even if it has more parameters than the number of observations. The cross-validation error in Fig 7B and the validation error in Fig 7F both indicate that the model is not overfitting, because the model is able to generalize to never-seen-before synthetic communities.

102Besides the leave-SynCom-out cross-validation experiment, we also performed a leave-biological-replicate-out cross-validation experiment. For this experiment, we used the data from 1 biological replicate to train each of the 3 models and used those trained models to predict the other biological replicate. The mean cross-validation error is shown in Table 5.

103Neural network methods are often regarded as “black-box” methods. However, Neural networks are not totally opaque and uninterpretable. There are many recent studies aimed to reveal and extract information from Neural networks. In order to extract biological insights from the NN, we performed a perturbation-based approach to quantify and understand the influence of specific features under different contexts.

104As our inputs are all binary features, we define sensitivity as the difference in the output between the “on” (1) and “off” (0) state of a particular feature. The sensitivity of the prediction with respect to each input context is where x(xi = 1) means change the ith feature in x to 1 and keep other features fixed.

105In a linear model, ρi(x) = wi. Hence, the sensitivity in a linear model is independent of the input context. In a linear model with interaction features, we have

106Therefore, a linear model with interaction features has an input context–dependent sensitivity. In NNs, the exact form of ρi(x) is complicated and hard to express in a closed form. In this paper, we calculate sensitivity ρi(x) numerically. Fig 7C shows the sensitivity of different models. For context-specific models, we consider all possible inputs with synthetic community of size |S| = 2.

107The code to estimate sensitivity of the different models is available at https://github.com/clingsz/wheelPi.

108We sought to identify cases in which replacing (swapping) one of the bacterial blocks (B) for a different block (C) in a reference synthetic community S1 = {A, B}, thereby producing a different synthetic community S2 = {A, C}, would produce a significant increase in shoot Pi content under a certain phosphate condition (pretreatment p phosphate-level q). We can use a trained model to estimate the mean and variance of the plant phenotypic output in every possible synthetic community. The expected Pi content for any input of interest can be calculated as f(xb,p,q,S). The worst-case variance estimate is used in our prediction, in which the largest residual variance (difference between observed value and predicted value) related to a bacterial block is transferred from the training data to all related conditions. Specifically,

109Six predicted samples are generated for every design:

110Given any 2 synthetic communities S1, S2, under a certain phosphate condition, p and q, the mean difference is and the p-value is calculated from a 2-sample t test on and.

111The code to generate candidate block replacements is available at https://github.com/clingsz/wheelPi.