Edinburgh Explorer Assembly of Microbial Communities in Replicate Nutrient-Cycling Model Ecosystems Follows Divergent Trajectories, Leading to Alternate Stable States

Summary We studied in detail the reproducibility of community development in replicate nutrient-cycling microbial microcosms that were set up identically and allowed to develop under the same environmental conditions. Multiple replicate closed microcosms were con-structed using pond sediment and water, enriched with cellulose and sulphate, and allowed to develop over several months under constant environmental conditions, after which their microbial communities were characterized using 16S rRNA gene sequencing. Our results show that initially similar microbial communities can follow alternative – yet stable – trajectories, diverging in time in a system size-dependent manner. The divergence between replicate communities increased in time and decreased with larger system size. In particular, notable differences emerged in the heterotrophic degrader communities in our microcosms; one group of steady state communities was enriched with Firmicutes , while the other was enriched with Bacteroidetes . The communities dominated by these two phyla also contained distinct populations of sulphate-reducing bacteria. This biomodality in community composition appeared to arise during recovery from a low-diversity state that followed initial cellulose degradation and sulphate reduction.


Introduction
Microbial communities mediate many of the key steps in global carbon, sulphur and nitrogen cycling (Falkowski et al., 2008), and are exploited by humans for applications including waste decomposition (Martins et al., 2004) and the production of bioenergy (Zhou et al., 2013). Microbial gut communities are also essential to human and animal health (Walter and Ley, 2011;Hanning and Diaz-Sanchez, 2015). However, the processes that drive the establishment, maintenance and function of microbial communities remain poorly understood Widder et al., 2016). In particular, it remains unclear whether microbial ecosystem development is inherently predictable, being controlled by deterministic factors, or unpredictable, being controlled by stochastic factors (Langenheder and Szekely, 2011;Pagaling et al., 2014). Establishing a clear picture of the factors controlling predictability in microbial ecosystem development would have significant environmental, economic and medical implications, including improved reliability of wastewater treatment, improved energy production or yield from bioreactors and better treatment of human gut conditions (Ju and Zhang, 2015;Widder et al., 2016).
Deterministic factors influencing microbial community development include selection by environmental conditions and species interactions. These factors are expected to act as an ecological filter, imposing 'species sorting' and driving microbial community composition towards a predictable stable state (Martiny et al., 2006;Langenheder and Szekely, 2011). However, ecosystem dynamics can also be influenced by stochastic processes, such as immigration and emigration events, and random births and deaths (Volkov et al., 2003;Khatri et al., 2012;Pagaling et al., 2014), which lead to variability in community development. The relative importance of stochastic and deterministic factors in shaping community composition can be dependent on time (Zhou et al., 2013;Zhou et al., 2014), space (Stegen et al., 2012) and initial conditions (Pagaling et al., 2014). This is further complicated by the fact that nonlinearities in microbial population dynamics can lead to chaotic ecosystem development, even in the absence of stochastic factors (Beninca et al., 2008;Beninca et al., 2009;Graham et al., 2007;Bush et al., 2015), although this is not always the case (Kato et al., 2005;Kato et al., 2008;Falk et al., 2009).
While in many cases, community assembly leads to a single stable compositional state, in other cases, multiple stable states may be possible, arising from feedbacks between species and/or species-environmental interactions (Scheffer et al., 2001). In these cases, stochastic factors can determine which of the possible states is observed. Alternative stable states have been observed in populations of macroscopic organisms in various environments, including southern boreal forest communities in Canada (Jasinski and Payette, 2005), forest and savanna communities in Africa (Staver et al., 2011), eutrophic lakes (Blindow et al., 1993;Scheffer et al., 1993;Ibelings et al., 2007), marine communities (Petraitis and Dudgeon, 2004) and pond mesocosms (Chase, 2003). In contrast, fewer examples of alternative stable states in microbial communities have been reported, possibly due to the difficulty in defining a stable state (as opposed to a transient state) in microbial populations (Allison and Martiny, 2008). However, it has been suggested that alternative microbial community states may occur after a stochastic event that affects the environmental conditions, for example in the human gut following small bowel transplantation (Hartman et al., 2009), in the mucus layers of coral following increased temperature (Mao-Jones et al., 2010), or in aquatic ecosystems in pitcher plants (Sirota et al., 2013). In bioreactors, alternative stable states may also arise stochastically after reactor start-up, for example in a microbial electrolysis cell (Zhou et al., 2013), a batch reactor for nitrogen removal from urine (Burgmann et al., 2011) and in wastewater treatment processes (Martins et al., 2004).
Microbial microcosm communities constitute an important resource for the testing and development of ecological and evolutionary theory (Horner-Devine et al., 2004;Jessup et al., 2004;Prosser et al., 2007), as they develop rapidly and are easy to manipulate in the lab (Widder et al., 2016). This makes it possible to perform replicate experiments on a scale that would be impossible for animals or plants. Here, we exploited these attributes to make a direct test of the reproducibility of microbial ecosystem development, using nutrient-cycling microcosms (Winogradsky columns) made from freshwater pond sediment and water (Rundell et al., 2014;Abbasian et al., 2015). The microbiology of the Winogradsky column is wellunderstood, with the key stages in its initial development being driven by anaerobic cellulose degradation and subsequent sulphate reduction (Rogan et al., 2005). In the presence of sulphate, much of the flux of fermentation products in the sediment is diverted from methanogenesis to sulphate reduction, which is more energetically favourable (Muyzer and Stams, 2008). This leads to a very low redox potential in the microcosms after 8-9 weeks of development, and a reduction in species diversity compared to the inoculating sediment (Pagaling et al., 2014). Thereafter, the redox potential increases in the upper layers of the microcosm as the produced sulphide is reoxidized by anaerobic sulphide oxidizers and phototrophic sulphur bacteria (Supporting Information Fig. S1). Dependent on the levels of sulphide present, oxygenic photosynthesis by cyanobacteria and eukaryotic algae may also emerge in the upper layers, and spatial environmental variation within and between microcosms can occur (Pagaling et al., 2014;Esteban et al., 2015;Widder et al., 2016). The different subcommunities of microorganisms which carry out these diverse functions interact with each other and with the abiotic environment to generate an emergent microcosm system; we refer to the total assemblage of organisms in the system as the 'microcosm community'.
In our previous work on this system (Pagaling et al., 2014), we showed that triplicate microcosms seeded with the same source sediment developed different microbial communities, although a signal of the source (or history) of the transplanted microbial communities was retained in the final community composition. In this study, we analysed the microbial communities in many tens of initially similar, parallel microcosms using 16S rRNA gene V4 region amplicon sequencing. We observe that replicate experiments result in the development of two alternative community groups at long times, one characterized by elevated levels of Firmicutes and the other characterized by elevated levels of Bacteroidetes. Additionally, we observe that the sulphate-reducing communities also vary between these two groups of microcosms, while phototrophs do not exhibit the same co-variation. Divergence between replicate communities increases in time and is dependent on the system size. Our results suggest that divergence of ecosystem trajectories into possible alternative stable states occurs following a low-diversity state associated with the redox potential minimum.

