Introduction

Plant leaf surfaces, commonly termed the phyllosphere, harbor a wide diversity of microorganisms [1]. These endophytic and epiphytic communities can influence plant health and fitness through a variety of means, including protection against pathogens [2, 3], plant growth promotion [4], primary productivity enhancement [5], protection against abiotic conditions including frost [6], and fixation of atmospheric nitrogen (N) [7]. Plant hosts can exert some control over the abundance and composition of their microbiome members by virtue of the differing chemical and physical features of resources provided on their surfaces [8], but also through immune activity, molecular signaling, and barrier formation [9,10,11,12,13,14,15]. This filtering effect can give rise to predictable differences in microbiome composition among hosts [10, 16,17,18], a phenomenon referred to as species identity (or genotype) effects. Evidence for such effects comes from phylogenetic clustering of associated microbial taxa [19, 20], deviation from null or neutral expectations [21, 22], changes consequent to host genetic manipulation [10], or compositional differences explained by species or genotype as a factor [16,17,18, 20, 23]. While species identity effects suggest the importance of host control over microbiota, such effects are often weak or variable when tested in broader environmental or ecological contexts [24, 25]. This raises the question of whether and how host effects can be swamped by environmental factors in shaping the microbiome.

The neighboring plant community constitutes a major component of a plant’s environmental and ecological context. Neighborhood effects, also known as associational effects, have been extensively studied for pathogen and herbivore transmission [26,27,28], revealing patterns of transmission that depend on the proximity to the nearest conspecific neighbor (i.e., conspecific negative density dependence) [29,30,31,32] as well as species frequency-dependent patterns of host fitness [33,34,35]. Much less effort has focused on the role of neighborhood effects for non-pathogenic plant-associated microorganisms [36]. Given that nearby vegetation, site, and soil have been identified as sources shaping phyllosphere microbial communities [23, 37,38,39,40,41], both neighbor identity and proximity are likely to be important factors shaping epiphytic microbial communities. Moreover, it has been shown both theoretically [42] and empirically [43] that in the presence of high dispersal rates, community members can persist even in the face of strong selection against them (e.g., as a result of plant filtering effects), a phenomenon termed mass effects. As such, differences in microbiome composition that arise between species could be diminished when inter-host dispersal is high. Indeed this has been shown in zebrafish, where differences in bacterial community composition among host variants were dramatically reduced when inter-host dispersal was allowed [44].

Recent observational research in tree communities has revealed detectable neighborhood effects on epiphytic communities [38], but many open questions remain. For instance, it is unclear whether neighborhood effects are general and causative, a crucial gap in knowledge if such effects are to be incorporated into agricultural practice. It is also unclear what role neighbor or focal plant identity and age/biomass play in microbiome assembly. We address these knowledge gaps using field- and greenhouse-based experiments involving tomato, pepper, and bean host plants. By manipulating both the focal plant species identity and the local neighborhood composition in the field, we were able to directly test the relative importance of plant host filtering versus local dispersal sources in shaping microbiome composition on leaves. We then performed a controlled cross-inoculation study in the greenhouse to directly examine the effects of host filtering and inoculum source on microbiome assembly of the three plant species involved. We hypothesized that: (1) neighborhood effects would increase as neighboring plants increase in age and biomass; (2) that neighborhood effects would depend on both neighbor and focal plant identity due to host filtering effects; (3) that host species effects would be diminished in the presence of neighboring plants; and (4) that experimental transplantation of microbiomes across hosts would result in compositional change as a result of host filtering, but that there would remain a detectable signal of the past host.

Methods

Experimental design: neighborhood study

To test for the relative influence of host species and neighborhood effects on foliar microbial communities, we implemented a fully factorial, randomized block design at the Oxford Tract, a research farm near the University of California, Berkeley. The study included three plant species: tomato (Solanum lycopersicum var. Moneymaker), pepper (Capsicum anuum var. Early Cal Wonder), and bean (Phaseolus vulgaris var. Bush Blue Lake 274). Plant neighborhoods were established in which a single 5-week-old tomato or pepper plant, or a cluster of 6, 20-day-old beans, was planted as the focal individual in the middle of a circle of eight neighborhood plants, each planted 0.61 m from the focal plant, with fully reciprocal combinations of focal and neighborhood plants (Fig. 1A) in a randomized complete block design with 6 replicates. Neighborhood plots were established in 3 zones (2 blocks per zone) spaced 1.22 m apart separated by at least 0.91 m from any other plants at the experimental site to minimize edge effects. Each neighborhood plot was 1.22 × 1.22 m, and separated by 1.22 m from adjacent plots. No neighbor control plots, in which focal plants had no neighbors encircling them, were included for each species. Thus, the experiment contained nine neighborhood comparisons and three no neighbor comparisons. Two weeks prior to planting, the soil was tilled for weed management and drip lines were installed underneath plastic sheeting to provide irrigation. The plastic sheeting prevented the growth of weeds and minimized the dispersal of soil and irrigation water onto plants. Plants were planted through small holes made through plastic sheeting. Individual tomato and pepper plants were grown to a height of about 20 cm before transplantation into the field, while beans were grown to a height of about 20 cm before transplantation. All greenhouse plants were watered using drip irrigation to minimize the wetting of the leaves, and thus the development of large epiphytic bacterial community sizes. All plants, including focal individuals and neighbors, were transplanted in the field on June 1, 2019.

Fig. 1: Experimental design of the field trial.
figure 1

A Experimental neighborhoods were constructed by planting a focal plant in the center of a ring of neighbor plants (or none), with fully factorial combinations of focal and neighbor plant species. Each block contained 12 comparisons, with focal and neighbor plant abbreviated (T = tomato, P = pepper, B = bean, and N = no neighbor), respectively. Focal plants were harvested and replaced each month, while the neighborhoods were left to continue growing. Shown are a tomato focal plant with tomato neighbors (top) and a pepper focal plant with bean neighbors (bottom) at the time of initial planting. B A pepper neighborhood surrounding a bean focal plant at the time of harvest number 3. C The biomass of neighborhood plants (g) increased to varying degrees with time for each plant species by harvests 1–3.

Full size image

