Synthetic microbiota reveal priority effects and keystone strains in the Arabidopsis phyllosphere.

Multicellular organisms, including plants, are colonized by microorganisms, some of which are beneficial to growth and health. The assembly rules for establishing plant microbiota are not well understood, and neither is the extent to which their members interact. We conducted drop-out and late introduction experiments by inoculating Arabidopsis thaliana with synthetic communities from a resource of 62 native bacterial strains to test how arrival order shapes community structure. As a read-out we tracked the relative abundance of all strains in the phyllosphere of individual plants. Our results showed that community assembly is historically contingent and subject to priority effects. Missing strains could, to various degrees, invade an already established microbiota, which was itself resistant and remained largely unaffected by latecomers. Additionally, our results indicate that individual strains of Proteobacteria (Sphingomonas, Rhizobium) and Actinobacteria (Microbacterium, Rhodococcus) have the greatest potential to affect community structure as keystone species.

Introduction

Higher organisms, including animals and plants, are hosts to a diverse microbiota. Part of such microbiota may exert beneficial or harmful effects on their host and— more generally— its phenotype. In consequence, the rules governing microbiota assembly may reveal fundamental ecological principles and are a pre-requisite to engineering microbiomes. Although studying the native microbiota composition of hosts provides valuable information and correlations, such data do not reveal causal relationships[1,2]. Probing microbial interactions outside of the host (i.e. on artificial media) may help dissect causality and interactions[3,4]; however, such systems may insufficiently capture in situ conditions (e.g. shown through proteomics[5,6]), where spatial structure and dynamic host responses can critically affect the microbiota[2,7].

Synthetic microbial community analyses in gnotobiotic systems are a valuable approach to create reproducible conditions to experimentally test microbial interactions in situ [2]. Such systems have been developed for animal[8-11] and plant models including the well-studied plant Arabidopsis thaliana [12-15]. Recently, an extensive bacterial culture collection[16] that captures a large part of the natural microbial diversity of healthy A. thaliana plants[17,18] was established. Importantly, upon in planta reconstitution, the assembled synthetic community comprising of 224 leaf-derived genome-sequenced isolates closely resembled the natural microbiota[16].

Apart from representing a formidable playground for testing fundamental ecological principles[19], the phyllosphere, which refers to the leaf-dominated aerial parts of plants, is a vast environment for ecosystem functions, with its estimated 1 billion km2 surface area worldwide[18,20]. Most of the inhabiting microbes are commensals or beneficial symbionts that affect plant fitness, for example by providing pathogen protection[21-23]. In spite of the phyllosphere's ephemeral nature, similar bacterial communities— dominated by the phyla Proteobacteria, and followed by the Actinobacteria, Bacteroidetes, and Firmicutes— are found in leaves year after year[5,24]. Such a phylogenetic conservation of community composition suggests that community assembly is not a random process but is instead governed by structuring principles, which are currently only poorly understood[2,25].

Aside from geography, plant species, genotype, and seasonal variation[12,24,26-31], ecological forces such as assembly history, i.e. the timing and order in which species arrive, may play a significant role in community structure. Several groups have shown that priority effects, the imprint of arrival order on community structure[32,33], play a role in shaping microbial communities of nectar yeasts[34], fungi[35-37], and cyanobacteria[38]. Additionally, arrival order can determine whether plant-protective bacteria or fungi are effective against pathogens[39-41]. For phyllosphere communities of greenhouse A. thaliana plants, a high variation in the relative abundance of various taxa was observed, which might indicate that stochastic colonisation and historical contingency may also play a role in structuring microbial communities in plants[42].

Here, we used a phylogenetically diverse, 62-leaf bacterial strain, gnotobiotic A. thaliana model system to conduct drop-out and late introduction experiments (Fig. 1) to shed light on the principles determining community structure in planta. Specifically, we aimed to understand the role of priority effects, i.e. how the order in which strains are introduced onto the plant affects the final outcome, and whether and to what extent single strains could drive community assembly.

Fig. 1

Concept of Drop-Out and Late Introduction Experiments Using Synthetic Microbiota

a) For the control community, axenic plants are inoculated with a 62-strain synthetic microbiota and 16S rRNA gene amplicon sequencing is used as a read-out for relative abundance. b) In order to understand how a particular group or strain ("A", "B", "C", or "D") affects other groups/strains, one group ("A") is removed to analyse the effects of its absence on the rest of the community. Subsequently, this group is introduced later in order to evaluate the extent to which late-arriving strains can invade a pre-established microbiota as well as the overall effects of arrival order on community assembly.

Results

Synthetic bacterial community structure

We designed a 62-strain community that covers essentially all OTUs previously isolated from the A. thaliana phyllosphere[16] and in which each single strain could be individually identified through 16S ribosomal RNA gene (16S rDNA) amplicon sequencing (Supplementary Table 1). This 62-strain "all" input community was made up of 32 Proteobacteria, 20 Actinobacteria, 6 Bacteroidetes, and 4 Firmicutes (Supplementary Table 2). Strains were mixed and used to inoculate A. thaliana plants in a gnotobiotic system (for inoculum validation see Supplementary Fig. 1). The relative composition and intrinsic variability of the synthetic community that established in the phyllosphere was examined three and five weeks post-inoculation. As expected[12], we confirmed that the community was stable between these timepoints (Supplementary Fig. 2; p ≥ 0.01). For initial overall community analysis, we thus combined leaf samples harvested at both timepoints resulting in a total of 48 replicates, each representing the community composition of one plant (Fig. 2).