Community development in multiple replicate microcosms
To investigate variability in microbial community development, 100 replicate sediment-water microcosms (Winogradsky columns) were set up identically. Previous data from multiple experiments in our lab and elsewhere has shown that development of mature Winogradsky column microcosm communities requires 3-4 months, with very few changes in community composition and microcosm function occurring after 12 weeks, while redox potential at all levels in the microcosm declines considerably in the first 10 weeks as cellulose is degraded, then subsequently recovers in the upper water layer while remaining low in the sediment (Pagaling et al., 2014;Esteban et al., 2015;Supporting Information Fig. S1A). To compare variability during the early and late stages of development in our 100 replicate microcosms, 50 microcosms were destructively sampled after 8 weeks, and the remaining microcosms were destructively sampled after 16 weeks. The bacterial communities from the microcosms at these two time points were visualized by denaturing gradient gel electrophoresis (DGGE) of 16S rRNA genes (Supporting Information Fig.  S2A), which revealed that the bacterial communities were very similar in composition at 8 weeks, but had diverged after 16 weeks.
An NMDS plot of the bacterial communities in the 8-and 16-week replicate microcosms, as determined by DGGE, showed that the microcosm communities clustered according to the time of sampling (Supporting Information Fig  S2B). PERMDISP analysis showed that the 16-week microcosm communities were significantly more dispersed than the 8-week microcosm communities (F 5 37.77, p 5 0.001), while PERMANOVA analysis indicated that the 8-and 16-week communities also differed significantly in composition (Pseudo-F 5 26.75, p 5 0.001). Interestingly, the bacterial communities in the 16-week microcosms appeared to cluster into two subgroups: one group was more similar to the 8-week microcosm communities (designated subgroup 16-1) and the other was less similar to the 8-week microcosm communities (designated subgroup 16-2) (Supporting Information Fig. S2B). These two subgroups were significantly different in composition (PERMANOVA Pseudo-F 5 3.68, p 5 0.001), but not in dispersion (PERMDISP F 5 3.16, p 5 0.126). Significant divergence of microbial communities in triplicate Winogradsky column microcosms made from the same source sediment has been demonstrated previously (Pagaling et al., 2014). The larger number of replicates used in this study provides stronger evidence that starting communities taken from the same, homogenized source can develop differently under identical external environmental conditions. The existence of two distinct groups of 16-week microcosms also hints at the possible existence of alternative stable states.
Microbial species present in the source sediment will be subject to strong selection in the microcosm environment due to the addition of organic matter (cellulose) and sulphate and incubation under non-native conditions (constant illumination at 258C). To analyse the change in species composition during this selection, the source sediment plus 18 of the 8-week microcosms and 20 of the 16-week microcosms selected at random were subject to metataxonomic analysis by Illumina sequencing of 16S rRNA genes (Supporting Information Fig. S3A). In comparison to the source sediment, which was dominated by the phyla Proteobacteria (57%) and Bacteroidetes (22%), the 8-week microcosms were dominated by Firmicutes (79 6 7%) in all but 2 samples (all values mean 6 standard deviation across the replicate groups). After 16 weeks the composition had again changed strikingly, with a resurgence in Proteobacteria (27 6 9%) and Bacteroidetes (10 6 7%), the emergence of Spirochaetes (10 6 8%) and variable persistence of Firmicutes (35 6 19%). These observations are similar to those of previous studies on Winogradsky column ecology which showed that Proteobacteria, Bacteroidetes and Firmicutes were the dominant phyla (Rundell et al., 2014;Abbasian et al., 2015), though the overall community composition was source sedimentdependent (Pagaling et al., 2014;Rundell et al., 2014).
NMDS plots of community composition based on Bray-Curtis similarity at the OTU level illustrated the drastic change in composition between inoculating sediment and microcosms, separate clustering of the 8-week and 16week microcosms, and the distinct composition of the two low-Firmicutes 8-week microcosms (Fig. 1). These results did not depend on our specific choice of similarity measure: similar results were obtained using a Jaccard presence-absence index (data not shown). We observed a significantly larger spread in community composition for the 16-week data points compared to the 8-week data points (PERMDISP: F 5 10.10, p 5 0.01), as in our DGGE analysis. Analysis by PERMANOVA also showed that the 8-and 16-week communities differed significantly in composition (Pseudo-F 5 26.75, p 5 0.001), in agreement with our DGGE analysis. Also in agreement with the fingerprint data, the sequencing data showed that our 16-week microcosm communities fell into two distinct subgroups ( Fig. 1). Super-position of the group classifications from the sequencing data onto the DGGE-based ordination showed that the subgroups were consistent between the two methods (Supporting Information Fig. S2C). The two out-lying 8-week data points (designated subgroup 8-2) fell between the two 16-week subgroups, suggesting that their microbial communities were intermediate in composition (Fig. 1).
Rarefaction analysis showed that there was a drastic loss of species diversity in the 8-week microcosms compared to the starting sediment (Supporting Information Fig.  S4). Analysis based on the Shannon diversity index showed that the richness and evenness of the two subgroups of 16-week microcosms were higher than those of the 8-week microcosms, indicating that there was some recovery of diversity in the later stages of community development (Supporting Information Fig. S4). Welch's unequal variances t-test showed that this difference was significant (8-1 vs. 16-1: F 5 11.48, p 5 0.004; 8-1 vs. 16-2: F 5 46.42, p < 0.001). In contrast to the majority of 8-week microcosms, the rarefaction curves for the two 8week outlier microcosm communities showed that they contained the highest microbial diversity of all, supporting the hypothesis that the communities in these outliers are distinct from the subgroups 8-1, 16-1 and 16-2.