Focal plants were harvested and replaced at 30-day intervals (3 successive plantings) while the neighborhoods were retained and continued to grow throughout the study (Fig. 1B, C). Focal plants were all the same age upon planting as the original cohorts to allow for direct comparisons across cohorts. Thus, the bacterial community composition of the focal plants was assessed on 3 separate occasions for each of 72 focal points in the neighborhoods (totaling 216 focal samples). Only one focal plant (pepper with no neighbors) died prior to sampling. Additionally, at each round of planting, bacterial community composition, but not abundance, was assessed on five plants of each species at the time of transplantation to identify taxa that had been established on plants in the greenhouse.

Neighborhood plant attribute measurements

Several attributes of neighborhood plants were measured before each monthly harvest of focal plants in order to determine how neighbors might impact phyllosphere communities of the focal plant. These attributes included average neighbor height, the distance of the focal plants to the nearest neighbor, the number of neighbors touching the focal plant (if any), the total number of flowers on the neighborhood plants, and whether the neighborhood plants had signs of herbivory, infestation, or disease (yes or no). Further, the biomass of the neighborhood was estimated without harvesting the plants by fitting a linear model relating the height and weight of the focal plants and extrapolating to that of the neighbor plant weights based on their height. Separate linear models were fit for each plant species (Supplementary Table 1), at the following ages: tomatoes were roughly 9 weeks old (flowering and early fruiting stage), peppers were roughly 9 weeks old (flowering and early fruiting stage), and beans were roughly 7 weeks old (flowering and fruiting stage).

Sample processing

Immediately before harvesting the focal plants, their height was measured. Focal plants were then excised at their base using ethanol-sterilized scissors, transferred to 1-gallon (3.79 L) sterile plastic bags, and transported in a chilled cooler to the laboratory. Plant weight was recorded and then plants were subsampled, collected, and re-weighed. Subsampling was necessary to reduce biomass differences among samples and to enable a more efficient collection of epiphytic bacteria by sonication. Foliar bacteria were collected from plant subsamples (range: 3.84–621.14 g, median 40.84 g) by adding 180 ml of sterile 10 mM MgCl2 to the sample bags and sonicating for 10 min in a sonicating water bath (Branson model 5800). Leaf wash was then filtered through an autoclaved coffee filter and distributed across four 50 ml conical tubes, which were then centrifuged at 4000 rcf at 10 °C for 10 min to pellet microbial cells. The supernatant was then decanted from each tube and the pellets were resuspended in 1.8 ml King’s broth (KB) and combined into a single tube. 600 μl of this resuspended pellet was frozen at −80 °C for subsequent DNA extraction, while the remaining two 600 μl aliquots were each mixed with 400 ul 1:1 KB:Glycerol and frozen at −80 °C for subsequent experimentation.

Experimental design—follow-up transplant study

To further test the importance of host filtering, inoculum source, and dispersal history, we conducted a follow-up greenhouse study in which bacterial communities recovered from field plants at harvest time point 2 were inoculated onto new plants of these same species under controlled conditions. Cryopreserved phyllosphere communities from a single focal field plant were transferred to either the same plant species from which they were isolated or onto the plant species that had previously surrounded that focal plant when it was in the field. For instance, the microbiome from a tomato that was surrounded by beans was applied equally to new tomato and bean plants. This was done for all combinations. Experimental blocks from the field trial were treated as experimental blocks in the greenhouse trial, using blocks 2–6 from the field (5 replicates per treatment). We deliberately did not equalize inoculum density, as we anticipated that bacterial abundances would vary according to plant species and thereby constitute an important component of species identity effects. The biomass of every donor plant, however, was recorded for downstream statistical analysis. Additionally, for each plant species, we included 5 replicate blank inoculum controls in which the same volume of sterile 10 mM MgCl2 that was used to resuspend inoculum was sprayed onto plants. We further included replicate heat-killed controls, in which field-derived leaf wash was autoclaved for 40 min before being applied to each of three plant hosts, in the same manner as the experimental inocula.

Inocula were prepared by thawing the freezer stock, centrifuging at 4000 rcf and 10 °C for 10 min to pellet cells, decanting the supernatant, and re-suspending cells in 7 ml 10 mM MgCl2 then splitting in half to make two 3.5 mL inocula. Twenty-two samples were inoculated each day (block) such that each block contained every comparison, and this was repeated for 5 days. Inocula were sprayed onto the adaxial (top) and abaxial (bottom) sides of leaves using ethanol- and UV-sterilized misting caps. After inoculation, the moist, sprayed plants were placed in a chamber maintaining ca. 100% relative humidity for 20 h in order to maintain leaf moistness, thus encouraging microbial growth, before being transferred to a greenhouse where their placement was randomized across 3 benches. After 7 days, the plants were returned to the humid chamber for 20 h immediately before harvest in order to facilitate further microbial multiplication on leaves and thus allow for maximal host filtering via differential growth of phyllosphere members. Plants were then harvested and processed as in the field study.

DNA extraction, PCR, library preparation, and sequencing

One-sixth of the total leaf surface microbial extraction per plant was used for DNA extraction with DNeasy Powersoil Kits (Qiagen). The sample order was randomized to avoid batch effects, and a blank (no sample) control was included in every round of DNA extraction. DNA concentration of each sample was quantified using the Qubit dsDNA HS Assay Kit. Sample DNA (10 μl) was used as template and polymerase chain reaction (PCR) amplified for 35 cycles at the University of California—Davis Host–Microbe Systems Biology Core using the 799F (5′-AACMGGATTAGATACCCKG-3′)–1193R (5′-ACGTCATCCCCACCTTCC-3′) primer combination, which targets the V5–V7 region of the 16S rRNA gene, and was designed to minimize chloroplast amplification [45, 46]. To further minimize host mitochondrial and chloroplast amplification, peptide nucleic acid clamps were added to each reaction [47]. The resulting amplicons were diluted 8:1 and were further amplified for 9 cycles to add sample-specific barcodes, then quantified using Qubit, pooled in equal amounts, cleaned with magnetic beads, and size selected via electrophoresis on a Pippin Prep gel (Sage Science, USA). The resultant library was then sequenced on the MiSeq (paired-end 300) platform (Illumina, USA).

Sequence processing

