The plant NADPH oxidase RBOHD is required for microbiota homeostasis in leaves.

The plant microbiota consists of a multitude of microorganisms that can affect plant health and fitness. However, it is currently unclear how the plant shapes its leaf microbiota and what role the plant immune system plays in this process. Here, we evaluated Arabidopsis thaliana mutants with defects in different parts of the immune system for an altered bacterial community assembly using a gnotobiotic system. While higher-order mutants in receptors that recognize microbial features and in defence hormone signalling showed substantial microbial community alterations, the absence of the plant NADPH oxidase RBOHD caused the most pronounced change in the composition of the leaf microbiota. The rbohD knockout resulted in an enrichment of specific bacteria. Among these, we identified Xanthomonas strains as opportunistic pathogens that colonized wild-type plants asymptomatically but caused disease in rbohD knockout plants. Strain dropout experiments revealed that the lack of RBOHD unlocks the pathogenicity of individual microbiota members driving dysbiosis in rbohD knockout plants. For full protection, healthy plants require both a functional immune system and a microbial community. Our results show that the NADPH oxidase RBOHD is essential for microbiota homeostasis and emphasizes the importance of the plant immune system in controlling the leaf microbiota.

In nature, plants are colonized by a diverse microbiota. The phyllosphere constitutes a vast microbial habitat in which microorganisms populate plant surfaces as epiphytes and the intercellular space inside leaves as endophytes[1,2]. The functional repertoire of the microbiome expands the capacity of the plant to adapt to its environment, as beneficial microbes can promote plant growth and increase tolerance to abiotic and biotic stress[3-5]. Plants must be prepared to recognize pathogens and mount defense responses against them[6,7] while being colonized by a complex microbiota that can potentially also trigger immunity and provide indirect plant protection[8,9]. Therefore, plant immunity must be balanced to allow accommodation of a microbiota[4,10]. Emerging evidence suggests that the plant immune system has a significant impact on the bacterial microbiota in the rhizosphere[9,11,12] and phyllosphere[13-15]. However, our current understanding of specific host factors involved in the perception and structuring of the microbiota and the molecular mechanisms underlying them is limited.

In recent years, experiments with synthetic communities (SynCom) have been introduced to identify host genetic factors that drive microbiota assembly under controlled environmental conditions[16]. The availability of representative collections of strains isolated from Arabidopsis thaliana and the reproducible reconstitution of the microbiota to similar phylogenetic structures as found in nature[11,17,18] provide opportunities to quantify community changes and probe available genetic resources to establish gene-phenotype causal relationships.

Here, we used a reductionist approach to link host immune components to the assembly of the bacterial microbiota. To test the impact of the plant immune system on the microbiota, we used single and higher order mutants of well-described key immunity factors, and probed changes in the microbiota composition compared to the one that establishes on wild type plants. To identify specific plant-microbiota interactions we adopted a SynCom approach to achieve strain-level resolution. We found a number of A. thaliana immunity mutants with altered bacterial community composition and demonstrate a critical function of an immunity-related NADPH oxidase in shaping the leaf microbiota. Ultimately, we explore the interplay between the host and individual community members and their roles in plant health and community composition.

Results

Bacterial community in the A. thaliana phyllo- and endosphere

For microbiota reconstitution experiments, we used an established gnotobiotic growth system consisting of A. thaliana and a bacterial synthetic community based on the At-LSPHERE collection of genome-sequenced strains[17,19]. The microbiota inoculum communities consisted of either 222 or 223 strains (termed SynCom-222 or SynCom-223), which represent the maximal genomic diversity of the At-LSPHERE and included strains with identical amplicon sequence variant (ASV) but different genomes. In follow-up experiments, we chose a subset of 137 strains (SynCom-137) in which each strain represent an ASV with a unique 16S rRNA gene (rDNA) amplicon sequence (Supplementary Table; Methods). Germ-free 10-day-old A. thaliana seedlings were inoculated with the SynCom and the aboveground part of the plants were harvested 5.5 weeks after germination to examine the bacterial communities by 16S rDNA amplicon sequencing (Fig. 1a).

Figure 1

Characterization of the At-LSPHERE model community in the phyllosphere and endosphere of Arabidopsis thaliana.

a) Schematic workflow of SynCom experiments. Two types of SynCom (SynCom-222 and SynCom-137) with similar taxonomic compositions but different numbers of strains as indicated below the stacked bar plot were used (see Supplementary Table). Colors represent bacterial phyla and classes for Proteobacteria. In the experiments, germ-free Arabidopsis thaliana seedlings were inoculated after 10 days, and aboveground plant tissue was harvested 5.5 weeks after germination. Plant phenotype was characterized by scoring disease symptoms and fresh weight measurements. Bacterial community profiles were determined for the entire phyllosphere community and the endosphere-enriched community by 16S rDNA amplicon sequencing. b) Community composition in the A. thaliana phyllosphere (“phyllo”) and endosphere (“endo”) after inoculation with SynCom-222 on phylum and class level for Proteobacteria. Numbers above bar stacks indicate a 1000-fold difference between bacterial loads in both compartments. c) Overall bacterial load on SynCom-inoculated plants measured by colony forming units per gram plant fresh weight (n=5) and by qPCR of bacterial 16S rDNA relative to plant gene abundance in the phyllosphere and endosphere (n=3). d) Volcano plot shows log2-fold changes in strain abundances of SynCom-222 between phyllosphere and endosphere samples (n=12). Dot size represents relative abundance in endosphere. Statistical significance of differential strain abundance is expressed with p-values determined by two-sided Wald test and Benjamini-Hochberg adjusted. Black labels highlight strains that were affected in all three SynCom experiments (see Source Data Figure 1). e) Relative abundance of strains (or ASVs indicated by superscript circles) determined by 16S rDNA amplicon sequencing of the SynCom-222 endosphere. Detected strains with relative abundance >0.01% are displayed (see entire graph in Extended Data Figure 1). Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. “nd” and dot size indicate the number of plant replicates where a given strain was not detected (n=12). Colors represent strain phylogeny.

The relative abundance of phyllosphere strains (Fig. 1b) was similar to the natural microbiota of A. thaliana at the phylum level and was consistent with SynCom experiments from previous studies[17,19]. The community structure in the A. thaliana Col-0 phyllosphere comprised Alphaproteobacteria (36%), Betaproteobacteria (30%), Gammaproteobacteria (2%), Actinobacteria (19%), Bacteroidetes (12%) and Firmicutes (0.0005%) (Fig. 1b). In addition to these phyllosphere samples, we analyzed the endosphere-enriched bacterial fraction (Fig. 1a; Methods). Comparing the community composition, we observed a relative enrichment of Gammaproteobacteria in the endosphere compared to the entire phyllosphere (Fig. 1b and Extended Data Fig. 1). The quantification of bacterial abundance by qPCR using 16S rDNA copies per plant gene revealed reduced colonization of the endosphere. Similarly, colony forming units (CFU) per gram of fresh weight recovered from endosphere samples were approximately 1000-fold lower than those recovered from phyllosphere samples (Fig. 1c), as expected[14,20].

Extended Data Fig. 1

SynCom-222 composition of inoculum, phyllosphere and endosphere.

Relative abundance of strains (or ASVs indicated by superscript circle) determined by 16S rDNA amplicon sequencing of samples from SynCom-222 inoculum mix, phyllosphere and endosphere samples of A. thaliana Col-0. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. "Undetected" and dot size indicate the number of plant replicates where a given strain was not detected (n=12). Colors represent strain phylogeny.

