Unexpected Diversity During Community Succession in the Apple Flower Microbiome

Shade et al. Apple flower microbiome

Supplemental Text for

Unexpected diversity during community succession in the apple flower microbiome

Ashley Shade, Patricia S McManus, and Jo Handelsman

Supplemental Methods

Nucleic acid extraction

Prior to nucleic acid extraction, flowers were stored frozen at -80 C. Samples of flowers were removed from the storage freezer, flash frozen in liquid nitrogen, and massed. Microbial cells were separated from flower material before nucleic acid extraction. Flowers were washed in 35-mL cold 1x PBS-0.15% Tween solution with five to seven 3-mm glass beads as a mild abrasive. Tubes were secured horizontally in a polystyrene foam container filled with ice, shaken at maximum speed on a rotator at 4 °C for 20 min, and sonicated for 5 min in a water bath (model 2210, Branson Ultrasonics Corporation, Danbury, CT). Large flower debris was removed by filtering the liquid through sterile cheesecloth. The remaining flower debris was pelleted by centrifugation at 1500 rpm (272 x g) and 4 °C for 5 min (model Avanti J-E, Beckman-Coulter, Brea, CA). Supernatant was transferred to a clean tube, and centrifuged for 30 min at 6500 rpm (5111 x g), 4 °C to pellet the microbial cells. The pellet was re-suspended with buffer CLS-TC from the FastDNA spin kit (MP Biomedicals, Solon, OH). DNA was extracted following manufacturer’s instructions with minor modifications. Samples were homogenized for 1 min using a vortex adapter. We performed the optional 5-min incubation at 55 °C to maximize DNA recovery prior to the final spin. DNA was quantified using a Nanodrop spectrofluorometer (ThermoScientific NanoDrop Products, Wilmington, DE). Samples were stored at -20 °C until amplification. Extracted DNA concentration was between 19 and 70 ng/µL.

PCR was performed according to Research and Testing Laboratories standard protocols, and included a 30 cycle PCR with a mixture of HotStart and HotStar high fidelity Taq polymerases. Tag-encoded FLX amplicon pyrosequencing utilized Roche 454 FLX instrument with Titanium reagents and procedures. The tag-sequences used and sample descriptions ("mapping file") are given in Table S1A. The average product length was 335 nucleotides (See Supplemental Methods).

Sequence analyses