Amplicon sequences were processed using the DADA2 pipeline [48] implemented in the R statistical environment [49], including the packages ShortRead [50], Biostrings [51], and Phyloseq [52]. Forward and reverse reads were truncated at 260 and 160 bp, respectively, and quality filtered using the function ‘filterAndTrim’ with default settings (i.e., maxN = 0, maxEE = c(2, 2), and truncQ = 2). Error rates for forward and reverse reads were determined using the ‘learnErrors’ function, and then applied to remove sequencing errors from reads and assign them to amplicon sequence variants (ASVs) using the ‘dada’ function. Paired reads were merged, converted into a sequence table, and then chimeric sequences were removed from the sequence table. Taxonomy was assigned to the remaining ASVs using the ‘assignTaxonomy’ function, which implements the RDP Naïve Bayesian Classifier algorithm with kmer size 8 and 100 bootstrap replicates [53]. This taxonomic classification used the Silva (version 138) SSU taxonomic training dataset formatted for DADA2 [54]. Chloroplast and mitochondrial sequences were filtered from the ASV table by removing any ASVs with a taxonomic assignment of ‘Chloroplast’ at the Order level or ‘Mitochondria’ at the Family level, respectively. Lastly, we applied the ‘isContaminant’ function (method = prevalence) from the package ‘decontam’ [55] to our samples using our blank (no sample) DNA extractions to identify and remove putative contaminants introduced during DNA extraction.

Bacterial quantification using droplet digital PCR (ddPCR)

In order to estimate foliar bacterial abundances of each plant sample from the field, we used droplet digital PCR (ddPCR) using the Bio-Rad QX200 system (Bio-Rad, USA) on the same bacterial DNA extracted from leaf wash described above. Comprehensive ddPCR methods are described elsewhere [56], but   briefly, we targeted the V5-V7 region of the 16S rRNA gene in sample DNA using the chloroplast-excluding 799 F (5′-AACMGGATTAGATACCCKG-3′)–1389R (5′-ACGGGCGGTGTGTRC-3′) primer combination. Five microlitres of 1:10 diluted DNA template were combined with 11 μl of 2× EvaGreen Supermix (Bio-Rad, USA) and 0.22 μl of each primer, and 5.56 μl of molecular grade water to a total volume of 20 μl. Reaction mixes were then loaded into the QX200 droplet generator with 70 μl of droplet generation oil and then transferred to a PCR plate. Thirty-nine cycles of PCR were performed under the following conditions: 95 °C for 10 min, 95 °C for 30 s, 55 °C for 30 s, 72 °C for 2 min, with steps 2–4 repeated 39 times, 4 °C for 5 min, and 90 °C for 5 min. EvaGreen signal was measured on the QX200 droplet reader, cutoff thresholds were set for each column based on background fluorescence in no-template controls, and concentrations were determined using the associated QuantaSoft software. Abundances are reported as 16S rRNA gene copies per g plant material as well as estimates of 16S rRNA gene copies per individual plant by taking into account the proportion of the total plant sample that was used for sample processing.

Statistical analysis

All statistical analyses were performed using R version 4.0.3 [49]. Community matrices were rarefied to 6400 counts per sample ten times and averaged in order to account for differences in sampling extent across samples. The same rarefying procedure was performed for the greenhouse study samples, but to a depth of 10800 counts per sample. ASV rarefaction curves were generated using the ‘rarecurve’ function in the vegan package in R [57]. Bray Curtis bacterial community dissimilarities were calculated between samples using the ‘vegdist’ function. Community structure differences among host species identity, neighbor species identity and experimental block were assessed using a PERMANOVA on Bray-Curtis distances using the ‘adonis’ function (also in the vegan package), which performs a sequential test of terms and uses the algorithm presented in [58]. The order of terms in the model for the field trial was: host, neighbor, and block, and the order of terms for the greenhouse experiment was host, previous host, previous neighbor, donor plant biomass, and block, although this did not impact our qualitative conclusions. To assess the change in the relative strength of these factors through time, the PERMANOVA was performed for each of the three harvesting time points separately. Since not all samples were successfully sequenced, generating slight differences in sample numbers among harvests, we adjusted R2 values to take into account sample numbers and degrees of freedom using the ‘RsquareAdj’ function in the vegan package. In instances of hypothesis testing on subsetted data, p values were adjusted for multiple comparisons using the Benjamini–Hochberg procedure with the ‘p.adj’ function (method = ‘hochberg’). Indicator taxa analysis was performed using the ‘multipatt’ function in the indicspecies package [59].

In order to assess the unique contribution of plant species identity for each neighborhood type, we used variation partitioning on Hellinger-transformed community matrices to partition out the effects of space. A geographic distance matrix was calculated for all experimental plots based on plot GPS coordinates using the program Geographic Distance Matrix Generator [60], and then pairwise distances were decomposed into principal coordinates using the ‘pcnm’ function in the vegan package. For each combination of plants, significant principal coordinates were selected using forward and backward (direction = “both”) model selection with the ‘ordistep’ function in the vegan package. The unique contribution of host species identity was then calculated after partitioning out the significant spatial PCNMs that were selected using the ‘varpart’ function in vegan. We further assessed the contribution of other host factors including plant height and weight by performing the above-mentioned model selection and variation partitioning. The statistical significance of various fractions was then tested by performing redundancy analysis ordination and declaring the non-focal factors as conditions.

In order to estimate the impact of neutral processes (including chance and dispersal) on microbiome assembly, we performed neutral modeling of phyllosphere communities using the ‘fit_sncm’ function in the package reltools [61]. This package fits the neutral model from [62], as implemented by [21]. In order to assess phylogenetic patterns in the phyllosphere communities, we constructed a phylogenetic tree of all ASVs with greater than 20 counts in the community matrix, which included 7949 ASVs. Sequences were aligned using the ‘AlignSeqs’ function in the DECIPHER package [63] using default settings. Next, pairwise distances between sequences were calculated using the ‘dist.ml’ function in the phangorn package version 2.5.5 [64]. These distances were then used to construct a neighbor-joining tree using the ‘NJ’ function in phangorn. Lastly, the neighbor-joining tree was used as a starting point to create a generalized time-reversible with gamma rate variation (GTR + G + I) maximum likelihood tree using the ‘pml’, ‘update’, and ‘optim.pml’ functions in the phangorn package. Lastly, we calculated the mean pairwise distance (MPD) of taxa in each sample and compared to the MPD of a null model to calculate the standardized effect size (SES) using the ‘ses.mpd’ function in the picante package [65]. We used the community randomization null model (null.model = species.pool, iterations = 999) whereby a randomized community matrix is constructed by drawing from the total species pool with equal probability. Using such a procedure, a Z-score (the SES of MPD versus the null community) below zero can be interpreted as phylogenetically clustered whereby taxa co-occurring in a sample are more closely related than the same number of taxa drawn at random from the species pool. By contrast, samples above zero are interpreted as phylogenetically overdispersed, i.e., the phylogenetic distance among co-occurring taxa is greater than the above-stated null expectation. It is worth noting that the use of the 16S rRNA gene for such phylogenetic analyses may have limited ability to capture functional and/or genomic differences among taxa, especially traits that have been acquired through horizontal gene transfer.

