Tree phyllosphere bacterial communities: exploring the magnitude of intra- and inter-individual variation among host species

Study site & host-tree species

The two study sites are located in a natural temperate forest stand in Gatineau (45°44′50″N; 75°17′57″W) and Sutton (45°6′46″N; 72°32′28″W) Quebec, Canada. These sites are characterized by a cold and humid continental climate with temperate summer. A total of six individuals (three at each site) from each of five tree species common to temperate forests and dominant in the canopy were sampled to provide representatives of both angiosperms and gymnosperms: Abies balsamea (Balsam fir), Acer rubrum (Red maple), Acer saccharum (Sugar maple), Betula papyrifera (Paper birch) and Picea glauca (White spruce).

Bacterial community collection

We sampled phyllosphere communities from trees on August 29, 2013 as part of another experiment (Laforest-Lapointe, Messier & Kembel, 2016). Sampling was carried out one week after the last rainfall event. We defined three strata within the canopy: bottom-canopy (1–2 m height), mid-canopy (2–4 m height), and top-canopy (4–6 m height). 30 individuals were randomly selected by picking random geographic coordinates and finding the closest individual at this location. For the first tree sampled from each species, we clipped 50–100 g of leaves at the four cardinal points at mid-canopy height, plus a single sample at bottom-canopy and top-canopy heights, into sterile roll bags with surface-sterilized shears. We also sampled bottom-canopy leaves from two other randomly chosen trees from each species. For bacterial community collection and amplification we used the protocols described by Kembel et al. (2014). We collected microbial communities from the leaf surface by five minutes of horizontal mechanical agitation of the samples in a diluted Redford buffer solution. We resuspended cells in 500 µL of PowerSoil bead solution (MoBio, Carlsbad, California). We extracted DNA from isolated cells using the PowerSoil kit according to the manufacturer’s instructions and stored at −80°C.

DNA library preparation and sequencing

We used a two-step PCR approach to prepare amplicon libraries for the high-throughput Illumina sequencing platform. The use of combinatorial primers for paired-end Illumina sequencing of amplicons reduced the number of primers while maintaining the diversity of unique identifiers (Gloor et al., 2010). First, we amplified the V5–V6 region of the bacterial 16S rRNA gene using chloroplast-excluding primers in order to eliminate contamination by host plant DNA (16S primers 799F-1115R (Redford et al., 2010; Chelius & Triplett, 2001)) following protocols described by Kembel et al. (2014). We cleaned the resulting product using a MoBio UltraClean PCR cleanup kit. We isolated a ∼445-bp fragment by electrophoresis in a 2% agarose gel, and recovered DNA with the MoBio GelSpin kit. We prepared multiplexed 16S libraries by mixing equimolar concentrations of DNA, and sequenced the DNA library using Illumina MiSeq 250-bp paired-end sequencing at Genome Quebec.

We processed the raw sequence data with PEAR (Zhang et al., 2014) and QIIME (Caporaso et al., 2010) software to merge paired-end sequences to a single sequence of length of 350 bp, eliminate low quality sequences (mean quality score <30 or with any series of 5 bases with a quality score <30), and de-multiplex sequences into samples. We eliminated chimeric sequences using the Uclust and Usearch algorithms (Edgar, 2010). Then, we binned the remaining sequences into operational taxonomic units (OTUs) at a 97% sequence similarity cutoff using the Uclust algorithm (Edgar, 2010) and determined the taxonomic identity of each OTU using the BLAST algorithm (Greengenes reference set) as implemented in QIIME (Caporaso et al., 2010). The number of sequences per sample ranged from 6,256 to 75,412. From these 1,499,777 sequences, we rarefied each sample to 5,000 sequences and repeated analyses on 100 random rarefactions. Re-analysis did not quantitatively change results and so we report only the result of the analysis of a single random rarefaction. We included the resulting 275,000 sequences in all subsequent analyses.

Statistical analyses

We created a database excluding OTUs represented fewer than 3 times to minimize the presence of spurious OTUs caused by PCR and sequencing errors (Acinas et al., 2005). We identified the OTUs that were present on all samples to define the “core microbiome” (Shade & Handelsman, 2012). Then we tested for significant associations between bacterial taxa and host species, and canopy location using the Linear Discriminant Analysis Effect Size (LEfSe) algorithm (Segata et al., 2011). This analysis allows the recognition of significant individual host-microbe associations and evaluates the strength of associations between organisms from different groups (Segata et al., 2011).

We performed analyses with the ape (Paradis et al., 2004), picante (Kembel et al., 2010), and vegan (Oksanen et al., 2007) packages in R (R Development Core Team, 2013) and ggplot2 (Wickham, 2009) for data visualization. We quantified the taxonomic variation in bacterial community structure among samples with the Bray–Curtis dissimilarity. To illustrate patterns of bacterial community structure, we performed a nonmetric multidimensional scaling (NMDS) ordination of Bray–Curtis dissimilarity. We identified relationships between bacterial community structure, host species identity, and sample canopy location by conducting a permutational multivariate analysis of variance (PERMANOVA, Anderson, 2001) on the community matrix. We employed a blocking randomization to account for the non-independence of observations among sites. To decompose the total variation in the community matrix explained by host species identity and canopy location, we performed a partial redundancy analysis (RDA; Legendre & Legendre, 1998). This technique measures the amount of variation that can be attributed exclusively to each set of explanatory variables. We performed three permutational tests of multivariate homogeneity of group dispersions (Levene’s test for variances’ homogeneity multivariate equivalent; Anderson, 2006; Anderson, Ellingsen & McArdle, 2006): one to test if variance in intra-individual canopy bacterial communities was equal between individuals (30 samples from five trees sampled at six canopy locations); a second to compare interspecific variation between species (30 bottom-canopy samples from 30 different trees); and finally a third to test per-species intra- and inter-individual variation (all 55 samples). We estimated phyllosphere bacterial alpha-diversity using the Shannon index calculated from OTU relative abundances for each community. We performed an analysis of variance (ANOVA) and subsequent post-hoc Tukey’s tests to compare differences in diversity across species. The authors declare that the experiment comply with the current laws of the country in which the experiment was performed.