Fig. 2

Relative Abundance in Control Community

Relative abundance of the 62-strain control community (t1 and t2 combined, n = 48, 1 independent replicate). Violin plot with swarm plot overlay and pie chart are coloured by strain phylogeny. For each strain, the horizontal line represents the median and points represent individual samples. Along the bottom, the number of replicates where a given strain was not detected are indicated by circle size and count.

The community was dominated by two strains, L68-Rhizobium and L203-Microbacterium, followed by L420-Devosia, L231-Sphingomonas, L177-Burkholderia, and L233-Rhodococcus (Fig. 2 and Supplementary Data 1, 2), all of which were detected in each of the 48 replicates. Notably, the relative abundance of a strain in the initial inoculum was not predictive of its colonisation success (Fig. 2, Supplementary Figs. 1 and 3, Supplementary Data 1). Of the 62 strains in the community, five remained undetected in all 48 samples, three of which were Bacillus and close relatives. At higher taxonomic levels, Proteobacteria was the most abundant phylum, followed by Actinobacteria, Bacteroidetes, and Firmicutes (Fig. 2). These data are in agreement with previous data on the Arabidopsis phyllosphere[16,18,29,43,44].

Drop-out and late introduction effects

In order to investigate if and to what extent priority effects influence community assembly, and to test whether late-arriving groups and strains could invade a pre-established microbiota, the most abundant phylum found in the phyllosphere— the Proteobacteria— and each of its subclasses (Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria) were individually dropped-out from the initial inoculum and introduced at a later timepoint (Fig. 1b, 3a). Due to the common evolutionary history of the proteobacterial classes, these drop-outs reflect a natural choice to allow for the monitoring of effects within and among confined bacterial groups. In order to evaluate both the effects of drop-outs on overall community structure as well as to identify the individual strains affected by the various removals, we used PERMANOVA and DESeq2, respectively.

Fig. 3

Proteobacteria Class Drop-Out and Late Introduction Experiment

a) Schematic of drop-out and late introduction experiments showing the various analysed comparisons ("I", "II", and "III") among the control, the late arrival, and the group absent communities (n = 15-24 per condition; 1 independent replicate). b) Overall effects of class drop-outs and late introduction (PERMANOVA). Rows are labeled with the test groups being compared according to panel (a). Circle size and shade represent effect size and stars denote significance. c) Strains affected by class drop-outs and late introductions (DESeq2). Column labels refer to the test groups being compared according to panel (a) for each drop-out/ late introduction condition. For mock spray control, see Supp Fig. 2b. Strains are phylogenetically clustered and heatmap colours represent the log2 fold change of the test condition/control community. p-values (Benjamini-Hochberg corrected): * = 0.05, ** = 0.01, *** = 0.001, and **** = 0.0001.

For each of the drop-out and late introduction experiments, we analysed the synthetic communities for potential effects on the rest of the community (i.e. all strains except the missing group) when a group was missing altogether (Fig. 3a, comparison "I") and when the group was introduced late (Fig. 3a, comparison "II"). We also asked whether the rest of the community established differently if a particular group was introduced early (i.e. present in the initial inoculum) or late (Fig. 3a, comparison "III"). In addition, we investigated if the late-arriving group as well as the community as a whole established differently as a result of the arrival order (early vs. late) of the group in question (Fig. 3a; comparison "III"). Because we introduced the missing strains in a MgCl2 solution, we tested whether spraying the control community with MgCl2 affected its structure, but the effect was only minor (Supplementary Fig. 2).

We found that the drop-outs and late introductions of entire groups resulted in significant (p ≤ 0.01) effect sizes based on PERMANOVA analysis (Fig. 3b and Supplementary Fig. 4). The largest effect sizes were found for the impact of arrival time (early vs. late) on the group that was introduced late (i.e. the invading group: Fig. 3a, comparison "III"; Fig. 3b, and Supplementary Fig. 4; effect sizes 19 - 62 %; p ≤ 0.002). Additionally, when considering all 62 strains, late introduction of various groups (Fig. 3a, comparison "III") produced communities that were different from those in which all strains were present from the beginning (with large effect sizes of 12 - 28 %, p = 0.0001), except in the case of the Gammaproteobacteria. When focusing only on the founding population (i.e. excluding the late-arriving group) and its response to early vs. late arrival of a drop-out group (Fig. 3a and b, comparison "III"), only small effect sizes (5 - 6 %, p ≤ 0.01) were found.