For univariate data such as ASV-level richness, MPD SES, and ddPCR-based abundance data, a three-way ANOVA was fit to test for significant effects of the host, neighbor, harvest time point, and block, with interactions therein. The appropriateness of this procedure was verified by checking for a normal distribution of residuals on the model. In order to calculate the proportion of variance explained by each variable in the ANOVA model used for bacterial abundance, we divided the sum of squares (SS) of each variable by the total SS (i.e., the sum of the SS of every variable, including the residual SS) and report this as an R2 value.

Results

Experimental manipulation of plant neighborhood in the field

We first compared the phyllosphere microbiome of plants that were surrounded by no neighbors, conspecific (same species) neighbors, or heterospecific (different species) neighbors. Because the field trial was conducted over the course of 3 months, with focal plants being replaced with a plant of the same species but at the original age of planting each month, we were also able to compare neighborhood age/biomass effects on microbiome assembly. After processing, 175 of the 216 focal plant samples from the field yielded high-quality sequencing reads (Tomato n = 63, Pepper n = 52, Bean n = 60). Of these 175, 64 were from harvest 1, 58 were from harvest 2, and 53 were from harvest 3. The 41 excluded samples each had less than 24 total reads and thus failed either to amplify or to yield sequences. The greenhouse control plants that were included to assess the bacterial communities established prior to transplantation into the field yielded very few reads (28 of the 45 samples had less than 50 reads), indicating that bacterial colonization prior to transplantation was minimal. Of the 17 greenhouse control samples with detectable sequences, communities were dominated by members affiliated with the Enterobacterales, Corynebacterales, Burkholderiales, and Pseudomonadales.

The field trial dataset contained 5,414,393 observations of 19,818 ASVs, 13,455 of which had >10 occurrences, and 2253 of which had >100 occurrences. Within-host ASV-level richness ranged from 22 to 769 ASVs across all treatments and hosts (Supplementary Fig. 1A–C). Richness levels varied significantly by harvest time (F2,174 = 24.21, p < 0.001), declining throughout the season, and varying by host identity (F2,174 = 4.96, p < 0.01), with beans harboring a greater richness, especially at the first harvest. The neighborhood did not impact bacterial richness, however, the total number of flowers on the neighborhood plant species at the time of focal plant harvest was positively correlated with bacterial richness on the focal plants (R2adj = 0.044, p = 0.01). Similar qualitative trends were observed for the Shannon index, except that host plant weight and height were also positively correlated with diversity (R2adj = 0.022, p = 0.03 and R2adj = 0.044, p < 0.01, respectively).

Bacterial abundance per plant varied significantly by host identity (F2,174 = 27.62, R2 = 0.196, p < 0.001), neighborhood (F3,174 = 8.04, R2 = 0.086, p < 0.001), and harvest time point (F2,174 = 3.89, R2 = 0.014, p = 0.051, Fig. 2A), with no significant interactions among variables and no block effect. Tomato and bean plants tended to harbor higher bacterial abundances than pepper plants (p < 0.001). As expected, abundance per plant was positively correlated with plant weight (R2adj = 0.187, p < 0.001) and plant height (R2adj = 0.117, p < 0.001). Several neighborhood attributes had associations with bacterial abundance on focal hosts. Specifically, estimated neighborhood biomass (R2adj = 0.146, p < 0.001), average neighbor height (R2adj = 0.123, p < 0.001), and total number of flowers (R2adj = 0.022, p = 0.038) were all negatively correlated with bacterial abundance per plant. Bacterial abundance per focal plant was negatively associated with community richness (R2adj = 0.016, p = 0.05) and the Shannon index (R2adj = 0.03, p < 0.01). Lastly, if we normalize bacterial abundances by the plant material weight used for processing, we see weaker effects of host identity (F2,174 = 3.86, R2 = 0.035, p = 0.023) and neighbor (F3,174 = 2.59, R2 = 0.034, p = 0.054), but a stronger effect of harvest time point (F2,174 = 5.61, R2 = 0.025, p = 0.019, Supplementary Fig. 2) as well as a block effect (F5,174 = 3.87. R2 = 0.025, p < 0.01).

Fig. 2: Bacterial abundance and composition vary across host species and harvest time.
figure 2

A The abundance (log10 16S rRNA gene copies measured using ddPCR of leaf washes, y-axis) for individual focal plant species (x-axis) surrounded by different neighbor plant species (box color) and at different successive harvest times (panels 1–3). B Relative abundance of the nine most abundant bacterial orders distinguished by host species and harvest time. All other less abundant or ambiguously assigned orders are grouped under ‘Other’.

Full size image

Overall, phyllosphere communities were dominated by taxa affiliated with the phyla Proteobacteria, Firmicutes, and Actinobacteriota. The most abundant bacterial orders were the Bacillales, Burkholderiales, Enterobacterales, Lactobacillales, Micrococcales, Rhizobiales, Sphingomonadales, and Xanthomonadales (Fig. 2B). Bray Curtis dissimilarities among samples were driven by harvest time point (R2 = 0.063, p = 0.003), host species (R2 = 0.055, p = 0.003), neighbor (R2 = 0.023, p = 0.048), and block (R2 = 0.036, p = 0.06), with a significant interaction between host and harvest (R2 = 0.034, p = 0.033), and a trending interaction between neighbor and harvest (R2 = 0.041, p = 0.072).

The effects of host identity on bacterial community composition decrease through time while neighborhood effects increase through time and vary by host identity

