Contrasting patterns of microbial dominance in the Arabidopsis thaliana phyllosphere

We first analyzed a previous microbiome dataset from wild A. thaliana rosettes, for which we had used plant DNA as a scaling factor to calculate the bacterial load of each amplicon sequence variant (ASV) in the V4 hypervariable region of the 16S rRNA gene (rDNA) (10). After classifying ASVs at the genus level, Sphingomonas and Pseudomonas were approximately equal in bacterial load overall, and twice as abundant as the third most common genus (Pedobacter). However, while Pseudomonas loads varied widely between samples, Sphingomonas was among the most consistent genera (Fig. 1A). Impressively, the consistency of Sphingomonas colonization also extended to individual Sphingomonas ASVs that independently colonized most plants (Fig. 1B). We named the most abundant Sphingomonas and Pseudomonas ASVs “SphASV1-V4” and “PseASV1-V4”, respectively.

We next reanalyzed slightly longer rDNA sequences in the V3V4 hypervariable regions from a different set of wild A. thaliana rosettes in ref. 19, for which we had partitioned A. thaliana leaf epiphytes from endophytes. We found a single V3V4 ASV each that matched SphASV1-V4 and PseASV1-V4, and we refer to these longer sequences henceforth simply as SphASV1 and PseASV1. PseASV1 exactly matched the representative sequence previously used to define the abundant Pv-ATUE5 lineage (19). SphASV1 was more abundant than PseASV1 in both the epiphytic and endophytic fractions of the leaf, and both PseASV1 and SphASV1 were well represented in the endophytic compartment (Fig. 1C). SphASV1 was in fact the most abundant single ASV in the dataset. The combination of numerical dominance, endophytic enrichment, and high consistency of colonization drove our interests to characterize genetic and ecological features of this ASV and local Sphingomonas more generally.

An abandoned train station in Eyach, Germany supports a population of hundreds to thousands of A. thaliana plants that are predominantly winter annuals, and has been a source of material for several recent studies of A. thaliana-associated microbiomes (6, 8, 10, 19). We revisited this site and collected A. thaliana and the most common neighboring plants in spring, when most A. thaliana plants were green, mature, and flowering. We also collected common non-A. thaliana plants and bulk soil samples in late summer after all local A. thaliana vegetation had dried up and died, and the next cohort had not yet germinated. We reasoned that by late summer, microbes that had resided on spring A. thaliana must have long since migrated to survive on other plants or habitats (13, 37, 38) or died along with their host, and thus at this later time point, isolates collected from other plant species or soil could not be recent migrants from A. thaliana.

We chose each non-A. thaliana plant species solely based on its abundance at the site and not by any prior preference for certain plant species, and we therefore chose some plants that we were initially unable to identify, despite being able to confidently classify and recognize different individuals by their morphology. For each plant species, we pooled at least one entire leaf from at least six independent individuals per sample, collecting seven independent samples per species. Although this pooling strategy that spanned many plants per sample precluded analysis of variation between individual plants, it maximized our ability to broadly survey the site (SI Appendix, Fig. S1). After bringing leaf samples back to the lab, we surface-sanitized them in 75% ethanol for 45 to 60 s (39) and macerated them in phosphate-buffered saline (PBS) buffer with a mortar and pestle. Soil samples were directly mixed with PBS. We then mixed the fresh lysates with glycerol to make −80 °C freezer stocks (Methods).