Removing any of the groups from the inoculation altered the community structure (i.e. the relative abundance contribution of each strain) of the rest of the community compared to the condition where all strains were present from the beginning (Fig. 3a, comparison "I"), with the removal of the Alphaproteobacteria (effect size = 6.4 %, p = 0.0002) and the entire Proteobacteria (effect size = 6.6 %, p = 0.0001) having a larger effect compared to the removal of the Beta- or Gammaproteobacteria (effect size < 3 %, p ≤ 0.01; Fig. 3b and Supplementary Fig. 4). Late introduction of the missing group (vs. mock spray with MgCl2; Fig. 3a, comparison "II") had no significant effect on the rest of the community (Fig. 3b and Supplementary Fig. 4).

We next assessed the contribution of individual strains to the observed community changes described above (Fig. 3c, Supplementary Fig. 5). This breakdown describes the directionality of the observed effects, i.e. positive or negative consequences (Fig. 1b). Additionally, it indicates whether strains within one phylogenetic group (and thus longer shared evolutionary history) react similarly to changes in the community or upon invasion at a later time point. The absence of Alphaproteobacteria in the early stages of colonisation had exclusively beneficial effects for most Betaproteobacteria and several Actinobacteria (Fig. 3a and c, comparison "I"). Congruently, the late introduction of Alphaproteobacteria had significant negative effects on most of the same strains, when compared to the effect of adding back the Alphaproteobacteria vs. mock spray (comparison "II"). These reciprocal results confirm negative effects between members of the Alphaproteobacteria and individual strains of the Betaproteobacteria and Actinobacteria. Interestingly, most Alphaproteobacteria were largely able to invade a pre-existing microbiota and could establish similarly whether they were present at the early stages vs. arriving late (comparison "III").

The absence of Betaproteobacteria in the early stages of colonisation (comparison "I") had only a small effect on the rest of the community, as it affected the relative abundance of only three strains (p ≤ 0.05), and their late introduction was inconsequential to the rest of the community (comparison "II"). Notably, however, the Betaproteobacteria themselves were rather sensitive to the time of their introduction (comparison "III") with an overall effect size of 62 % (Fig. 3b). Three out of eight strains (L160-Acidovorax, L220-Variovorax, and L416-Methylophilus) benefited significantly from late arrival while one strain (L177-Burkholderia) was negatively impacted (p ≤ 0.0001). The absence or late introduction of Gammaproteobacteria (comparison "I" or "III", respectively) had no significant consequences (p ≥ 0.01), with the exception that the Betaproteobacterium L416-Methylophilus benefited when Gammaproteobacteria were introduced late.

The absence of the entire Proteobacteria phylum (comparison "I") caused changes in various Actinobacteria and Bacteroidetes; however, the affected strains were not a direct sum of the strains impacted by the absence of Alpha-, Beta-, and Gammaproteobacteria, respectively. On the other hand, upon late introduction, the Proteobacteria subclasses seemed to largely react as they did in their individual class drop-outs (comparison "III"). None of the Firmicutes were affected by any of the drop-outs or late introductions; similarly, only two strains of the genus Flavobacterium in the phylum Bacteroidetes were responsive to any perturbations.

Single strain drop-out effects

Next, we tested whether the removal of individual strains from the initial 62-strain inoculum could significantly alter community assembly by acting as keystone species. In total, 25 strains (Supplementary Table 2) spanning a median relative abundance of 0.01 to 27 % were individually dropped-out from the initial inoculum and community structure was examined five weeks post-inoculation. Of all the strains tested, about one third (i.e. eight strains: L160-Acidovorax, L203-Microbacterium, L231-Sphingomonas, L233-Rhodococcus, L262-Rhizobium, L265-Pseudorhodoferax, L405-Chryseobacterium, and L416-Methylophilus) had a significant impact on community structure (Fig. 4a and Supplementary Fig. 6; effect sizes = 4 - 13%, p ≤ 0.001). With regard to effect size, the absence of four of these strains, all of which had a median relative abundance > 1%, caused large, significant effect sizes, e.g. L203-Microbacterium (13%), L231-Sphingomonas (12%), L233-Rhodococcus (10 %), and L262-Rhizobium (9 %) (p-values = 0.0001; Fig. 4b). We then asked whether all strains that showed a higher detected relative abundance also had a larger effect size when removed from the rest of the community, but the correlation was only moderate (Supplementary Fig. 7, Pearson correlation = 0.443, p = 0.027).

Fig. 4

Single Strain Drop-Outs

a) Overall effects of single strain drop-outs (PERMANOVA). Data from two independent replicates (n = 11-12 per condition per replicate) was analysed separately and together. Strains are phylogenetically clustered and coloured. Circle size and shade represent effect size and stars denote significance. p-values (Benjamini-Hochberg corrected): * = 0.05, ** = 0.01, *** = 0.001, and **** = 0.0001. b) Exemplary PCA plots of the L203-Microbacterium, L231-Sphingomonas, L233-Rhodococcus, and L262-Rhizobium drop-outs generated from the two independent replicates as well as from the combined dataset. R1: replicate 1, R2: replicate 2. c) Strains affected by the single strain drop-outs (p ≤ 0.01). Effector strains (left side) and affected strains (right side) are clustered and coloured by phylogeny. Blue lines and red lines represent negative and positive effects, respectively, inferred from the single strain drop-out experiments. Line thickness correlates linearly with fold changes.