We next examined the relative influence of host species identity and neighborhood on focal plant microbiome structure at each time point during the field experiment by performing a PERMANOVA on Bray Curtis dissimilarities using host species identity (i.e., tomato, pepper, or bean), neighborhood (i.e., tomato, pepper, bean, or no neighbor), and experimental block (1–6) as independent variables. The effect of host identity was significant, but diminished in size over the three time points (Harvest 1: Adj. R2 = 0.096, p < 0.001; Harvest 2: Adj. R2 = 0.068, p < 0.001; Harvest 3: Adj. R2 = 0.027, p < 0.001, Fig. 3A, see Table 1 for pre-adjusted R2 values). In contrast, the effect of neighborhood status was initially not statistically significant, but increased in size over the three time points (Harvest 1: Adj. R2 = −0.001, p = 0.208; Harvest 2: Adj. R2 = 0.017, p < 0.014; Harvest 3: Adj. R2 = 0.032, p = 0.003, Table 1, Fig. 3A). Block effects were not statistically significant and tended to decrease throughout the experiment (Harvest 1: Adj. R2 = 0.011, p = 0.063, Harvest 2: Adj. R2 = 0.009, p = 0.06, Harvest 3: Adj. R2 = 0.009, p = 0.063, Fig. 3A). No significant interactions among variables were observed at individual harvests.

Fig. 3: The effects of host identity on bacterial community composition decrease through time while neighborhood effects increase and vary by host plant species.
figure 3

A Adjusted R2 values (y-axis) are the result of PERMANOVA analyses on Bray Curtis dissimilarities for each harvest, accounting for sample number and degrees of freedom from slight differences in sample number. See Table 1 for pre-adjusted R2 values. The effect of host identity (solid maroon line), the effect of the neighborhood (blue dot-dash line), and the effect of the experimental block (green dotted line) are shown. The harvest time point is shown on x-axis. Filled circles indicate statistical significance (p < 0.05), while open circles represent statistically insignificant effects (p > 0.05). B Host plant species experience neighborhood effects on phyllosphere bacterial communities differently through time. Adjusted R2 values (y-axis) and harvest time point (x-axis) are as described for plot A. Tomato hosts (solid red line), pepper hosts (light green dot-dash line), and bean hosts (light blue dotted line) are shown. Filled circles indicate statistical significance (p < 0.05), while open circles represent statistically insignificant effects (p > 0.05).

Full size image
Table 1 Results of a PERMANOVA on phyllosphere bacterial community Bray Curtis dissimilarities for harvest time points 1–3 in the field trial.
Full size table

By excluding the ‘no neighbor’ controls, we then tested for an effect of neighbor identity by treating neighbor type (i.e., tomato, bean, or pepper) as an independent variable. On this subset of plants we see similar trends through time: host identity effects diminish (Harvest 1: Adj. R2 = 0.114, p = 0.002; Harvest 2: Adj. R2 = 0.050, p = 0.002; Harvest 3: Adj. R2 = 0.024, p = 0.004) and neighbor identity effects increase (Harvest 1: Adj. R2 = −0.003, p = 0.319; Harvest 2: Adj. R2 = 0.010, p = 0.044; Harvest 3: Adj. R2 = 0.018, p = 0.044). A significant block effect was only observed at time point 3 (Adj. R2 = 0.032, p = 0.03).

We further asked whether a closer approximation of bacterial taxon absolute abundances might impact our conclusions. The relative abundance of each taxon in each sample was multiplied by the total 16S rRNA copies per 10 μl DNA (the same volume used for sequencing library preparation) and then divided by the total g plant material processed for each sample to yield an estimate of each taxon’s absolute abundance (quasi-absolute abundance) per g plant material. In this new dataset, sample dissimilarities were modeled using the same PERMANOVA procedure as above. This generated the same qualitative findings as the relative abundance data, but with slightly stronger effect sizes for the influence of neighborhood plant species (Supplementary Table 2). One new result revealed by this approach, however, was a host by neighborhood interaction, which increased from harvest 2 (Adj. R2 = 0.006, p = 0.06) to harvest 3 (Adj. R2 = 0.022, p = 0.057, see Supplementary Table 2 for pre-adjusted R2 values). In other words, the effect of neighborhood depended on the host’s species identity, and this effect became stronger over time.

To further assess whether the three host plants species differed in their susceptibility to neighborhood effects, we subset the data by plant species and assessed the strength of neighborhood effects separately over the three time points. Similar to the combined data, no plant species exhibited a detectable neighborhood effect of microbiome composition at harvest 1. The tomato focal plants only exhibited a detectable neighborhood effect at harvest 3 (R2adj = 0.086, p = 0.009, Fig. 3B), and no block effects at any harvest. For the pepper focal plants, there was only a neighborhood effect at harvest 2 (R2adj = 0.108, p = 0.048, Fig. 3B). Lastly for the bean focal plants, a neighborhood effect was only detected at harvest 3 (R2adj = 0.069, p = 0.039, Fig. 3B). No significant block effects were observed for pepper or bean at any harvest.

Neighbor status and identity diminish effects of host species identity on phyllosphere bacterial communities

To test the hypothesis that focal plants with experimental neighbors will have less influence from host filtering effects than plants with no immediate neighbor, we assessed the unique contribution of host species identity to focal plant microbiomes for each neighbor type by statistically excluding the effects of space (i.e., the spatial proximity of neighborhoods). We find supporting evidence for our hypothesis in harvests 2 and 3, but not harvest 1 (Fig. 4), indicating a dependency on the age structure of neighbors. For harvests 2 and 3, the unique contribution of host species identity was stronger for the plants having no neighbors than the plants having either tomato, pepper, or bean as neighbors. In fact, at harvest 3 focal plants that were surrounded by bean or tomato neighbors had no detectable effect of host species identity after separating out the effect of space. At harvest 1, focal plants with neighbors had higher effects of host species identity than the no neighbor controls, and this was especially the case for plants with tomato or bean neighbors. For the focal plants surrounded by peppers, host species effects followed a hump-shaped relationship (increasing at first, then decreasing), rather than the monotonic decrease observed for tomato- or bean-surrounded plants. It thus appears that host-mediated selection of epiphytic bacterial communities becomes subordinate to the effect of immigrant inoculum from neighboring plants as the biomass of neighboring plants increases.

Fig. 4: The unique contribution (adjusted R2) of host plant identity to phyllosphere bacterial structure decreases over time for plants with neighbors.
figure 4

Unique contribution was calculated by partitioning out spatial principal coordinates using RDA-based variation partitioning. The order of depiction of the neighbor plant species is by estimated neighborhood biomass from lowest to highest. Boxes 1–3 represent different harvest time points. In cases where Radj = 0, host species identity did not significantly explain variation in phyllosphere bacterial composition.