We extracted metagenomic DNA directly from aliquots of the frozen lysates and prepared 16S rDNA amplicon libraries spanning the V3 and V4 hypervariable regions (Methods) using peptide nucleic acids (PNAs) to reduce organelle amplification (40). The residual chloroplast sequences in each sample were consistent within each plant pool, and allowed us to narrow down the identity of the unknown plants (SI Appendix, Fig. S2). Unfortunately, we could not obtain sufficient bacterial reads from dandelion and thistle due to a natural mutation in their chloroplast sequences that made our PNAs ineffective (41), leading to an overabundance of chloroplast sequences. To contrast bacterial communities across the sampled plant hosts, we binned the ASVs into bacterial families, and clustered samples by their pairwise Bray–Curtis dissimilarity (Fig. 1D). This revealed substantial similarities in bacterial family membership between groups of samples, with clear plant clustering both by sampling season and also by some plant taxa, consistent with previous publications linking plant genotype and seasonal effects to bacterial community composition (14, 42). The family Pseudomonadaceae, of which 92.9% was the genus Pseudomonas, reached high relative abundances not only in the A. thaliana phyllosphere, but also in some other hosts—particularly the other Brassicaceae Draba verna and Cardamine hirsuta (Fig. 1E). Sphingomonadaceae, of which 80.7% corresponded to the genus Sphingomonas, was ubiquitous and abundant across all plant hosts (Fig. 1E and SI Appendix, Fig. S3). Across all plants, PseASV1 and SphASV1 accounted for a substantial fraction of the reads from Pseudomonas and Sphingomonas (44.7% and 39.4%, respectively) (Fig. 1 F and G).

Pseudomonas isolates from local A. thaliana populations in Germany, the majority of which has the PseASV1 16S rDNA sequence, have been previously characterized (19). To investigate Sphingomonas genome diversity across host species, and likewise to characterize the genomic features associated with highly abundant Sphingomonas groups, we cultured Sphingomonas from frozen plant lysates using both remaining A. thaliana lysates from ref. (19) as well as from A. thaliana plus diverse plant hosts in the present study. We enriched for Sphingomonas using Luria Broth (LB) media supplemented with cycloheximide and streptomycin, an antibiotic to which most Sphingomonas are resistant due to a natural mutation in the rpsL gene (43), and isolated 410 Sphingomonas colonies. Using LB supplemented with cycloheximide and nitrofurantoin as previously described (19), we also isolated an additional 103 Pseudomonads. We generated draft genome assemblies, annotated open reading frames, and extracted 16S rDNA sequences (Methods). We first analyzed the V3 and V4 regions of the 16S rDNA, which allowed us to match the genomes to our existing culture-independent ASVs. Critically, for both Pseudomonas and Sphingomonas, we recovered isolates in relative abundances consistent with those from culture-independent surveys (Fig. 1 F and G), suggesting that culturing did not bias recovery rates and captured broad patterns of diversity on leaves.

Previously we observed low genomic diversity among PseASV1/Pv-ATUE5 strains isolated from A. thaliana (19). To similarly evaluate genomic diversity for an analogous set of Sphingomonas, we selected all the isolates in our collection that both came from A. thaliana and had the SphASV1 16S rDNA sequence (representing 174 with SphASV1 out of 340 total A. thaliana isolates). We compared the genomes with each other using the Mash algorithm (44), which decomposes genomes into k-mers and calculates a distance based on the fraction of shared k-mers. As a comparator, we also included a representative set of  99 diverse PseASV1-associated genomes isolated from A. thaliana, comprising 82 strains from ref. (19) that were at least 0.1% different in Mash distance from all others across the core genome, and an additional group of 17 PseASV1-associated strains from this study. We converted the Mash distances to a similarity score between 0 (least similar) and 100 (identical), which closely corresponds to average nucleotide identity (ANI) (Methods) (45, 46). When PseASV1-associated genomes were compared with each other, all pairs had Mash similarities > 96. However, pairwise comparisons between SphASV1-associated Sphingomonas genomes had Mash similarities as low as 81 and averaging 89, indicating that the V3V4 region of the 16S rDNA sequence was a relatively poor predictor of genome similarity (Fig. 2 A and B). Generally, the longer full 16S rDNA sequence (as opposed to one or two variable regions) provides increased discriminatory power between strains (47), and we therefore clustered the full-length 16S rDNA sequences extracted from the genomes into operational taxonomic units (OTUs) at 99.5% identity. This yielded no additional subgroups for PseASV1-associated strains, but partitioned SphASV1-associated strains into six subgroups, five of which included more than one strain and could be compared (Fig. 2A). Genomes within these five subgroups were more similar, with average intragroup Mash similarities of 90, 91, 91, 96, and 97. All of these values were lower than for Pseudomonas (Fig. 2B). We further explored the ability of the gyrase B (gyrB) gene, a commonly used high-resolution phylogenetic marker, to distinguish strains, and extracted gyrB sequences from the assembled genomes (Fig. 2A). gyrB outperformed 16S rDNA as a strain-specific marker, and each Sphingomonas shared its gyrB sequence with on average 1.6 other SphASV1-associated strains, while each Pseudomonas shared its gyrB sequence with on average 4.3 other PseASV1-associated strains.