The analysis of the relative abundance of strains revealed that a few bacteria were significantly enriched or depleted in the endosphere across three SynCom experiments (Fig. 1d). Xanthomonas Leaf131 in particular but also Pseudomonas Leaf59 were enriched in the endosphere samples, thus explaining the higher relative abundance of Gammaproteobacteria (Fig. 1b). In contrast, Methylophilus Leaf416, Aeromicrobium Leaf350 and Pseudorhodoferax Leaf265 were depleted, indicating that these strains or ASVs have a greater relative abundance on the leaf surface compared to the community in the apoplast (Fig. 1d). The experiment was repeated with SynCom-223 and SynCom-137 and revealed similar results (Source Data Fig. 1d), which shows that the strains in SynCom-137 are representative for the ASVs.

Because the differences in relative abundance suggest a distinct adaptation to an endophytic or epiphytic lifestyle of some strains, we conducted leaf infiltration experiments to test bacterial growth in the apoplast. Indeed, the endophytic-enriched strains Xanthomonas Leaf131 and Pseudomonas Leaf59 were able to grow to a greater extent (37-fold and 3.5-fold increase, respectively) after low-dose infiltration into leaves of soil-grown Col-0 plants, in contrast to the depleted strain Pseudorhodoferax Leaf265 (2.5-fold increase) (Extended Data Fig. 2), confirming the potential of the former, especially Xanthomonas Leaf131, to thrive as endophytes. Importantly, no disease symptoms were observed after leaf infiltration, and bacterial titers were lower than those of the model pathogen Pseudomonas syringae pv. tomato DC3000 (Pst) (Extended Data Fig. 2).

Extended Data Fig. 2

A. thaliana Col-0 leaves after infiltration with bacterial strains.

a) Bacterial load measured as colony forming units (CFU)/cm2 at 0 and 5 days post infiltration (dpi). Xanthomonas Leaf131, Pseudomonas Leaf59, Pseudorhodoferax Leaf265, Pseudomonas syringae pv. tomato DC3000 (Pst) hrcC- and wild-type Pst were infiltrated at OD=0.002 (~106 CFU/ml) into leaves of soil-grown, four-week-old Col-0 plants. Infiltrated leaves were harvested, surface sterilized in 70% ethanol and homogenized before serial dilution plating. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. Statistical differences were calculated with two-tailed Mann-Whitney U-test (0 dpi, n=4; 5 dpi, n=8; ns, not significant; * p<0.05, ** p<0.01). b) Photographs of leaves from Col-0 five days after treatments as described above. White bar indicates one cm.

Screening A. thaliana mutants for an altered bacterial community

To identify host genetic factors that impact the bacterial community composition in the phyllosphere, we tested A. thaliana mutants with defects in distinct parts of the immune system. These included mutants for immune signaling or regulators that are involved in pattern-triggered immunity (PTI)[21], effector-triggered immunity (ETI)[22] and homeostasis of defense-associated plant hormones[6] (Supplementary Table).

A. thaliana Col-0 wild-type and different mutants were inoculated with SynCom-222, and the community composition was determined by 16S rDNA amplicon sequencing. We examined the impact of the host genotype on the community composition in a principal component analysis (PCA), which quantifies how much variation in the strain (or ASV) abundances can be explained by individual genotypes, and determined the effect size and statistical significance by PERMANOVA. Bacteria inhabiting the apoplast are in closer contact with host cells and, as such, may be subject to stronger host control than bacteria on the leaf surface. Therefore, we primarily focused our analysis on the bacterial community in the endosphere compartment.

The largest impact on the endophytic community was observed in the A. thaliana mutant rbohD[23], which lacks the plasma membrane-localized NADPH oxidase Respiratory Burst Oxidase Homologue (RBOH) D, with a 25% effect size (Fig. 2a-c). Among the genotypes with defects in PTI signaling, the effect of rbohD was larger than that of other genotypes with defects in pattern recognition receptors (e.g., fls2efrcerk1 (fec)[24], lym3lym1[25], pepr1pepr2[26]) or coreceptors (e.g., bak1-5bkk1cerk1 (bbc)[27]). The genotypes compromised in plant hormone biosynthesis or perception had an effect size between 17% by jar1 (JA signaling)[28] and 7% by ein2 (ethylene signaling)[29]. Notably, the higher-order mutants with defects in multiple hormone systems, jar1ein2npr1 (jen)[30] and dde2ein2pad4sid2 (deps)[31], did not have a stronger impact on the community composition than the single mutants tested (Fig. 2a). The salicylic acid-dependent immunoregulatory mutant npr1[32] and the stomatal opening regulator ost1[33] had a minor impact on the endophytic community composition. The quadruple mutants min7fec (mfec) and min7bbc (mbbc)[27] also affected community assembly, corroborating a recent study reporting that these mutants are involved in bacterial microbiota homeostasis[14]. The hyperimmune genotypes cpr5[34] and cpr6[35], which have constitutively high salicylic acid (SA) levels, also had an impact on the community.

Figure 2

Effect of plant genotype on the leaf endosphere community.

a) The endosphere bacterial community (SynCom-222) of each plant genotype was compared to Col-0 wild-type plants in principal component analysis (PCA, n=12) followed by PERMANOVA (permutations=10,000), and the effect size of the genotype was plotted in decreasing order. b) PCA of bacterial endosphere communities in Col-0 (blue), bbc (orange), jar1 (pink) and rbohD (green). PC1 and PC2 are principal components PC1 and PC2 with their explained variance (%). c) Exemplary PCA of endosphere communities in rbohD, jar1 and bbc. Effect size represents variance explained by genotype, and statistical significance is expressed with p-values determined by PERMANOVA (Benjamini-Hochberg adjusted, n=12). PCA results for all genotypes can be found in Source Data of Fig. 2a. d) Relative abundance of phyla (and classes for Proteobacteria) of endosphere bacteria on the indicated plant genotypes. Genotypes are ordered by decreasing abundance of Gammaproteobacteria. Asterisks (or hash tags for Firmicutes) denote significant differences in taxa on a genotype compared to Col-0 in a two-sided t-test (p<0.05, Benjamini-Hochberg adjusted, n=12). e) Community spread of genotypes rbohD, jar1 and bbc relative to Col-0 in PCA. The community spread was calculated as the Euclidean distance of data points to the centroid in PCA. Distances of genotypes were normalized to the median distance of Col-0 as the z-score. Community spread was calculated for SynCom-222, SynCom-223 and SynCom-137. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. Statistical differences were determined by two-sided t-test (n=12; ns, not significant; *p<0.05; **p<0.01; ***p<0.001).

In fact, all genotypes with large effects on the endosphere community also affected the entire phyllosphere community (Extended Data Fig. 3a). The malectin-like receptor-like kinase Feronia knockout fer[36] had the highest effect size (28%) on the phyllosphere community but was not tested for the endosphere community due to poor germination and impaired development[37]. In addition to the genotype effect on data variance in a PCA, we also analyzed the variability of community composition upon colonization of immunity mutants (Fig. 2e). The microbiota variability in rbohD plants was higher than that in wild-type plants, while it was equally variable in jar1 and wild-type plants in all SynCom experiments.

Extended Data Fig. 3

Genotype effect on bacterial community in phyllosphere.

a) Effect of plant genotype on phyllosphere community. The bacterial community of each genotype was compared to Col-0 wild-type in principal component analysis (PCA, n=12) followed by PERMANOVA (permutations=10,000), and the effect size of the genotype was plotted in decreasing order. Effect size represents variance explained by genotype (p-value<0.05, Benjamini-Hochberg adjusted; n=12). b) Relative abundance of phyllosphere bacteria on the indicated plant genotypes. Asterisks (or hash tags for Firmicutes) denote significant differences in taxa on genotypes compared to Col-0 in a two-sided t-test (p<0.05, Benjamini-Hochberg adjusted, n=12).