Generally, drop-outs of single strains affected the relative abundance of other strains (Fig. 4c and Supplementary Fig. 8 and 9). Interestingly, although the removal of a given strain impacted other strains positively or negatively, the affected strains themselves reacted to all tested drop-outs in only one direction (i.e. positively or negatively, but not both). The only exception was L145-Arthrobacter, which responded positively to one drop-out but negatively to all others (Supplementary Fig. 8).

Next, we identified the strains that were affected by the various single strain drop-outs (Supplementary Fig. 10). Consistent responses (i.e. similar, significant changes in the same affected strains) were found upon removal of the strains that had the largest, significant (p ≤ 0.01) effects sizes based on the PERMANOVA analysis (Fig. 4a, b; L203-Microbacterium, L231-Sphingomonas, L233-Rhodococcus, L262-Rhizobium). Notably, these strains were found to cluster together when the experiments were combined for the DESeq2 analysis (Supplementary Fig. 8), as their absence produced similar changes in the community, predominantly by increasing the relative abundance of strains that were often undetected or had a very low median relative abundance (≤ 0.1 %; e.g. L70-Stenotrophomonas, L129-Pseudomonas, L182-Brevibacillus, L187-Exiguobacterium, L220-Variovorax, L222-Agromyces, L359-Flavobacterium). The most prominent sensitive strain across all drop-out conditions was L416-Methylophilus, with 15 of 25 drop-outs causing a significant (p ≤ 0.01) increase in its relative abundance. Because these sensitive strains had low relative abundance, we also tested whether there was a correlation between the relative abundance of a strain and the likelihood that it was impacted by a drop-out, but this was not the case (Supplementary Fig. 11; Pearson correlation = -0.180, p = 0.161).

The single-strain drop-out data provided a unique opportunity to create a causal network, as it is known which strain (i.e. the dropped-out strain) affected other strains (Fig. 5). The network consists of effector and affected nodes[45] (here strains) connected by edges (i.e. direct or indirect interactions), which are illustrated by directional arrows. The network based on the combined dataset was strictly hierarchical; i.e. no bidirectional edges were observed for the tested strains (Fig. 5). As follows from the nature of the effects described above, the network was predominantly inhibitory and the keystone strains which were identified as having a reproducible, significant impact on the community based on the PERMANOVA analysis (Fig. 4A; strains L203-Microbacterium, L231-Sphingomonas, L233-Rhodococcus, L262-Rhizobium) affected many strains (≥ 10 strains; p ≤ 0.01). However, with the exception of the positive influences of L262-Rhizobium and L145-Arthrobacter on L203-Microbacterium, none of these keystone strains were themselves affected by any of the other tested dropped-out strains.

Fig. 5

Causal network (p ≤ 0.01) based on single-strain drop-outs.

Nodes (strains) are coloured by phylogeny and labels refer to strain names. Rectangles represent effector strains that were individually removed, while ovals are strains that were affected by various drop-outs. Red and blue arrows (depicting positive and negative effects, respectively) represent edges (direct or indirect interactions).

We then asked whether there was a correlation between outgoing node degree (the number of strains significantly affected by a node strain) and effect size upon a node's removal and found a strong positive correlation (Pearson correlation = 0.904, p = 5.8 x 10-10; Fig. 6). Thus, we conclude that keystone strains, directly or indirectly, interact with more strains.

Fig. 6

Correlation between node out degree and effect size upon node removal

Pearson correlation (r = 0.904, p = 5.8 x 10-10) between node out degree and effect size upon node removal. Labels refer to strain names. See Fig. 2 or Supplementary Table 2 for strain names and phylogeny.

Discussion

The recent development of higher-complexity gnotobiotic systems[9,16,46] has provided the opportunity to better simulate native-like conditions in the laboratory to understand the principles governing microbial community assembly and dynamics[2,23]. Though co-occurrence networks have emerged as a useful tool for exploring bacterial interactions through correlation analysis of relative abundance data[45-49], this approach is insufficient for conclusively establishing causality or identifying keystone species[50]. Here, we moved towards this goal by conducting drop-out experiments and methodically removing groups or single strains in an effort to understand the nature and directionality of bacterial interactions (direct or indirect) in a complex community using a gnotobiotic Arabidopsis model system (Fig. 1b).

The examination of variation among 48 control samples provided first insights on community composition and intrinsic variability. The overall community structure was broadly consistent with previous results in terms of phylum and class distribution[16], despite considerable plant to plant variation (Fig. 2). As has previously been shown with the human microbiome, where even the most abundant OTU can vary nearly 5,000-fold across individual samples[51], our finding suggests that stochastic variation has the potential to play a large role in microbiome structure at the level of individual plants. However, in contrast to the stochasticity seen in human microbiome samples, where variation is expected to arise due to the intrinsic heterogeneity of human lifestyles and genetics[52,53], the compositional variation we observe in our experiments arises even when controlling for growth conditions and plant genotype. Though not investigated in detail, the stochastic variations may sometimes involve functional substitutions, e.g. by strains competing for similar substrates. For example, in the rare cases where L90-Methylobacterium or L86-Methylobacterium were the most abundant strains (n = 2 and n = 6, respectively), the other Methylobacterium and Methylophilus strains were below average relative abundance (Supplementary Data 2), with one exception. Since all of these strains presumably utilise methanol, they may compete for this carbon and energy source.