We initially suspected that the relatively high overall genome diversity of SphASV1-associated isolates compared with PseASV1-associated isolates might be due to high variation in the accessory genomes—specifically the differential presence of plasmids. After closing 12 of our genomes with long-read sequencing, we detected in them a total of 16 plasmids with up to three per genome and comprising up to 14% of the total genome size (Fig. 2D). To see if plasmids and mobile elements might be responsible for low Sphingomonas genome similarity, we considered only conserved sequences in genomes by first creating a “soft” core genome using panX (48) to identify open reading frames present in at least 70% of genomes, resulting in a set of 274 shared genes. Mash comparisons of core genomes of SphASV1-associated isolates (Fig. 2C) were still more diverse than even whole genomes of PseASV1-associated isolates. Thus, despite the pervasiveness of the SphASV1 sequence in our dataset, one cannot reliably extrapolate what this ASV sequence means at the level of genomic content.

To explore relatedness of SphASV1-associated strains, and how they compare with other Sphingomonas, we calculated maximum-likelihood (ML) “soft” core genome phylogenies (48) from all 410 local Sphingomonas isolates (340 from A. thaliana and 70 from other local plant hosts), along with 70 sequence-related isolates from NCBI RefSeq (49). Gene presence or absence in the accessory genome (Fig. 2E) tracked well with differences in the core genome, as tested by correlating pairwise Jaccard distances calculated on the presence/absence matrix of accessory genes to Mash distances calculated on core genes (R² = 0.86, SI Appendix, Fig. S4). We compiled a short list of Sphingomonas genes that had a high likelihood to improve survival among competing bacteria, or to facilitate interaction with a plant host, and made a custom database of Sphingomonas protein sequences from RefSeq (Dataset S1). We searched for these features in our genomes by aligning the assemblies using BLASTX (50). While nearly all isolates have diagnostic genes for the type 2 secretion system, and genes for the type 4 and type 6 secretion systems are common among some clades including those of SphASVI, only a handful of less abundant strains seems to potentially have a type 3 secretion system (Fig. 2E). Flagellar motility is common. A group of SphASV1 genomes has a full suite of genes for anoxygenic photosynthesis, a fascinating feature that can supplement heterotrophic energy production and likely improves survival in well-illuminated and nutrient-poor conditions such as the phyllosphere (51, –54). All genomes are rife with toxin–antitoxin systems, likely to stabilize plasmids or superintegrons (55). Indeed, we found widespread evidence of plasmids in our draft genomes, with 250 (60%) of our isolates showing signatures of one or more of the 16 plasmids identified in the closed genomes (Fig. 2F).

To increase confidence in the preceding gene searches in our draft genomes, we repeated the analysis comparing the 12 closed genomes to their draft genome counterparts. Both draft and complete versions of each genome showed essentially the same presence/absence patterns (SI Appendix, Fig. S5), demonstrating that draft genomes were sufficient for analysis of gene content at this level of detail.