Analysis of the community composition at higher taxonomic levels, such as phylum (or class for Proteobacteria), between individual genotypes revealed that several taxa were significantly altered in their relative abundance in rbohD, mbbc, bbc, mfec, cpr5, jar1 and fer (fer tested only for phyllosphere) compared to wild-type plants (Fig. 2d and Extended Data Fig. 3b).

Hierarchical clustering of the differential strain abundances between individual genotypes and Col-0 of endosphere-enriched strains grouped rbohD separately from the other genotypes, but no other significant clusters were observed (Extended Data Fig. 4a). Cluster analysis of the phyllosphere community profiles revealed clustering of all genotypes but fer and cpr5 separately (Extended Data Fig. 4b).

Extended Data Fig. 4

Overview and clustering of community profiles on genotypes versus Col-0.

a) Strain changes in endosphere communities are displayed in a heatmap as log2 fold-change of strains in the endosphere of the individual genotypes versus Col-0 (columns). Strains or ASVs (indicated by superscript circle) of SynCom-222 are ordered and colored by phylogeny. Hierarchical clustering (R command hclust, method "single") of genotypes was performed on an Euclidean distance matrix of log2 fold-changes between test conditions and controls. b) Strain changes in phyllosphere communities shown in the heatmap with genotype clustering as described above. Differential strain abundance was calculated using DESeq2, and statistical significance was expressed with p-values (two-sided Wald test, Benjamini-Hochberg adjusted): the black cell rectangle highlights significant changes p<0.05.

Overall, we found striking genotype effects on the leaf bacterial community in both endosphere-enriched and total phyllosphere samples.

rbohD assembles a dysbiotic microbiota

Because rbohD showed the largest impact on the endosphere community with an effect size of 25% (Fig. 2a), compared to 17% and lower for other immunity mutants, we examined the impact of RBOHD on the microbiota more closely. RBOHD is part of a family of NADPH oxidases that generate extracellular reactive oxygen species (ROS)[38] and function in diverse physiological processes, including biotic and abiotic stress signaling[39,40]. In particular, RBOHD is regulated through activated immune receptors and is responsible for inducible apoplastic ROS production during PTI and ETI[23,41,42], playing critical roles in microbe-associated molecular pattern (MAMP)-induced stomatal closure[43], cell wall damage-induced lignification[44,45] and resistance to pathogens[23,46,47].

In addition to assembling an altered endophytic community, as described above, rbohD plants developed disease symptoms after inoculation with a SynCom (Fig. 3a), while axenic rbohD plants were indistinguishable from the wild-type. SynCom-induced disease was also reflected in the reduced average plant weight (Fig. 3b) and points towards a dysbiosis phenotype[48]. Although symptoms could be observed in most rbohD plants, the disease severity varied (Fig. 3a). We classified SynCom-induced disease into five categories (1, healthy, to 5, dead), and increasing disease severity matched with decreasing plant fresh weight (Extended Data Fig. 5b). Next, we examined the dysbiotic microbiota of rbohD at the strain level by using the SynCom-137. A defined set of bacteria was affected in their relative abundance in the endosphere and phyllosphere compared to Col-0 (Fig. 3c). Mostly Gammaproteobacteria (Xanthomonas Leaf131, Serratia spp. Leaf50 and Leaf51, Pseudomonas spp. Leaf83, Leaf58, Leaf127, Leaf15 and Leaf434), an Alphaproteobacterium (Sphingobium Leaf26) and two Actinobacteria (Sanguibacter Leaf3 and Microbacterium Leaf161) showed a significant increase in relative abundance in rbohD compared to wild-type plants (Fig. 3c) and were among the top colonizers in the endosphere of rbohD (Fig. 3d). Notably, the same strains were consistently enriched in rbohD in three SynCom experiments, which underlines the robustness of the observed effect. SynCom-137 recapitulated the results of the full community, thus showing that single strains were representative of the ASV with respect to their behavior in rbohD. Furthermore, most of the affected strains were increased in both the endosphere and phyllosphere of rbohD (Fig. 3c).

Figure 3

The plant mutant rbohD shows a dysbiosis phenotype and assembles a microbiota enriched in Gammaproteobacteria.

a) Representative pictures of five-week-old Col-0 and rbohD inoculated with SynCom-137. White bar indicates one cm. b) Fresh weight of aboveground plant tissue of Col-0 and rbohD inoculated with SynCom-137. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range (n= 12; two-tailed Mann-Whitney U-test; ***, p<0.001). c) Heatmap shows log2 fold-changes of strains in phylogenetic order in rbohD compared to Col-0 wild-type plants. Columns show phyllosphere and endosphere samples from plants inoculated with either SynCom-222, SynCom-223 or SynCom-137. Black rectangles show significant changes, p<0.05 (n=12; Wald test, Benjamini-Hochberg adjusted). Strains (or ASVs indicated by superscript circles) are colored according to phylogeny. Red dots highlight enriched strains in rbohD across multiple experiments. d) Volcano plot shows the relative abundance of rbohD-enriched strains of SynCom-137 in rbohD endosphere (log2-fold-changes in rbohD compared to Col-0 with adjusted p-value <0.05).

Extended Data Figure 5

Microbiota-induced disease in rbohD is linked to the enrichment of Xanthomonadaceae.

a) Screening of plant phenotypes and disease symptoms in rbohD after colonization with individual bacterial strains. Germ-free, 10-day-old rbohD seedlings were inoculated with a bacterial suspension (OD600 ranging from 0.02 to 0.08), and the plant phenotype was examined after 3.5 weeks (n=8). Inoculated plants showed a variety of phenotypes, e.g., stunted plants, necrotic lesions, curled leaves or dead plants. Phenotype-inducing strains were selected based on symptoms and are highlighted by red rectangles in the phylogeny. b) Disease index was assessed in rbohD plants: 1, healthy; 2, mild symptoms on individual leaves; 3, stronger symptoms on multiple leaves; 4, strong symptoms on whole plant; 5, severe symptoms or dead plant. The graph displays the plant fresh weight (mg) of SynCom-137-inoculated rbohD plants (n=20) per disease category. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. c) Bacterial load in SynCom-137-inoculated Col-0 and rbohD plants measured as colony forming units (CFU) per gram of plant fresh weight isolated from the endosphere and phyllosphere (n=16; two-tailed Mann-Whitney U-test; ns, not significant; **, p<0.01). Bacterial load in Col-0 and rbohD plants inoculated with the indicated SynCom represented by qPCR of the bacterial 16S rDNA gene relative to plant gene. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. Statistical significance was calculated by two-tailed Mann-Whitney U-test (n=3, each pool of four DNA samples; ns, not significant). qPCR data for SynCom-222 in Col-0 are the same as in Figure 1c. d) Correlation analysis between the relative abundance of Xanthomonadaceae in phyllosphere samples of rbohD inoculated with the indicated SynCom and plant fresh weight. The Spearman coefficient ρ (p-value<0.01) was calculated using ggscatter command (ggpubr, R), grey area shows 95% confidence interval of regression line in green. Correlation data on endosphere samples was not possible due to bulk surface sterilization of plants. e) ROS accumulation was measured in leaf discs with a luminol-based assay after treatment with extracts from heat-killed bacteria. ROS production was recorded for 45 min, and luminescence counts were integrated over time. ROS triggered by individual treatments were normalized to ROS production by 10 nM flg22. Normalized ROS accumulation is shown for each bacterial strain. Barplots show mean and error bars show standard deviation (n=16; combined data from two independent experiments). Red dots indicate rbohD-enriched strains.