Potential heterotrophic degraders drive community variation
The drastic reduction in species diversity in most 8-week microcosms, followed by a partial recovery by 16 weeks, suggests strong selection for certain highly-specialized functions in the initial stages of microcosm development, with later limited functional re-diversification. The key selective factors in our microcosms are likely to be the high levels of cellulose and sulphate added uniformly throughout the sediment. Of the major phyla identified within our microcosms (Supporting Information Fig. S3A), we expect to find the ability to degrade cellulose and its products within the Firmicutes, Bacteroidetes and Fibrobacteres; members of the Actinobacteria, which have also been shown to degrade cellulose (Lynd et al., 2002), were not abundant in the microcosms. Members of the Spirochaetes have been found to associate with Bacteroidetes to enhance cellulose degradation in the rumen (Kudo et al., 1987), and with the Firmicutes species Clostridium thermocellum for cellulose degradation in co-culture (Pohlschroeder et al., 1994) due to their fermentation of the mono-and di-saccharide products of cellulose breakdown. Moreover, Spirochaetes of the genus Treponema were selected in other studies of Winogradsky columns supplemented with cellulose (Esteban et al., 2015). We therefore examined the temporal variation of OTUs mapping to the four phyla Firmicutes, Bacteroidetes, Spirochaetes and Fibrobacteres during microcosm development.
The five most abundant OTUs in each of our 39 samples belonging to any of the phyla Firmicutes, Bacteroidetes, Spirochaetes and Fibrobacteres were extracted from our sequence data set, giving 24 OTUs in total since some of these OTUs were abundant in several samples. Emphasizing the strong selection for individual OTUs in the microcosms, these 24 OTUs accounted for 69-86% of the total abundance of these phyla in each microcosm sample, and 45% of their abundance in the source sediment. A heat-map was used to monitor the abundance of these OTUs in the source sediment and microcosm samples (Fig. 2). Those OTUs which were abundant in the source sediment (mostly Bacteroidetes) were rare in the 8-week microcosms, although some reappeared in 16-week microcosms. The 8-1 group of microcosms was dominated by a group of 10 OTUs from the phylum Firmicutes, suggesting a key role in the initial breakdown of cellulose. In the 8-2 outlier microcosms, these OTUs were less abundant, while two OTUs assigned to Treponema (phylum Spirochaetes; denovo50143) and Fibrobacteres (denovo6269) were dominant. The same two Treponema OTUs were abundant in the 16-week microcosms (especially group 16-2), which also exhibited increased abundance of a number of Bacteroidetes OTUs and one from the Firmicutes (denovo16758). Firmicutes OTUs were more abundant in group 16-1 than group 16-2 microcosms, which were conversely more dominated by Bacteroidetes (Fig. 2). This suggests a changing selection for different potential heterotrophic degraders during microcosm development, and a different final population of such organisms in the two subgroups of 16-week microcosms.