In both the single-strain and class drop-outs, there were significant changes in community structure as a result of removing strains, presumably due to changes in the interactions of the strains present. Overall, the frequency of effects suggests that interactions among strains in phyllosphere communities are common. A major conclusion is that once established, an initial community is robust, i.e. difficult to perturb (Fig. 3b). In each drop-out experiment, the initial community established without the drop-out group was not significantly altered upon introduction of the drop-out group (except perhaps Alphaproteobacteria, p = 0.03; Fig. 3a and b, comparison "II"). As previously reported in other studies[34,35], this indicates that an initially established microbiome is likely to be resistant to perturbations associated with the introduction of new species. Hence, founding taxa, i.e. microbes that colonise a plant early, are likely to have a persistent, long-term influence on community structure[54] even if the community is later exposed to other microbial populations. An interesting caveat is that, in many cases, new strains seem to be able to slip in without perturbing the pre-existing microbiota, suggesting that although the initial population is itself resistant to alteration (i.e. the founding population remains stable when new strains are introduced), it is not resistant in the sense of impenetrable (i.e. late-arriving strains can still colonise the host). Similar results were obtained when investigating invasion and exclusion patterns among co-existing Pseudomonas syringae strain pairs on leaves[55], where about half of the tested late-arriving strains could colonise the host after pre-inoculation with another strain. In our experiments, many strains that were introduced late established with a relative abundance similar to that obtained when present from the beginning (Fig. 3a and c, comparison "III"). However, as was previously reported for some Pseudomonas syringae strains[55], though common, the ability to invade a pre-established microbiota is not universal, as some strains were significantly inhibited by a late arrival. Combining these observations, our results point toward an initial competition for resources and a process of becoming established that has long-lasting effects; nevertheless, it appears that there are also unsaturated niches for newcomers to colonise even if they are not able to significantly change the pre-existing community. Analogously, late-arriving strains have also been shown to effectively colonise mice caecal communities, where it was also hypothesised that many niches were unsaturated[54].

It is interesting to note that three of the Betaproteobacteria significantly benefited from late arrival (Fig. 3a and c; comparison "III"), suggesting they can grow best after other strains have had a chance to colonise the host. Whether late-arrivers benefit due to niche modification either indirectly via changes induced in the plant or directly due to a dependency on metabolites, nutrients, or other substances originating from other bacteria is unclear. However, benefiting from late arrival appears to be the exception, as most strains either showed no significant impact from late arrival, or were inhibited by it. These findings are in line with previous research showing that inhibitory (rather than beneficial) priority effects are also predominant in other systems[34,35,38,54], and is of practical relevance upon use of biocontrol agents to prevent, for example, pathogen infection[40,41] or frost damage[56]. Overall, the class drop-out and late introduction experiments suggest that timing of bacterial arrival is an important determinant of phyllosphere community structure for incoming bacteria, and small changes in the order and timing of arrival likely play a role in shaping the microbiota.

While the class drop-out experiments revealed that community structure was sensitive to the presence or absence of many strains, the single strain drop-outs showed that a) there are single keystone taxa that play a large role in shaping phyllosphere community structure and b) keystone taxa seem to influence community structure mostly by affecting strains that have very low relative abundance (≤ 0.1%). These results are in line with previous findings highlighting the importance of the "rare biosphere" (i.e. microbes occurring at low densities[57-59]) in microbial communities, where low abundance microorganisms are thought to provide a robust reservoir of ecological function and resiliency[59] and can become more (or less) prominent in response to environmental changes[60].

About half of the strains (52 %) were not affected by any of the single strain drop-outs tested (Fig. 4c and Supplementary Fig. 8). This may reflect the diversity of the selected strains (≤ 97 % OTU similarity), with some bacteria engaging in competition or cooperation and others largely occupying niches that are relatively specialised and distinct so as to largely avoid interaction with other taxa. Alternatively, this apparent robustness might reflect a high dynamic adaptability of these strains. To distinguish between these scenarios, dedicated follow-up experiments will be required. Interestingly, strains that were significantly affected by various drop-outs always responded either positively or negatively to all the drop-outs, but not both (the exception being L145-Arthrobacter). This points towards particular strains having, for the most part, overlapping niches or requirements and more rarely relying on other strains for metabolites or nutrients. L416-Methylophilus reacted positively to the removal of most strains, which may suggest that this strain has a high niche overlap with most other strains. In either case, L416-Methylophilus itself is a key, albeit rather stochastic, component of the community, as mentioned above, as its presence/absence and relative abundance has a significant effect on community structure (Fig. 4a). Notably, the relative abundance of Methylophilus in natural A. thaliana populations has also been observed to be highly variable[16].

Several interesting observations about the nature of microbe-microbe interactions in the phyllosphere can be made based on the causal network (Fig. 5) produced through the single-strain drop-outs. For example, although incomplete, the network is predominantly inhibitory (~ 75 %), with only two strains (L262-Rhizobium and L145-Arthrobacter) being responsible for half of the positive effects. Additionally, we note that the network is hierarchical and without observed bidirectional effects. Both of these observations support the hypothesis that competition (rather than cooperation) is prevalent in phyllosphere bacterial communities.