To investigate whether disease in rbohD is accompanied by an increased bacterial load, we assessed the bacterial cell numbers in the phyllosphere and endosphere. Overall bacterial abundance measured by CFU and qPCR normalized to plant weight or plant genomes, respectively, was not significantly increased in rbohD compared to Col-0 across multiple experiments (Extended Data Fig. 5c). This finding is in line with previous plant pathology studies that did not observe a correlation between disease severity and pathogen titer[49-51], although a correlation had been described in other studies[20,52-54]. In fact, the outcome might depend on the specific plant-microbe interaction, the microbiota background, and environmental conditions. To test whether certain bacteria are linked to disease in rbohD, we correlated plant weight and the relative abundance of bacterial families. Indeed, we observed a strong negative correlation for Xanthomonadaceae (Spearman coefficient ρ=-0.82 and ρ=-0.83 for experiments with SynCom-222 and SynCom-137, respectively) in the phyllosphere samples of rbohD (Extended Data Fig. 5d).

As mentioned above, in our experiments, rbohD had a larger effect on the endosphere community profile and on individual strains than other PTI mutants, such as fec, bbc, mfec, and mbbc (Fig. 2a, Extended Data Fig. 6a). For example, beside Xanthomonas Leaf131 that was more abundant also in other PTI mutants, additional Pseudomonas strains, Leaf83, Leaf58, Leaf127, Leaf15, Leaf434, showed increased abundance only in rbohD. In addition to RBOHD, RBOHF also plays a role in ROS production and disease resistance, although to a lesser extent[23,47,55]. We tested the double knockout mutant rbohDrbohF[23] and found that it had an even stronger impact on SynCom-137 than rbohD alone. For example multiple strains were specifically enriched in the endosphere of rbohDrbohF but not in rbohD (Extended Data Fig. 6b), such as Rhizobium Leaf311 and Brevundimonas Leaf168 (Alphaproteobacteria); Pseudorhodoferax Leaf265 and Acidovorax Leaf160 (Betaproteobacteria); Serratia Leaf50, Erwinia Leaf53 and Pseudomonas Leaf59 (Gammaproteobacteria); Arthrobacter strains Leaf145 and Leaf141, and Plantibacter Leaf314 (Actinobacteria).

Extended Data Fig. 6

The community assembly of rbohD and rbohDrbohF substantially differs from that of other PTI genotypes.

The heatmap shows log2 fold-changes in strain abundance on different genotypes compared to Col-0 wild-type. Columns show endosphere and phyllosphere samples from plants inoculated with either a) SynCom-222 or b) SynCom-223 or SynCom-137. Strains (or ASVs indicated by superscript circles) are ordered and colored by phylogeny. Statistical significance of differential strain abundances was calculated with the two-sided Wald test (DESeq2 package, R) and highlighted with a black cell rectangle for p<0.05 (n=12; Benjamini-Hochberg adjusted).

The substantial impact of rbohD and rbohDrbohF compared to other PTI-associated genotypes in our experiments (Fig. 2a and Extended Data Fig. 6) could be explained by the central roles of NADPH oxidases in converging signaling pathways in immunity and beyond[40,56]. Overall, our results suggest that dysbiosis of rbohD is caused by a shift in community structure with an increased abundance of certain taxa dominating the community.

The bacterial leaf community contains opportunistic pathogens

Many plant-associated bacteria carry immunogenic MAMPs and can potentially trigger immunity[8]. We selected phylogenetically diverse species, including rbohD-enriched strains, and tested whether they can elicit an apoplastic ROS burst. Extracts of most rbohD-enriched strains and other strains triggered plant ROS, suggesting that ROS-eliciting potential is common, especially among Beta- and Gammaproteobacteria (Extended Data Fig. 5e).

To identify the causal agent(s) of SynCom-induced disease in rbohD, we tested the pathogenicity of the individual SynCom-137 strains. We inoculated single strains onto germ-free rbohD seedlings and monitored the plant phenotype. Based on this screening we identified 28 phenotype-inducing strains that caused plant symptoms, such as necrotic lesions or curled leaves, indicating a strong host-microbe interaction (Extended Data Fig. 5a), with Xanthomonas Leaf131 triggering the most severe disease symptoms when inoculated alone (Fig. 4a and 4b). The At-LSPHERE strain collection contains another Xanthomonas strain, Leaf148, which shares the same 16S rDNA amplicon sequence but has only an 87% average nucleotide identity (ANI) with Leaf131. We tested both Xanthomonas strains for their phenotype on Col-0, rbohD and rbohD/RBOHD complementation lines. While rbohD plants became diseased or died after inoculation with Leaf131 or Leaf148, Col-0 and the complemented rbohD/RBOHD looked healthier and had higher plant weights (Fig. 4a and 4b). In addition, bacterial extracts of both Xanthomonas strains triggered a plant ROS burst, indicating that the plant is able to perceive these bacteria (Fig. 4c). Moreover, infiltration of Xanthomonas Leaf131 into leaves of soil-grown Col-0 wild-type and rbohD plants resulted in mild disease symptoms only in rbohD (Fig. 4d). We concluded that Leaf131 and Leaf148 are opportunistic pathogens that develop their deleterious potential on immunocompromised rbohD plants.

Figure 4

Xanthomonas Leaf131 and Leaf148 are opportunistic pathogens in immunocompromised rbohD.

a) Pictures of Col-0, rbohD and rbohD/RBOHD plants inoculated with Xanthomonas Leaf131 or Leaf148. Germ-free 10-day-old seedlings were inoculated with bacterial solution at OD=0.02 (~107 CFU/ml). White bar indicates one cm. b) PTI-associated phosphorylation sites of RBOHD are involved in resistance to opportunistic pathogens and in microbiota homeostasis. Fresh weight of phyllosphere plant tissue from Col-0, rbohD, rbohD/RBOHD, rbohD/RBOHD-S39AS339AS343A, rbohDrbohF/RBOHD, rbohDrbohF/RBOHD-S343A-S347A and rbohDrbohF after inoculation with either single strain Xanthomonas Leaf131 or Leaf148 or SynCom-137. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. Brackets denote comparison groups, and the number above indicates the p-value of two-tailed Mann-Whitney U-test and the number below the fold-change in plant weight (n=10). c) Apoplastic ROS production was measured over time in leaf discs with a luminol-based assay after treatment with water, 10 nM flg22 or extracts from heat-killed Xanthomonas Leaf131 and Leaf148 adjusted to the indicated total protein concentrations. ROS production is displayed as relative light units (RLU) and curves show the mean (n=8) and error bars the standard error. d) Photographs of leaves from Col-0 and rbohD five days after infiltration with 10 mM MgCl2 or Xanthomonas Leaf131 at OD=0.002 (~106 CFU/ml). White bar indicates one cm.