Besides reaching similarly high abundances in the same leaves (Fig. 1 A and B), both Sphingomonas and Pseudomonas grow on many similar substrates in vitro (35), making it likely they interact in the phyllosphere. Previous work revealed that certain strains of Sphingomonas, in particular Sphingomonas melonis Fr1, can ameliorate symptoms in A. thaliana leaves caused by pathogenic Pseudomonas and Xanthomonas in a gnotobiotic system (24, 35, 56). More recently, a strain of S. melonis was demonstrated to protect rice against the bacterial seedling blight pathogen Burkholderia plantarii (22). We first sought to screen some of our local isolates for potential plant-protective activity against virulent Pv-ATUE5 strains in a gnotobiotic system as we used previously (19).

We germinated an A. thaliana accession endemic to our field site, Ey15-2, on half-strength MS solid media in the presence of each of 19 diverse strains of local Sphingomonas, as well as S. melonis Fr1. Following 10 d of cocultivation, we challenged the plants with Pseudomonas, either the virulent local P. viridiflava strain Pv-ATUE5:p25c2 (19) or the model pathogen P. syringae pv. tomato (Pst) DC3000, and monitored plant health over the next 6 d by measuring green pixels (Methods and Dataset S2). Pst DC3000 slowed plant growth compared with buffer control, while Pv-ATUE5:p25c2 killed plant tissues (Fig. 3 A–C and SI Appendix, Fig. S6). Surprisingly, seedlings germinated in the presence of S. melonis Fr1 were consistently stunted compared with all other plants (Fig. 3B and SI Appendix, Fig. S6 and Discussion S1). However, despite this negative effect on growth, S. melonis Fr1 protected plants from the worst effects of Pv-ATUE5:p25c2, with infected plants not dying but instead growing more slowly, not significantly different from plants treated with the less lethal Pst DC3000 across two replicates (FDR-adjusted Mann–Whitney U test, P > 0.05). Our local Sphingomonas strain SphATUE:S139H133 also protected plants, reducing Pv-ATUE5:p25c2 virulence such that it was no worse than Pst DC3000 across two independent experiments (Fig. 3 A and B and SI Appendix, Fig. S6) (FDR-adjusted Mann–Whitney U test, P > 0.05). Unlike S. melonis Fr1, SphATUE:S139H133 did not stunt growth.

To test whether protective effects might extend to more natural conditions, we grew plants on potting soil, and set up a similar experiment in which 2-wk-old Ey15-2 seedlings were sprayed first with one of four Sphingomonas strains or a boiled Sphingomonas control, and 4 d later were challenged with Pseudomonas sprayed at high concentrations (O.D.600 = 1.0) (57). We also included enhanced disease susceptibility 1 (eds1-1, Ws-0 background) as an infection control because this mutant, defective for numerous defense responses mediated by salicylic acid (SA), is hypersusceptible to Pst DC3000 (58). At 5 d post infection (dpi), we observed classic Pst DC3000 symptoms on most plants including chlorotic leaves, stunted growth, and increased anthocyanin at the apical meristem (57), especially on eds1-1 plants. However Pv-ATUE5:p25c2, which was consistently deadly on agar plates, did not produce any obvious symptoms on soil-grown A. thaliana and did not greatly affect rosette size (Fig. 3D and SI Appendix, Discussions S2 and Fig. S8).

In the protection experiment on Ey15-2 and eds1-12 plants, we did not observe a consistent protective effect of any local Sphingomonas strain against Pst DC3000 symptoms on soil. However, for both Ey15-2 and eds1-1 plants, pre-treatment with S. melonis Fr1 did result in larger plants than pretreatment with boiled Sphingomonas or with buffer (FDR-adjusted Mann–Whitney U test, P < 0.05, SI Appendix, Fig. S7). In contrast to the agar system, S. melonis Fr1-treated Ey15-2 plants grown on soil were not stunted (SI Appendix, Fig. S7). To confirm that Sphingomonas was still present at the end of the experiment, and if pathogen titers had been affected by Sphingomonas, we quantified bacterial communities in the leaves of eight plants per condition at 5 dpi using hamPCR (59). We observed the inoculated Sphingomonas 16S rDNA sequences enriched in each end point sample, as expected (Fig. 3E). Samples pre-treated with S. melonis Fr1 supported less Pst DC3000 proliferation than those pretreated with other Sphingomonas or buffer in both the Ey15-2 and eds1-1 backgrounds, although Fr1 pretreatment only significantly differed from buffer in eds1-1 (Mann–Whitney U test, P < 0.05). However, we found no evidence that other isolates, including the local Sph139H133, which was protective on agar, reduced pathogen titers in plants grown on soil in these strong infection conditions.