Biomarkers of the 16-week sample groups
To identify particular microbial taxa characteristic of the sediment and each of the subgroups of microcosms, a metagenomic biomarker discovery tool, LEfSe (Segata et al., 2011) was used to test the differential abundance of all OTUs from our four selected phyla across the sample groups. This analysis showed that the differences between the sediment and microcosm subgroups could be explained largely by differences in abundance of the OTUs examined above (Table 1). Of the 19 OTUs with LDA scores of >4.0, only two, denovo45589 (Firmicutes) and denovo43414 (Spirochaetes) were not represented in the set of 24 most abundant OTUs from these phyla.
LEfSe was also used to test for phyla with differential abundances between the two 16-week subgroups of microcosms, considering all phyla (Table 1). Abundances of the Firmicutes (group 16-1) and Spirochaetes and Bacteroidetes (group 16-2) were strongly discriminatory (LDA scores between 4.5 and 5.0; p values between 0.0002 and 0.0003), while Chloroflexi and Proteobacteria were also biomarkers for these groups, but to a lesser extent. This emphasizes that the phyla which we have identified as potential heterotrophic degraders drive the variation between these groups of microcosms. When the abundances of these heterotrophic degraders were plotted as vectors on a PCoA plot, Firmicutes were again found to correlate with subgroup 16-1, while Bacteroidetes, Spirochaetes, Verrucomicrobia and Proteobacteria were found to correlate strongly with subgroup 16-2 (Fig. 3). However, the vectors showed that the Bacteroidetes and Firmicutes varied strongly along the major PCO1 axis (correlation coefficients are 20.942 and 0.973 respectively), and therefore, accounted for most of the variation between the two 16-week subgroups.
To identify patterns of co-occurrence between taxa in our microcosm communities, we used co-occurrence network analysis (Barberan et al., 2012;Berry and Widder, 2014;Widder et al., 2014). Although not conclusive evidence of a direct interaction between two taxa, (which may show apparent correlations for indirect reasons even if they are spatially separated within the microcosms), such analysis can assist in identifying possible direct or indirect interactions. Applying this analysis to the bacterial communities for our 16-week microcosms revealed the presence of two distinct clusters of OTUs, which we denoted cluster 1 (purple box in Fig. 4A) and cluster 2 (green box in Fig.  4A). OTUs within the same cluster tended to show similar occurrence patterns across microcosm samples, while OTUs in different clusters tended to show different occurrence patterns. To analyse whether these two clusters of OTUs were characteristic of our two distinct subgroups of 16-week microcosms, the fractional abundance of OTUs that fell into clusters 1 and 2 were computed for each microcosm community (Fig. 4B). OTUs within cluster 1 were enriched in microcosms in subgroup 16-1 (purple bars in Fig. 4B), while OTUs within cluster 2 were enriched in microcosms in subgroup 16-2 (green bars in Fig. 4B). Repeating the network analysis at the Class level showed clearly that cluster 2 contained Bacteroidetes, Spirochaetes and Verrucomicrobia (Supporting Information Fig. S5), while cluster 1 contained Firmicutes. Interestingly, this analysis also revealed a new association; cluster 1 also contained Anaerolineae (within the phylum Chloroflexi), suggesting that they were co-selected in our 16-week microcosms with the Firmicutes (Supporting Information Fig. S5). Verrucomicrobia isolates have been found to contain genes that are associated with cellulose degradation, though this has not been demonstrated in vitro (Martinez-Garcia et al., 2012;Wertz et al., 2012). Relatively little is known about the Anaerolineae; however, a cultivated member of the genus Anaerolinea (A. thermolimosa) is a strict anaerobe that can utilize a number of sugars as substrates for growth (Yamada et al., 2006). It is possible that it plays a similar role to Treponema, i.e., it uses cellobiose and other cellulose degradation intermediates as substrates.
All the microcosms in this study were set up identically and maintained under identical conditions, without immigration. Thus, the split into two distinct groups observed in the 16-week microcosms cannot be driven by external factors. The sequencing analysis does not point to a specific mechanism driving microcosms down one of these alternative developmental trajectories, although it is intriguing that prior to the split (i.e., at 8 weeks), we observe strong, consistent selection of a few Firmicutes taxa. A relevant factor that could be driving the community splitting is species interactions within the microcosm communities (e.g., Lima-Mendez et al., 2015). Moreover, minor differences in the relative abundance of key Firmicutes and/or Bacteroidetes taxa, which are rare in the starting inoculum (Fig. 2), may become amplified during development and result in major differences in the mature microcosms. It is also very likely that coupling between the species composition of the emergent microbial community and its chemical environment, which is driven by cellulose degradation, sulphate reduction and subsequent phototrophic nutrient cycling, plays an important role.

