|Home | About | Journals | Submit | Contact Us | Français|
NB – NewLeaf Symbiotics, St. Louis, MO 63132 USA
LdL -- Department of Plant and Soil Sciences, University of Kentucky, Lexington, KY 40546 USA.
Plants live in biogeochemically diverse soils that harbor extraordinarily diverse microbiota. Plant organs associate intimately with a subset of these microbes; this community’s structure can be altered by soil nutrient content. Plant-associated microbes can compete with the plant and with each other for nutrients; they can also provide traits that increase plant productivity. It is unknown how the plant immune system coordinates microbial recognition with nutritional cues during microbiome assembly. We establish that a genetic network controlling phosphate stress response influences root microbiome community structure, even under non-stress phosphate conditions. We define a molecular mechanism regulating coordination between nutrition and defense in the presence of a synthetic bacterial community. We demonstrate that the master transcriptional regulators of phosphate stress response in Arabidopsis also directly repress defense, consistent with plant prioritization of nutritional stress over defense. Our work will impact efforts to define and deploy useful microbes to enhance plant performance.
Plant organs create distinct physical and chemical environments that are colonized by specific microbial taxa1. These can be modulated by the plant immune system2 and by soil nutrient composition3. Phosphorus is present in the biosphere at high concentrations, but plants can only absorb orthophosphate (Pi), a form also essential for microbial proliferation4, 5 and scarce in soil6. Thus, plants possess adaptive phosphate starvation responses (PSR) to manage low Pi availability that typically occurs in the presence of plant-associated microbes. Common strategies for increasing Pi uptake capacity include rapid extension of lateral roots foraging into topsoil where Pi accumulates7 and establishment of beneficial relationships with some soil microorganisms8, 9. For example, the capacity of a specific mutualistic fungus to colonize Arabidopsis roots is modulated by plant phosphate status implying coordination between the PSR and the immune system8, 10. Descriptions of pathogen exploitation of PSR-immune system coordination are emerging11, 12.
We demonstrate that Arabidopsis mutants with altered phosphate starvation responses (PSR) assemble atypical microbiomes, either in phosphate-replete wild soil, or during in vitro colonization with a synthetic bacterial community (SynCom). This SynCom competes for phosphate with the plant and induces PSR in limiting phosphate. PSR in these conditions requires the master transcriptional regulator PHR1 and its weakly redundant paralog, PHL1. The severely reduced PSR observed in phr1 phl1 mutants is accompanied by transcriptional changes in plant defense leading to enhanced immune function. Negative regulation of immune system components by PHR1 is direct, as measured by target gene promoter occupancy, and functional, as validated by pathology phenotypes. Thus, PHR1 directly activates microbiome-enhanced response to phosphate limitation while repressing microbially-driven plant immune system outputs.
We linked PSR to the root microbiome by contrasting the root bacterial community of wild-type Arabidopsis Col-0 with three types of PSR mutants (Fig. 1a, b, Supplementary Text 1; Extended Data Fig. 1; Supplementary Table 1). PSR, historically defined in axenic seedlings and measured by Pi concentration in the plant shoot, is variable across these mutants. In replete Pi and axenic conditions, phr1 plants accumulate less free Pi than wild-type13; pht1;1, pht1;1 pht1;4 and phf1 accumulate very low Pi levels and express constitutive PSR14, 15; and pho2, nla and spx1 spx2 express diverse magnitudes of Pi hyper-accumulation16– 18. We grew plants in a previously characterized wild soil19 that is not overtly phosphate deficient (Extended Data Fig. 2). Generally, the Pi concentration of PSR mutants grown in this wild soil recapitulated those defined in axenic conditions, except for phf1 and nla which displayed the opposite phenotype to that observed in axenic agar, and phr1 which accumulated the same Pi concentration as Col-0 (Fig. 1b). These results suggest that complex chemical conditions, soil microbes, or a combination of these can alter Pi metabolism in these mutants.
Bacterial root endophytic (EC) community profiles were consistent with previous studies2, 19. Constrained ordination revealed significant differences between bacterial communities across the Pi accumulation gradient represented by these PSR mutants [5.3 % constrained variance, canonical analysis of principal coordinates (CAP)] (Fig. 1c). Additionally, CAP confirmed that phr1 and spx1 spx2 carried different communities, as evidenced by their separation on the first three ordination axes, and that phf1 was the most affected of Pi-transport mutant (Fig. 1c). Specific bacterial taxa had differential abundances inside the roots of mutant plants compared to wild-type. Mutants from the same PSR type had a similar effect on the root microbiome at a low taxonomic level [97 % identity Operational Taxonomic Unit (OTU)] (Fig. 1d), while they had no overlapping effect at a higher taxonomic level (Family, Extended Data Fig. 1g). This suggests that closely related groups of bacteria have differential colonization patterns on the same host genotypes. Importantly, we found that the enrichment and depletion profiles were better explained by PSR mutant signaling type rather than the mutant’s capacity for Pi accumulation: all of the Pi-transport-related mutants had a similar effect on the root microbiome, and the antagonistic PSR regulators phr1 and spx1 spx2 each exhibited unique patterns (Fig. 1a, d and Extended Data Fig. 1f, g). Our results indicate that PSR components influence root microbiome composition in plants grown in a phosphate-replete wild soil, leading to alteration of the abundance of specific microbes across diverse levels of Pi accumulation representing diverse magnitudes of PSR.
Our observations in a wild soil suggested complex interplay between PSR and the presence of a microbial community. Thus, we deployed a tractable but complex bacterial synthetic community (SynCom) of 35 taxonomically diverse, genome-sequenced bacteria isolated from the roots of Brassicaceae (nearly all from Arabidopsis) and two wild soils. This SynCom approximates the phylum level distribution observed in wild-type root endophytic compartments (Extended Data Fig. 3, Supplementary Table 1, Supplementary Table 2). We inoculated seedlings of Col-0, phf1 and the double mutant phr1 phl1 (a redundant paralogue of phr113) grown on agar plates in low or high Pi (Supplementary Text 2). Twelve days later, we noted that the SynCom had a negative effect on shoot Pi accumulation of Col-0 plants grown on low Pi, but not on plants grown on replete phosphate (Fig. 2a). As expected, both PSR mutants accumulated less Pi than Col-0; the SynCom did not rescue this defect. Thus, in this microcosm, plant-associated microbes drive a context-dependent competition with the plant for Pi.
We sought to establish whether PSR was activated by the SynCom. We generated a literature-based core set of 193 PSR transcriptional markers and explored their expression in transcriptomic experiments (Extended Data Fig. 4a, b, Supplementary Table 3). In axenic low Pi conditions, only the constitutive Pi-stressed mutant phf1 exhibited induction of these PSR markers. By contrast, Col-0 plants expressed only a marginal induction of PSR markers compared to those plants grown at high Pi (Fig. 2b). This is explained by the purposeful absence of sucrose, a key component for the PSR induction in vitro20 (Supplementary Text 2; Extended Data Fig. 5) that cannot be used in combination with bacterial SynCom colonization protocols. Remarkably, the SynCom greatly enhanced the canonical transcriptional response to Pi starvation in Col-0 (Fig. 2b); this was dependent on PHR1 and PHL1 (Fig. 2b; Extended Data Fig. 4b). Various controls validated these conclusions (Supplementary Text 2; Extended Data Fig. 4–Fig. 6). Importantly, shoots of plants pre-colonized with SynCom on 0 or 50 µM Pi, but not on 650 µM Pi, accumulated 20–40 times more Pi than shoots from similarly treated non-colonized plants when subsequently transferred to full Pi conditions in the absence of additional bacteria (Fig. 2c and Supplementary Table 4). This demonstrates functional PSR activation by the SynCom. We thus propose that the transcriptional response to low Pi induced by our SynCom reflects an integral microbial element of normal PSR in complex biotic environments.
We evaluated agar- and root-associated microbiomes of plants grown with the SynCom (Supplementary Text 3; Fig. 2d, e; Extended Data Fig. 7e, f; Supplementary Table 5). In line with results from plants grown in wild soil, we found that PSR mutants failed to assemble a wild-type SynCom microbiome (Fig. 2f). Some strains were differentially abundant across PSR mutants phf1 and phr1 phl1 (Fig. 2e, f; Extended Data Fig. 7c), Pi concentration (Fig. 2g; Extended Data Fig. 7d), or sample fraction (Extended Data Fig. 7b, e, f). These results established a microcosm reconstitution system to study plant PSR under chronic competition with plant-associated microbes and allowed us to confirm that the tested PSR mutants influence root microbiome membership.
We noted that phr1 phl1 and phf1 differentially activated transcriptional PSR in the presence of our SynCom (Fig. 2b). Therefore, we investigated the transcriptomes of plants growing in the SynCom to understand how these microbes activate PHR1-dependent PSR. We identified differentially expressed genes (DEGs) that responded to either low Pi, presence of the SynCom, or the combination of both (hereafter PSR-SynCom DEGs) (Supplementary Text 4; Extended Data Fig. 8a, b; Supplementary Table 6). Hierarchical clustering (Fig. 3a, Supplementary Table 7) revealed gene sets (c1, c2, c7 and c10) that were more strongly activated in Col-0 than in phr1 or phr1 phl1. These clusters contained most of the core PSR markers regulated by PHR1 (Fig. 3b). They were also enriched in PHR1 direct targets identified in an independent ChIP-seq experiment (Fig. 3c, Supplementary Table 8), PHR1 promoter binding motifs (Extended Data Fig. 4c), and genes involved in biological processes related to PSR (Fig. 3d and Supplementary Table 9). PHR1 surprisingly contributed to transcriptional regulation of plant immunity. Five of the twelve clusters (Fig. 3a; c3, c6, c7, c8 and c11) were enriched in genes related to plant immune system output; four of these were over-represented for jasmonic acid (JA) and/or salicylic acid (SA) pathway markers (Fig. 3d, c3, c7, c8, and c11; Supplementary Table 9) and three of these four were enriched for PHR1 direct targets (Fig. 3c). SA and JA are plant hormone regulators of immunity and at least SA modulates Arabidopsis root microbiome composition2.
To further explore PHR1 function in the regulation of plant immunity, we generated transcriptomic time course data for treatment-matched Col-0 seedlings following application of Methyl Jasmonate (MeJA) or the SA analogue Benzothiadiazole (BTH; Supplementary Table 10). We found a considerable over-representation of SA- and JA-activated genes among the PSR-SynCom DEGs (468 versus 251 predicted for SA and 165 vs. 80 predicted for JA; p<0.0001, hypergeometric test) (Extended Data Fig. 8c–h, Supplementary Table 7). A large proportion of SA-responsive genes were more strongly expressed in phr1 and phr1 phl1 than in Col-0; these were strongly enriched for classical SA-dependent defense genes (Extended Data Fig. 8d, e). A second group of SA-responsive genes that were less expressed in phr1 and phr1 phl1 than in Col-0 lacked classical SA-dependent defense genes and were weakly enriched for genes likely contributing to PSR (Extended Data Fig. 8d). By contrast, most JA-responsive genes exhibited lower expression in phr1 and phr1 phl1 (Extended Data Fig. 8g, h), including a subset of 18 of 46 genes known or predicted to mediate biosynthesis of defense-related glucosinolates (Extended Data Fig. 8i)21. This agrees with the recent observation that phr1 exhibited decreased glucosinolate levels during Pi starvation22. Analyses of SA- and JA- up-regulated genes revealed enrichment of direct PHR1 targets (Fig. 3e), consistent with the converse observation that some PHR1-regulated clusters enriched in direct targets were also enriched in defense genes (Fig. 3c, d). Many of the SA- and JA- responsive genes were PSR-SynCom DEGs (Fig. 3f; Extended Data Fig. 8c–h, Supplementary Table 7). Thus, PHR1 directly regulates an unexpected proportion of the plant immune system during PSR triggered by our SynCom.
We tested whether PHR1 also controls the expression of plant defense genes under conditions typically used to study PSR (axenic growth, sucrose present, no microbiota involved). We performed RNA-seq in response to low Pi in these conditions and identified 1482 DEGs in Col-0 and 1161 DEGs in phr1 phl1 (Fig. 4a, b; Extended Data Fig. 9, Supplementary Table 11). A significant number of our BTH/SA-activated genes were also up-regulated in phr1 phl1, but not in Col-0 in response to low Pi (Fig. 4a, b; Supplementary Table 12). A large number of these overlapped with the defense genes induced in phr1 phl1 by our SymCom (Fig. 4c; red ellipse, 113/337 = 33%; clusters c3 and c8 from Fig. 3a). At least 14/113 are direct PHR1 targets (Supplementary Table 12).
To underscore the role of PHR1 in the regulation of response to microbes, we analyzed the transcriptional profile of Col-0 and phr1 phl1 plants exposed to the flagellin peptide flg22. We chose a chronic exposure to flg22 to mimic the condition of plants in contact with a microbiome. We found that phr1 phl1 plants displayed higher expression of flg22-responsive genes23 than Col-0, independent of phosphate status (Supplementary Text 5; Fig. 4d; Extended Data Fig. 9a, b; Supplementary Table 11; Supplementary Table 13). This indicates that PHR1 negatively regulates the immune response triggered by flg22.
We hypothesized that phr1 phl1 would express an altered response to pathogen infection. The phr1 and phr1 phl1 mutants exhibited enhanced disease resistance against both Hyaloperonospora arabidopsidis isolate Noco2, and Pseudomonas syringae DC3000 (Fig. 4e, f). Collectively, these results confirm the role of PHR1 as a direct integrator of PSR and the plant immune system.
Plant responses to phosphate stress are inextricably linked to life in microbe-rich soil. We demonstrate that genes controlling PSR contribute to assembly of a normal root microbiome. Surprisingly, our SynCom enhanced the activity of PHR1, the master regulator of the PSR, in plants grown under limited phosphate. This led to our discovery that PHR1 is a direct regulator of a functionally relevant set of plant immune system genes. Despite being required for the activation of JA-responsive genes during PSR24, we found that PHR1 is unlikely to be a general regulator of this response (Extended Data Fig. 9c–e, Supplementary Table 12), but rather may fine-tune JA response in specific biological contexts.
We demonstrate that PSR and immune system outputs are directly integrated by PHR1 (and, likely, PHL1). We provide a mechanistic explanation for previous disparate observations that PSR and defense regulation are coordinated and implications that PHR1 is the key regulator8, 11, 12, 24. We provide new insight into the intersection of plant nutritional stress response, immune system function, and microbiome assembly and maintenance; systems that must act simultaneously and coordinately in natural and agricultural settings. Our findings will drive investigations aimed at enhancing phosphate use efficiency using microbes.
Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper.
For experiments in wild soil, we collected the top-soil (approx. 20 cm) from a site free of pesticide and fertilizer at Mason Farm (MF; North Carolina, USA; +35°53′30.40′′, −79°1′5.37′′) 19. Soil was dried, crushed and sifted to remove debris. To improve drainage, soil was mixed 2:1 volume with autoclaved sand. Square pots (2 × 2 inch square) were filled with the soil mixture and used to grow plants. Soil micronutrient analysis is published by Lundberg, D.S. et al. 19.
All Arabidopsis thaliana mutants used in this study were in the Columbia (Col-0) background (Supplementary Table 16). All seeds were surface-sterilized with 70 % bleach, 0.2 % Tween-20 for 8 minutes, and rinsed 3X with sterile distilled water to eliminate any seed-borne microbes on the seed surface. Seeds were stratified at 4 °C in the dark for 2 days.
To determine the role of phosphate starvation response in controlling microbiome composition, we analyzed five mutants related to the Pi-transport system (pht1;1, pht1;1 pht1;4, phf1, nla, and pho2) and two mutants directly involved in the transcriptional regulation of the Pi-starvation response (phr1 and spx1 spx2). All these genes are expressed in roots 13– 18.
Seeds were germinated in sterile square pots filled with MF soil prepared as described above. We also used pots without plants as “bulk soil” controls. All pots, including controls, were watered from the top with non-sterile distilled water to avoid chlorine and other tap water additives 2 times a week. Plants were grown in growth chambers with a 16-h dark/8-h light regime at 21 °C day 18 °C night for 7 weeks. In all experiments, pots with plants of different genotypes were randomly placed in trays according to true random numbers derived from atmospheric noise; we obtained those numbers from www.random.org. We positioned trays in the growth chamber without paying attention to the pots they contained, and we periodically reshuffled them without paying attention to the pot labels.
Plants and bulk soil controls were harvested and their endophytic compartment (EC) microbial communities isolated as described19. DNA extraction was performed using 96-well format MoBio PowerSoil Kit (MOBIO Laboratories) following the manufacturer’s instruction.
The method of Ames25 was used to determine the phosphate concentration in the shoots of seedlings grown on different Pi regimens and treatments. Main root length elongation was measured using ImageJ software26 and for shoot area and number of lateral roots WinRhizo software27 was used.
For wild soil experiment 16S sequencing, we processed libraries according to Caporaso, et al.28. Three sets of index primers were used to amplify the V4 (515F-806R) region of the 16S rRNA gene of each sample. In each case, the reverse primer had a unique molecular barcode for each sample28. PCR reactions with ~20 ng template were performed with 5 Prime Hot Master Mix in triplicate using plates 2, 4 and 5 from the 16S rRNA Amplification Protocol28. PCR blockers mPNA and pPNA29 were used to reduce contamination by plant host plastid and mitochondrial 16S amplicon. The PCR program used was:
Reactions were purified using AMPure XP magnetic beads and quantified with Quant IT Picogreen. Amplicons were pooled in equal amounts and then diluted to 5.5 pM for sequencing. Samples were sequenced on an Illumina MiSeq machine at UNC, using a 500-cycle V2 chemistry kit. The library was spiked with 25 % PhiX control to increase sequence diversity. The raw data for the wild soil experiments is available in the EBI Sequence Read Archive (accession PRJEB15671).
For SynCom experiment 16S library, we 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, D.S. et al.29. Basically, the molecule-tagging step was changed to an exponential amplification to account for low DNA yields with the following reaction:
|5 µL||Kapa Enhancer|
|5 µL||Kapa Buffer A|
|1.25 µL||5 µM 338F|
|1.25 µL||5 µM 806R|
|0.375 µL||mixed PNAs (1:1 mix of 100 µM pPNA and 100 µM mPNA)|
|0.5 µL||Kapa dNTPs|
|0.2 µL||Kapa Robust Taq|
Following 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, D.S. et al.29. Sequencing was performed at UNC on an Illumina MiSeq instrument using a 600-cycle V3 chemistry kit. The raw data for the SynCom experiments is available in the EBI Sequence Read Archive accession PRJEB15671.
For wild soil census analysis, sequences from each experiment were pre-processed following standard method pipelines from2, 19. Briefly, sequence pairs were merged, quality-filtered and de-multiplexed according to their barcodes. The resulting sequences were then clustered into Operational Taxonomic Unit (OTUs) using UPARSE30 implemented with USEARCH7.1090, at 97 % percent identity. Representative OTU sequences (Supplementary Dataset 1) were taxonomically annotated with the RDP classifier31 trained on the Greengenes database (4/February/2011; Supplemental Dataset 1). We used a custom script (https://github.com/surh/pbi/blob/master/census/1.filter_contaminants.r) to remove organellar OTUs, and OTUs that had no more than a kingdom-level classification, and an OTU count table was generated (Supplementary Table 1, Supplementary Dataset 1).
SynCom sequencing data were processed with MT-Toolbox32. 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 Q-score under 20, and trimmed from the 5’ end to a final length of 270 bp. The resulting sequences were matched to a reference set of the strains in the SynCom generated from Sanger sequences, the sequence from a contaminant strain (47Yellow) that grew in the plate from strain 47 (Supplementary Table 2) and Arabidopsis organellar sequences. Sequence mapping was done with USEARCH7.1090 with the option –”usearch_global” at a 98 % identity threshold. 90 % of sequences matched an expected isolate, and those sequence mapping results were used to produce an isolate abundance table. The remaining unmapped sequences were clustered into OTUs with the same settings used for the census experiment, the vast majority of those OTUs belonged to the same families as isolates in the SynCom, and were probably unmapped due to PCR and/or sequencing errors. We combined the isolate and OTU count tables into a single master table. The resulting table was processed and analyzed with the code at (https://github.com/surh/pbi/blob/master/syncom/7.syncomP_16S.r). Matches to Arabidopsis organelles were discarded. PCR blanks were included in the sequencing and the average counts per strain observed on those blanks were subtracted from the rest of the samples following33. Extended Data Fig. 7a shows the number of usable reads across samples, and the remaining number after subtracting sterile controls (blanks).
For physiological, transcriptional analysis or pathology experiments, we used phr1, phr1 phl1, phf1, and coi1–16, sid2-1 mutants, which are all in the Col-0 genetic background (Supplementary Table 16). For all physiological and transcriptional analysis in vitro, Arabidopsis seedlings were grown 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] solidified with 1 % bacto-agar (BD, Difco). Media were supplemented with Pi (KH2PO4) at distinct concentrations depending on the experiment; 1 mM Pi was used for complete medium and approximately 5 µM Pi (traces of Pi in the agar) was the Pi concentration in the medium not supplemented with Pi. Unless otherwise stated, plants were grown in a growth chamber in a 15-h dark/9-h light regime (21 °C day /18 °C night).
For Synthetic Community experiments, plants were germinated on Johnson medium containing 0.5 % sucrose, with 1 mM Pi, 5 µM Pi or supplemented with KH2PO3 (phosphite) at 1 mM for 7 d in a vertical position, then transferred to 50 µM Pi or 625 µM Pi media (without sucrose) alone or with the Synthetic Community at 105 c.f.u/mL, for another 12 d. For the heat-killed SynCom experiments, plants were grown as above. Heat-killed SynComs were obtained by heating different concentrations of bacteria: 105 c.f.u/mL, 106 c.f.u/mL and 107 c.f.u/mL at 95 °C for 2 h in an oven. The whole content of the heat-killed SynCom solutions were added to the media.
For the functional activation of the PSR by the SynCom, plants were germinated on Johnson medium containing 0.5 % sucrose, 1 mM Pi for 7 d in a vertical position, then transferred to 0, 10, 30, 50 and 625 µM Pi alone, or to 0, 50 and 625 µM Pi with the Synthetic Community at 105 c.f.u / mL, for another 12 d. At this point, we harvested our time zero (3 replica per conditions, each replica was 5 shoots harvested across all plates used). The remaining plants were transferred again to 1mM Pi to evaluate the capacity of the plants for Pi accumulation in a time series analysis. We harvested plant shoots every 24 h for 3 days and Pi-concentration was determined. Pi increase was calculated as: (Pi concentration day n – Pi concentration day 0) / Pi-concentration day 0. Relative increase in Pi concentration is plotted in Fig. 2c. Both relative and absolute Pi concentration values are provided in Supplementary Table 4.
We repeated this experiment twice. For the first experiment, we used 6 plates with 10 plants per condition (48 plates and 480 plants in total). We harvested three replicas per time point with 5 shoots each. In all cases, shoots were harvested across all plates used. For the second experiment, we used 11 plates with 10 plants per condition (88 plates and 880 plants). In this case, we harvested 6 replicas for 1, 2 and 3 days after the re-feeding with Pi, and three replicas for time zero. Each replica contains 5 shoots harvested across all the plates used.
For the demonstration that sucrose is required for the induction of PSR in sterile conditions, plants overexpressing the PSR reporter construct IPS1:GUS13 were grown in Johnson medium containing 1 mM Pi or 5 µM Pi supplemented with different concentrations of sucrose. After 12 days, the expression of the reporter constructs IPS1:GUS, highly induced by low Pi, was followed by GUS staining. Plants were grown in a growth chamber in a 15-h light/9-h dark regime (21 °C day /18 °C night).
For the ChIP-seq experiment, phr1 harboring the PromPHR1:PHR1-MYC construct18 and Col-0 seedlings were grown on Johnson medium 1 mM Pi, 1 % sucrose for 7 days and then transferred to a media not supplemented with Pi for another 5 days. Plants were grown in a growth chamber in a 15-h light/9-h dark regime (21 °C day /18 °C night). A total of 2364 genes were identified as regulated by PHR1. The ChIP-seq data will be fully presented in de Lorenzo, et al.34.
For the transcriptional analysis under conditions typically used to study PSR (axenic growth with sucrose present; no microbiota involved), with Methyl Jasmonate (MeJA) and the 22-amino acid flagellin peptide (flg22), plants were germinated on Johnson medium (1 % sucrose) containing 1 mM Pi for 7 d in a vertical position and then transferred to 1 mM Pi and 5 µM Pi media containing 1 % sucrose either alone or supplemented with 10 µM MeJA (Sigma) or 1 µM flg22 (Sigma) for 12 d.
For growth inhibition assays, seedlings were grown on Johnson medium (1 % Sucrose) in 1 mM Pi and 5 µM Pi conditions for 5 d, transferred to 1 mM Pi and 5 µM Pi media supplemented or not with 10 µM MeJA for 5 d. Main root length was then measured using ImageJ software26.
For Synthetic Community experiments, we selected 35 diverse bacterial strains. 32 of them were isolated from roots of Arabidopsis and other Brassicaceae species grown in two wild soils19. Two strains came from Mason Farm unplanted soil19, and Escherichia coli DH5α was included as a control (Supplementary Table 2). More than half (19/35) of the strains belonged to families enriched in the EC of plants grown in Mason Farm soil (Supplementary Table 2)2, 19. The strains were chosen from a larger isolate collection in a way that maximizes SynCom diversity while retaining enough differences in their 16S rRNA gene to allow for easy and unambiguous identification.
A single colony of bacteria to be tested was inoculated in 4 mL of 2XYT medium (16 g/L Tryptone, 10 g/L Yeast Extract, 5 g/L NaCL, ~5.5 mM Pi) in a test tube. Bacterial cultures were grown while shaking at 28°C overnight. At this point, the Pi concentration was reduced to by dilution to 5 mM Pi average in the supernatants (10 cultures used for the quantification). Cultures were then rinsed with a sterile solution of 10 mM MgCl2 followed by a centrifugation step at 2600 g for 8 min. This process was repeated twice. The concentration of Pi in the supernatant after the first wash with MgCL2 was 0.06 mM Pi and after the second wash it was reduced to 0.005 mM Pi. In the suspension of SynCom member cells in MgCl2, the average concentration of Pi was 0.08 mM. The OD600nm was measured and assuming that 1 OD600nm unit is equal to 109 c.f.u/mL we equalized individual bacterium concentration to a final value of 105 c.f.u/mL of medium. The concentration of Pi in the final SynCom was 0.09 µM Pi. Thus, based on these results, we were not Pi fertilizing the plant by adding the SynCom. Medium was cooled down (to 40–44°C) near the solidification point and then the bacteria mix was added to the medium with agitation. We monitored the pH in the media after adding 1, 5, 10 mL of 10 mM MgCl2 which represents almost ten times the volume we used to add the SynCom. After adding MgCl2 the pH in the media remained stable. We also analyzed the pH after adding the SynCom at 105, 106 and 107 c.f.u/ml of media and found no pH changes. Therefore, we considered that the MES buffer we used was appropriate for this experiment.
To isolate and quantify bacteria from plant roots in the SynCom experiment, plant roots were harvested, and rinsed 3 times with sterile distilled water to remove agar particles and weakly associated microbes. Plant material was then freeze-dried. Root pulverization and DNA extraction was conducted as described above.
To isolate and quantify bacteria from agar samples, a freeze and squeeze protocol was used. Syringes with a square of sterilized miracloth at the bottom were completely packed with agar 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 min and most of the supernatant discarded. The remaining 1–2 mL of supernatant, containing the pellet, was moved into clean microfuge tubes. Samples were centrifuged again, supernatant was removed, and pellets were used for DNA extraction. DNA extraction was performed using 96-well format MoBio PowerSoil Kit (MOBIO Laboratories).
For oomycete pathology studies, Hyaloperonospora arabidopsidis (Hpa) isolate Noco2 was propagated on the susceptible Arabidopsis ecotype Col-0. Spores of Hpa were suspended in deionized sterile water at a concentration of 5×104 spores/mL. The solution containing spores was spray-inoculated onto 10-d-old seedlings of Arabidopsis grown in fertilized potting soil. Inoculated plants were grown at 21 °C under a 9-h light regime. Asexual sporangiophores were counted 5 d post-inoculation on at least 100 cotyledons for each genotype.
For bacterial pathology studies, Pseudomonas syringae pv. tomato DC3000 was suspended in 10 mM MgCl2 to a final concentration of 105 c.f.u/mL. 35–40-d-old plants of Arabidopsis grown on soil were hand-infiltrated using a needle-less syringe on the abaxial leaf surface. Leaf discs (10 mm diameter) were collected after 1 h and 3 d post inoculation, and bacterial growth was measured as described35.
We performed 3 different sets of RNA-seq experiments in this study. (I) The first set (Fig. 2b, Fig. 3 and Extended Data Fig. 4b) evaluated the effect of the SynCom on the phosphate starvation response of Arabidopsis seedlings. In addition to wild-type Col-0 (4 replicates), phf1 (4 replicates) and phr1 phl1 (4 replicates) were included in the experiment shown in Fig. 2b, whereas Col-0 (10 replicates), phr1 (10 replicates) and phr1 phl1 (6 replicates) were used in the experiment shown in Fig. 3 and Extended Data Fig. 4b. (II) The second experiment (Extended Data Fig. 6a, b) is an expansion of the first and was designed to evaluate whether different pre-treatments (1 mM Pi, 5 µM Pi, 1 mM Phosphite [Phi]) influence the phosphate starvation response triggered by the SynCom. We used Col-0 (4 replicates), phf1 (4 replicates) and phr1 phl1 (4 replicates) in this experiment. (III) Finally, the third experiment evaluated the effect of MeJA and flg22 on the phosphate starvation response (Fig. 4 and Extended Data Fig. 9) of Arabidopsis seedlings. The genotypes Col-0 (6 replicates) and phr1 phl1 (6 replicates) were used. The experiments listed above were repeated between two and five independent times and each repetition (defined as “batch” in the generalized linear model, see RNA-seq data analysis, below) included two biological replicates per genotype per condition. Supplementary Table 15 contains the metadata information of all RNA-seq experiments. Raw reads and read counts are available at the NCBI Gene Expression Omnibus under accession number GSE87339.
Total RNA was extracted from roots of Arabidopsis according to Logemann, et al.36. 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 (20000 g, 10 min) for phase separation. The aqueous phase was transferred to a new 1.5 ml tube and 0.05 volumes of 1N acetic acid and 0.7 volumes 96 % ethanol were added. The RNA was precipitated at −20 °C overnight. Following centrifugation, (20000 g, 10 min, 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.
Illumina-based mRNA-seq libraries were prepared from 1000 ng 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 min. 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.
Initial quality assessment of the Illumina RNA-seq reads was performed using the FASTX-Toolkit. Cutadapt37 was used to identify and discard reads containing the Illumina adapter sequence. The resulting high-quality reads were then mapped against the TAIR10 Arabidopsis reference genome using Tophat38, with parameters set to allow only one mismatch and discard any read that mapped to multiple positions in the reference. The Python package HTSeq39 was used to count reads that mapped to each one of the 27,206 nuclear protein-coding genes. Extended Data Fig. 10 shows a summary of the uniquely mapped read counts per library. Raw sequencing data and read counts are available at the NCBI Gene Expression Omnibus accession number GSE87339.
Differential gene expression analyses were performed using the generalized linear model (glm) approach40 implemented in the edgeR package41. This software was specifically developed and optimized to deal with over-dispersed count data, which is produced by RNA-seq. Normalization was performed using the trimmed mean of M-values method (TMM42; function calcNormFactors in edgeR). The glmFit function was used to fit the counts in a negative binomial generalized linear model with a log link function40. For the SynCom experiment (Fig. 3), the model includes the covariates: phosphate content (High or Low), bacteria (present or absent) and batch effect. A term for the interaction between Phosphate and Bacteria was included as represented below:
The model used to analyze the effect of MeJA and flg22 (Fig. 4) included the following covariates: phosphate content (High or Low), MeJA (present or absent), flg22 (present or absent) and batch effect.
In each model, the term “Batch” refers to independent repetitions of the experiment (see the Genome-wide gene expression analyses section). Data from the different genotypes were fitted independently with the same model variables. The Benjamini-Hochberg method (False Discovery Rate; FDR)43 was applied to correct the p-values after performing multiple comparisons. Genes with FDR below or equal to 0.01 and fold-change variation of at least 1.5X were considered differentially expressed.
Transcriptional activation of the phosphate starvation response was studied using a literature-curated set of phosphate starvation marker genes (Extended Data Fig. 4a, Supplementary Table 3). This core set consists of 193 genes that were up-regulated by phosphate starvation stress across four different gene expression experiments13, 44– 46. The RPKM (Reads Per Kilobase of transcript per Million mapped reads) expression values of these 193 genes were z-score transformed and used to generate box and whiskers plots to show the distribution of the expression values of this gene set.
Hierarchical clustering analyses were performed with the heatmap.2 function in R from the gplots package47 using the sets of differentially expressed genes identified in each experiment. Genes were clustered based on the Euclidean distance and with the complete-linkage method. Genes belonging to each cluster were submitted to Gene Ontology (GO) enrichment analyses on the PlantGSEA platform48 in order to identify over-represented biological processes.
Genes whose transcription is induced by MeJA (672 genes), BTH/SA (2096 genes) or both hormones (261 genes) were used as markers of the activation of these immune response output sectors in Arabidopsis (Supplementary Table 10)49. These gene sets were defined using two-week old Col-0 seedlings grown on potting soil and sprayed with MeJA (50 µM; Sigma), BTH (300 µM; Actigard 50WG) or a mock solution (0.02 % Silwet, 0.1 % ethanol). Samples were harvested 1 h, 5 h and 8 h after the treatment in two independent experiments. Total RNA was extracted with the RNeasy Plant Mini kit (Qiagen) and then used to prepare Illumina mRNA-seq libraries. The bioinformatics pipeline to generate count tables and the criteria used to define differentially expressed genes between conditions (Hormone treatment vs. Mock treatment) was the same as described above. Raw sequencing data are available at the NCBI Gene Expression Omnibus under the accession number GSE90077.
Most statistical analyses were performed in the R statistical environment50 and follow methods previously described2. As described in the following subsections, a number of packages were used, and many were called through AMOR-0.0–1451, which is based on code from Lebeis, et al.2. All scripts and knitr52 output from R scripts are available upon request. Most plots are ggplot253 objects generated with functions in AMOR51. For all linear modeling analyses (ANOVA, ZINB, GLM), terms for batch and biological replicate were included whenever appropriate. Code for both census and SynCom analysis is available at https://github.com/surh/pbi.
For wild soil and SynCom experiments, the number of samples per genotype and treatment was determined based on our previously published work, which showed that seven and five samples are enough to detect differences in wild soils and SynCom experiments, respectively.2, 19. For RNA-seq experiments, we used at least four replicates per condition, which is sufficient for parameter estimation with the edgeR software41.
Alpha and beta diversity were calculated on count tables that were rarefied to 1000 reads. Samples with less than this number of usable reads (i.e. high quality non-organellar reads) were discarded. Alpha diversity (Shannon index, richness) metrics were calculated using the “diversity” function in vegan54, and differences between groups were tested with ANOVA (Extended Data Fig. 1a). Site diversity (Extended Data Fig. 1b) was calculated with the “sitediv” function in AMOR51. Unconstrained ordination was performed with vegan (Bray-Curtis), and Principal Coordinate Analysis (PCoA) was performed with AMOR (Extended Data Fig. 1d)51. Canonical Analysis of Principal Coordinates (CAP) is a form of constrained ordination55 and was performed using the “capscale” function of the vegan package in R54. CAP was performed on the full counts of the EC samples only, using the “Cao” distance. Constraining was done separately on plant genotype while conditioning on sequencing depth and biological replicate. This approach allowed us to focus on the portion of variation that is associated with plant genotype, conditionally, independent of other factors.
For the SynCom experiments, richness was directly calculated in R. Principal Coordinate Analysis was performed with the “PCO” function of AMOR51 using the “Cao” distance which was calculated with vegan54 on an abundance table rarefied to 1500 reads per sample. Canonical Analysis of Principal Coordinates (CAP) was performed using the “capscale” function of the vegan package54 in R. CAP was performed on the full counts of the root samples only, using the “Cao” distance. Constraining was done separately on Fraction, Pi level and plant genotype while conditioning on sequencing depth and the other covariates.
Differentially abundant bacterial taxa across fraction and genotype in the wild soil experiments were identified using the same approach as in Lebeis, et al.2. Briefly, we used a Zero-Inflated Negative Binomial (ZINB) framework that allowed us to test for the effect of specific variables, while both controlling for the other covariates and accounting for the excess of zero entries in the abundance tables. These zero-entries likely represented under-sampling and not true absences. The same analysis was performed at the family and OTU-level on the measurable OTUs (taxa that have an abundance of at least 25 counts in at least five samples)19. Results are in Extended Data Fig. 1e–h, Supplementary Table 1. Extended Data Fig. 1h shows the distribution of significant genotypic effects on bacterial abundances at both taxonomic levels; in both cases the behavior is similar, indicating small and even effects of all genotypes.
For the comparison of enrichment profiles between genotypes, we followed the same Monte-Carlo approach described in Lebeis, et al.2. Briefly we looked at the enrichment/depletion profile of bacterial taxa for each mutant compared to wild-type Col-0, and asked, for each pair of mutants, if they were more similar than expected by chance and assed significance by random permutation. Results are in Fig. 1d, Extended Data Fig. 1g.
To define differentially abundant strains in SynCom experiments, we found that a Negative Binomial GLM approach gave more stable results than the ZINB approach. We used the edgeR package41 to fit a quasi Negative Binomial GLM model with the glmQLFit function, and significance was tested with the glmQLFtest function56. Results of all relevant pairwise comparisons are in Extended Data Fig. 7 and Supplementary Table 5.
For the definition of robust colonizers in synthetic community experiments, we calculated the average relative abundance of E. coli on all root samples and counted, for each strain, how many times it was more abundant than E. coli’s average on the same set of root samples. Then we used a one-sided binomial test to ask if the probability of a given strain to be more abundant than the average E.coli was significantly higher than a coin toss (50%). Strains that passed the test were labeled as robust-colonizers, the rest of the strains were labeled as Sporadic or Non-Colonizers. The results are indicated in Fig. 2e and Supplementary Table 2.
All data generated from this project is publicly available. Raw sequences from soil census and SynCom colonization are available at the EBI Sequence Read Archive under accession PRJEB15671. Count tables, metadata, taxonomic annotations and OTU representative sequences from the Mason Farm census and Syncom experiments are available as Supplementary Datasets 1 and Supplementary Datasets 2 respectively. Custom scripts used for statistical analysis and plotting are available at (https://github.com/surh/pbi). Raw sequences from transcriptomic experiments are available at the NCBI Gene Expression Omnibus under the accession number GSE87339. The corresponding metadata information is provided in Supplementary Table 15. All code is available upon request.
Support by NSF INSPIRE grant IOS--1343020 and DOE-USDA Feedstock Award DE-SC001043 to J.L.D. S.H.P was supported by NIH Training Grant T32 GM067553--06 and is a Howard Hughes Medical Institute International Student Research Fellow. P.J.P.L.T was supported by The Pew Latin American Fellows Program in the Biomedical Sciences. J.L.D is an Investigator of the Howard Hughes Medical Institute, supported by the HHMI and the Gordon and Betty Moore Foundation (GBMF3030). M.E.F and O.M.F are supported by NIH NRSA Fellowships F32-GM112345-02 and F32-GM117758-01, respectively. N.W.B was supported by NIH NRSA Fellowship F32--GM103156. J.P-A is funded by the Spanish Ministry of Economy and Competitiveness (MINECO BIO2014-60453-R and EUI2008-03748). We thank Stratton Barth and Emily Getzen for technical assistance, the Dangl lab microbiome group for useful discussions and Drs. Sarah Grant, Derek Lundberg, Farid El Kasmi, Paul Schulze-Lefert and his colleagues for critical comments on the manuscript. Supplement contains additional data. Raw sequence data is available at the EBI Sequence Read Archive accession PRJEB15671 for microbiome 16S profiling, and at the Gene Expression Omnibus accessions GSE87339 for transcriptomic experiments. J.L.D is a co-founder of, and shareholder in, and S.H.P collaborates with, AgBiome LLC, a corporation whose goal is to use plant-associated microbes to improve plant productivity.
Supplementary Information is linked to the online version of the paper at www.nature.com/nature
Author contributionsG.C., P.J.P.L.T., S.H.P and J.L.D. designed the project, G.C., S.H.P, T.F.L and M.E.F. set up the experiments, harvested samples and organized construction of 16S sequencing libraries. G.C and T.F.L performed control experiments related with PSR induced by the SynCom. G.C, N.W.B, M.E.F and T.F.L set up the experiments, harvested samples and isolated RNA. P.J.P.L.T. organized, performed construction of RNA-seq libraries and analysed RNA-seq data. S.H.P. analysed 16S sequencing data. S.H.P and P.J.P.L.T oversaw data deposition. G.C, T.F.L. and P.J.P.L.T performed pathology experiments. G.C., P.J.P.L.T, S.H.P, T.F.L, O.M.F and J.L.D. analysed data and created figures. L.d.L. performed the ChIP-seq experiment. C.D.J and P.M. advised on sequencing process and statistical methods. G.C., P.J.P.L.T., S.H.P., and J.L.D. wrote the manuscript with input from O.M.F., C.D.J., and J.P.-A.
All data generated from this project is publicly available. Raw sequences from soil census and SynCom colonization are available at the EBI Sequence Read Archive accession PRJEB15671. Count tables, metadata, taxonomic annotations and OTU representative sequences from the Mason Farm census and Syncom experiments are available as Supplementary Datasets 1 and Supplementary Datasets 2 respectively. Custom scripts for analysis and visualization are at https://github.com/surh/pbi. Raw sequences and counts from RNA-seq experiments are available at the NCBI Gene Expression Omnibus under accession number GSE87339. The corresponding metadata information is provided in Supplementary Table 15. All code is available upon request. Reprint and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interest. Readers are welcome to comment on the online version of this article at www.nature.com/nature.