Sequences, OTUs and taxonomy

High-throughput Illumina sequencing of the bacterial 16S rRNA gene (Claesson et al., 2010) identified 5,005 bacterial operational taxonomic units (OTUs, sequences binned at 97% similarity) in the phyllosphere of five temperate tree species, an average of 1,055 ± 57 OTUs (mean ± SE) per tree sampled. Most of these bacterial taxa were relatively common across samples, with only 3.4% of OTUs occurring on a single tree and 0.8% of OTUs occurring on all trees. The OTUs present on all samples represent the “core microbiome”: the microbial taxa shared among multiple communities sampled from the same habitat (Shade & Handelsman, 2012). In this study, the core microbiome consisted of 42 OTUs (Table 1) representing 61% of all sequences, of which 72% were Alphaproteobacteria, 9% Cytophagia, 7.8% Betaproteobacteria, 5% Acidobacteria, 2% Gammaproteobacteria and 2% Actinobacteria. The most abundant order was Rhizobiales (49%) from which 77% of sequences were assigned to the family Methylocystaceae. While there was some variation in the most abundant classes both across the five tree species and among canopy locations (Figs. 1 and 2), the class Alphaproteobacteria was always the dominant taxon, with relative abundances ranging from 42% on P. glauca to 84% on B. papyrifera (Fig. 1).

Intra-individual vs. inter-individual and interspecific variation

Host species identity and individual identity effects could not be distinguished statistically due to the fact that analyses of intra-individual variation were based on a single individual per species. This host species/individual effect explained 65% of variation in phyllosphere bacterial taxonomic community structure while the impact of canopy location was not statistically significant (PERMANOVA on Bray–Curtis dissimilarities; Table 2). We then tested whether canopy position had an effect on community structure after accounting for the variation explained by host species/individual using a partial redundancy analysis (RDA) on bacterial community structure constrained by host species identity. The RDA showed that when differences in bacterial community structure driven by host species identity were accounted for, sample canopy location explained 22% of the remaining variation in community structure. In comparison, in the dataset with 30 different individuals, host species identity explained only 47% of variation in phyllosphere bacterial community structure (PERMANOVA on Bray–Curtis dissimilarities; Table 2). When considering intra-individual and inter-individual samples, host species identity (R² = 47%) was the strongest driver of variation in phyllosphere bacterial community structure closely followed by individual identity (R² = 32%) and finally by canopy location (R² = 6%; PERMANOVA on Bray–Curtis dissimilarities; Table 2). Community composition of samples clustered based both on the individual (Fig. 3A) and species (Fig. 3B) from which they were collected (non-metric multidimensional scaling (NMDS) based on Bray–Curtis distances among samples).

The first permutational multivariate test of variance homogeneity (an analogue of Levene’s test of homogeneity of variances) on intra-individual phyllosphere communities indicated a significant difference between P. glauca and B. papyrifera (Tukey’s post hoc test; P = 0.03). The second test of the homogeneity of inter-individual variance between host species showed that P. glauca’s variance in community structure (mean distance to centroid = 0.34) was higher than A. saccharum (0.25; P < 0.01) and A. rubrum (0.26; P < 0.05) while all other comparisons were not significant. Finally, the third test between per species intra-individual and inter-individual variation indicated one significant difference in variation for B. papyrifera (P = 0.005; Fig. 4).

The alpha-diversity of leaf bacterial community differed significantly across host species identity but not across canopy locations. Post-hoc Tukey honestly significant differences tests confirmed that Shannon alpha-diversity is higher on conifer species (4.9 ± standard error (SE) of 0.04 for A. balsamea and 5.3 ± SE 0.04 for P. glauca) than on angiosperm species (3.7 ± SE 0.06 for A. rubrum, 4.1 ± SE 0.05 for A. saccharum and 3.6 ± SE 0.09 for B. papyrifera).

Bacterial indicator taxa

The LEfSe analysis successfully identified indicator taxonomic groups associated with different host species, but not across different canopy locations (Table 3). The conifers, A. balsamea and P. glauca, had the highest number of associated bacterial indicator taxa (46 and 188 respectively). The strongest bio-indicators of A. balsamea were the Frankiaceae family and multiple taxonomic levels of the phylum Acidobacteria: Acidobacteria, Acidobacteriales and Acidobacteriaceae. For P. glauca, the strongest bioindicators were multiple taxa from the Bacteroidetes phylum (Cytophagia, Cytophagales, Cytophagaceae, Spirosoma and Saprospirae, Saprospirales, Chitinophagaceae), and from the Actinobacteria, Chloroflexi, and Deltaproteobacteria. In contrast, B. papyrifera showed an overrepresentation of 24 bacterial taxa including the phylum Proteobacteria, the class Alphaproteobacteria and several of its orders (Rhodospiralles, Rickettsiales, Caulobacterales). Finally, the two Acer species (A. rubrum and A. saccharum) were associated with 19 and 32 indicators respectively, including the order Rhizobiales: A. rubrum being associated with the family Methylocystaceae and A. saccharum with the order Methylobacteriaceae.