RBOHD activation in A. thaliana occurs through N-terminal phosphorylation by receptor-like cytoplasmic kinases, such as Botrytis-Induced Kinase 1 (BIK1)[21,42]. Multiple complementary and overlapping RBOHD phosphorylation sites have been reported to be critical for PTI signaling[56-60]. We probed whether RBOHD phosphorylation sites can be linked to microbiota-induced disease. After validation of impaired MAMP-induced ROS production in RBOHD mutants (Extended Data Fig. 7a)[57,58], we tested whether RBOHD with mutated phosphorylation sites, i.e., rbohD/RBOHD-S39A-S339A-S343A and rbohDrbohF/RBOHD-S343A-S347A, could complement the rbohD disease phenotype after inoculation with SynCom-137 or single strains of Xanthomonas Leaf131 and Leaf148. The presence of SynCom-137 caused a significant weight reduction of rbohD and rbohDrbohF and, to a lesser extent, of the phosphorylation site mutants compared to their respective controls (Fig. 4b). The intermediate disease phenotype of the RBOHD phosphorylation site mutants compared to the knockout lines rbohD and rbohDrbohF could be due to regulatory redundancy of multiple posttranslational modifications resulting in tunable activation of RBOHD and ROS production[56-60]. Inoculation with Xanthomonas Leaf131 and Leaf148 caused disease symptoms and reduced the average plant weight of rbohD and the two phosphorylation site mutants compared to their control lines Col-0, rbohD/RBOHD, and rbohDrbohF/RBOHD, respectively (Fig. 4b and Extended Data Fig. 7b). In contrast, inoculation with Pst reduced the plant weight to a similar level on all genotypes irrespective of RBOHD presence, while colonization by the nonvirulent Pst hrcC- did not affect the plant weight of the different genotypes (Extended Data Fig. 7b).

Extended Data Fig. 7

PTI-associated phosphorylation sites of RBOHD are involved in resistance to opportunistic pathogens.

a) ROS production of A. thaliana Col-0, rbohD, rbohF, rbohD/RBOHD, rbohD/RBOHD-S39AS339AS343A, rbohD/RBOHD-S343A-S347A, rbohD/RBOHD-S343A-S347A, rbohDrbohF/RBOHD, rbohDrbohF/RBOHD-S343A-S347A and rbohDrbohF, after treatment with 100 nM flg22. ROS production was measured in leaf discs from soil-grown plants with a luminol-based assay and expressed as integrated luminescence over 45 min (AU, arbitrary units). Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. b) Fresh weight of Col-0, rbohD, rbohD/RBOHD, and rbohD/RBOHD-S343A-S347A inoculated with individual strains of Xanthomonas spp. Leaf131, Leaf148, Pst hrcC-, and Pst wild-type or mock-inoculated with 10 mM MgCl2. Germ-free 10-day-old seedlings were inoculated with OD=0.02 of the respective strains. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. Significant differences were calculated between the mutants and their respective controls (n=20; two-tailed Mann-Whitney U-test). Brackets above bar plots indicate comparison groups with p-values displayed above and fold-change below. Data from two independent experiments are shown in separate graphs.

In summary, our data highlight the importance of RBOHD and its regulation via immunity-associated phosphorylation sites for plant health during colonization with opportunistic Xanthomonas pathogens and for microbiota homeostasis.

Dysbiosis of rbohD is driven by opportunistic pathogens

Due to the enrichment of Xanthomonas Leaf 131 in rbohD (Fig. 3c and 3d) and the phenotype it caused in mono-association (Fig. 4a), we next wondered whether the opportunistic pathogen was sufficient to explain the phenotype in a community context and speculated on the relative contribution of plant immunity and the microbiota to plant health. We designed several synthetic communities by removing either the opportunistic pathogen Xanthomonas Leaf131 or all rbohD-enriched and plant phenotype-inducing (REPI) strains (Supplementary Table). We inoculated germ-free Col-0, rbohD and rbohD/RBOHD with the selected bacterial mixes. The full community control, i.e., SynCom-137, reduced plant health and weight compared to axenic control plants only in rbohD, not in Col-0 and rbohD/RBOHD (Fig. 5), confirming our previous observation that SynCom causes disease in rbohD (Fig. 3a). Xanthomonas Leaf131 single inoculation was detrimental to rbohD and reduced plant weight in Col-0 and rbohD/RBOHD, corroborating its pathogenic potential (Fig. 4). We then mimicked the dysbiotic community observed in rbohD plants by assembling SynCom-REPI, consisting of the 32 strains for which we observed phenotypes or enrichment in rbohD (Fig. 3c and Extended Data Fig. 5a). SynCom-REPI severely reduced the plant weight of rbohD and had a mild effect on Col-0. Remarkably, when Leaf131 or all SynCom-REPI members were removed from SynCom-137, the plant weight difference between Col-0 and rbohD disappeared and both genotypes had phenotypes that were not distinguishable from axenic control plants (Fig. 5). Because Xanthomonas Leaf131 is also part of the REPI strains, we assumed that it is responsible for the negative effect on rbohD. Thus, we repeated the dropout SynCom experiment with an additional condition consisting of Leaf131 dropout from SynCom-REPI. Again, we observed that whenever the SynCom did not contain Leaf131, plant health and weight were increased in Col-0, rbohD and rbohD/RBOHD (Extended Data Fig. 8). This shows that Xanthomonas Leaf131 is responsible and sufficient for SynCom-induced disease in rbohD. Importantly, a diverse bacterial community can alleviate the detrimental effects of the opportunistic pathogen Xanthomonas Leaf131 on immunocompromised rbohD; however, a healthy plant requires functional RBOHD (Fig. 5).

Figure 5

Xanthomonas Leaf131 causes dysbiosis in rbohD.

a) Fresh weight of Col-0, rbohD and rbohD/RBOHD plants inoculated with either 10 mM MgCl2 (axenic), SynCom-137, single strain Leaf131, SynCom-REPI (containing 32 rbohD-enriched and phenotype-inducing strains), SynCom-137 without (w/o) Leaf131, SynCom-137 w/o REPI. Germ-free 10-day-old seedlings were inoculated with OD=0.02 of SynCom-137, and other inocula with lower strain numbers were diluted with 10 mM MgCl2 to obtain equal amounts of cells from each strain in the inoculum. Box plots show the median with upper and lower quartiles and whiskers present 1.5x interquartile range. Significant differences were calculated with ANOVA and Tukey’s HSD post-hoc test (n=20, letters indicate significance groups, α=0.05). b) Pictures of plants treated as indicated above. Genotypes Col-0, rbohD and rbohD/RBOHD were grown in the same microbox in rows, and different genotypes are highlighted by colored arrows (Col-0, blue; rbohD/RBOHD, light blue; rbohD, green). White bar indicates one centimeter. The experiment was repeated twice with additional treatments (see Extended Data Figure 8).

Extended Data Fig. 8

Leaf131 is required and sufficient for disease in rbohD.

Two independent replicate experiments with dropout synthetic communities as presented in Figure 5. In addition, SynCom-REPI without (w/o) Leaf131 was tested and compared to SynCom-REPI. Significant differences were calculated with ANOVA and Tukey's HSD post-hoc test (n=20, letters indicate significance groups, α=0.05).

Discussion

The assembly and maintenance of a healthy microbiome is crucial for host fitness[10,48]. Prior knowledge about the function of key plant immunity genes makes it now possible to test how these genes influence not only pathogens, but the entire plant microbiota and interactions within these. By using a reverse genetic screen, we found multiple plant immunity mutants that impact bacterial leaf community assembly. In particular, we uncovered the pivotal role of the NADPH oxidase RBOHD to prevent microbiota dysbiosis in the phyllosphere. NADPH oxidases are evolutionarily conserved in multicellular eukaryotes with functions in defense, development and redox-dependent signaling[61,62]. Strikingly, NADPH oxidases play a major role in gut microbe homeostasis in mammals[63,64], zebrafish[65] and insects[66,67], highlighting the importance of understanding the interplay between ROS production, host immunity and microbiota.