Other functional groups in the microcosm communities
As sulphur cycling and phototrophy are key processes in the Winogradsky column, it is of interest to examine whether microorganisms which are predicted to contribute to these functions also vary among replicate microcosms. Sulphate reduction can be relatively reliably assigned to particular known bacterial groups based on phylogeny (Martiny et al., 2013;Muyzer and Stams, 2008). In our microcosms, predicted sulphate-reducing bacteria (SRB) were dominated by the genera Sulfurospirillum, Desulfovibrio and Desulfolobus, while Geobacter and Desulfolobus species were present in the sediment inoculum but at lower levels (Supporting Information Fig. S6). Again, strong selection of individual OTUs was evident. LEfSe analysis showed that a Sulfurospirillum species was elevated in group 16-1 (LDA score 5 5.12, p 5 0.0109), while particular Desulfovibrio, Desulfobulbus, Desulfobulbaceae Desulfobacteraceae and Desulfarculaceae OTUs were elevated in 16-2 (Desulfovibrio LDA score 5 4.63, p 5 0.0168; see Supporting Information Table S2). This phenomenon was only observed at the OTU level (not at class level), suggesting that this variation is fine-scale. When the similarity between predicted SRB OTUs in different microcosms was visualized by NMDS (Fig. 5), the community composition was found to be significantly different in the two 16-week microcosm subgroups 16-1 and 16-2 (PERMANOVA Pseudo-F 5 3.33, p 5 0.008). It is reasonable to suppose that alternative states of the cellulose degrader community might also result in alternative sulphate reducer communities as SRB rely on degraders to carry out the initial breakdown of complex biopolymers (Boyer, 1986;Muyzer and Stams, 2008). Anaerobic cellulose degradation produces various organic acids (Leschine, 1995), which can be used as substrates for sulphate reducing bacteria (Muyzer and Stams, 2008). It is possible that differences in the levels of these various organic compounds produced by different groups of degraders are influencing the sulphate reducing community composition in the microcosms. Alternatively, specific variations in the activities of different SRB could feedback onto the composition of the degrader community at later stages of microcosm development, e.g., via sulphide toxicity. Variability between replicate microcosms in the abundance of phototrophic sulphur bacteria (corresponding to blooms of the family Chlorobiaceae) and sulphur oxidizers (corresponding to blooms of the epsilonproteobacterium Sulfuricurvum) was observed in the 16-week microcosm communities (Supporting Information Fig. S3A), but these variations were independent of the subgroups of microcosms driven by differences in the heterotroph and SRB communities. In each case a single, highly-abundant OTU which was rare in the sediment inoculum was responsible. These blooms probably resulted from the production of sulphide by the SRB described above in the anaerobic environment created by the degradation of cellulose [Chlorobi are obligate anaerobes (Bryant and Frigaard, 2006) and Sulfuricurvum can only tolerate micro-aerobic conditions (Kodama and Watanabe, 2004), possibly growing at the interface between sulphide and oxygen]. Sulfuricurvum kujiense is the only cultivated member of this genus, which was isolated from groundwater (Kodama and Watanabe, 2004), so like the Chlorobiaceae, tends to be found in aquatic environments. Due to their exclusively aquatic habitat, sporadic blooms of these microorganisms may contribute to their patchy distribution in the microcosms.

