Phytobiomes are compositionally nested from the ground up

Sampling

We selected our plant species and sites to minimize environmental heterogeneity. S. taccada is an evergreen shrub, native to Hawai’i, which produces fruits and flowers across seasons. It is also one of the most common littoral plants on Oahu’s East (windward) shore, enabling sampling at regular spatial intervals. The coastal habitat was selected because local climates at this elevation are minimally influenced by land topography, and because all coastlines are accessible to the public by virtue of state law. All samples were collected on January 26, 2017. A single S. taccada individual was sampled from each of 10 locations along the windward shore of O’ahu at 6 km intervals (Fig. 1). Generally, there were few individuals containing both reproductive parts (fruits and flowers) at a given site, and if more than one was present the sampled individual was selected haphazardly. Within the Kailua site, nine additional individuals were selected haphazardly within a 50 m² plot. Plant locations were recorded with a Garmin Rhino GPS. A single mature, disease-symptom free leaf, flower, fruit, stem and axil was sampled by swabbing with a sterile cotton swab that had been moistened with an alkaline Tris extraction buffer (containing EDTA and KCL; comparable to Extract-N-Amp solution; Sigma Aldrich, St. Louis, MO, USA). At two locations we were unable to locate and sample flowers. Root (rhizosphere) samples were collected by selecting an area of the root with diameter 0.5–1.7 cm, and buried 2–10 cm below the surface. Soil was collected on swabs taken within 10 cm of the plant base, at a depth of 2–5 cm. Because soil was likely the most heterogeneous component of our study, bulk soil samples were also collected for soil description; including soil taxonomy and pH. Finally, sterile swabs dipped in extraction solution were exposed to the air for approximately 20 s as extraction negative controls. After sampling, swabs were immersed in 200 µL extraction buffer and stored in a cooler until DNA extraction later in the day. Although soil is not technically a plant part, we refer to it as such when differentiating sample types below.

Nucleic acid extraction and library preparation

Swabs were incubated in 50 µL of the extraction buffer at 65 °C for 10 min and at 95 °C for an additional 10 min. After brief centrifugation, swabs were removed from the buffer, and 200 µL of 3% BSA was added to each sample. DNA extractions were purified and concentrated using a carboxylated paramagnetic bead solution. Briefly, 100 µL of gDNA was mixed with 80 µL of bead solution, immobilized with magnets, washed twice with 70% ethanol, and suspended in 30 µL of 0.1x TE buffer.

DNA library preparation and sequencing