Sphingomonas and Pseudomonas are generalists, with no known exclusive hosts of any given strain (13). As a first step toward determining the breadth of hosts in our study area, we compared Sphingomonas and Pseudomonas ASVs across local plant species. While SphASV1 was more abundant overall in the spring collection, it was consistently detectable on most plant species in both seasons. Plant taxa colonized by SphASV1 were also frequently colonized to appreciable levels by multiple other Sphingomonas ASVs (Fig. 4A). In contrast, PseASV1 was relatively more abundant on A. thaliana, C. hirsuta and D. verna – all from the family Brassicaceae – and less abundant, though still easily detectable, on other taxa (Fig. 4B). PseASV1 was especially enriched in the spring, which matches our previous finding of more Pv-ATUE5 isolates from spring vs. winter collected A. thaliana (19).

To test for strain sharing across hosts, we compared PseASV1 and SphASV1 genomes of isolates from A. thaliana, T. officinale, T. repens, C. hirsuta, and I. glandulifera (Fig. 4 C and D). These species were chosen because we had been able to easily isolate both Pseudomonas and Sphingomonas from the same plant pool lysates prepared for each of these species. We compared the genomes’ ANI(45, 46). While 29 out of the 86 SphASV1 isolates from the 2018 harvest shared at least 99.9% ANI with at least one other isolate, with only one exception these highest similarity isolates came from the same pool of plants, implying that there was little evidence for the same clonal Sphingomonas lineage independently colonizing different plant individuals. To exclude potential clones from the same plant, we recalculated all pairwise genome similarities between SphASV1 isolates from different lysates. We compared similarity within A. thaliana isolates (intrahost isolate similarity) and between A. thaliana isolates and those from any other plant species (interhost isolate similarity). A higher intrahost similarity would be evidence of host-specificity, possibly due to unique selective forces within A. thaliana, or easier migration between individuals. We found that the distribution of intrahost ANI values closely matched that of interhost ANI values (Fig. 4E), with the intrahost ANI being only marginally higher (0.9%, Mann–Whitney U test, P < 0.001).

PseASV1 genomes were more similar to each other than SphASV1 genomes (Fig. 2A). Every isolate shared at least 97.5% ANI with at least one isolate from a different plant species, and 47 of the 50 isolates shared at least 99.9% ANI with at least one other isolate (Fig. 4D), although as with SphASV1, these highest similarity isolates tended to come from the same plant pool. However, three groups of highest similarity isolates were shared across different plant pools, different species, and even different seasons (Fig. 4D). This is consistent with previous observations that Pv-ATUE5 strains are common and persistent on diverse A. thaliana populations (19), and demonstrates that other local plant species host these strains as well. As with SphASV1 isolates, the intrahost ANI values followed a similar distribution to interhost ANI values, with intrahost isolate ANI values again being marginally higher (0.4%, Mann–Whitney U test, P < 0.001, Fig. 4F).

From our genome-sequenced isolates, we had found that closely related Sphingomonas and Pseudomonas strains could colonize diverse hosts. To broaden our survey and extend these observations, we adopted a time and cost-efficient approach to enrich each genus in bulk from plant lysates and sequence the enriched pool as a metagenome. We plated glycerol stock from each lysate on either selective Pseudomonas or Sphingomonas medium and grew colonies en masse (Fig. 5A and SI Appendix, Fig. S9). As a control for lysate viability and a rough reference of bacterial load, we also cultured bacteria from lysates on nonselective rich LB medium. We counted the colonies on each plate after 2 d at room temperature for Pseudomonas or mixed bacteria on LB, or after 7 d for slower growing Sphingomonas (Fig. 5B) (Methods).

