The effects of soil phosphorous content on microbiota are driven by the plant phosphate starvation response

Phosphate starvation response (PSR) in non-mycorrhizal plants comprises transcriptional reprogramming resulting in severe physiological changes to the roots and shoots and repression of plant immunity. Thus, plant-colonizing microorganisms – the plant microbiota – are exposed to direct influence by the soil’s phosphorous (P) content itself, as well as to the indirect effects of soil P on the microbial niches shaped by the plant. The individual contribution of these factors to plant microbiota assembly remains unknown. To disentangle these direct and indirect effects, we planted PSR-deficient Arabidopsis mutants in a long-term managed soil P gradient, and compared the composition of their shoot and root microbiota to wild type plants across different P concentrations. PSR-deficiency had a larger effect on the composition of both bacterial and fungal plant-associated microbiota composition than P concentrations in both roots and shoots. The fungal microbiota was more sensitive to P concentrations per se than bacteria, and less depended on the soil community composition. Using a 185-member bacterial synthetic community (SynCom) across a wide P concentration gradient in an agar matrix, we demonstrated a shift in the effect of bacteria on the plant from a neutral or positive interaction to a negative one, as measured by rosette size. This phenotypic shift is accompanied by changes in microbiota composition: the genus Burkholderia is specifically enriched in plant tissue under P starvation. Through a community drop-out experiment, we demonstrate that in the absence of Burkholderia from the SynCom, plant shoots accumulate higher phosphate levels than shoots colonized with the full SynCom, only under P starvation, but not under P-replete conditions. Therefore, P-stressed plants allow colonization by latent opportunistic competitors found within their microbiome, thus exacerbating the plant’s P starvation.


Introduction 22
Plant-derived carbon is the primary energy source for terrestrial heterotrophs, most of which are 23 microbial. The interaction of these microbial heterotrophs with plants ranges between the extremes of 24 mutualistic symbiosis [1] and pathogenesis [2,3]. However, the vast majority of plant-associated microbial 25 diversity, the plant microbiota, lies between these two extremes, inducing more subtle, context-26 dependent effects on plant health [4][5][6]. The microbiota consumes plant photosynthate [7][8][9], and it 27 provides benefits via protection from pathogens [10][11][12][13][14] or abiotic stress [15,16] or by increasing nutrient 28 bioavailability [4,17,18]. 29 components of the plant immune system, which can lead to enhanced pathogen susceptibility but also to 48 the alteration of the plant's microbiota under phosphate starvation [4]. Arabidopsis microbiota are altered 49 in phr1 phl1 and phf1 mutants [4,36] in experiments using both natural and synthetic microbial 50 communities [4]. 51 Here, we examined (i) the effect of soil phosphorus (P) content on plant microbiota composition; (ii) how 52 PSR modulates the plant microbiota and (iii) the interplay between PSR and soil P content in shaping the 53 plant microbiota composition. We used a combination of greenhouse experiments with differentially P-54 fertilized soils, Arabidopsis PSR mutants and laboratory microcosms utilizing tractable synthetic bacterial 55 communities. Using PSR mutants planted in P-amended soil, we demonstrate that the plant PSR regulators 56 have a profound effect on the composition of root and shoot microbiota, overshadowing the effect of the 57 soil P content. We constructed an ecologically tractable system utilizing a complex bacterial synthetic 58 community (SynCom) as a model of the plant root microbiome and used this system to study the 59 interactions between microbiota assembly and abiotic stress. We demonstrate deterministic responses 60