Full size image

Additionally, in certain cases, host identity combined with host height, weight, or both height and weight in a way that increased the explanatory power of host factors. While in several instances this boosted the explanatory power of host factors (e.g., at harvest 1), our qualitative conclusions remain the same (Supplementary Fig. 3). In other words, the effects of host identity are weaker for all plants that have neighbors at both harvests 2 and 3.

Phylogenetic clustering and neutral model fit vary by host and through time

To better understand the predominance of deterministic processes in shaping phyllosphere community membership and determine whether the three plant species might be influenced by different assembly processes, we tested for patterns of phylogenetic clustering. Evidence of phylogenetic clustering within a host species would suggest that phylogenetically-conserved traits are being selected for in a host-specific way [19, 20]. We tested this idea using the standardized effect size (SES) of the mean pairwise distance (MPD) of bacterial ASVs in each sample. MPD SES was significantly influenced by host species identity (F2,175 = 219.86, p < 0.001), harvest time point (F2,175 = 12.33, p < 0.001), a host by harvest interaction (F4,175 = 43.84, p < 0.001), and neighborhood (F3,175 = 3.38, p < 0.05, Fig. 5A). Tomato- and pepper-associated communities were more phylogenetically clustered than would be expected by chance (as indicated by negative SES values), and this was the case for all three-time points for both plants. In contrast, bean-associated communities showed evidence of phylogenetic overdispersion (as indicated by highly positive SES values). Bean SES values tended to decrease (i.e., tend towards clustering) over time, but were highly variable. One exception was at harvest 3, where beans with no neighbors had high variability in MPD SES, beans with conspecific neighbors were phylogenetically clustered, and beans with tomato neighbors were overdispersed (Tukey’s HSD conspecific neighbor vs. tomato neighbor p = 0.02).

Fig. 5: Bacterial leaf surface community assembly processes differ between plant hosts, suggesting differences in host filtering.
figure 5

A Standardized effect size (SES) of mean pairwise distance (MPD) of phyllosphere communities split by the host (x-axis), neighbor (box color), and harvest time point (panels 1, 2, or 3). SES = (MPDobs − MPDnull)/SD(MPDnull), whereby values below 0 suggest phylogenetic clustering. B The fit of a neutral model declines through time but differs strongly by host identity. Neutral model goodness-of-fit values (R2, y-axis) at each harvest (x-axis) for tomato (solid red line), pepper (green dot-dashed line), and bean (blue dotted line). Filled circles indicate statistical significance, open circles indicate not significant (negative or 0 goodness-of-fit values).

Full size image

We next asked how well the occupancy–abundance relationships within each host species could be fit by a neutral model, whereby passive dispersal and ecological drift are the primary drivers of establishment, and then asked whether the fit to that model changed over time. Of the three hosts, bean-associated communities had the highest goodness-of-fit values followed by tomato-associated communities, suggesting differences among hosts in the role of neutral processes in shaping community structure (Fig. 5B). Both bean and tomato hosts showed a decline in the fit of a neutral model from harvest 1 to 2 (Bean: harvest 1 R2 = 0.483, harvest 2 R2 = 0.074; Tomato: harvest 1 R2 = 0.214, harvest 2 R2 = 0.052), and the neutral model failed to fit either set of plants for harvest 3 (as indicated by a negative goodness-of-fit). At all three time points, pepper hosts were never fit by a neutral model (indicated by negative goodness-of-fit values).

Experimental greenhouse transplantation of field study-derived inocula replicates host filtering and reveals effects of inoculum source

The subsequent greenhouse study allowed us to more closely examine the effects of inter-host dispersal of bacterial taxa, as phyllosphere bacterial communities that were recovered from the field were transplanted onto either the same plant species from which they were collected or onto the plant species that had previously surrounded the source plant. From the 105 total reciprocal inoculations, 103 samples yielded sufficient high-quality sequences for analysis. The resultant dataset contained 2,640,588 observations of 1734 ASVs, 1379 of which had greater than 10 observations, and 621 of which had over 100 observations.

We observed a linear relationship between the ASV-level richness of the sample from which the inoculum was derived and the number of inoculum ASVs that were detectable in the greenhouse samples (R2adj = 0.105, p = 0.004, Supplementary Fig. 4). The number of overlapping ASVs between the inoculum and the experimental plants was significantly related to the host species identity (F2,88 = 8.503, p < 0.001), the previous host species from which the inoculum was sampled (F2,88 = 3.871, p = 0.028), and interactions between the host and previous host (F4,88 = 2.598, p = 0.049) as well as between the host and previous neighbor (F4,88 = 2.968 p = 0.03).

Similar to the field study, phyllosphere communities were dominated by members affiliated with the phyla Proteobacteria and Firmicutes. The most abundant bacterial orders were the Bacillales, Burkholderiales, Enterobacterales, and Pseudomonadales (Supplementary Fig. 5). The bacterial community structure on treated plants differed significantly from that of control plants to which only sterile buffer had been applied (PERMANOVA R2 = 0.0156, p = 0.048).

We used the heat-killed inoculum as a control to gain insights into the level of host selection and subsequent inoculum establishment across our experimental plants. To do so, we asked how strongly phyllosphere microbiomes were differentiated by host species identity. Plants that received a heat-killed inoculum were slightly less differentiable than plants that received a live inoculum (PERMANOVA Adj. R2 = 0.100, p = 0.04 vs. Adj. R2 = 0.103, p < 0.001, respectively). The plants that received the sterile buffer as control were even less differentiable by host species identity (Adj. R2 = 0.056, p = 0.07). Interestingly, if we subset samples based on whether experimental plants received an inoculum from heterospecific (different species) or conspecific (same species) hosts, we see that heterospecific transplants resulted in more differentiable hosts than conspecific transplants (Adj. R2 = 0.123, p < 0.001 vs. Adj. R2 = 0.109, p < 0.001, respectively). It thus appears that the treatment plants receiving live cells more efficiently filtered communities, driving differentiation of host species and that heterospecific inoculum sources further bolstered host differentiation.