Research and Testing Laboratories performed sequence denoising and chimera-checking, as per their standard analysis protocols ( Additional quality control and processing of sequences was performed using the default QIIME 1.3 workflows. There were 225,243 raw input sequences. Sequences with ambiguous bases, a mean quality score below 25 within a 50-bp window, more than six homopolymers, any primer mismatches, missing quality scores, or not between 200 and 1000 bp in length were removed. From these quality control steps, 171,996 sequences passed, resulting in a minimum 1838, maximum 10776, and mean 5733 sequences per sample. Rarefaction to 1838 sequences was conducted for all samples prior to community analyses.

It is common to have a wide range of amplicons per sample in 454 sequencing because of sample to sample variability (e.g., as in (1)), which is why datasets are often rarefied to an equal number of sequences per sample prior to analysis (2). For this dataset, we noticed an interesting trend of an increased number of amplicons recovered and the age of the blossom, and subsequently showed a significant, direct correlation between these (see Main Text, and Supplemental Results below). Thus, we suggested that the gradual increase in the number of amplicons is the result of increase microbial load on each flower through time, as studies of Erwinia amylovora infection have shown that bacterial growth rate increases with apple flower age (3, 4). Note that the dataset was subsampled (rarefied) to an equal number of sequences per sample prior to any community analyses.

Operational taxonomic units (OTUs) were set at 97% sequence identity with uclust (5) and taxonomy was assigned using Ribosomal Database Project (RDP, (6)) with a minimum confidence of 0.8. To be conservative with diversity estimates, singleton OTUs (i.e., those detected only once over the entire dataset) were removed from the rarefied dataset. Sequences were aligned with PyNAST using the green genes template (7), and a phylogenetic tree was built with FastTree (8) to calculate UniFrac distances between samples (singletons omitted, rarefied to 1838 sequences), a measure of beta-diversity (9-11). Weighted UniFrac distances ultimately were chosen over unweighted to prevent over-emphasis of rare taxa in an uneven microbial dataset, and because of its greater explanatory value in describing beta-diversity (principal components analysis of weighted UniFrac axis 1 = 32.9% of variation explained, unweighted = 12.2%). Using QIIME, rarefied richness (no. OTUs) and rarefied Faith’s phylogenetic diversity metric (12) were calculated as measures of alpha (within-sample) diversity.

Among 225,243 raw sequences, 171,996 met our quality control criteria and were retained for analysis. The average sequence length was 335 nucleotides. After operational taxonomic unit (OTU) picking and clustering at 97% sequence identity, 4,531 OTUs were identified. Greater than half of these OTUs occurred only once throughout the dataset (singletons). To estimate diversity conservatively, singletons were removed after re-sampling (rarefaction) to 1,838 sequences per sample, the minimum number of sequences observed within a sample. Omitting singletons did not affect beta diversity (community differences across samples), as assessed with Mantel test (Pearson’s R between full dataset and dataset omitting singletons = 0.997, p < 0.001), and Procrustes rotation (p = 0.001). After singleton removal, 1,677 OTUs from 50,865 tag-sequences remained in the analysis.

Evaluating differences in microbial community structure

There are two multivariate ways that distinguish whether communities have different structures, by their centroid (mean) or by their spread (variability (13-15)). Different membership is visualized by different centroids in an ordination, representing different mean composition across groups (Figure S1Aa andS1Ab). Communities may also distinct because of differences in the variability, or spread (Figure S1Ab andS1Ac). Difference in variability is visualized by dispersion from a centroid. Thus, we asked whether there were differences in either centroid or spread across trees, as assessed by weighted UniFrac distance. Weighted UniFrac distance accounts for differences due to the number, relative abundance, composition, and phylogenetic representation of OTUs (16).

Hierarchical clustering

Clustering is an alternative to ordination for uncovering patterns across species or samples. All 1,677 OTUs from the 30 samples were included in the analysis. To cluster the OTUs by similar patterns and not by similar abundances, each OTU was relativized to its total number of occurrences (Figure S3A). Then, pair-wise Bray-Curtis similarities were calculated between all OTUs. The similarities were used to begin the hierarchical clustering based on the complete linkage algorithm of the hclust function in R. The algorithm first combined the most similar OTUs to form clusters, and then combined the most similar clusters iteratively until one large cluster was achieved. We explored different levels of clusters for underlying temporal patterns. The first level below the largest cluster contained six clusters (Figure 4a), and we found that these six clusters were comprised of groups of temporally cohesive OTUs. Finally, before creating bar charts, the dataset was relativized by the total number of sequences observed on each day.

MultiCoLA

To understand how omitting rare taxa from the dataset affected beta diversity, we performed a multivariate cutoff level analysis (MultiCoLA; (17)) using two standard metrics of community similarity, Bray-Curtis and Sørenson, as well as weighted UniFrac distance (Table S3). In MultiCoLA, the least abundant taxa are removed from the dataset in increasing proportions, and then each reduced dataset is correlated to the full dataset to determine how the omission of rare taxa influences beta diversity.

Supplemental Results

Temporal analysis of apple flower communities

Because time was a significant descriptor of microbial phylogenetic diversity and variability (and because we did not detect differences across trees), we aggregated the six trees each day to explore further temporal patterns, and found that sequencing and sampling efforts were sufficient but not exhaustive for all days (Figure S4A), likely because of the high frequency of rare OTUs. Temporal patterns were not apparent at the phylum level (Figure S1). For example, no time points were distinguished by the presence or absence of particular phyla. Proteobacteria decreased in relative abundance when the flowers opened, while Archaea and Chloroflexi increased and maintained the same proportion until the final collection. TM7 also increased in relative abundance at flower open. Verrucomicrobia peaked on 01 May, early flowering, although it was relatively abundant at all time points.

There was no difference in Pielou’s evenness (equitability of taxon representation) through time, except between 01 May and 29 April, and 01 May and 02 May (analysis of variance between days, F = 3.09 on 4 d.f., p = 0.03; Tukey’s HSD p = 0.04 for both above comparisons, p > 0.05 for all others). In both cases, the 01 May community had lower evenness and higher variability. A similar pattern was observed for richness (analysis of variance F = 4.97 on 4 d.f., p < 0.05). There were no differences in the number of OTUs across days, except for 01 May, which was significantly different from both 29 April and 02 May (Tukey HSD p < 0.5 for both). However, there was no general trend in evenness or richness through time.

It is expected that the number of sequences will vary across samples due to differences in sequencing reactions, which is the rational for rarefying datasets to an equal sequencing depth (2). However, this variability should not be systematic (e.g., through time, as in Figure S2). We found that the number of sequences recovered, potentially serving as a very rough approximation for the microbial load on flowers (18), peaked in the middle of the time series, correlating with flower biomass (Pearson’s r = 0.42, p = 0.02, Figure S2). We speculate that the increase in the number of amplicons is because of increased microbial load on each flower through time, as studies of Erwinia amylovora infection have shown that bacterial growth rate increases with apple flower age (3, 4). The correlation with flower biomass may be due to increased flower niche spaces available for microbial colonization, such as on expanding petals before peak flowering. However, because of known technical biases in DNA extraction, PCR efficiency, primer specificity, and 454 noise, this result should be interpreted with caution.

Diversity and representation of OTUs across phyla and successional groups

To visualize the occurrences of OTUs across phyla, a tree was built omitting unidentified Bacteria (Figure 6, main text). From OTUs identified to the phylum level, there were 24 that had relative abundances between 0.006 and 0.10; 18 of these were Deinococcus-Thermus or TM7. The other OTUs were evenly represented, and had relative abundances less than 0.006. Though proportionally less than Deinococcus-Thermus and TM7, Proteobacteria, Bacteroidetes, and Actinobacteria were prevalent in apple flower communities.

We also tested for differences in evenness and richness across successional groups. There was evidence for modest differences in richness (analysis of variance p = 0.08) and clear differences in evenness (Table S5). Specifically, the Mid successional group had lower evenness than the other groups (all p < 0.001). Indeed, the dominant OTU in the dataset, affiliated with Deinococcus-Thermus, was a member of the Mid successional group (main text Figure 6 pink bars). Additionally, the Climax group had significantly higher evenness than the others (all p < 0.05), except the Early-Succession group (p = 0.20, main text Figure 6 dark green bars). All other successional groups did not differ in evenness. These results suggest that, with a few exceptions, successional groups contained similar and diverse representation of phylogenetic lineages, and that the dynamics of any one phylum, order, or family did not characterize overall successional patterns.

Normalizing the dataset to flower biomass

We asked whether normalizing the community dataset to flower biomass would influence our interpretation of community patterns, and found that the biomass-normalized dataset exhibited the same community patterns as the original dataset (Protest R = 0.998, p = 0.001). Therefore, we continued analysis using the dataset that was not normalized to flower biomass.

More and interesting network analysis results

For all significant associations, positive (collinear or direct) associations were distinguished from negative (inverse or indirect) correlations. Positive associations are representative, for instance, of mutualisms or shared resources, while negative associations are representative, for instance, of antagonism. However, the biological processes underpinning the direction of these associations cannot be informed by the LSA. There was a similar number of positive and negative associations (772 positive associations, 817 negative), but the majority of associations were time-delayed (1233 delayed associations, 322 unilateral).

Associations additionally could be described by whether they were unilateral (occurring on the same day) or delayed in time (lagged). For example, a time-delayed association may occur if the growth of one taxon alters environmental conditions, and, after time, the altered environment impacts the growth of a neighboring taxon. The proportion of delayed and unilateral associations was consistent across successional groups (Figure S5A). Generalist taxa had a greater proportion of negative associations than others (Figure S5A), supporting the possibility that many Generalist OTUs, though persistent, perhaps were not as competitive as the members of the other successional groups.

Among the 20% most prevalent OTUs for each tree, 72-84% of OTUs had at least one significant association with another OTU (Table 3). However, among the 20% most prevalent OTUs for the full network, 52% had associations. In general, the networks of individual trees had comparable diameters, geodesic distance (a metric of the "small world effect"- the average distance from any one OTU to any other), clustering coefficients (a metric of the number of OTUs that are tightly associated together in small cliques). An exception was Tree 4, which had a smaller diameter than the others. Furthermore the average number of associations (mean degree) varied from 15.20 to 28.2 among trees. One interesting emergent property of the full network, though, that was not manifested at the individual tree level, was the clustering coefficient, C, which was much higher than those observed in individual trees. This suggests that the full network is not simply the sum or average of the individual tree networks, but that using the trees as biological replicates informs identification of clusters of associated OTUs that otherwise may not have be apparent. These basic network properties suggest that, like many other networks (e.g.,(19)), the apple flower microbiome has highly clustered OTUs with a clustering coefficient of 0.32, but also maintains the "small world" property of minimal linkages between disparate nodes (here, represented as OTUs). Notably, the random network generated with 175 nodes and 1539 edges (to be comparable to the full network), had comparable geodesic distance and mean degree, but a higher diameter and much lower clustering coefficient than the full network.

Detection of Erwinia affiliated taxa

One taxon out of the full dataset (singletons included), OTU 1138, was affiliated with the Erwinia genus at 97% sequence identity, but we were unable to determine whether OTU 1138 represented Erwinia amylovora (the causal agent of fire blight). OTU 1138 was rare in the dataset and occurred 81 times over the time series, and peaked in abundance before flowers opened and before streptomycin was sprayed. However, there were no differences in OTU 1138 occurrence between trees sprayed with streptomycin and trees that were not sprayed (p = 0.71). Notably, close relatives of E. amylovora represent both plant pathogens and commensals. Some of these close relatives may be protective against E. amylovora (e.g., Pantoeavagans C9-1). We cannot distinguish these closely related, non-pathogenic species from OTU 1138.

Supplemental Discussion

Ecology of successional groups on apple flowers

The ecology of successional groups can be informed by the timing of each group’s peak occurrence, as well as by the identity of prevalent members. For example, Late successional members peaked in abundance just prior to petal fall. This group included a high abundance of Lactobacillus and Acetobacter taxa, whose occurrences align with conditions of flower decomposition by yeast (20-22). Additionally, the Climax group had highest evenness and phylogenetic diversity compared to the other groups, suggesting that resources were abundant and competition was relatively low at petal fall.

The Generalist taxa persisted on the flowers, but because they never were the most abundant, these taxa likely were not as competitive or as fast-growing as members of the other successional groups. However, these relatively stable, persistently detected Generalist members may inform as to the core microbiome of apple flowers at its most conservative definition (23). These Generalists may be targeted for follow-up studies, as their stability and relative abundance indicates that they may maintain intimate relationships with the host, with potential utility for crop management.

Caveats: Why didn't we observe an effect of streptomycin?