Effects of system size on microcosm variability
We hypothesized that the low levels of community diversity which we observe in our 8 week microcosms might be associated with the subsequent divergence in community composition between replicate microcosms. Stochastic elimination of certain low-abundance taxa during the lowdiversity state at 8 weeks could influence the subsequent community trajectory, resulting in divergent communities at 16 weeks. To test this idea, we varied the system size as a method to manipulate population size at the 8 week 'bottleneck' (Wertz et al., 2007). Triplicate microcosms of different volumes ranging from 5 ml to 1 l were set up and incubated for 21 weeks. Illumina sequencing of these mature microcosms showed that the most dominant group was Chlorobi (12%), followed by Chloroflexi (10%), Firmicutes (9.9%), Proteobacteria (9.6%) and Bacteroidetes (4.6%) (Supporting Information Fig. S3B) -phyla that were similarly found in the multiple replicate microcosms. However, these microcosms were found to contain larger blooms of anoxygenic phototrophs and did not contain such large populations of Spirochaetes as the multiple replicate microcosms. Again, individual OTUs which were rare in the sediment inoculum were strongly selected in the microcosms, with a preponderance of likely cellulose degraders, sulfate reducers, sulfide oxidizers and anoxygenic phototrophs (Supporting Information Fig. S7).
To assess how variability between replicate microcosms changed with system size, the mean Bray-Curtis similarity values between fully-developed (21 week) triplicate microcosms were calculated and plotted against microcosm size (Supporting Information Fig. S8). Strikingly, smaller microcosms were more variable; we observed a sharp increase in the average Bray-Curtis similarity between replicates for microcosms larger than 100 ml. We further checked whether the variability in the 16-week multiple replicates experiment corresponded with the variability in the microcosms of a similar size in the system size experiment, and found that they were consistent [Bray-Curtis values are 43.2 6 15.2 (mean 6 SD) and (37.4 6 19.9) respectively].
An NMDS plot based on the similarity values showed a large change in composition between sediment and microcosm communities and clear groupings according to microcosm size (Fig. 6A). Different-sized communities were significantly different in composition (PERMANOVA Pseudo-F 5 2.51, p < 0.001); this may in large part be due to the change in the surface area-volume ratio as the size of the microcosm increases. To focus on community variability, the data were adjusted such that the points corresponding to each set of triplicate microcosms had a common centroid (Supporting Information Fig. S9) (Pagaling et al., 2014). This showed significantly higher divergence between replicates for the smaller microcosms (PERMDISP F 5 41.6, p < 0.001, comparing between 'small' microcosms with volumes less than 100 ml and 'large' microcosms of 100 ml and above; Supporting Information Fig. S8). Further analysis showed that the composition of and variability in the potential cellulosedegrading phyla also changed with microcosm size. A correlation between variability in the Firmicutes-Bacteroidetes ratio and the size of the microcosm was observed, with this ratio being more variable in smaller microcosms (Fig.  6B). Taken together, these data show that system size can be an important factor in microbial community variability, and that at smaller system sizes, microcosms can exist in states of extreme domination by either Bacteroidetes or Firmicutes, while a more consistent Firmicutes-rich state is observed for larger microcosms.