In spite of our experimental system bearing similarity to a natural system, it is not clear to what degree priority effects dominate community assembly in natura. In our drop-out and late introduction experiments, we set up a strict inoculation timing cutoff, where the initial population arrived at t0 and the late population arrived at t1 (three weeks later). This artificial timing cutoff allowed us to study and identify priority effects, but it is likely that historical contingency is much less pronounced in nature, where the timing of strain arrival is more fluid. Nonetheless, the fact that we do see priority effects in the experimental system used here points towards historical contingency being relevant in natura. Overall, both historical contingency and deterministic factors seem to play a role in shaping phyllosphere communities. In nature, however, removal of one strain would probably have a milder effect, as it is likely that there would be another very similar strain ready to fill the vacated niche, i.e. they would be functionally redundant. Additional drop-out experiments utilising many strains within the same OTU could address this question.

Future work analysing the role of single strains in this and more complex communities needs to be done to further elucidate microbe-microbe interactions in planta. However, it is clear that synthetic community studies are a powerful tool for understanding microbiota structure and dynamics in the phyllosphere and other host-microbe systems. The observed priority effects, and in some cases resistance to late arrival, suggests that practical manipulations of microbiomes could be more successful if applied early in a host’s life cycle when the microbial community is still developing[40,41,56]. Still, further work is needed to translate knowledge of microbiome interactions into a predictable science that could be harnessed. By gaining a better understanding of microbiota assembly and dynamics, it may be possible to selectively manipulate microbiomes in ways that provide advantages to the host, for example by providing disease protection or enhancing crop yields.

Methods

Plant growth conditions

Plants were grown using calcined clay in round, gamma-sterilised microboxes (model: O118/80 + OD118) fitted with an XXL+ filter (Combiness, Nazareth, Belgium)[13]. Calcined clay (Diamond Pro Calcined Clay Drying Agent, Arlington, USA), used as an inert soil substitute, was washed with dH2O until runoff was almost clear (10 washes). Extra water was then poured off, and the clay was dried (open trays) and subsequently heat sterilised (closed trays) each for 3 hours at 200 °C. 175 mL clay and 75 mL 0.5x Murashige and Skoog (MS) medium (Duchefa, Haarlem, Netherlands; pH 7, including vitamins, without sucrose) were added to each microbox, mixed, and softly flattened. A. thaliana Col-0 seeds were surface sterilised [21] and stratified for 4 days at 4 °C. In order to ensure that each box contained four plants at the time of inoculation, extra seeds were placed in each box and removed prior to inoculation. Plants were watered weekly with 8 mL, 0.1x MS medium (Duchefa, Haarlem, Netherlands; pH 7, including vitamins, without sucrose) per microbox.

Growth chambers (CU-41L4, Percival, Perry, USA) were fitted with a mixture of full spectrum lights (Philips Master TL-D 18W/950 Graphica) and lights emitting some UVA and UVB light (Sylvania Reptistar F18W/6500K) in order to simulate natural conditions. Combined light intensity was set to 190-200 µmol/m2/s for wavelengths of 400-700 nm (PAR) and 4-5 µmol/m2/s for wavelengths of 280-400 nm (UV). Growth chambers were set to 22 °C, 54 % humidity, and 11 hour photoperiod.

62-strain synthetic community selection

The V5-V7 region of the 16S rDNA gene amplified with primers 799F[61] and 1193R[44,62] was used to cluster 203 strains previously isolated from native A. thaliana [16] at ≥ 97 % similarity using mothur (v. 1.36.1)[63]. From each cluster, one representative strain was chosen based on the following criteria: a) priority was given to the location where most isolates originated[16] (Brugg > Hoengg > Tuebingen > Seebach > ETH), b) for each location, the leaf from which most isolates originated was given priority[16], and c) the strain's ability to colonise the phyllosphere should be ≥ 105 CFUs / g fresh weight. However, if the colonisation value was ≤ 105 CFUs / g fresh weight, the strain with the highest colonisation value was chosen regardless of location or isolation leaf. Strains that were suspected to be opportunistic pathogens based on pre-test experiments in an agar-based growth system were excluded (two Serratia strains, one Erwinia, one Acinetobacter). Additionally, a Deinococcus strain, which was the only representative in its phylum, was also excluded. See Supplementary Table 2 for detailed information on the selected 62 strains.

Proteobacteria class drop-out and late introduction