Likewise to dysbiosis in the human intestine, which has been characterized by an increase of relative abundance of Gammaproteobacteria[68] and dysbiotic communities being more diverse than communities in healthy patients[69], we also found an enrichment of multiple Gammaproteobacteria (Fig. 3) and a higher variability among the microbiota between individual rbohD plants compared to wild-type (Fig. 2e). However, in contrast to observations made in the gut, we identified Xanthomonas spp. as the most conspicuous, acting as opportunistic pathogens when inoculated on plants alone (Fig. 4).

Xanthomonas spp. are common plant pathogens and are ubiquitously found in the phyllosphere microbiota on a variety of plants[70-74]. The type-3 secretion system (T3SS) represents a major virulence determinant in pathogenic Xanthomonas[74]. A preliminary genome analysis of Xanthomonas Leaf131 and Leaf148 revealed the absence of a T3SS17 which is in accordance with their opportunistic lifestyle. Note also that the opportunistic Xanthomonas pathogens identified here have been isolated from phenotypically asymptomatic A. thaliana plants in their natural environment[17]. Correspondingly, numerous ecological studies have detected pathogenic species, including Xanthomonas, as members of the plant microbiota that are asymptomatic for the host in a community context but cause disease during mono-association[53,75-77]. Our data support the concept of conditional pathogens[78], which defines species with pathogenic potential as an integral part of the microbiota that are kept in check by other community members and the host, and underlines the importance of both the genotype and the microbiota to plant protection (Fig. 5).

The mechanisms by which the microbiota reduces the prevalence of opportunistic pathogens might be manifold and remain to be elucidated. Screening of bipartite interactions between At-LSPHERE strains in vitro did not reveal growth inhibition of Xanthomonas Leaf131 and Leaf148 by other community members[79]; nonetheless, the plant could be protected through a concerted action of multiple commensals or indirectly via stimulation of plant immunity[80-84]. Furthermore, Xanthomonas WCS2014-23, a close relative to the opportunistic pathogen Xanthomonas Leaf148 (97.5% ANI), has been isolated from A. thaliana roots and interacts with other bacteria to promote plant growth and to induce disease resistance[85].

Apoplastic ROS production by RBOHD is a typical response during plant-pathogen interactions. As the community shift in rbohD could be observed in endosphere and phyllosphere samples, the internal and external bacterial populations might form a continuum[1]. The absence of RBOHD could facilitate endophytic growth of certain taxa, which then relocate to the leaf surface. RBOHD-produced ROS can act directly as antimicrobials as well as in local and systemic signaling[23,38,40]. In our experiments, many bacteria were able to induce a ROS burst, indicating that microbiota members have the capacity to trigger PTI. Future experiments will help to disentangle whether RBOHD selectively impacts individual bacteria, for example through antimicrobial ROS, cell wall modification via cross linking or downstream immune responses[23,40,58,44], and whether pathogenic Xanthomonas act as keystone species by affecting other microbes either directly or indirectly via the plant.

Recently, Chen et al. proposed a genetic network to control the microbiota and prevent dysbiosis[14]. Our data indicate that RBOHD and RBOHF could be central parts of such a network and integrate signals from converging pathways, as rbohD and rbohDrbohF knockout lines had a more pronounced impact on the endophytic community than other immunity mutants under our experimental conditions. While Gammaproteobacteria generally benefited from gene knockouts in bbc mbbc, mfec, rbohD and rbohDrbohF, we could see overlap but also differences at the strain level between the genotypes (Fig. 2d and Extended Data Fig. 6). Notably, RBOHD has been identified as a central link between the immune signaling pathways PTI and ETI[86,87], which could explain its extraordinary effect on the microbiota. In addition, RBOHD could be a convergence point of biotic and abiotic stress signaling[39,40] that jointly regulate microbiota homeostasis. It remains to be determined how the different host components interact to control the microbiota. The impact of RBOHD on selected members of the bacterial microbiota potentially has strong implications for plant defense against pathogens as these bacteria may prime plant immunity and convey health-protective properties.

In summary, we identified a pivotal function of the NADPH oxidases RBOHD and RBOHF in maintaining microbiota homeostasis in the phyllosphere to prevent dysbiosis. Our findings emphasize the role of the commensal microbiota together with a functional immune system in the control of opportunistic pathogens.

Methods

Plant growth conditions

Gnotobiotic plants were grown in calcined clay (Diamond Pro Calcined Clay Drying Agent, Arlington, Texas, USA) supplemented with 0.5x Murashige and Skoog (MS) medium including vitamins, pH 5.8 (M0222.0050; Duchefa, Haarlem, Netherlands), in round gamma-irradiated microboxes (No. O118/80+OD118 with green filter lid; Saco2, Deinze, Belgium) as described previously[19].

A. thaliana seeds were surface sterilized using 70% EtOH for two minutes and 7% NaOCl containing 0.2% Triton X-100 for eight minutes[88]. Seeds were washed with water six times and stratified for 4 days in the dark at 4°C. Several seeds were individually placed in five positions per microbox to account for unsuccessful seed germination. Additional seedlings were removed one day before inoculation to keep five plants per microbox. Microboxes only contained plants of one genotype and treatment for the initial screening experiment with SynCom-222. Each plant was watered with 200 μl of 0.5x MS, pH 5.8, at 9, 24 and 38 days after germination. A list of plant genotypes and their references can be found in Supplementary Table.

Growth chambers (CU-41L4; Percival, Perry, Iowa, USA) had air flow from the bottom of each shelf and were equipped with full-spectrum lights (Philips Master TL-D 18 W/950 Graphica; Philips, Amsterdam, Netherlands) and UVA-/UVB-emitting light (Sylvania Reptistar F18 W/6500 K; Sylvania, Etobicoke, Ontario, Canada) to simulate natural conditions. The combined light intensity was set to 220 μmol m-2 s-1 for wavelengths of 400-700 nm (measured with PAR quantum sensor SKP215; Skye Instruments, Powys, UK) and 5.4 μmol m-2 s-1 for wavelengths of 300-400 nm (UV light measured with MU-250; Apogee, Utah, USA). The growth chambers were set to 22°C and 54% relative humidity and run under an 11-hour light cycle. Due to light activity during the day, the temperature inside microboxes increased to 24°C, and water condensates appeared on sidewalls.

Synthetic community selection

SynCom composition was selected based on the At-LSPHERE strain collection17. SynCom-222 contained 222 strains (Supplementary Table) and covered maximal phylogenetic diversity to investigate changes in the community assembly on A. thaliana genotypes. SynCom-223 was used in a follow-up experiment and contained 223 strains (SynCom-223 additionally included Serratia Leaf50, which was left out in SynCom-222 due to concerns of being a strong pathogen, which turned out not to be the case). The V5-V7 region of the 16S rDNA gene amplified with the primers 799F89 and 1193R90 allows us to distinguish 137 of these strains with 100% sequence identity representing ASVs. To be able to attribute changes in strain abundance to individual strains, we selected SynCom-137 containing one strain as a representative for each ASV. To distinguish whether the data shown in the Figures belongs to individual strains or an ASV containing multiple strains (as indicated in Supplementary Table), the ASVs with more than one strain are marked by a superscript circle after the strain name (i.e., Xanthomonas Leaf131°) in the axis labels. Experiments with SynCom-137 were always done with a single strain for each ASV.

Synthetic community mixing and plant inoculation