To harvest the bacteria, we scraped all colonies from each plate and prepared whole metagenome shotgun libraries from the pools (Fig. 5A). After sequencing the metagenomes, we mapped the reads to a comprehensive reference database including all of our local genomes plus selected publicly available genomes (49, 60) (Methods and SI Appendix, Figs. S10–S12). Because some colonies did not belong to our targeted genera, we also included in our reference database “decoy” genomes of common plant-associated bacteria to capture reads from these “contaminant” strains (SI Appendix, Fig. S11). A total of 165 metagenomes (91 Sphingomonas + 74 Pseudomonas) passed our quality thresholds. For both bacterial genera, there was a clear shift in dominant strains across seasons, including for C. hirsuta and T. repens, which were alive and were sampled at both timepoints (Fig. 5 C and D). This seasonal shift was also apparent in the Sphingomonas and Pseudomonas amplicon data (Fig. 4 A and B). Finally, many metagenome reads from diverse plant samples mapped to the same strains, in agreement with them being widely shared across plants. In particular, this could be observed for strains associated with the most abundant 16S rDNA sequences, SphASV1 and PseASV1. The relative representation of both SphASV1 and PseASV1 isolates in the summer soil samples was low compared with that in plants, suggesting that as the new cohort of A. thaliana germinates in the following fall, it may be more likely that A. thaliana seedlings acquire these strains from nearby plants rather than from the soil.

We consistently noticed a slower growth rate for Sphingomonas than for Pseudomonas, but SphASV1 seemed to be as successful as PseASV1 in establishing substantial population sizes in leaves of wild plants. We hypothesized that the LB medium might advantage Pseudomonas, and growth of the bacteria in the more complex nutrient and chemical milieu supplied by plants might result in more equal performance between the genera. To test this, we collected wild leaves of both A. thaliana and Brassica napus from a site in Kusterdingen, Germany in May 2021. For both A. thaliana and B. napus leaves, PseASV1 and SphASV1 were present and highly abundant in all starting leaf material. We macerated a subset of the leaves with a mortar and pestle, and then compared the growth of bacteria over the next 2 d within macerated leaves and within detached but unmacerated control leaves (Fig. 6), with all samples kept in a plant growth chamber. At each timepoint, we quantified bacterial relative abundances with 16S rDNA amplicon sequencing to detect all genera, and we cultivated and counted Pseudomonas CFUs to estimate changes in absolute bacterial abundances.

PseASV1 proliferated to represent a strong majority of bacteria in most A. thaliana and B. napus leaf macerates (Fig. 6 A and C). While PseASV1 also increased its relative abundance in detached whole leaves of both plant species, the magnitude was much less pronounced than in macerates. We also plated liquid collected from A. thaliana leaf macerate onto solid LB agar and incubated it in the same conditions and examined the bacterial population after 2 d (Fig. 6B). As in macerated plant tissue, PseASV1 strains on LB markedly increased their abundances in the community after 2 d, while the relative abundance of SphASV1 strains and indeed all Sphingomonadaceae decreased to essentially zero (Fig. 6B). These results strongly suggest that PseASV1 thrives in dead or dying leaves and furthermore on rich nonliving substrates, with the intriguing corollary being that Sphingomonas depends on healthy leaves to maintain its competitive edge.

D.S.L., R.d.P.J., and P.P.N.A. designed research; D.S.L., R.d.P.J., P.P.N.A., T.L.K., O.S., K.P., A.B.-G., and D.W. performed research; D.S.L., W.D., and I.B. contributed new reagents/analytic tools; D.S.L., R.d.P.J., and P.P.N.A. analyzed data; and D.S.L., R.d.P.J., P.P.N.A., and D.W. wrote the paper.

The authors declare competing interest. The authors have organizational affiliations to disclose, D.W. consults for KWS SE, a plant breeder and seed producer. The authors have stock ownership to disclose, D.W. holds equity in Computomics, which advises plant breeders.