All 62 strains were streaked onto R-2A agar (Sigma-Aldrich) supplemented with 0.5 % v/v methanol and allowed to grow for six days at 22 °C. A sterile, 1 µL plastic loop was used to scratch off a "loop full" of each strain and re-suspend it in 1 mL of 10 mM MgCl2. Tubes containing the re-suspended strains were then vortexed for 10 minutes, and strains that formed aggregates were filtered through a 10 µm filter (CellTrics®, Sysmex Suisse AG, Switzerland) prior to mixing. Strains were mixed into their respective phylogenetic groups (Alphaproteobacteria, Betaproteobacteria, Gammaproteobacteria, and "others") in equal ratio (v/v; 750 µL each). Drop-out inocula were then prepared by mixing these groups according to the experimental design and adjusting the total volume so that all mixes added up to the same total (6.2 mL). For example, for the "No Alphaproteobacteria" inoculum, 0.8 mL of Betaproteobacteria (8 strains), 0.3 mL Gammaproteobacteria (3 strains), 3 mL "others" (30 strains), and 2.1 mL of 10 mM MgCl2 (to account for the volume of the 21 missing Alphaproteobacteria strains) were added. Each mixture was adjusted to an OD600 of 0.007 and one week old seedlings were inoculated via pipetting (1 mL per seedling). Because of how these mixtures were prepared, the inoculum density of individual strains in the class drop-outs was higher (maximum by a factor of two) than that in the "all" community control. This change was considered tolerable, i.e. not expected to cause observable effects to be linked to a difference in bacterial density, as it has previously been shown that the synthetic community established with the full At-LSPHERE collection was robust against an even stronger imbalance in the inoculum (factor of four)[16]. Axenic control plants were mock inoculated with 10 mM MgCl2.

Phyllosphere samples were harvested (n = 15 to n = 18 per condition) three weeks post-inoculation to examine microbial community composition (see harvest procedure below). Half of these harvested plants were harvested as-is, while the other half were sprayed with 10 mM MgCl2 immediately before harvesting to assess the effect of the spray on the community (mock spray control). Other plants (n = 15 to n = 16 per condition) were sprayed with the missing group of strains (e.g. adding back Alphaproteobacteria to the condition where they were originally lacking) or mock sprayed with 10 mM MgCl2 (n = 15 to n = 24 per condition) and allowed to grow for another two weeks.

For the late introduction of missing strains, all 62 strains were grown on R-2A plates as described above. Each strain in the Proteobacteria phylum and the "others" (not Proteobacteria) were mixed equally (v/v; 750 µL each) into its respective group (Alphaproteobacteria, Betaproteobaceria, Gammaproteobacteria, and "others"). For each group, a proportionate volume (per number of strains in the mixture) was added to a 50 mL Falcon tube and diluted with 10 mM MgCl2 to the same final volume as the "all" combination (e.g., if adding 1 mL of each strain, the "all" total would be 62 mL; in the Alphaproteobacteria mixture, one would add 21 mL Alphaproteobacteria for 21 Alphaproteobacteria strains and 41 mL of 10 mM of MgCl2 to dilute each strain to the same concentration as in the "all" community). For the Proteobacteria mix, the Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria groups were mixed proportionately and then adjusted to the final "all" strains total volume. OD600 of the "all" mixture was measured and diluted to correspond to an OD600 of 0.02 in 50 mL. The same dilution factor was then used to dilute each group mixture (Alphaproteobacteria, Betaproteobacteria, Gammaproteobacteria, and all Proteobacteria) containing a smaller subset of strains. Each microbox was sprayed six times with a chromatographic TLC reagent sprayer (Faust Laborbedarf AG, Schaffhausen, Switzerland). Axenic plants or control plants without the late introduction strains were mock sprayed with 10 mM MgCl2.

Single strain drop-out

In total, 25 strains were chosen to be individually dropped-out. Of these, 10 out of the 15 most abundant strains (based on the "all" control community; Fig. 2) were randomly selected in addition to 15 out of the remaining 47 strains. Strains and boxes were prepared as described above. 26 mixtures were prepared— an "all" control and 25 mixtures each lacking one of the following strains: L33, L61, L68, L83, L86, L88, L90, L145, L160, L203, L231, L233, L262, L265, L285, L288, L289, L334, L344, L359, L371, L391, L405, L416, and L420 (see Supplementary Table 2 for strain details and phylogeny). A core mix of 37 strains that was common to all mixtures was pre-mixed (1:1 v/v; 750 µL each). For the "all" mix, 740 µL (37 strains x 20 µL each) of the core strains were added, together with 20 µL of each of the remaining 25 strains. For the drop-out mixtures, 740 µL of the core mixture were added together with 20 µL of 24 other strains and 20 µL of 10 mM MgCl2 to compensate for the volume of one missing strain. The total volume was 1240 µL. The OD600 of the "all" 62-strain mix was measured and diluted to correspond to OD600 0.007 in 20 mL of 10 mM MgCl2 solution. The same dilution factor was used for the other mixtures. Inoculation was done as described above.

Plant sampling for bacterial gDNA

For the Proteobacteria class drop-out/late introduction experiment, plants (n = 15 to n = 24 per condition; see above) were harvested at t1 = four weeks and t2 = six weeks old, corresponding to three weeks post-inoculation and two weeks after late introduction. For the single strain drop-out, plants (n = 11 to n = 12 per condition) were harvested only at six weeks old in each of two independent replicates. In one of the independent replicates, one strain (33-Sphingomonas) was accidentally omitted; however, as this strain was often undetected or had a very low relative abundance, its absence did not lead to significant differences between the replicates (Supplementary Fig. 12).