Phosphate stress-induced changes in the root microbiome. 241
The shifting role of the SynCom from increasing shoot size under replete Pi to decreasing shoot size and 242 PSR induction under Pi limitation (Fig 5A and S5B Fig)  SynCom composition in wt plants, we measured alpha and beta diversity along our Pi gradient (0, 10, 30, 247 50, 100, 1000 µM KH2PO4) in roots, shoots and agar. We observed a positive correlation between alpha 248 diversity and Pi concentrations, resembling a partial ecological diversity-productivity relationship -the 249 prediction/observation of a bell-shaped response of ecological diversity to environmental productivity 250 [49,50] -in roots and shoots, but not in the surrounding agar ( Fig 5E). As for beta diversity, the 251 composition of the SynCom shifted significantly along the Pi concentration gradient ( for the plant, to a net-negative one, as measured by shoot size (Fig 5A). 254

Burkholderia respond to Pi stress-induced changes in the plant 255
In a previous publication [4], we demonstrated that PHR1 negatively regulates defense-related genes 256 under low-Pi conditions. Suppression of plant defense and consequent alterations in colonization could 257 account for some of the shift we observed from a beneficial to a detrimental community. We thus aimed 258 to identify bacteria that respond to Pi stress-induced changes in the plant, rather than the Pi concentration 259 itself. To do so, we searched for USeqs that displayed a strong Pi:fraction (shoot, root, agar) interaction in 260 our GLM (S7A Fig, S11 Table and Materials  to plants inoculated with a SynCom excluding all five Burkholderia isolates. We also included a SynCom 270 11 excluding all members of the neighboring Ralstonia clade (Fig 4A), which didn't display any discernible Pi-271 response. We measured Pi concentrations in the shoots (a proxy for PSR) of plants grown in high (1000 272 µM) and low (50 µM) Pi with the different SynComs. In addition, we measured shoot Pi in a re-feeding 273 treatment with SynCom-inoculated plants grown in low (50 µM) Pi and then transferred to high Pi (1000 274 µM) conditions. All SynCom treatments decreased shoot Pi content in the low Pi conditions compared to 275 the uninoculated plants but recovered to a higher shoot Pi level than the uninoculated treatments upon 276 transferring to high Pi conditions, reproducing our previous report [4] (Fig 6B). Among  Despite the fact that phosphate is a critical nutrient for plants and their microbiota, differences in 287 phosphate content have relatively subtle effects on plant and soil microbiome compositions compared to 288 abiotic factors like pH or drought, which cause pronounced, phylogenetically consistent changes in 289 community configurations [19,27,29]. Several studies link host physiological response to the soil 290 phosphate status with the bacterial [4,51] and fungal [34,36] microbiome. A recent report of Arabidopsis 291 planted in a 60-year-long annual phosphorus fertilization gradient (the same soil used in the current study) 292 showed a modest P effect on plant microbiome composition [41]. Previously, we showed that PSR mutants 293 in Arabidopsis have subtly different bacterial microbiomes in Pi replete [4] conditions and a recent 294 publication showed that PSR mutants had a slightly altered fungal microbiome in Pi replete but not in Pi 295 depleted conditions [36]. 296 Here, we analyzed fungi and bacteria side by side and demonstrated a pronounced effect of PSR 297 impairment on both bacterial and fungal components of the plant microbiota. We noted an intriguing 298 difference that emerged in the patterns of niche sorting between bacteria and fungi (S3A S3B, S3D and 299 S3E Fig). The bacterial microbiota composition is strongly dependent on the soil bacterial community 300 composition, whereas changes to the fungal microbiota are uncoupled from changes to the soil fungal 12 community composition. This indicates that the plant is markedly more selective as to the fungi allowed 302 to proliferate in its tissue than it is with bacteria. Similar to [36], amendment of the soil with P at the time 303 of the experiment (low P vs low+P) caused a shift in the microbial community, albeit weaker than the 304 effect of knocking-out PSR genes. Our results show that impairment of PSR genes profoundly affects the 305 composition of the plant microbiota, under a range of P conditions, and that observed shifts in root-306 derived microbial communities may not be a result of sensitivity to P concentrations, but rather a response 307 to PSR regulation in the hosts. This raises the alternative hypotheses that PSR-regulated shifts in 308 microbiota composition are either adaptive to the plant, or reflect opportunistic strategies by bacteria, 309 exploiting the repression of immunity by PSR regulation [4,17]. Under the former hypothesis microbes 310 recruited by the plant under Pi stress provide the plants with an advantage vis a vis coping with this stress, 311 whereas under the latter, opportunistic microbes might be making a bad situation worse for the plant. In 312 the case of Burkholderia in our SynCom, results support the latter hypothesis. Burkholderia contribute to 313 depletion of shoot Pi stores, only under Pi-limiting conditions. However, plant-adaptive microbial 314 recruitment under low Pi has been shown to occur as well [17]. The fact that in soil bacteria responding 315 to PSR genes are not a monophyletic group indicates that multiple mechanisms are involved. It is likely 316 that these mechanisms encompass both plant-adaptive and opportunistic strategies. 317 The genus Burkholderia emerges as a PSR-responsive taxon. We examined the effect of Burkholderia on 318 shoot Pi accumulation, which we've shown to be a reliable marker for PSR (S1F Fig). We compared the 319 effect of Burkholderia on shoot Pi accumulation from within a full SynCom (a realistic proxy for the 320 bacterial community) to that of the full SynCom lacking Burkholderia, a strategy akin to knocking-out a 321 gene of interest, also recently applied in [52]. The control treatment for this type of approach is the full 322 SynCom, while in a plant-bacterium binary association experiment it would typically be sterile conditions. 323 As both sterile conditions and binary association are strong deviations from conditions that may be 324 encountered in the field, the results of binary association experiments may be correspondingly distorted. 325 Using the drop-out approach, we expect to see more subtle differences, as the microbial load on the plant 326 doesn't change much, but also that these differences be more relevant to the field; an expectation that is 327 yet to be empirically tested. Our observation that dropping Burkholederia out of the SynCom increased 328 shoot Pi in Pi limiting conditions (50 µM Pi) but not in Pi replete conditions (1000 µM Pi) suggests that 329 strains in this genus shift their relationship with the plant from a seeming commensal to a 330 competitor/pathogen. It is likely that this shift is related to the repression of plant immune function by 331 key regulators of PSR in low Pi [4], suggesting a specific plant-dependent trade off during PSR. 332 13 This study shows that despite 60 years of differential fertilization, differences in PSR and in microbiome 333 composition between the low P and high P soils are subtle, possibly because Pi status is highly buffered 334 by the plant ionomic regulatory network [53]. Only when comparing the low P vs the P-supplemented 335 low+P samples, is there a discernible difference in PSR (Fig 1A), which correlates to a stronger effect on 336 microbiota composition. This suggests that bioavailable Pi added to the soil is quickly consumed, and 337 short-term amendments are needed in order to detect changes. It is easier to produce Pi-limiting 338 conditions in vitro using defined media, as evidenced from shifts in both PSR gene expression and 339 microbiome composition in the microcosm system introduced here. Our SynCom, comprising 185 340 genome-sequenced endophytic bacterial isolates, was designed to resemble a natural bacterial 341 community ( Fig 4B). The community assembly patterns shown for this system are highly reproducible, 342 demonstrating that microbiome assembly is largely a deterministic process. The reproducibility and 343 editability of this system can be used for detailed mechanistic study of the processes that determine 344 community assembly and its influence on plant phenotype and fitness. 345 Materials and Methods 11°59′33.3′′E) [40,54]. Soil cores (10 cm diameter × 15 cm depth) were taken from 18 6 X 5 m unplanted 351 plots, belonging to two strips. These plots represent three P fertilization regimens: low, medium and high 352 P (0, 15 and 45 kg P ha −1 year −1 , respectively). Strips were harvested independently in the middle of March 353 (strip 1) and beginning of April (strip 2). Approximately 2 cm of the topsoil was discarded and the 354 remaining lower 13 cm of soil was stored at 4 °C until use. Soils were homogenized with a mesh sieve wire 355 (5 × 5 m 2 ) and about 300 g of soil were added to each pot (7 × 7 × 7 cm 3 ). 356

b. Experimental design 357
Each of the three Arabidopsis genotypes was grown in soil from all 18 plots (6 plots per P treatment). In 358 addition, a 4 th P regimen designated 'Low+P' was created by adding additional P to a set of pots with low 359 P. The amount of P added to these pots is based on the difference in total P between Low and High P 360 14 plots. The average difference between Low and High P over all the plots is 42 mg P per kg soil [41]. Per 361 pot, this is 12.6 mg P (accounting for 300 g soil per pot). Thus, a 10 ml solution consisting of 4.2 mg P in 362 the form of 20% K2HPO4 and 80% KH2PO4 was added to the pots in 3 applications (Week 2, 4 and 6) before 363 watering (in order to distribute the P through the soil). 364 Thus, the experiment included four soil treatments (low P, medium P, high P, low+P) and three genotypes 365 After eight weeks of growth, pots were photographed, and shoot size was quantified using WinRhizo 375 software (Regent instruments Inc. Québec, Canada). For DNA extraction, two roots, two shoots and soil 376 from each pot were harvested separately. Roots and shoots were rinsed in sterile water to remove soil 377 particles, placed in 2 ml Eppendorf tubes with 3 sterile glass beads, then washed three times with sterile 378 distilled water to remove soil particles and weakly associated microbes. Root and shoot tissue were then 379 pulverized using a tissue homogenizer (TissueLyser II; Qiagen) and stored at ˗80 ˚C until processing. Five 380 ml of soil from each pot was suspended in 20 ml of sterile distilled water. The resulting slurry was sieved 381 through a 100 µm sterile cell strainer (Fisher Scientific) and the flow-through was centrifuged twice at 382 maximum speed for 20 minutes, removing the supernatant both times. The resulting pellet was stored at 383 ˗80 ˚C until processing. For RNA extraction, one root system and one shoot were taken from three 384 replicates of each treatment, washed lightly to remove soil particles, placed in 2 ml Eppendorf tubes with 385 three glass beads and flash frozen with liquid nitrogen. Tubes were stored at ˗80 ˚C until processing. which are more compatible with our experimental system, in our SynCom. A detailed description of this 407 collection and isolation procedures can be found in [47]. One week prior to each experiment, bacteria 408 were inoculated from glycerol stocks into 600 µL KB medium in a 96 deep well plate. Bacterial cultures 409 were grown at 28°C, shaking at 250 rpm. After five days of growth, cultures were inoculated into fresh 410 media and returned to the incubator for an additional 48 hours, resulting in two copies of each culture, 7 411 days old and 48 hours old. We adopted this procedure to account for variable growth rates of different 412 SynCom members and to ensure that non-stationary cells from each strain were included in the inoculum. 413 16 After growth, 48-hour old and 7-day old plates were combined and optical density (OD) of the culture was 414 measured at 600 nm using an Infinite M200 Pro plate reader (TECAN, Switzerland). All cultures were then 415 pooled while normalizing the volume of each culture according to the OD (we took a proportionally higher 416 volume of culture from cultures with low OD). The mixed culture was then washed twice with 10 mM 417 MgCl2 to remove spent media and cell debris and vortexed vigorously with sterile glass beads to break up 418 aggregates. OD of the mixed, washed culture was then measured and normalized to OD=0.2. 100 µL of 419 this SynCom inoculum was spread on each agar plate prior to transferring seedlings. We amplified the V3-V4 regions of the bacterial 16S rRNA gene using primers 338F (5′-ACTCCTACGGGAG 488 GCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). Two barcodes and 6 frameshifts were added 489 to the 5' end of 338F and 6 frameshifts were added to the 806R primers, based on the protocol in [58].  product was indexed using 96 indexed 806R primers with the same reaction mix as above, and 9 cycles of 497 the cycling conditions described in [29]. PCR products were purified using AMPure XP magnetic beads 498 (Beckman Coulter) and quantified with a Qubit 2.0 fluorometer (Invitrogen). Amplicons were pooled in 499 equal amounts and then diluted to 10 pM for sequencing. Sequencing was performed on an Illumina 500 MiSeq instrument using a 600-cycle V3 chemistry kit. The raw data for the natural soil experiment is 501 available in the NCBI SRA Sequence Read Archive (accession XXXXXXX). The raw data for the SynCom 502 experiment is available in the NCBI SRA Sequence Read Archive (accession XXXXXXX). 503 b. Fungal/Oomycete ITS sequencing 504 We amplified the ITS1 region using primers ITS1-F (5′-CTTGGTCATTTAGAGGAAGTAA-3′; [59]) and ITS2 (5′-505 GCTGCGTTCTTCATCGATGC-3′; [60]). Samples were diluted to concentrations of 3. Illumina-based mRNA-Seq libraries were prepared from 1 μg RNA following [51]. mRNA was purified from 524 total RNA using Sera-mag oligo(dT) magnetic beads (GE Healthcare Life Sciences) and then fragmented in 525 the presence of divalent cations (Mg 2+ ) at 94 °C for 6 minutes. The resulting fragmented mRNA was used 526 for first-strand cDNA synthesis using random hexamers and reverse transcriptase, followed by second-527 strand cDNA synthesis using DNA Polymerase I and RNAseH. Double-stranded cDNA was end-repaired 528 using T4 DNA polymerase, T4 polynucleotide kinase, and Klenow polymerase. The DNA fragments were 529 then adenylated using Klenow exo-polymerase to allow the ligation of Illumina Truseq HT adapters (D501-530 D508 and D701-D712). All enzymes were purchased from Enzymatics. Following library preparation, 531 quality control and quantification were performed using a 2100 Bioanalyzer instrument (Agilent) and the 532 Quant-iT PicoGreen dsDNA Reagent (Invitrogen), respectively. Libraries were sequenced using Illumina 533 HiSeq4000 sequencers to generate 50-bp single-end reads. 534

Data processing and Statistical analyses 535 a. Quantification of plant phenotypes -soil experiment 536
To measure correlation between all measured plant phenotypes (shoot Pi, shoot weight, shoot size) we 537 applied hierarchical clustering based on the all vs all pairwise correlation coefficients between all the 538 phenotypes measured. We used the R package corrplot v.0.84 [61] to visualize correlations. To compare 539 shoot Pi accumulation, we performed a paired t-test between low P and P-supplemented low P (low+P) 540 samples, within each plant genotype independently (α < 0.05). 541 b. Amplicon sequence data processing -soil experiments 542 Bacterial sequencing data was processed with MT-Toolbox [62]. Usable read output from MT-Toolbox 543 (that is, reads with 100% correct primer and primer sequences that successfully merged with their pair) 544 were quality filtered using Sickle [63] by not allowing any window with a Q-score under 20. After quality 545 filtering, samples with < 3000 reads, amounting to 51 samples, all soil samples, were discarded. The 546 resulting sequences were collapsed into amplicon sequence variants (ASVs) using the R package DADA2 547 v1.8.1 [64]. Taxonomic assignment of each ASV was performed using the naïve Bayes kmer method 548 implemented in the DADA2 package using the Silva 132 database as training reference [64]. 549 21 Fungal sequencing data was processed as previously described [41]. Briefly, a combination of QIIME [65] 550 and USEARCH [66] pipelines were used to cluster the fungal reads into 97% OTUs. Filtering of non-fungal 551 OTUs was performed by aligning each representative against a dedicated ITS database. Finally, taxonomic 552 assignment of each OTU was performed using the WarCup fungal ITS training set (2016) [67]. 553 The resulting bacterial and fungal count tables were deposited at https://github.com/isaisg/hallepi 554 c. Community analyses -soil experiments 555 The resulting bacterial and fungal count tables were processed and analyzed with functions from the 556 ohchibi package [68]. Both tables were rarefied to 3000 reads per sample. An alpha diversity metric 557 (Shannon diversity) was calculated using the diversity function from the vegan package v2.5-3 [69]. We 558 used ANOVA to test for differences in Shannon Diversity indices between groups. Tukey's HSD post-hoc 559 tests here and elsewhere were performed using the cld function from the emmeans R package [70]. Beta 560 diversity analyses (Principal coordinate analysis, and canonical analysis of principal coordinates) were 561 based on Bray-Curtis dissimilarity calculated from the rarefied abundance tables. We utilized the capscale 562 function from the vegan R package v.2.5-3 [69] to compute a constrained analysis of principal coordinates 563 (CAP). To analyze the full dataset (all fraction, all genotypes all phosphorus treatments), we constrained 564 by fraction, plant genotype and phosphorus fertilization treatment, while conditioning for the plot effect. 565 We performed the Genotype: phosphorus interaction analysis over each fraction independently, 566 constraining for the plant genotype and phosphorus fertilization treatment while conditioning for the plot 567 effect. In addition to CAP, we performed Permutational Multivariate Analysis of Variance (PERMANOVA) 568 over the two datasets described above using the adonis function from the vegan package v2.5-3 [69]. 569 Finally, we used the function chibi.permanova from the ohchibi package to plot the R 2 values for each 570 significant term in the PERMANOVA model tested. 571 The relative abundance of bacterial phyla and fungal taxa were depicted using the stacked bar 572 representation encoded in the function chibi.phylogram from the ohchibi package. 573 We used the R package DESeq2 v1.22.1 [71] to compute the enrichment profiles for both bacterial ASVs 574 and fungal OTUs. For the full dataset model, we estimated main effects for each variable tested (Fraction, 575 Plant Genotype, and Phosphorus fertilization) using the following design: 576 Abundance ~ Fraction + Genotype + Phosphorus Treatment 577 22 We delimited ASV/OTU fraction enrichments using the following contrasts: Soil vs Root, Soil vs Shoot and 578 Root vs Shoot. An ASV/OTU was considered statistically significant if it had q-value < 0.1. 579 We implemented a second statistical model in order to identify ASVs and OTUs that exhibited statistically 580 significant differential abundances depending on genotype. For this analysis we utilized only root-derived 581 low P and P-supplemented low P (low+P) treatments. We utilized a group design framework to facilitate 582 the construction of specific contrasts. In the group variable we created, we merged the genotype and 583 phosphate levels per sample (e.g. Col-0_lowP, phf1_low+P or phr1 phl1_lowP). We controlled the paired 584 structure of our design by adding a plot variable, resulting in the following model design: 585 Abundance ~ Plot + group 586 We delimited 6 sets (S1, S2, S3, S4, S5, S6) of statistically significant (q-value < 0.1) ASVs/OTUs from our 587 model using the following contrasts: 588 S1 = {Samples from Col-0, higher abundance in low treatment in comparison to low+P treatment} 589 S2 = {Samples from phf1, higher abundance in low treatment in comparison to low+P treatment} 590 S3 = {Samples from phr1 phl1, higher abundance in low treatment in comparison to low+P treatment} 591 S4 = {Samples from Col-0, higher abundance in low+P treatment in comparison to low treatment} 592 S5 = {Samples from phf1, higher abundance in low+P treatment in comparison to low treatment} 593 S6 = {Samples from phr1 phl1, higher abundance in low+P treatment in comparison to low treatment} 594 The six sets described above were used to populate Figures 3c-f. 595 The interactive visualization of the enrichment profiles was performed by converting the taxonomic 596 assignment of each ASV/OTU into a cladogram with equidistant branch lengths using R. We used the 597 interactive tree of life interface (iTOL) [72] to visualize this tree jointly with metadata files derived from 598 the output of the statistical models described above. The cladograms for both bacteria and fungi can be 599 downloaded from the links described above or via the iTOL user hallepi. 600 In order to compare beta diversity patterns across samples, we only used samples coming from pots 601 where sequence data from all three fractions (soil root and shoot) passed quality filtering. Then, for each 602 fraction we estimated a distance structure between samples inside that fraction using the Bray Curtis 603 dissimilarity metric. Finally, we computed Mantel [73] correlations between pairs of distance objects (e.g. 604 samples from Root or samples from Shoot) using the vegan package v2.5-3 [69]  (that is, reads with 100% correct primer and primer sequences that successfully merged with their pair) 611 were quality filtered using Sickle [63] by not allowing any window with Q-score under 20. The resulting 612 sequences were globally aligned to a reference set of 16S rRNA gene sequences extracted from genome 613 assemblies of SynCom member strains. For strains that did not have an intact 16S rRNA gene sequence in 614 their assembly, we generated the 16S rRNA gene using Sanger sequencing. The reference database also 615 included sequences from known bacterial contaminants and Arabidopsis organellar 16S sequences (S12 616 The resulting count table was processed and analyzed with functions from the ohchibi package. Samples 629 were rarefied to 1000 reads per sample. An alpha diversity metric (Shannon diversity) was calculated using 630 the diversity function from the vegan package v2.5-3 [69]. We used ANOVA to test for differences in alpha 631 diversity between groups. Beta diversity analyses (Principal coordinate analysis, and canonical analysis of 632 principal coordinates) were based on were based on Bray-Curtis dissimilarity calculated from the rarefied 633 abundance tables. We used the capscale function from the vegan R package v.2.5-3 [69] to compute the 634 canonical analysis of principal coordinates (CAP). To analyze the full dataset (all fraction, all phosphate 635 24 treatments), we constrained by fraction and phosphate concentration while conditioning for the replicate 636 effect. We performed the Fraction:Phosphate interaction analysis within each fraction independently, 637 constraining for the phosphate concentration while conditioning for the rep effect. In addition to CAP, we 638 performed Permutational Multivariate Analysis of Variance (PERMANOVA) analysis over the two datasets 639 described above using the adonis function from the vegan package v2.5-3 [69]. Finally, we used the 640 function chibi.permanova from the ohchibi package to plot the R 2 values for each significant term in the 641 PERMANOVA model tested. 642 We visualized the relative abundance of the bacterial phyla present in the SynCom using the stacked bar 643 representation encoded in the chibi.phylogram from the ohchibi package. 644 We used the package DESeq2 v1. 22 We calculated the USeqs/OTUs fraction enrichments using the following contrasts: Agar vs Root, Agar vs 649 Shoot and Root vs Shoot. A USeq/OTU was considered statistically significant if it had q-value < 0.1. In 650 order to populate the heatmaps shown in Figure 5C, we grouped the Fraction and Phosphate treatment 651 variable into a new group variable that allowed us to fit the following model: 652 Abundance ~ Replicate + group 653 We used the fitted model to estimate the fraction effect inside each particular phosphate level (e.g. Root 654 vs Agar at 0Pi, or Shoot vs Agar at 1000Pi). 655 Additionally, we utilized a third model for the identification of USeqs/OTUs that exhibited a significant 656 Fraction:Phosphate interaction between the planted agar samples and the plant fractions (Root and 657 Shoot). Based on the beta diversity and alpha diversity results, we only used samples that were treated 658 with 0, 10, 100 and 1000 µM of phosphate. We grouped the samples into two categories based on their 659 phosphate concentration, low (0 µM and 10 µM) and high (100 µM and 1000 µM). Then we used the 660 following model specification to derive the desired interaction effect: 661 Abundance ~ Fraction + Category + Fraction:Category + Replicate 662 Finally, we subset USeqs that exhibited a significant interaction (Fraction:Category, q-value < 0.1) in the 663 following two contrasts (Planted Agar vs Root) and (Planted Agar vs Shoot). 664 25 In order to compare beta diversity patterns across samples, we only used samples coming from pots 665 where sequence data from all three fractions (soil root and shoot) passed quality filtering. Then, for each 666 fraction we estimated a distance structure between samples inside that fraction using the Bray Curtis 667 dissimilarity metric. Finally, we computed Mantel [73] correlations between pairs of distance objects (e.g. 668 samples from Root or samples from Shoot) using the vegan package v2.5-3 [69] implementation of the 669

Mantel test. 670
For the drop-out experiment, we ran an ANOVA model inside each of the phosphate treatments tested 671 (50 µM Pi, 1000 µM Pi and 501000 µM Pi). We visualized the results of the ANOVA models using the 672 compact letter display encoded in the CLD function from the emmeans package. 673 All scripts necessary to reproduce the synthetic community analyses are deposited in the following GitHub 674 repository: https://github.com/isaisg/hallepi. 675

e. Phylogenetic inference of the SynCom Isolates 676
To build the phylogenetic tree of the SynCom isolates, we utilized the super matrix approach previously 677 described in [47]. Briefly, we scanned 120 previously defined marker genes across the 185 isolate genomes 678 from the SynCom utilizing the hmmsearch tool from the hmmer v3.1b2 [77]. Then, we selected 47 markers 679 that were present as single copy genes in 100% of our isolates. Next, we aligned each individual marker 680 using MAFFT [78] and filtered low quality columns in the alignment using trimAl [79]. Afterwards, we 681 concatenated all filtered alignments into a super alignment. Finally FastTree v2.1 [80] was used to infer 682 the phylogeny utilizing the WAG model of evolution. 683 We utilized the inferred phylogeny along with the fraction fold change results of the main effect model to 684 compute the phylogenetic signal (Pagel's λ) [ Initial quality assessment of the Illumina RNA-seq reads was performed using FastQC v0.11.7 [84]. 690 Trimmomatic v0.36 [85] was used to identify and discard reads containing the Illumina adaptor sequence. 691 The resulting high-quality reads were then mapped against the TAIR10 [86] Arabidopsis reference genome 692 using HISAT2 v2.1.0 [87]with default parameters. The featureCounts function from the Subread package 693 [88] was then used to count reads that mapped to each one of the 27,206 nuclear protein-coding genes. to define differentially expressed genes (DEGs) using the raw count table described above (Materials and 699 Methods 4f). We used only samples from low P and P-supplemented low P (low+P) treatments along the 700 three genotypes tested (Col-0, phf1 and phr1 phl1). We combined the Genotype and Phosphorus 701 Treatment variables into a new group variable (e.g. Col-0_lowP or phf1_low+P). Because we were 702 interested in identifying DEGs among any pair of levels (6 levels) of the group variable (e.g. Col-0_lowP vs 703 Col-0_low+P) we performed a likelihood ratio test (LRT) between a model containing the group variable 704 and a reduced model containing only the intercept. Next, we defined DEGs as genes that had a q-value < 705 0.1. 706 For visualization purposes, we applied a variance stabilizing transformation to the raw count gene matrix. 707 We then standardized (z-score) each gene along the samples. We subset DEGs from this standardized 708 matrix and for each gene calculated the mean z-score expression value in a particular level of the group 709 variable (e.g. Col-0_lowP); this resulted in a matrix of DEGs across the six levels in our design. Next, we 710 created a dendrogram of DEGs by applying hierarchical clustering (method ward.D2, hclust R-base [90]) 711 to a distance object based on the correlation (dissimilarity) of the expression profiles of the genes across 712 the six levels in our design. Finally, we delimited the cluster of DEGs by cutting the output dendogram into 713 five groups using the R-base cutree function [90]. Gene ontology enrichment was performed for each 714 cluster of DEGs using the R package clusterProfiler [91]. 715 For the PSR marker gene analysis we downloaded the ID of 193 genes defined in [4]. Then, we subset 716 these genes from our standardized matrix and computed for each gene the mean z-score expression value 717 27 in a particular level of the group variable. Finally, we visualized the average expression of this PSR regulon 718 across our groups of interest utilizing the function chibi.boxplot from the ohchibi package. 719 All scripts necessary to reproduce the RNA-Seq analyses are deposited in the following GitHub repository: 720 https://github.com/isaisg/hallepi.  (Materials and methods 4b). The arrows on the bo�om of the panel denote the direc�on of the enrichment rela�ve to the name of the contrast tested, the up arrow signifies enrichment in the le� frac�on of the contrast, whereas the down arrow signifies enrichment in the right frac�on of the contrast (e.g. RootvsSoil, up arrow enriched in root rela�ve to soil, bo�om arrow enriched in soil rela�ve to root). A detailed interac�ve visualiza�on of the fungal enrichment pa�erns across the mul�ple taxonomic levels can be found at (h�ps://itol.embl.de/tree/1522316254174721551987262).