Bacteria were individually streaked on R-2A-agar (Sigma-Aldrich) plates containing 0.5% v/v methanol (R2A-MeOH), grown for 4 days at 22°C, restreaked on fresh R2A-MeOH plates and grown for 3 days. To prepare the SynCom inoculum, a “loopfull” of biomass from each strain was scraped off with a sterile 1-μl plastic loop and resuspended individually in 1 ml of 10 mM MgCl2. Tubes containing the resuspended cells were vortexed for 5 min to disperse cell aggregates. Strains that were not fully resuspended after vortexing were treated with TissueLyzer II (Qiagen) at 30 Hz for 2x 45 sec. The SynCom inoculum was obtained by mixing each strain in an equal volume ratio. The SynCom inoculum was adjusted to an OD600 of 0.02, and each plant was inoculated with 200 μl of bacterial solution or 10 mM MgCl2.

Aliquots of 1 ml from the SynCom mixtures were taken for DNA extraction to determine the community composition of the inoculum. Almost all strains could be detected in the SynCom inoculum mixes, and undetected strains were most likely due to insufficient sequencing depths (Extended Data Fig. 1 and Supplementary Table). In addition, a dilution series was made from all resuspended strains, and 4 μl was spotted on R2A-MeOH plates to count colony forming units (CFU) and verify the viability of each strain in the inoculum.

Dropout SynCom inocula was prepared as described previously[19]. The resuspended strains were mixed in an equal ratio for SynCom-137 without (w/o) REPI containing the least number of strains. For inoculum preparation, 10 mM MgCl2 was added to account for “missing strains” compared to SynCom-137 containing most strains. The remaining strains were added at an equal ratio to an aliquot of SynCom-137 w/o REPI to obtain SynCom-137 w/o Leaf131, and MgCl2 was added to account for missing Leaf131 in the inoculum. Similarly, SynCom-137 was assembled by proportionally adding Leaf131 to an aliquot of SynCom-137 w/o Leaf131. SynCom-137 was adjusted to an OD600 of 0.02, and all dropout SynComs as well as single strain inoculum Xanthomonas Leaf131 were diluted accordingly to ensure an equal number of cells from each strain in every inoculum mix. Different plant genotypes (n=20) were grown in the same microbox to minimize box variations, and treatments were blinded after inoculum preparation to avoid unconscious biases.

Plant sampling and surface sterilization

Plants were inoculated at 10 days after germination and harvested 5-5.5 weeks after germination. Four plants per microbox were harvested for DNA extraction (n=12-16 plants per condition). The aboveground part of the plant was cut with a scalpel under sterile conditions, the cotyledons were removed with sterilized forceps, and the plant was placed into a 2-ml screw cap sampling tube (FastDNA SPIN Kit for Soil, MP Biomedicals). Plant fresh weight was determined, and phyllosphere samples were immediately frozen in liquid nitrogen. To obtain endosphere samples, plants were placed in regular 2-ml tubes, weighed, and endophytic bacteria were enriched by surface sterilization (in bulk for each genotype) using 70% ethanol for 30 sec and 3% NaOCl for 2 min to lyse epiphytic bacterial cells and degenerate their DNA[75,91]. Plants were washed three times in sterile water, placed in 2-ml screw cap sampling tubes, frozen in liquid nitrogen and stored at -80°C until further processing. Samples were randomized after harvest to avoid batch effects during 16S rDNA gene amplification and sequencing library preparation.

To determine bacterial abundance, plants were harvested as described above, weighed and placed in 2-ml tubes filled with 200 μl of phosphate buffer, pH 7, and a metal bead (5 mm diameter). Plants were homogenized in a TissueLyzer II (Qiagen) at 30 Hz for 2x 45 sec, and a dilution series was plated on R2A-MeOH plates to determine CFU. Mock treated plants using 10 mM MgCl2 were harvested to verify axenic growth conditions.

DNA extraction and 16S rDNA amplicon sequencing

Samples were freeze-dried at -40°C and 0.12 mbar for 16 h (Alpha 2-4 LD Plus, Christ, Osterode am Harz, Germany) and pulverized at 30 Hz for 2x 45 sec (TissueLyzer II, Qiagen). DNA was extracted using the FastDNA SPIN Kit for Soil (MP Biomedicals) following the manufacturer’s instructions. The samples were transferred to DNA low-binding 96-well plates (Frame Star 96, semiskirted), and the DNA concentration was quantified using dsDNA QuantiFluor (Promega) and normalized to 1 ng/μl. The 16S rDNA amplicon library was generated as previously described[17]. PCR amplification was performed with primers 799F and 1193R using Taq polymerase (Bioron DSF Taq) in triplicate and the following program: 94°C for 2 min, followed by 25 cycles of 94°C for 30 sec, 55°C for 30 sec, and 72°C for 60 sec, and a final extension at 72°C for 10 min. Technical triplicates were combined, and 5 μl of PCR run on a 2% agarose gel to visualize the amplification product. The PCR product was cleaned up using Antarctic phosphatase (New England Biolabs, NEB) and exonuclease I (NEB) to digest unused primers. The reaction mix was incubated at 37°C for 30 min and at 85°C for 15 min to inactivate enzymes. Sample 96-well plates were spun down at 4000 g for 10 min to pellet denatured enzymes. Barcoding of 16S rDNA amplicons was performed using barcoded forward and reverse primers as previously described[17]. Barcoding PCR was done in triplicate with the following program: 94°C for 2 min, followed by 10 cycles of 94°C for 30 sec, 55°C for 30 sec, 72°C for 60 sec, and a final extension at 72°C for 10 min. Triplicate reaction mixes were pooled, and the PCR product was visualized on a 2% agarose gel. To assemble the final library, individual samples were quantified according to the band intensity of the bacterial 16S rDNA amplicon on the agarose gel picture using ImageJ to calculate equal amounts of barcoded DNA for each sample. To account for differences between agarose gels, amplicon samples from each agarose gel were first added into separate subpools. DNA from each subpool was cleaned using commercial AMPure or homemade SPRI magnetic beads according to the manufacturer’s protocol. Subpools were run on a 2% agarose gel, the bacterial 16S rDNA amplicon band was excised from the gel, and DNA was extracted using a QIAquick Gel Extraction Kit (Qiagen, Venlo, Netherlands). The DNA concentration of each subpool was measured and combined in volumes according to respective sample numbers. The pooled 16S rDNA amplicon library was cleaned twice with AMPure magnetic beads using a bead-to-DNA ratio of 0.7 to remove small DNA fragments. The amplicon length distribution of the library was assessed on a 2200 TapeStation using HS D1000 (Agilent, Santa Clara, California, USA), resulting in an effective library size of 554-643 bp. Sequencing was performed on a MiSeq desktop sequencer (Illumina) at the Genetic Diversity Centre Zurich using the MiSeq reagent kit v.3 (paired end, 2x 300 bp). Denaturation, dilution and addition of 15% PhiX to the DNA library were performed according to the manufacturer’s instructions. Custom sequencing primers were used as described previously[17].

16S rDNA amplicon data processing

16S rDNA reference sequences of SynCom strains were extracted as described previously[19] and each strain is represented by one ASV reference. Paired-end sequencing reads were merged using the USEARCH v.11.0.667-i86 linux64 command fastq_mergepairs[92], with a minimum overlap of 16 bp and a minimum identity of 90%. Merged reads were quality filtered using -fastq_filter with a maximum expected error of 1 and a minimum length of 200 bp. The reads were classified and counted using -otutab with 100% identity to 16S rDNA reference sequences to generate an initial ASV table with a count for each reference per sample. Unclassified sequences of each sample (potentially originating from additional 16S rRNA gene copies of SynCom strains, or from sample or plant contamination, or from DNA artifacts of PCR amplification or sequencing) were dereplicated (-fastx_uniques) and clustered using the UPARSE algorithm (-cluster_otus[93], with a minimum cluster size of 1 and a fixed identity threshold of 97%). De novo OTU clusters were annotated using the SILVA SSU Ref NR database (release 132) and added to the initial ASV table. The unclassified reads that could not be assigned to the reference sequences were summed and included in the final ASV table as additional variables (Supplementary Table). 16S rDNA sequencing samples of axenic plants and water controls were used to determine the origin of OTUs from unclassified reads.