PCRs targeted the v4 region of the 16S rRNA gene using the primers 515f (GTGCCAGCMGCCGCGGTAA) and 806r (GGACTACHVGGGTWTCTAAT) with overhangs complementary to the 3′ region of a construct containing dual 8bp indices and the Illumina i5 and i7 adaptors following the standard Illumina 16S Metagenomic Sequencing Library Preparation protocol (https://support.illumina.com/content/dam/illumina-support/documents/documentation/chemistry_documentation/16s/16s-metagenomic-library-prep-guide-15044223-b.pdf; Illumina, San Diego, CA), using KAPA3G plant PCR kit reagents. Briefly the first PCR was carried out in a 25 µL reaction consisting of 25 cycles. Amplicons were diluted 1:10 and used as a template for the second “Index PCR” in which sample-specific illumina indices and sequencing adapters were added using an additional 8 cycles.

Negative PCR controls (extraction and field blanks), were added to the library preparations. PCR products were purified and normalized to equimolar concentrations using the just-a-plate 96 PCR Purification and Normalization Kit (Charm Biotech, San Diego, California, USA). Normalized PCR products were pooled and sequenced at the Hawai’i Institute of Marine Biology CORE Laboratory using the Illumina MiSeq V3, 600 cycle, paired-end protocol. PhiX (15%) was added to the sequencing reaction to increase template complexity. A region of the flow-cell was deemed defective after its use, resulting in fewer high-confidence sequences data than are typically recovered using this technology. Sequence data were deposited to the SRA as PRJNA385181.

Data processing

Average reverse sequence quality scores declined to <25 after ca 150 bp and R2 sequences were, therefore, not further considered. Chimeras were detected and removed from 16S sequences using VSEARCH (Rognes et al., 2016). Sequences were filtered by length (75 bp min) and quality score (mean score 25) and demultiplexed within the QIIME environment (Caporaso et al., 2010). OTU binning at 97% identity and taxonomic assignments were conducted in QIIME’s “open reference” workflow using default parameters (Caporaso et al., 2010) and the Greengenes 13_8 reference database. Taxa that were not assigned to bacteria (chloroplasts, and two mitochondria sequences) were removed from datasets, as were OTUs represented by ten or fewer sequences. None of the extraction or field negatives produced sequence data.

Because of uneven sequencing depth across our dataset, data were down-sampled to 1,000 reads per sample, a cutoff that included the majority of samples. Although there is no consensus about how best to treat uneven sampling depth, down-sampling to a common depth is generally robust for multivariate comparisons of community compositions (Weiss et al., 2017). Eighty-five samples containing greater than 1,000 reads were retained for subsequent analysis.

Data analyses

To determine which factors (site and/or plant part) predicted community composition, binary Jaccard distances were used in a single PERMANOVA model (function “adonis” in the vegan package (Oksanen et al., 2013)) considering marginal values of both factors and their interactions. In order to balance the number of samples among sites, a single individual plant was selected for the Kailua site. This analysis was repeated, including all samples from the Kailua site with qualitatively similar results. Because the analysis can be sensitive to dispersion of beta-diversity values, these were evaluated separately for plant parts. Following a globally significant RDA analysis showing that community dissimilarity was predicted by site, plant part and their interaction, a variance partitioning approach was used to determine the proportion of variance explained by each.

To address whether soil or plant parts might serve as a microbial reservoir for self-inoculation, we used two measures of nestedness, both of which indicate the extent to which OTUs present in species-poorer plant parts are contained in those with higher species richness. A bipartite matrix was calculated by summing all reads within each plant part and then downsampling to the lowest sample sum. We calculate the nestedness temperature (Atmar & Patterson, 1993), a popular, though older, method calculating nestedness as a range from 0 (perfectly nested) to 100 (perfectly random) based on hierarchical entropy order and presence absence data. We also calculated nestedness based on overlap and decreasing fill (NODF; Almeida-Neto et al., 2008) scaled from 0 to 1. Both were compared with a distribution of randomized null communities generated using the quasiswap method in vegan. Differences in community membership among habitats were decomposed into nestedness and turnover components (Baselga & Orme, 2012) as implemented in the nestedbetajac function in vegan.

The extent to which certain plant parts “select for” microbes, was calculated using the network-wide H2 index (Blüthgen, Menzel & Blüthgen, 2006), which returns a value from 0 (generalized) to 1 (specialized) based on potential associations given OTU abundance totals. Significant deviations from null expectations were quantified using a distribution of null community matrixes calculated using the quantitative r2dtable method in vegan in which marginal totals are maintained. We calculated the related d’ statistic (Blüthgen, Menzel & Blüthgen, 2006) to assess selectivity of individual plant parts.

To test for dissimilarity by distance patterns, binary Jaccard community dissimilarities were compared with geographic distances between samples using Mantel correlations (tests using Bray Curtis indices were qualitatively identical). So that the intensively sampled Kailua site did not skew detected spatial structure, we randomly selected a single Kailua individual for analysis (Mantel correlations including all samples showed similar effect sizes, except for above ground plant parts, which no longer showed a significant pattern; Table S2).

Indicator species analysis was calculated using the indicspecies package (De Caceres & Jansen, 2016) on single groups (multi-group combinations were suppressed due to the high number of OTUs).

Sequencing results

Illumina sequencing of the V4 region of bacteria 16S gene generated 15,354,523 reads. Following demultimplexing, quality score filtering and chimera detection 8,901,385 bacteria reads remained. Removal of OTUs identified as chloroplasts and mitochondria, and OTUs present in read abundances <10 resulted in 4,643,809 reads and 2,192 bacteria OTUs retained for subsequent analysis.

Specialization and nestedness of habitats

Bacteria were significantly more nested than null expectations (Table 1; Fig. 1) and nestedness contributed 9% of total between-sample dissimilarity. Nestedness hierarchy followed the vertical structure of the plant: below-ground samples contained the highest species richness, above-ground vegetative structures were intermediate, and reproductive structures were lowest in the hierarchy (Fig. 2).

Microbes were significantly specialized on particular plant parts. We examine the degree and nature of specialization across plant parts using a bipartite network architecture index, and also within plant part using an indicator species analysis. Network-wide degrees of specialization were high and 762 “indicator species” were statistically associated with plant parts (Table S3). The d’ statistic, which measures degree of specialization for individual habitats, showed that soil and rhizosphere were the most specialized parts.

Plant part and abiotic drivers of community composition

A PERMANOVA model demonstrated that both plant part and location shape S. taccada phytobiomes (Table 2), but that there was no interaction between those terms. Further examination of community partitioning demonstrates that plant part is a stronger determinant of bacteria community than location. Partitioning of variance, based on RDA analysis, demonstrated differences in the relative importance contributed by these two variables. Plant part accounted for 88% of the explained variance (RDA, DF = 50, F = 1.33, P = 0.043), although residual variance accounted for 81% of the total, suggesting that other unmeasured variables may be more predictive of community composition.

NMDS ordination (Fig. 3) demonstrates that bacteria communities assorted into three principle groups associated with plant parts: belowground (soil and rhizosphere), mid-plant vegetative (stem and axil), and reproductive (fruit and flower). Leaf samples are spread among the other above-ground plant parts. Dispersion of bacteria beta diversity (bacteria betadispersion ANOVA, DF = 6, P = 0.08; Fig. S1) did not differ significantly between any sample types, suggesting that results of the PERMANOVA analysis reflect differences in multidimensional community “location” rather than spread. Bacteria were strongly partitioned at the order level by plant parts (Fig. 4). Bacillales were ubiquitous and abundant in soil and leaf samples, whereas Rhizobiales and Actinomycetales were abundant across root and vegetative plant parts.

While site location was a strong predictor of compositional differences, tests of dissimilarity by distance for microbial community compositions showed mixed results. There was a strong distance-decay relationship for below-ground samples and the relationship for above-ground plant parts and all samples was weaker, though significant (Table 3). Although geographic distance was generally a poor predictor of microbial composition, site location strongly impacted community composition (Table 2). This may be due to a deterministic influence of environmental variables that were not accounted for in our study, including the influence of the surrounding vegetation.

Distribution of OTUs across the plant

We show that the surface phytobiome composition of S. taccada is strongly shaped by plant part. Both community similarity and richness patterns are structured by plant parts: stem, axil, and reproductive parts were most similar to each other, rhizosphere samples were most similar to soil samples, and leaves were allied with above-ground parts more generally. Segregation of microbial communities by plant part (Junker et al., 2011; Ottesen et al., 2013; Junker & Keller, 2015) and even within-individual part location (Leff et al., 2015) has been reported for bacteria previously. Notably, Leff et al. (2016) demonstrated that both fungi and bacteria are partitioned by plant part (seed, rhizosphere and root) among sunflower cultivars. By examining differences among plant parts within a network context, however, we gain insight into the potential sources and sinks of plant microbes.

Compositional differences between samples might be attributed to both species replacement and nestedness components, and the relative contributions of these variables can suggest causal mechanisms (Baselga, 2010). High turnover typically suggests stochastic assembly processes, whereas nestedness is more characteristic of deterministic factors such as environmental filtering (Si et al., 2016). We found that 9% of compositional differences was due to nestedness, a number that is consistent with the global mean of 0.093 ± 0.054 (SD), calculated in a recent meta-analysis of 99 studies, including microbes (Soininen, Heino & Wang, 2018). In contrast, the proportion of beta-diversity explained by turnover, 84%, was significantly higher than the global mean of 53% ± 18.4 (SD).

We show that soil may be an important source for plant-surface microbiomes, including above-ground parts. Our analysis demonstrates that plant microbes are nested with respect to plant parts, that nestedness accounts for a large proportion of the among-sample diversity, and that below-ground parts (soil and rhizosphere) are at the top of the hierarchy. Even putative “soil” dwelling microbes, such as Bacillales, are found throughout the plant. These results suggest that soils may serve as an important source for plant-surface microbes, either through initial inoculation as the plant emerges from the ground, or repeated introduction via wind and water dispersal. An alternative hypothesis, that soils and plants could both be sinks from an alternative source (such as airborne bacteria), with nestedness patterns resulting from ordered extinctions, could be tested via experimental manipulations.

Several other studies have noted nestedness patterns in other microbial systems, including mycorrhizal tree islands displaying strong species area relationships (Peay et al., 2007), and mouse guts whose bacterial diversity decreases, orderly, over generations (Sonnenburg et al., 2016). Notably, bacteria phyla from 97 independent studies spanning diverse geography and habitat were significantly nested (nestedness decreased with taxonomic rank; Thompson et al., 2017). Although this study did not consider plant parts per se, they also found significant nestedness patterning among plant samples, with rhizosphere being the most species rich category, followed by plant surfaces, and internal plant tissues.

Despite nestedness patterns overall, each plant part contained a distinct subset of indicator bacteria that were not shared. These patterns suggest that distinct assembly processes may govern different plant parts. Pollinators, for example, are known to be effective vectors for microbes found in nectar (Herrera et al., 2010; Aleklett, Hart & Shade, 2014). Less is known about the origins or mechanisms of inoculation of other plant parts.

Factors other than physiochemical plant traits also likely play a role in shaping phytobiomes. The strong division between below- and above-ground communities, for example, may result from dispersal limitation in addition to climatic differences. Furthermore, plant part longevity may play a role in compositional differences. Flowers and fruits, the plant parts with lowest species richness, are also the most ephemeral, present on a plant for days to weeks. S. taccada flowers studied on the French Polynesian Island of Moorea, for example, remain open for an average of fewer than four days (Liao, 2008), seemingly scant time to recruit microbes from the environment.

Plant microbial distributions over space

In microbial biogeography studies, distance between samples is generally a reliable predictor of microbial community dissimilarity, particularly among communities that are not attached to a macro-organism (Martiny et al., 2011; Zinger, Boetius & Ramette, 2014; Tedersoo et al., 2014), even at the scale of centimeters (Bell, 2010). The distance decay of plant-associated microbial similarity is less resolved compared with studies of “host-independent” substrates such as ocean, water, or soil. A positive pattern could be attributed to at least two factors and or their cumulative effects. First, dispersal limitation could lead to clinal dissimilarity among communities by enabling “ecological drift” over relatively short timescales (Martiny et al., 2006). Second, geographic distance could correlate with some environmental cline, which selects for a distinct adapted community. Factors such as rainfall (O’Rorke et al., 2017), soil pH (Fierer & Jackson, 2006), soil nutrients, or even host genotypic clines (Bálint et al., 2013) all contribute to community dissimilarity.

Among plant-associated microbes, evidence for the dissimilarity by distance pattern is mixed. Redford et al. (2010) found little evidence of geographic structure in pine bacteria phyllosphere communities spanning continents, and Meaden, Metcalf & Koskella (2016) found no evidence for dispersal limitation among Oak phyllosphere bacteria within a 20 ha plot. In contrast, Stone, Bacon & White (2000) found strong isolation by distance patterns among Magnolia phyllosphere bacteria located within 400 m of each other. Oono, Rasmussen & Lefèvre (2017) found that low abundance pine needle fungal endophytes were spatially structured over ∼100 km, but that high abundance communities were not.

In our study, both above-ground and belowground bacteria were correlated with geographic distance (Table 3), although the latter relationship was stronger. Differing community turnover rates among plant parts could be attributable to different dispersal rates since airborne bacteria would presumably travel farther and more quickly due to prevalent trade winds compared to subsurface soil bound communities. Although there were no obvious environmental clines along our transect, microbes may have responded differently to additional unmeasured variables. Multiple studies show that bacteria are particularly sensitive to soil properties such as pH. Substrates change rapidly over short distances from the high tide, and although pH did not covary in a linear fashion along the transect, pH was variable among sample sites. For example at a single site for this study, we sampled soils ranging from ph 5.28 to 8.65 (Table S1). Finally, differences in the distance decay relationship could be driven by interactions with the host plant. Obligacy or strong host selection of some community members, for example, could depress turnover if the microbe was required for growth and survival.