Conclusions
We measured the extent of divergence between trajectories of microbial community development in identically prepared multiple replicate Winogradsky column microcosms made from a single batch of homogenized freshwater pond sediment and water, and incubated under identical laboratory conditions. After 8 weeks' development, microcosm community diversity was very low compared to the inoculating sediment, and several OTUs belonging to the phylum Firmicutes, which were rare in that sediment had been strongly selected in the majority of microcosms. These Firmicutes species may have been involved in the degradation of the microcrystalline cellulose which was added uniformly to the sediment: in human gut ecosystems, Firmicutes are key early colonizers of insoluble fibres which may act as keystone species in their degrading communities (Berry, 2016). Two of the 8-week microcosms exhibited higher species diversity and a composition intermediate between the other 8-week microcosms and those which had developed for 16 weeks; it is possible that the communities in these microcosms were transitioning between the 8-week and 16week states.
After 16 weeks' development, diversity had partially recovered and microcosm communities had diverged into two apparent alternative stable states. These two community states differed strongly in their potential heterotroph communities, with one containing elevated abundances of Firmicutes OTUs (mostly shared with the 8-week microcosms) and the other containing elevated abundances of Bacteroidetes OTUs which had to some extent replaced the Firmicutes. Anaerolineae were co-selected with the Firmicutes, while Treponema (phylum Spirochaetes) and Verrucomicrobia were co-selected with the Bacteriodetes. It seems likely that the reorganization of the heterotroph community is a result of the exhaustion of the supplied cellulose and a re-diversification of heterotrophic metabolism in the microcosms in a stochastic manner. This stochasticity has apparent consequences for the SRB, which utilize organic acids from organic matter degradation.
We also observed that microcosm size played an important role in influencing community variability and composition, with smaller microcosms being more variable than larger ones and displaying more extreme variation in the ratio of Bacteroidetes to Firmicutes. As our microcosms were set up identically and maintained under identical, isolated conditions, the driver of this split in community composition and the consistent co-selection of certain phyla is likely to be an internal mechanism -such as stochastic variations in abundance, followed by nonlinear amplification through community and environmental feedbacks -which warrants further investigation. Divergence appears to occur following passage through a low-diversity bottleneck associated with a redox potential minimum at 8-9 weeks (Pagaling et al., 2014). Our observation that smaller microcosms are more variable suggests that they may be more susceptible to the stochastic elimination of rare taxa at the low diversity bottleneck, which may bias the subsequent recovery trajectory. This work has implications for other natural and engineered microbial communities as it suggests that divergence into alternative states (with possible consequences for function) may occur independently of external factors.

Experimental procedures
Sampling and microcosm set-up Sediment and water samples were collected from an urban, eutrophic freshwater pond, Blackford Pond, Edinburgh, Scotland (longitude 23.11 latitude 55.55). Sediment was taken from depths of <0.5 m with an ethanol-sterilized trowel; water samples were collected in sterile Duran bottles. All microcosms were set up following a standard Winogradsky Column recipe in which nutrients consisting of 2.5 g of microcrystalline cellulose, 5 g of CaSO 4 and 0.25 g of CaCO 3 were added per 100 g of sediment, which was then homogenized (Pagaling et al., 2014). One hundred replicate 15 ml microcosms were set up as follows: 7.5 g of homogenized sediment mixture were added to each 15 ml falcon tube and topped up with 5 ml pond water (so that sediment and water were in equal volumes), and incubated at 258C with 24 h Northlight illumination for 8 or 16 weeks. Microcosms were staggered in the incubator to avoid shading from adjacent tubes and the incubator was lined with a reflective surface to ensure even distribution of light to all the microcosms. No correlation between microbial community development and the position of the microcosms in the incubator has been observed (Pagaling et al., 2014;this study). To test the effects of system size, another 21 microcosms were set up with volumes ranging from 5 ml to 1000 ml (5, 20, 100, 130, 250, 500, 1000 ml), in triplicate. The dimensions of these microcosms are stated in Supporting Information Table S3 and a photo can be seen in Supporting Information Fig. S10. Sediment and water was added in equal volumes to each microcosm. These microcosms were incubated for 21 weeks.
Community DNA extraction, PCR amplifications, denaturing gradient gel electrophoresis (DGGE) and statistical analysis After incubation for 8 or 16 weeks (for the multiple replicates experiment), or 21 weeks (for the system size experiment), microcosms were destructively sampled by vortexing, 1 ml of the slurry was removed and DNA was extracted using an UltraClean Soil DNA Isolation Kit (MoBio). Vortexing prior to sampling removed spatial variability so that our sample was a representation of the total microcosm community. DNA was also extracted from the homogenized source sediment. Previous work has shown that variability between replicate DNA extractions from the same homogenized sediment (Supporting Information Methods) or mature microcosm (Pagaling et al., 2014) is much lower than that between individual microcosms. Community profiles were obtained from the microcosms by DGGE as previously described (Pagaling et al., 2014). Normalization of the DGGE gels and band-matching analysis was done in Bionumerics version 7.0 (Applied Maths, Sint Martens-Latem, Belgium). These data were exported and used for statistical analysis in Primer-E version 6.1.12 (Primer-E, Ivybridge, UK) to calculate Bray-Curtis similarity matrices, non-parametric multi-dimensional scaling (NMDS) plots, and PERMDISP and PERMANOVA tests.