16S rDNA amplicon data analysis and visualization

The final ASV table of each experiment was processed in R v.3.6.3 as described previously[19]. Briefly, the ASV table was arranged in compatible datasets, each consisting of samples from a test condition and control. To account for varying sequencing depths between samples, the ASV table was log-normalized and variance-stabilized by DESeq2 v.1.14.194. The data from each SynCom experiment were analyzed separately.

To examine the effect on individual strains between the test and control conditions, the output of DESeq2 provided log2 fold-change values and p-values (Wald test, Benjamini-Hochberg adjusted). The differential strain abundances between the test and control conditions were visualized as heatmaps or volcano plots.

To assess the overall effect on communities, PCA was performed with the transformed ASV table using the prcomp command. The effect size represents the variance explained by the compared factor and was calculated on Euclidean distances followed by a PERMANOVA to test for statistical significance using the adonis command of the package vegan v.2v.5-4.

To summarize the relative abundance of each strain or taxon in a sample, the relative abundance values were calculated by proportional normalization of each sample by its sequencing depth.

The following R packages were used during analysis and visualization: agricolae v.1.3-3[95], ape v.5.4[96], ggplot2 v.3.3.1[97], vegan v.2.5-4[98], DESeq2 v.1.14.1[94].

Phylogenetic analysis

The strain phylogeny was constructed based on full-length 16S rRNA gene sequences (see Supplementary Table for accession numbers). Sequences were aligned using SINA v.1.3.399 and SILVA SSU Ref NR database release 132. To build a maximum likelihood phylogeny from the alignment, PhyML v.3.3.2018214 was used with default parameters[100].

qPCR for bacterial abundance quantification

qPCR assays were run on a QuantStudio 7 Flex using FastStart Universal SYBR green Master mix (Roche Applied Science) in a 20 μl volume containing 300 nM of each primer and approximately 5 ng of DNA per reaction. Primer sequences have been described previously[13]. Bacterial 16S rDNA was amplified using primers 799F (5’-AACMGGATTAGATACCCKG) and 904R (5’-CCCCGTCAATTCITTTGAGTTTYAR) in a 3-step qPCR run with 95°C for 30 sec, 55°C for 60 sec, and 72°C for 60 sec for 40 cycles followed by a melting curve. The reference plant gene AT4G33380 was used to normalize the 16S rDNA copy number to the number of plant genomes. AT4G33380 was amplified using primers ExpF1 (5’-ATACAGAAACAACCACCCAAAAG) and ExpF2 (5’-ACGGTACTCCAATTTTCAAGACTC) in a 2-step qPCR run with 95°C for 30 sec and 72°C for 60 sec for 40 cycles followed by a melting curve. We tested the primers for background amplification by titrating increasing amounts of bacterial DNA to plant DNA and did not observe amplification of plant organelle 16S rDNA within the amplification range Ct12 to Ct32. In addition, linear DNA amplification was tested for both primer pairs with a dilution series of template DNA. Phyllosphere and endopshere samples from 16S rDNA amplicon sequencing were pooled (four technical replicate DNA samples per reaction) for the qPCR run. All samples were run in technical duplicates. The cycle number threshold (Ct) was set by evaluating the levels of background amplification in the water controls. The amplification efficiency (E) was determined by the following equation: E = -1+10(-1/slope). The 16S rDNA gene copy number was normalized to the plant gene and calculated as described[13]: 16S rRNA/plant gene = ECt plant/ECt 16S.

ROS burst measurements

ROS production after treatment with bacterial extracts was measured using a luminol-based assay as described previously[101]. Eight leaf discs (4 mm diameter) per treatment or A. thaliana genotype were collected, placed in white 96-well plates (Nr: 655075; Greiner Bio-One) and incubated in sterile water overnight in the dark at 22°C. The water was replaced with a solution containing 17 mg/mL luminol (Sigma‐Aldrich, St. Louis, Missouri, USA), 200μg/ml horseradish peroxidase (HRP, type VI-A; Sigma‐Aldrich) and either bacterial extracts (1:10 volume ratio of final reaction mix, i.e., 20 μg/ml) or 10 nM flg22 peptide (EZBiolab, Westfield, IN, USA). Luminescence was recorded as photon counts with a 0.5 millisecond exposure time in 90-sec intervals over a total period of 60 min in a Victor3 plate reader (Perkin Elmer, Waltham, MA, USA).

To obtain crude bacterial extracts, bacteria were grown on R2A-MeOH agar plates for three days at 22°C. Cells were resuspended in 10 mM MgCl2 with an OD600 between 0.2 and 0.8, and 100 μl of bacterial suspension was plated on R2A-MeOH agar plates. Plates were incubated at 22°C for 3-5 days until bacteria covered the plate in a continuous lawn. Bacteria were scraped off the plate and resuspended in sterile, deionized water, and the OD600 was adjusted to 10 in a 1-ml volume. Cells were lysed by boiling at 100°C for 10 min with intermediate vortexing and immediately placed on ice. The bacterial crude extract was obtained by centrifugation of the cell lysate for 10 min at 16,000 g at 4°C, and the supernatant was collected and stored in aliquots at -20°C for further analysis. Total protein concentration of bacterial extracts was determined in a BCA assay and extracts adjusted to 200 μg/ml.

Bacterial infections by leaf infiltration

A. thaliana seeds were sterilized using 70% ethanol for 1 min and stratified for 4 days at 4°C in the dark. Plants were grown in potting soil (Substrat 1, Klasmann-Deilmann, Germany) for 4-5 weeks in growth chambers (CU-41L4, Percival) fitted with full spectrum lights (Philips Master TL-D, 18 W/840) at 22°C and 60% relative humidity under an 11-hour light cycle. Bacteria were grown on R2A-MeOH agar plates for four days, restreaked and incubated for another three days at 22°C. Cells were resuspended in 10 mM MgCl2, and the OD600 was adjusted to 0.002 (~106 CFU/ml). Three to four leaves of four plants were infiltrated per treatment. Plants were covered with a plastic dome to increase humidity, and infection was monitored for 5 days.

Samples were taken from infiltrated leaves at the indicated time points. Leaves were surface sterilized by submerging them in 70% ethanol for 30 sec and subsequently washed in water three times. Leaf discs were taken with a cork borer (diameter 7.5 mm), and two leaf discs from the same plant were combined in tubes filled with 200 μl of 10 mM MgCl2 and two glass beads (3 mm diameter, Merck, Darmstadt, Germany). Leaf discs were homogenized with TissueLyzer (QIAGEN, Hilden, Germany) at 30 Hz for 2x 45 sec. The bacterial solution was plated as a dilution series on R2A-MeOH agar plates to count CFUs.

Gnotobiotic plant inoculation with single strains

Plants were grown in microboxes as described for SynCom experiments. Bacteria were grown on R2A-MeOH plates and bacterial suspension for inoculum was prepared as described above for the SynCom experiments. Germ-free plants were inoculated with 200 μl bacterial solution of OD600 of 0.02 in 10 mM MgCl2. For screening of bacteria that induce a plant phenotype during mono-association on rbohD, each inoculum was prepared from a 1:10 dilution of the resuspension of a “loopful” bacterial biomass in 1 ml of 10 mM MgCl2 (final inocula had an OD600 between 0.02 and 0.08). The plant phenotype was examined at 25 days after inoculation.