Indicator taxon analysis allowed us to examine the taxa enriched on each of the three host species in the field trial (at harvest 2) and the greenhouse trial. Of the 34 taxa that distinguished pepper plants in the field, 6 were found in the greenhouse dataset (Supplementary Table 3). The collective relative abundance of these taxa was significantly higher on pepper plants in the greenhouse than on tomatoes (Tukey’s HSD p < 0.001) or beans (Tukey’s HSD p < 0.001). Of the 65 taxa that distinguished tomato hosts in the field, 25 were detected in the greenhouse dataset (Supplementary Table 3). The collective relative abundance of these taxa was significantly higher on tomato and pepper plants relative to beans (Tukey’s HSD p < 0.05 and p < 0.01, respectively). Of the 300 taxa that distinguished bean, 27 were found in the greenhouse (Supplementary Table 3), but their collective relative abundance was not significantly different among hosts (p > 0.05).

The effect of donor plant biomass on recipient plant phyllosphere richness in the greenhouse depends on the origin of field inoculum

Within-host ASV-level richness of treatment plants ranged from 24 to 231 ASVs and varied significantly by host species identity (F2,74 = 12.41, p < 0.001), an interaction between host and previous host identities (F4,74 = 3.50, p < 0.05), and experimental block (F4,74 = 2.43, p = 0.058). Of the three plant species, peppers harbored significantly higher richness than tomatoes or beans (p < 0.001), which were indistinguishable from each other (p > 0.05). When the inoculation was a conspecific transfer (i.e., moving between two plants of the same species), a negative but weak relationship was observed between donor plant biomass and recipient plant richness (Adj. R2 = 0.09, p < 0.05, Fig. 6A). However, when the inoculation was a heterospecific transfer (i.e., between two different plant species), a positive and stronger relationship was observed between donor plant biomass and recipient plant richness (Adj. R2 = 0.21, p < 0.01, Fig. 6B). No significant differences in richness were observed between conspecific transplants and heterospecific transplants (p > 0.05).

Fig. 6: The effect of donor plant biomass (g) on phyllosphere community richness depends on the origin of the inoculum.
figure 6

A Conspecific (within the same species) transfers, where ASV-level richness (y-axis) is negatively, but weakly, correlated with the donor plant biomass (g, x-axis). Points are colored according to the recipient plant species. B Heterospecific (across different species) transfers, where ASV-level richness (y-axis) is positively correlated with the donor plant biomass (g, x-axis). Point colors correspond to host species and point shapes correspond to the donor plant species. For both plots, R2 and p values are derived from linear models.

Full size image

Experimental transplantation of phyllosphere communities reveals an influence of current and previous host identity, as well as donor plant neighbor

PERMANOVA on Bray Curtis dissimilarities revealed that host species identity, previous host identity (i.e., the species from which inoculum was derived), and experimental block all significantly contributed to phyllosphere community structure differences among experimental greenhouse plants (Host species R2 = 0.127, p < 0.001, Previous host species R2 = 0.0468, p < 0.001, Block R2 = 0.156, p < 0.001, Table 2). We also observed an effect of donor plant biomass on recipient plant phyllosphere community structure (R2 = 0.018, p = 0.05, Table 2). Moreover, if we interrogate the dataset by plant species, we see a “grandparent effect” in the phyllosphere community structure of greenhouse-grown pepper plants. That is, pepper phyllosphere communities were most strongly shaped by the previous host’s neighbor (R2 = 0.194, p < 0.003), i.e., the plant species that previously surrounded the donor plant. In contrast to what was observed on pepper, neither previous neighbor effects nor previous host effects were observed for tomato or bean subsets.

Table 2 Results of a PERMANOVA on phyllosphere bacterial community Bray Curtis dissimilarities for the greenhouse experiment.
Full size table

Discussion

Plant-microbe associations form in part through host filtering of microbiota that arrive via dispersal [66]. Our study experimentally manipulated neighbor presence, identity, and age in order to understand how these factors influence host filtering of phyllosphere communities. Over the course of the experiment, we found that host species identity effects on focal plant microbiomes decreased, while the effects of neighborhood increased (Fig. 3A). This finding builds on past studies showing that host species- or genotype-level differences in microbiota change over the growing season (e.g., [18, 67, 68]). However, an important distinction is that we experimentally held the host developmental stage constant throughout the experiment, thereby demonstrating that changes in host species identity effects over time are not simply due to host ontogeny, but hinge on characteristics of the neighboring plant community such as neighbor identity and biomass.

The increasing strength of neighborhood effects throughout the field experiment suggests that as neighboring plants grow and enrich for host species-specific microbes, they become larger sources of microbial propagules to their surrounding neighborhood, and thus alter the outcome of host filtering through compositional changes to the local species pool. Two other lines of evidence from our study further underscore the importance of neighbor identity for phyllosphere community assembly. First, when we directly control the directionality of dispersal in our greenhouse microbial transplant study, we see that the source of inoculum (i.e., the species identity of the donor plant) significantly contributed to the microbiome composition of the recipient plants (Table 2). Second, the field experiment uncovered strong differences in host species identity effects depending on the identity of neighbors (Fig. 4). For instance, at harvests 2 and 3, hosts that were surrounded by tomato or bean neighbors were substantially less differentiable in their phyllosphere community structure than plants surrounded by pepper or that had no neighbors. This suggests firstly that having a neighbor impacts the differentiation of hosts, but crucially that the outcome of neighborhood effects depends on neighbor identity. Interestingly, it has also been reported that inter-host dispersal among zebrafish greatly diminished genotype-level microbiome differences [44]. Our results not only reinforce this concept in plants but suggest that this effect depends on neighbor identity.