PCR amplification conditions for illumina sequencing
Variable (V4) regions of bacterial 16S rRNA genes were amplified for Illumina sequencing using the barcoded 515F-806R reverse primer series (Caporaso et al., 2012). PCR reactions contained 1 x Taq buffer plus additional MgCl 2 (final concentration 2.5 mM), 0.2 mM of each dNTP, 0.25 lM of each primer, 1 ll of template DNA and 1.25 U of Taq DNA polymerase (Roche Applied Science), with the volume made up to 25 ll with PCR grade water (Sigma). Thermal cycling conditions were: 948C for 3 min; followed by 35 cycles of 948C for 45 s, 508C for 1 min, 728C for 1.5 min; followed by 728C for 10 min. The efficiency of all PCRs was tested by electrophoresis of the products on 0.8% (w/v) TAE-agarose gels. PCR products were quantified using a Quant-iT PicoGreen ds DNA Assay Kit (Life Technologies) prior to pooling in equimolar quantities for Illumina sequencing.

Illumina sequencing and sequence data analysis
Pooled, barcoded amplicons were sequenced as part of a single run of an Illumina MiSeq v2 sequencer, yielding 250 bp paired-end reads. The reads were paired and analysed using QIIME 1.8.0 (Caporaso et al., 2010) as described in Supporting Information Methods, yielding 3.8 million (multiple replicates experiment) and 5.6 million (system size experiment) quality-controlled reads. To pick out taxa and predicted gene classes that were significantly elevated in different biological subgroups, OTU and phylum tables were uploaded to the LEfSE Galaxy server (http://huttenhower. sph.harvard.edu/galaxy), where the differential abundances were found and tested for biological significance (Segata et al., 2011). The output files of LDA scores and p values were exported and viewed in Excel to construct tables of differentially abundant features.

Co-occurrence network analysis
Co-occurrence network analysis was used to identify potential interactions between different taxa within our multiple replicate microcosm communities using OTUs picked at 97% identity from the 16-week sample data (Barberan et al., 2012;Berry and Widder, 2014;Widder et al., 2014). Full details are given in Supporting Information Methods.

Accession numbers
Raw sequence data were deposited into the NCBI Sequence Read Archive (SRA) under accession number SRP063557.

Supporting information
Additional Supporting Information may be found in the online version of this article at the publisher's web-site: Fig. S1. Mean redox potential (6 standard deviation) measured in the top (surface water), middle (sediment-water interface) and bottom (base of sediment) layers of replicate (n 5 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 to 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). . 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 subgroups 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 subgroups identified by sequencing (C). Microcosms whose grouping was not determined by sequencing are represented by yellow circles. Fig. S3. Phylum distributions for multiple replicate microcosms (A) and different system size microcosms (B) represented as fractional abundances. Fig. 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 subgroup; error bars represent standard deviations. Figures for the corresponding Shannon diversity (6 SD) at the maximum rarefaction depth are shown adjacent to the curves. Fig. 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 Supporting Information Table S4. Fig. S6. Heatmap of the abundances of OTUs corresponding to predicted sulphate-reducing bacteria (SRB) in the inoculating sediment (S) and the subgroups 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 subgroup of microcosms. Fig. 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. Fig. S8. Variation in the mean Bray-Curtis similarity values between microcosm communities (n 5 3) over different system sizes. Error bars represent standard deviations. Fig. S9. Centroid-adjusted NMDS plot of the similarity between bacterial communities in microcosms of different system sizes. Fig. S10. Photographs of the mature microcosms developed in glass bottles during the system size experiment. Table S1. Taxonomic assignments of the OTUs displayed in the heatmap in Fig 2. Table S2. LEfSe analysis of predicted sulphate-reducing OTUs in the 16-week microcosms. The Desulfovibrio and Desulfobacteraceae species listed are different 97% OTUs assigned to these taxa. Only OTUs with LDA scores higher than 4.0 are shown. Table S3. The dimensions of the microcosms used for investigating system size effects on microbial community development. Table S4. Taxonomy of sequences that were less than 1% abundant in the 16-week microcosms from Supporting Information Fig. S5.