Assembly of Microbial Communities in Replicate Nutrient-Cycling Model Ecosystems Follows Divergent Trajectories, Leading to Alternate Stable States
EULYN PAGALING, KRISTIN VASSILEVA, CATHERINE G. MILLS, TIMOTHY BUSH, RICHARD A. BLYTHE, JANA SCHWARZ-LINEK, FIONA STRATHDEE, ROSALIND J. ALLEN, ANDREW FREE
Supplementary Methods
Illumina Sequencing and Sequence Data Analysis
The 250 bp paired-end Illumina reads were paired and analysed using QIIME 1.8.0 (Caporaso et al, 2010), a minimum quality score of 20 and an overlap of 200 bp for maximum accuracy (Kozich et al, 2013). Paired reads of <250 bp were removed from the dataset using the Python script filter_short_reads.py (https://gist.github.com/walterst/7602058) and chimeric sequences were removed using UCHIME (Edgar et al, 2011) with version 13_8 of Greengenes (DeSantis et al, 2006) as the reference database. A total of 3,816,084 (multiple replicates experiment) and 5,578,717 (system size experiment) quality-controlled, non-chimeric sequences which remained were analysed in QIIME 1.8.0. De novo 97% operational taxonomic units (OTUs) were generated by clustering using Uclust, and the taxonomy of the representative sequence of each OTU was assigned using Uclust against Greengenes 13_8. Sequences that failed to align with the rest of the dataset using PyNast were excluded from the resulting OTU table; singleton OTUs were also removed. Relative abundance plots of assigned taxonomy and analysis of alpha diversity by rarefaction on the two sample sets were generated using QIIME. Welch’s unequal variances t-test was used to test the significance of the differences in alpha diversity.
Co-occurrence Network Analysis
For co-occurrence network analysis, “network nodes” were defined as OTUs picked at 97% similarity from the 16-week sample data, which had a non-zero population in half the samples or more. This stringent approach greatly reduced network complexity. The resulting species set contained 860 OTUs from a total of 16,159, preserving 94.5% of the total abundance. The Spearman rank correlation coefficient between each pair of OTUs was computed with the Scipy Stats package Spearmanr in Python, based on the list of relative abundances of each OTU across multiple samples, such that a pair of OTUs with high abundances in the same samples, and low abundances in the same samples would result in a positive coefficient. A “link” (or potential interaction) was assigned between a given pair of nodes (or OTUs) if the magnitude of the Spearman coefficient ρ was both greater than 0.55 and had a p value less than 0.01. Once a link had been assigned, its “weight” was assumed to be equal to ρ (which could be positive or negative). Including only those nodes with a non-zero 10 15 20 25 30 35 40
number of links resulted in a network containing 859 OTUs, or 94.4% of the total abundance.
The network was visualised in an adjacency plot. In this plot, nodes are listed along the horizontal and vertical axes. Links between two nodes are represented by a square at intersections of the two nodes, the colour of which denotes the weight of the link. Since our analysis produces a non-directed network in which links are symmetric, the adjacency plot is symmetric along the diagonal. The relative ordering of nodes along the axes of an adjacency plot is in principle arbitrary, but can be optimised to allow clusters of interactions between nodes to be visualised. To identify clusters of interactions between taxa in our data, we used a sorting algorithm to set the order of nodes along the axes. To carry out this algorithm, a measure of “disorder” of the network was used, which was defined as the weight of each positive edge, multiplied by its distance to the diagonal, summed over all positive edges, summed over all rows. The disorder would therefore be at a minimum if high positive weight edges were positioned close to the diagonal, which therefore meant that nodes with strong positive links with one another would be positioned adjacent to one another. To attain this optimum configuration, nodes were selected at random and then their positions were swapped if the act of swapping would result in the disorder across the network decreasing. Networks were explored with Networkx (Hagberg et al, 2008) and visualised with Matplotlib. Clusters were defined by eye.
Since only OTUs with non-zero abundance in half the samples or more were considered for Network Analysis, we applied the same filter to the multiple replicate microcosm data (but excluded the sediment sample, as many of its abundant OTUs will only occur in that single sample). The resulting NMDS plot showed that the patterns of microcosm similarity were preserved (data not shown).
Replicate DNA Extractions
To test the variability between replicate DNA extractions from homogenized source sediment, technical replicate samplings, DNA extractions, PCR amplifications and sequencing were conducted. The resulting microbial community compositions (determined by Illumina sequencing) showed a Bray-Curtis similarity 67.2. This is 45 50 55 60 65 70
replicates experiment (41.3 ± 18.0 across all the multiple replicate microcosms, 58.9 ± 13.7 across the 8-week microcosms and 48.3 ± 14.1 across the 16-week microcosms). We have also shown previously that replicate DNA extractions from mature
microcosms give indistinguishable DGGE fingerprints (Pagaling et al., 2014). 75
Supplementary Figure Legends
Figure S1. Mean redox potential (± standard deviation) measured in the top (surface water), middle (sediment-water interface) and bottom (base of sediment) layers of replicate (n=10) microcosms incubated for 66 weeks (A) measured using a platinum wire against an Ag/AgCl reference electrode. Bacterial community profiles obtained by DGGE of 16S rRNA gene V3 regions from replicate microcosms in the set analysed in panel A fortnightly from weeks 2-16 and after 53 weeks of incubation (B). NMDS plot of the bacterial community profiles from the DGGE in panel B. Solid arrows indicate trajectories of bacterial community profiles from week 2 to week 16. Dashed arrows indicate trajectories of bacterial community profiles from week 16 to 53 (C).
Figure S2. Bacterial community profiles obtained by DGGE from microcosms destructively sampled at 8 and 16 weeks (A). Fingerprints from multiple gels were aligned in Bionumerics 6.0 using internal marker lanes, and the normalized, aligned banding patterns used for band matching by this software are presented in the figure. For the 16-week microcosms, the individual fingerprints are colour-coded according to the clusters observed in panel B. NMDS plot of the bacterial community profiles obtained by DGGE (B). The communities at 8 weeks (diamonds) and 16 weeks (squares) are indicated. Two sub-groups and an outlier within the 16-week microcosms identified by eye are indicated by ovals and differential colouring. NMDS plot based on the DGGE fingerprints (identical to panel B) but coloured according to the sub-groups identified by sequencing (C). Microcosms whose grouping was not determined by sequencing are represented by yellow circles.
Figure S3. Phylum distributions for multiple replicate microcosms (A) and different system size microcosms (B) represented as fractional abundances.
Figure S4. Rarefaction curves of the sequencing data for the inoculating sediment and microcosms sampled at 8 and 16 weeks. Rarefaction curves were calculated for each sub-group; error bars represent standard deviations. Figures for the corresponding Shannon diversity (± SD) at the maximum rarefaction depth are shown adjacent to the 80 85 90 95 100 105 110
Figure S5. Network analysis of the bacterial community in the 16-week microcosms at the Class level. Two clusters of positive associations (black/blue clusters) and negative associations (pink/red clusters) can be identified. Sequences that are less than 1% abundant are listed in Table S4.
Figure S6. Heatmap of the abundances of OTUs corresponding to predicted sulphate-reducing bacteria (SRB) in the inoculating sediment (S) and the sub-groups of 8-week and 16-week microcosms. The five most abundant OTUs corresponding to SRB were selected from each sample, combined into an abundance table and then sorted by total abundance in the 8-1 sub-group of microcosms.
Figure S7. Heatmap of OTU abundances in the inoculating sediment (S) and mature microcosms of different sizes. The five most abundant OTUs were selected from each sample, combined into an abundance table and then sorted by total abundance in the 1000 ml microcosms.
Figure S8. Variation in the mean Bray-Curtis similarity values between microcosm communities (n=3) over different system sizes. Error bars represent standard deviations.
Figure S9. Centroid-adjusted NMDS plot of the similarity between bacterial communities in microcosms of different system sizes.
Figure S10. Photographs of the mature microcosms developed in glass bottles during the system size experiment.
115
120
125
130
References
Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. (2010). QIIME allows analysis of high-throughput community sequencing data. Nat Methods 7:335-336.
DeSantis TZ, Hugenholtz P, Larsen N, Rojas M, Brodie EL, Keller K, et al. (2006). Greengenes, a chimera-checked 16S rRNA gene database and workbench compatible with ARB. Appl Environ Microbiol 72:5069-5072.
Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. (2011). UCHIME improves sensitivity and speed of chimera detection. Bioinformatics 27:2194-2200.
Hagberg, AA, Schult DA, Swart PJ. (2008). Exploring network structure, dynamics, and function using NetworkX. Proceedings of the 7th Python in Science Conference (SciPy 2008), Pasadena, CA, USA.
Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. (2013). Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq Illumina sequencing platform. Appl Environ Microbiol 79:5112-5120.
Pagaling, E., Strathdee, F., Spears, B.M., Cates, M.E., Allen, R.J., and Free, A. (2014). Community history affects the predictability of microbial ecosystem development. ISME J 8:19-30. 140 145 150 155 160 165