One of the key differences among experimental neighborhoods was plant biomass, which varied by species (Fig. 1C). This is likely an important driver of the neighborhood effects observed in this study. Although we did not experimentally distinguish the effects of neighborhood age and biomass, larger biomass plants could harbor higher abundances of microorganisms by virtue of having more microbial habitat, thereby increasing the frequency of propagules dispersing onto focal plants. The effects of plant aging could shape the diversity of microbes, potentially filtering for more plant-associated or species-specific taxa. While we chose not to sample the epiphytic communities of the neighborhood plants in order to leave them undisturbed for the duration of the experiment, we observed a positive correlation between focal host biomass and epiphytic bacterial abundance. Thus, larger biomass neighborhoods could have diminished the strength of host filtering through mass effects, i.e., rescuing via dispersal the taxa that went locally extinct due to host selection. It is worth noting that neighborhood biomass and average neighbor height were negatively correlated with focal plant bacterial abundance. These negative trends could be driven by low rates of recruitment when specialized taxa are rare (e.g., in conspecific neighborhoods). We see evidence for the importance of neighbor biomass in several other results. First, as the higher-biomass tomato and bean neighborhoods grew, we see that focal plant species identity effects became weaker, so much so that by harvest 3 hosts were indistinguishable by their species identity if they were surrounded by tomato or bean (Fig. 4). Interestingly, at harvest 1, relative to no neighbor controls, focal hosts surrounded by tomatoes or beans exhibited higher species identity effects, suggesting that at an early stage, neighbor plants may bolster host filtering by providing a higher abundance and/or diversity of propagules. Moreover, the effects of the smaller pepper neighborhoods followed a similar but lagged trend, whereby host differentiation was highest at harvest time point 2, followed by more diminished host species effects at harvest 3. This could be driven by the observed slower growth of peppers relative to beans or tomatoes. In the greenhouse study, we see that the biomass of the donor plants significantly contributed to variation in both community composition (Table 2) and diversity (Fig. 6) of recipient plants. Here, the relationship between recipient plant phyllosphere richness and donor plant biomass was negative and weak if the transplant was conspecific (same species) but stronger and positive if it was heterospecific (between different species). Together these results underscore that biomass is not only an important component of neighborhood effects but that such influences of biomass may depend at least partly on the identity of the neighbors.

Our study also makes clear that differences in the strength of host filtering across species may impact susceptibility to neighborhood effects. We found that plant species exhibited neighborhood effects differently through time (Fig. 3B). Together with the observation of a host-by-neighborhood interaction effect, these results suggest that the relative impact of local neighborhood differs among focal host species, perhaps due to differences in the degree to which dispersal versus host filtering influence phyllosphere assembly. For instance, host-specific carrying capacities could be driving the observed differences in bacterial abundances across species (Fig. 2A). This could mean that species with lower abundances of bacteria (e.g., peppers) are more invasible, and hence taxa that are selected for may more quickly become outnumbered by immigrating taxa. This may explain why pepper plants exhibited neighborhood effects at harvest 2, while tomato and beans did not do so until harvest 3. Several lines of evidence also suggest that the plants may differ in their selective abilities. For instance, tomato- and pepper-associated communities consistently exhibited phylogenetic clustering (Fig. 5A), indicating that closely related taxa were often observed to co-occur on a single plant, perhaps due to host selection of shared traits. In contrast, bean microbiomes tended to be phylogenetically overdispersed and hence exhibited much less clustering than one would expect by chance. Overdispersion could indicate that traits under selection are not phylogenetically conserved, that competition is strong (e.g., between ecologically similar taxa), or that selection is relatively weak [19]. That we also see a strong fit of bean microbiomes to neutral models at harvests 1 and 2 (Fig. 5B) suggests that the bean phyllosphere may be less selective relative to the tomato and pepper phyllosphere, and therefore more susceptible to dispersal effects.

Interestingly in the greenhouse experiment, only pepper plant microbiomes exhibited a “grandparent effect” of inoculum, i.e., an effect of the donor plants’ previous neighbor. This result demonstrates that for certain species, microbiome composition not only reflects its contemporary host and its source history but also its previous dispersal history. This is analogous to a child horizontally acquiring a parent’s microbiome that carries with it traces of the parent’s former house, pet, or domestic partner. Results from the field trial may help shed light on why only pepper microbiomes contained detectable traces of dispersal history. The aforementioned patterns of phylogenetic clustering in pepper plant communities and the observation that neutral models failed to fit pepper plants (Fig. 5A, B) together suggest that the pepper phyllosphere may impose particularly strong selection on microbial communities. This strong filtering ability may not only select against taxa, but it could act to amplify taxa that previously dispersed from pepper neighbors, thus giving rise to effects of the previous neighbor. In other words, selective plants such as pepper may bolster pepper-associated microorganisms upon arrival, even if they have become rare through multiple dispersal events.

One major implication of our work is that it highlights the potential importance of local host frequencies in recruiting microbiota. At harvest 3, bean phyllosphere communities exhibited phylogenetic clustering only when the plants were surrounded by conspecific (same species) neighbors. As phylogenetic clustering could be interpreted as host filtering for phylogenetically conserved traits, this result suggests that for beans in particular, the frequency of conspecifics in the plant metacommunity may be an important determinant of host filtering efficacy. Interestingly, for pepper- or tomato-associated communities, the con- or hetero-specific status of neighborhoods had little influence over phylogenetic clustering, suggesting that this effect may depend on the strength of host filters. Similar findings have recently been reported for Acer saccharum trees [38], where the abundance of conspecific trees in the local metacommunity was positively correlated with the degree of host specialization in the phyllosphere. Moreover, in cacao trees, leaf litter of healthy conspecific hosts was shown to protect against pathogen damage [69]. While we did not test the fitness effects of host filtering for specialized microbial taxa, further investigation of this topic could bring about a new perspective on the conspecific negative density-dependence (i.e., Janzen–Connell) hypothesis, which posits that higher local densities of conspecifics may be disadvantageous due to the possibility of spreading specialized pests or pathogens [29, 30]. While there remain many examples of conspecific negative density-dependence [70], particularly in the tropics, meta-analyses seeking general trends have generated mixed results [71, 72]. It is therefore possible that recruitment of specialized beneficial taxa from nearby conspecific hosts could outweigh the negative effects of specialized pest/pathogen pressure in certain contexts, or for certain host species.

Overall, our work makes clear that local neighborhood identity and biomass are key components that shape the assembly of the phyllosphere microbiome. In both the field trial and greenhouse experiment, we find that although plants are able to select upon their microbial communities, the outcome of this selection is shaped by both neighbor identity and local biomass. Moving forward, this work has opened a number of critical questions regarding how neighborhood effects on the plant microbiome might shape plant health, fitness, and—in agricultural settings especially—yield. The work also raises questions about how invasive plant species might alter microbial dispersal within their communities, and potentially negatively feedback on native plant species’ fitness by reducing their ability to filter an optimal microbiome. In sum, our work demonstrates that host filtering and local dispersal are intimately intertwined and represent crucial considerations for the study of host–microbe associations.

Data accessibility

Raw sequence data can be accessed at the NCBI BioProject database (BioProject ID: PRJNA789311). Sequence data, sample meta-data, and reproducible R code are available in the following online repository: https://doi.org/10.6084/m9.figshare.16575194.v1.