The phyllosphere of each plant was harvested using sterilised tweezers and scalpels. Roots, cotyledons, and clay particles were carefully removed. Each harvested plant was then transferred directly into a lysing matrix E tube (FastDNA SPIN Kit for Soil, MP Biomedicals), frozen in liquid nitrogen, and stored at - 80 °C until further processing.

MiSeq library preparation and sequencing

Samples were freeze-dried (Alpha 2-4 LD Plus, Christ, Osterode am Harz, Germany) overnight. Lyophilised samples were homogenised (TissueLyser II, Qiagen) twice for 1.5 minutes at 25 Hz (rotating the plates once in between) and DNA was extracted using the FastDNA SPIN Kit for soil (MP Biomedicals) following the manufacturer's protocol. 16S rDNA amplicon library preparation was done as previously described[16] and sequencing was carried out at the Functional Genomics Centre Zurich using a MiSeq reagent kit v3 (PE 2 x 300 bp).

Generation of 16S rDNA-based community compositional data

Reference sequences were extracted from draft genomes using USEARCH v8.0.1623_win32 (Edgar 2010) by performing an in silico PCR (-search_pcr) with the 799F[61] and 1193R[44,62] primers and allowing for a maximum of two mismatches, a minimum/maximum amplicon length of 100 and 600 bp, respectively, and searching on both DNA strands. Sequences with ambiguous base pairs were then corrected using the available Sanger sequences [16].

Paired-end sequencing reads were merged using the USEARCH v10.0.240_i86linux64 [64] command fastq_mergepairs 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 0.1 and a minimum length of 100 bp. The reads were then classified and counted using -otutab with a minimum identity of 97 % to the 16S reference sequences, producing an initial OTU table with a count for each reference in each sample. Additional bash code totalled the number of unclassified sequences in each sample to add to the end of the OTU table.

If an experimental sample had been contaminated by a bacterium outside of the 62 strain community, we would expect its 16S sequence to also have been amplified and sequenced. To check for such a case and correctly quantify the relative abundances of the references, the unclassified reads were dereplicated (-fastx_uniques) and clustered using the UPARSE algorithm (-cluster_otus [65] with a minimum cluster size of 1 and a fixed identity threshold of 97 %). The representative sequences from these de novo OTUs were then added to the 16S reference sequences before all quality-filtered reads were re-classified and counted using -otutab with a minimum identity of 97 %. This final OTU table was checked against the initial OTU table and since there was no substantial change in the counts of the reference sequences nor the presence of high abundance de novo OTUs, we concluded that there had been no significant contamination.

Phylogenetics

The phylogeny of the strains was constructed based on their full length 16S rRNA gene sequences (see Supplementary Table 2 for accession numbers), which were aligned using SINA 1.3.3[66] and the SILVA SSU Ref NR database release 132, December 2017. PhyML 3.3.20180214[67] was then used with default parameters to build a maximum likelihood phylogeny from the alignment.

Data Analysis

For each experiment, the final OTU table was imported and processed in R 3.3.3. The data was organised into compatible data sets, each consisting of the samples from a single condition and the control samples. The OTU table was then log-normalised for sequence depth and variance-stabilised by DESeq2 1.14.1[68]. For the class drop-out experiment, samples were compared by treatment (“Effect of Mock Spray in Controls”, “Late Arrival vs. Mock (II)”), time point (“Effect of Time in Controls”), or initial condition (“Group Absent vs. Control (I)”, “Late Arrival vs. Control (III)”). For the single drop-out experiment, due to the noticeable batch effect between replicates one and two, replicate was added to the design, which otherwise compared samples by initial condition. All strains plus the sum of the unclassified OTUs were included during the normalisation procedure, and then later excluded from visualisation and results reports where they had been deliberately absent.

For the effects on individual strains, the output of DESeq2 provided log2 fold-change values and p-values based on Wald tests, adjusted by the Benjamini-Hochberg procedure. DESeq2 was designed for analysis of differential gene expression and assumes that the distribution of expression of a single gene can be fitted to a negative binomial distribution, and that most genes are not differentially expressed. For our data, these assumptions become that the distribution of relative abundance of a single strain can be fitted to a negative binomial distribution, and that most strains do not show differential abundance, which we consider reasonable.

For the visualisation of the overall effects on communities, principle component analysis (PCA) was applied to the transformed OTU table using the prcomp command. The effect size, which is the variance explained by the compared factor, and p-value of the comparison across a specific factor was calculated by PERMANOVA using the adonis function of the package vegan 2.5-4 (https://cran.r-project.org/web/packages/vegan/) with Euclidean distance. For the single drop-out experiment, the PERMANOVA was modified to account for the batch effect between replicates with the strata argument.

Data Visualisation

For the figures that summarise the relative abundance of each strain or the change in relative abundance in a comparison between samples, relative abundance values were calculated by proportionally normalising each sample by its sequencing depth, so that they would be more accessible to the reader. For the figures that correlated relative abundance against other factors, the median value of the DESeq2-transformed OTU table was used instead. For the pie chart of relative abundance, counts were summed across control samples.

The following R packages were also used during analysis and visualisation: corrplot, plotrix, circlize, ape, calibrate, igraph, bezier, beeswarm. Cytoscape (3.7.0) was used to illustrate the causal network.