Bacterial community assembly in activated sludge: mapping beta diversity across environmental variables

Abstract Effect of ecological variables on community assembly of heterotrophic bacteria at eight full‐scale and two pilot‐scale activated sludge wastewater treatment plants (AS‐WWTPs) were explored by pyrosequencing of 16S rRNA gene amplicons. In total, 39 samples covering a range of abiotic factors spread over space and time were analyzed. A core bacterial community of 24 families detected in at least six of the eight AS‐WWTPs was defined. In addition to the core families, plant‐specific families (observed at <50% AS‐WWTPs) were found to be also important in the community structure. Observed beta diversity was partitioned with respect to ecological variables. Specifically, the following variables were considered: influent wastewater characteristics, season (winter vs. summer), process operations (conventional, oxidation ditch, and sequence batch reactor), reactor sizes (pilot‐scale vs. full‐scale reactors), chemical stresses defined by ozonation of return activated sludge, interannual variation, and geographical locations. Among the assessed variables, influent wastewater characteristics and geographical locations contributed more in explaining the differences between AS‐WWTP bacterial communities with a maximum of approximately 26% of the observed variations. Partitioning of beta diversity is necessary to interpret the inherent variability in microbial community assembly and identify the driving forces at play in engineered microbial ecosystem.


O R I G I N A L R E S E A R C H
Bacterial community assembly in activated sludge: mapping beta diversity across environmental variables Siavash Isazadeh 1,2 | Shameem Jauffur 1 | Dominic Frigon 1

| INTRODUCTION
Activated sludge (AS) system employed in wastewater treatment plants (WWTPs) is among the world's largest biotechnological processes.
Over the last century, AS process has experienced continuous modifications in design and operation to improve its efficiency (Seviour & Nielsen, 2010). At the heart of the AS processes, consortia of heterotrophic microorganisms transform incoming organic compounds into biomass and CO 2 through specific metabolisms. Modifications of the process often aim at promoting control over specific eco-physiological characteristics of the microbial community. Thus, rationally developing further process modifications require the recognition and understanding of the key abiotic factors influencing community assembly within AS systems. As such, quantifying how the variations in the environmental or operational variables influence the microbial community composition remains a fundamental goal of wastewater microbiology.
Advancements in molecular biology techniques used to study the diversity of 16S rRNA or functional genes (e.g., FISH, T-RFLP, DGGE, and amplicon cloning) have provided new opportunities to better understand the complexity of microbial ecosystems. High-throughput | 1051 sequencing techniques have expanded considerably the depth of description of microbial diversities in wastewater treatments plants with AS process (AS-WWTPs) (Pinto & Raskin, 2012). Yet, the challenge remains to ascertain the effect of abiotic variations (e.g., operational, spatial, and temporal) on the structure of the assembly of bacterial communities.
Several lines of evidence link bacterial community compositions and performance of biological wastewater treatment processes to operational variations such as: influent composition (Akarsubasi et al., 2005;Lee, Kim, Hwang, O'Flaherty, & Hwang, 2009), reactor configuration (Rowan et al., 2003), temperature variation (Siripong & Rittmann, 2007), and solids retention time (SRT) (Akarsubasi, Eyice, Miskin, Head, & Curtis, 2009). In addition, ecosystem size has been found to influence the observed distances in community composition (called beta diversity over time (Soininen, 2010), and reactor scale can affect species selection based on cellular morphology (Martins, Pagilla, Heijnen, & van Loosdrecht, 2004). In a comprehensive laboratoryscale study, Pholchan, Baptista, Davenport, and Curtis (2010) found that the microbial community of AS reactors was affected by different operational conditions and reactor configurations; however, the relationship between the performance and community diversity was not explicitly associated. Zhang, Shao, and Ye (2012)  and then 4 years after to assess the relative stability of the community at one plant. Differences among the plants in influent characteristics, treatment processes (conventional, oxidation ditch, and sequence batch reactor [SBR]), and geographic locations (over a region of 68-km radius) provided the other environmental variables. In addition, environmental variables covering chemical stress, reactor scale, SRT, and reactor configuration (fully aerobic vs. anoxic/aerobic) and temporal drift (interannual) were investigated in more details using samples collected during a pilotscale experiment at one of the treatment plants (LaPrairie-WWTP) over a 2-year period. For one of the pilot-scale reactors, ozone was applied to the return activated sludge (RAS), which provided a chemical stress by enhancing bacterial mortality (decay) and modified the substrate composition in the system. In return, ozone effects on the RAS are likely to modify the AS community structure, a hypothesis that was tested by comparing the communities of two pilot-scale reactors: RAS-ozonated vs. nonozonated. The reactor-scale gradient was studied at the same location by comparing pilot vs. full-scale reactors. Interannual variations in the bacterial population assembly of full-scale reactor were studied over a period of 4 years. The results generated from the combination of these studies enabled us to explore a range of environmental variables, which could potentially explain observed bacterial population assemblies and their response to abiotic changes.

| Sampling at full-scale AS-WWTPs
The mixed liquor samples were obtained at eight full-scale AS-WWTPs located within a 68-km radius of the LaPrairie-WWTP ( Fig. S1). The selected plants differed in size, influent flow rates and characteristics, and treatment processes (Table S1). The sampling campaign was conducted over two time scales: samples were obtained during the same year in two seasons (summer and winter), and then one sample was obtained 4 years after to evaluate interannual stability in bacterial community composition. Samples were collected during August/September 2008 (summer of first year), February 2009 (winter of first year), February 2013 (winter 4 years after).
In addition, three mixed liquor samples were obtained on consecutive weeks in August 2008, at one of the plants (Granby-WWTP), to explore the population variation in weekly base. All mixed liquor samples, from each WWTP, after collection were transported to the laboratory on ice, and the solids were centrifuged and frozen at −80°C the same day until molecular biology analysis.

| Pilot-scale study with ozonation or return activated sludge
Ozonation of return activated sludge is a technique for minimization of waste biosolids production by chemical oxidation and enhanced solids degradation (Yasui & Shibata, 1994). Direct exposure of bacterial community to ozone has the potential to shape the community assemblies either through the high mortality induced by the ozonation or changes in the increased availability of organic substrates.
Previously, we investigated the ozone impact on the microbial community composition by comparing two pilot-scale reactors (a control and a RAS-ozonated test reactor) in oxic condition (Isazadeh, Ozcer, & Frigon, 2014). Our comparison of FISH versus pyrosequencing experiment showed that some phyla were underrepresented in pyrosequencing method. Therefore, in this study, we used different forward primer sets (3 primers) and looked for more detailed operational conditions to explore the ozone effect on microbial population.
Two experimental periods of 6-8 months were performed using pilot-scale AS systems receiving the same wastewater as the LaPrairie-WWTP; each periods were divided in phases of a few weeks with different AS operational conditions (SRT and conventional vs. anoxic/ aerobic process) and ozone dosages to determine the potential of RAS-ozonation effect on bacterial community (for additional details see supplementary materials and Table S1). Five and six biomass samples were collected, respectively, from the RAS-ozonated and control (non RAS-ozonated) reactor. Samples were collected at the end of each experimental phase (longer than 3 × SRT) to ensure the impact of the operation phase on the community composition.
Solids from mixed liquor samples were centrifuged in the plant, and immediately frozen at −80°C until molecular biology analysis. Details of experimental design and operational conditions along with influent characteristics were presented in supplementary materials ( Fig. S2 and Table S2) and elsewhere (Isazadeh, Feng, Urbina Rivas, & Frigon, 2014;Isazadeh, Urbana, Ozcer, & Frigon, 2015).
In parallel to the pilot-scale reactor sampling, seven samples were also collected from the full-scale LaPrairie AS-WWTP to investigate the effect of scale (control pilot-scale reactor samples vs. full-scale reactor samples) and temporal variation on the bacterial community assembly.

| Determination of microbial community composition
Genomic DNA was extracted from the mixed liquor suspended solids using a DNA extraction kit. Details of DNA extraction and PCR amplification are given in the supplementary material The V3-V4 region of the 16S rRNA genes was then PCR amplified using a mixture of three forward primers and a reverse primer (to increase the phyla coverage), all with additional pyrosequencing adaptor sequences (Pinto & Raskin, 2012). The PCR amplicons were purified using a PCR purification kit; and sequenced using a GS first year (Y1.PII) did not provide a reliable reads (less than 1000 read) and therefore has been removed from the downstream alpha and beta diversity analyses.
Postsequencing analysis was performed using the Qiime pipeline following the sequential analyses explained in the supplementary material. (Caporaso et al., 2010). We normalized the number of sequences prior to beta diversity to have the depth of coverage for even sampling.

| Data analyses and variation partitioning of beta diversity
Detailed exploratory data analyses were carried out in the R software (V3.0.1) using the Vegan (V.2.08) (Oksanen et al., 2011) and cluster packages for heatmap (Maechler, Rousseeuw, Struyf, Hubert, & Hornik, 2011). To measure the statistical significance of within plant seasonal and interannual variations we used the multivariate dispersion approach based in distances to centroids (Anderson, Ellingsen, & McArdle, 2006). In this analysis Betadisper function in vegan was used for multivariate, permutation-based hypothesis tests for differences in structure centroid and dispersion (beta diversity). Alpha diversities indices were compared based on a two-sample t test using nonparametric method in Qiime. The details of this section are provided in the supplementary materials. Our initial study (Isazadeh, Ozcer et al., 2014) and other reports (Gobet, Quince, & Ramette, 2010) showed no significant changes in beta diversity results before and after denoising, therefore we did not used the denoised data. Principal coordinate analyses (PCoA) were performed based on weighted Unifrac and Hellinger distances (Legendre & Legendre, 2012). Since both distances resulted in similar ordination plots, the Hellinger distance is presented herein. Partitioning of variations were performed using the varpart function of the Vegan Package (Oksanen et al., 2011). For LaPrairie-WWTP, two explanatory matrices (containing either environmental or temporal explanatory variables) were used to investigate the variation in beta diversity (Table   S3). The environmental matrix included reactors scale (control pilot-scale vs. full-scale), RAS-ozonation (control vs. RAS-ozonated pilot-scale reactors), and season (summer vs. winter); whereas temporal matrix included the interannual variations (2008-2009 vs. 2013).

| Sequence accession numbers
The 454 pyrosequencing reads have been deposited in the NCBI Sequence Read Archive, with the accession numbers: from SRS673933 to SRS673972. in winter season. This allowed the evaluation of seasonal and interannual differences in each treatment plant. From these samples, 83,248 of 16S rRNA gene amplicon sequence reads were obtained with the number of reads per sample ranging between 1,505 and 4,262 (average of 3,619 reads/sample). Unique OTUs were defined by grouping together sequence reads with ≥97% identities, which resulted in 5,610 OTUs. The average OTU richness within a sample was 643 with a range between 420 and 823 ( Table 1).

| Bacterial community assemblies in regional full-scale AS-WWTPs
Comparison of the observed OTU alpha diversities (i.e., within sites) between samples gotten in summer 2008 and winter 2009 (only a few months apart) showed significantly higher diversity in the winter than in the summer (paired student t test, p < .05; Table 1). In general, this result is attributable to both higher OTU richness and evenness. The only noticeable counter example is the case of Cowansville WWTP where the community diversity was observed to be much lower in the winter than in the summer; however, this may be due to a lower number of sequence reads recovered for the Cowansville winter 2009 sample (Table 1). On average, the temperature drop between the summer and winter season, at these WWTPs, was about 8°C (averages for summer being 23°C and winter being 15°C).

Differences in community compositions between samples obtained
from the eight AS-WWTPs were visualized using a PCoA of the Hellinger distances (Fig. 1A). In general, the differences in community composition between samples from different plants were relatively greater than the ones between samples of the same plant. This is true despite the fact that mixed liquor samples were taken in different seasons and 4 years apart for one set. Furthermore, the season did not seem to generate specific trends in the variations in the composition of the microbial community as the relative position of summer 2008 and winter 2009 T A B L E 1 Sequence reads, number of OTUs (total and shared), and biodiversity numbers in eight AS-WWTPs samples varied between plants (Fig. 1A). This low impact of the season on the overall microbial community composition can clearly be seen in the samples from Granby. In this plant, the differences between the summer and winter samples were less than the differences observed between samples obtained on three consecutive weeks (Fig.1A). Altogether, these data suggested a relatively slow community turnover within plants and a relatively low impact of the season on the microbial community composition compared to differences between plants.

AS-WWTPs
Although the relative differences in community composition between plants were much greater than ones within plants; it was still possible to identify a core set of phyla and families that were present in most samples.  S4a). Within these phyla, a group of common families could be defined as those present in at least six WWTPs based on the total number of reads that they represented ( Fig. S6). Some of the most abundant of these families (and Table S6) were as follows: the Flavobacteriaceae and Saprospiraceae (phylum Bacteroidetes), the Rhodobacteriaceae and Sphingomonadaceae (class Alphaproteobacteria), and the Comamonadaceae (class Betaproteobacteria). It is mainly differences between the abundances of these families that is depicted in Figure 1.
Beside the core group of families, a group of plant-specific and abundant (>6% reads from the plants) families could also be identified (
Principal coordinate analysis (PCoA) was performed to visualize the differences in community composition (i.e., beta diversity). In the full-scale reactor, the community composition did not vary much over the three sampling years as all the points appear close to each other on the PCoA plot ( Fig. 2A). To provide an idea of magnitude of differences, the points reported for the LaPrairie-WWTP in Fig. 1 (among WWTP differences) are equivalent to the points for the full-scale reactor in Consequently, the similarity of the microbial composition between the two pilot-scale reactors (i.e., the proximity of paired samples in Fig. 2A) remained relatively high throughout the study. Sphingobacteriales (phylum Bacteroidetes) (Fig. S5). These families explain for the most part the PCoA1 ordination ( Fig. 2A). Together, the abundances of the class Anaerolineae population were higher in the full-scale than in the pilot-scale reactors; and they decreased in both pilot-scale reactors (RAS-ozonated and control) in Year 2 of the pilot-scale experiment when the treatment systems were changed T A B L E 2 Sequence reads, OTUs (total and shared), and biodiversity numbers for LaPrairie AS-WWTP from anoxic/aerobic to completely aerobic. Reviewing the full-scale WWTP operation data suggested that the front-end of the aerobic plug-flow reactor was deficient in aeration, and that significant denitrification occurred (data not shown). Therefore, the presence of Anaerolineae seemed to be related to the denitrification in this system. Thus variations in treatment process conditions and aeration efficiencies as opposed to specifically reactor scale may be the source of differences in community composition between the full-scale and pilot-scale reactors.
The family Xanthomonadaceae also change significantly in abundance in both pilot-scale reactors through the Year 2 pilot-scale experiment. When the SRTs decreased from ~12 day to ~6 day, Xanthomonadaceae increased in abundance; an observation similar to the one obtained during a previous pilot-scale study at the same site during which the SRT was also 6 day (Isazadeh, Ozcer et al., 2014).
Thus, low SRTs (i.e., high growth and dilution rates) appear to favor the growth of Xanthomonadaceae in this system.

| Partitioning of the beta diversity
The  (Table 3).
However, most of the explanatory power of these two sets of state variables was shared. The environmental conditions (process types and seasons) could only add 5% to the explanation of the variance in community compositions (Table 3) (Table 3). Along with the complete full-scale dataset, this result suggests that a strong underlying drift in community composition is more important in determining the composition than the summer/winter seasonal effect.
The variance among the LaPrairie-WWTP community compositions (full-scale and both pilot-scale reactors) was partitioned over a set of variables describing the environmental conditions (namely scale of the reactor, process operation, and season), and another one describing the sampling year. Together, these two sets of state variables could only explain 13% of the variance in OTU community composition (Table 3).

| Core bacterial community of conventional activated sludge systems
The OTU structures of the AS-WWTP bacterial communities determined in the current investigation were in line with previously observed studies; the OTU distributions followed a strong power law rank-abundance model in which the top 10-20 OTUs (grouped in as many genera) accounting for 70-80% of the sequence reads, and a long tail of OTUs (~150 genera) accounting for the remaining reads (Hoffmann et al., 2007;Xia et al., 2010). The identities and abundances of observed phyla were also highly similar to the previous studies (Seviour & Nielsen, 2010;Zhang et al., 2012).
The resemblance of the community assemblies observed in this study were high over time and between WWTPs as determined at the levels of OTU (97%-identity cut-off), family and phylum (Table   S6). Such similarities in bacterial community have been reported by others when comparing bacterial community assemblies at AS-WWTPs in China and in North America (Zhang et al., 2012). In this study, common families (Table S6)  The differences in microbial communities revealed by the principal coordinate analysis were mainly related to variations in the abundance of these core microbial families (Fig. 1). Among the families describing the differences between plants are the ones belonging to the class Anaerolineae and the family Saprospiraceae (phylum Bacteroidetes), which are pointing upward along PCoA2, and the family Xanthomonadaceae, which is pointing downward along PCoA2 (Fig. 1). This pattern reveals a positive correlation in the abundances of the class Anaerolineae and the family Saprospiraceae, and a negative correlation with the abundances of the family Xanthomonadaceae. The same correlations were present within the data of the pilot-scale reactors At LaPrairie WWTP despite a much reduced diversity in operation conditions. Finally, the same correlations were also observed in other studies studying either differences between plants or temporal variations within a plant (Ju & Zhang, 2015;Vuono et al., 2015). Since the variations in operation conditions were very different in all these datasets, this suggests some biotic interactions (e.g., competition or cooperation) exist between the members of these families. While members of the family Saprospiraceae have been identified as strong protein-hydrolysers (Xia, Kong, & Nielsen, 2007), the role of members of the two other taxa (class Anaerolineae and family Xanthomonadaceae) remain to be defined.

| Plant-specific abundant families of conventional activated sludge systems
Despite the presence of a core set of families, this study identi- the importance of influent wastewater composition on shaping the bacterial community structure. Analysis of the influent wastewater from LaPrairie-WWTP revealed a high methanol and nitrate concentrations, which in turn explains the high abundance of the Methylotenera genus in this AS-WWTP (Isazadeh, Ozcer et al., 2014).
At other full-scale AS-WWTPs studied, a number of plant-specific abundant bacterial families were identified (Table S6). Among these, the genera from the families Gordoniaceae, Thiotrichaceae are wellknown for their implication in solid-liquid separation problems including bulking and foaming due to their filamentous morphology (Martins et al., 2004). The presence of these organisms has also been linked to influent compositions. For examples, members of the Gordoniaceae are present at high lipid loading rates (Frigon, Muyzer, van Loosdrecht, & Raskin, 2006); and members of the Thiotrichaceae are present at high volatile fatty acid concentration and in the presence of sulfides as they are capable of mixotrophic growth (Martins et al., 2004). Although the links between the influent composition and other environmental factors is not established for all the plantspecific highly abundant populations, these observations reported here demand further investigations of their role in activated sludge microbial communities.

| Level of explanation of variation in microbial community composition
Partitioning in the influent (Isazadeh, Ozcer et al., 2014). Consequently, wastewater treatment microbial ecologists need to go beyond the conventional influent characterization such as BOD 5 , COD, total phosphorus and ortho-phosphate, total Kjeldhal nitrogen, and volatile suspended solids to describe the composition of incoming wastewater with a high enough precision to understand the community structure and meaningfully explore the structure-functions relationships of heterotrophs.
Domestic wastewater contains: proteins (40%-60%), carbohydrates (25%-50%), fats and oils (10%), urea and a large number of organic compounds including pesticides and herbicides (Bitton, 2011 Consequently, both groups, who worked with rRNA-gene-targeted TRFLP or DGGE, concluded that higher resolution molecular biology techniques (i.e., deeper sampling of community diversity) would reveal the link between the community composition and operational changes. This study was performed by pyrosequencing and it reached the same results. Therefore, it did not support such speculation.
This being said, the results of community composition variance partitioning could be related to other mechanistic and theoretical reasons, namely that communities could be assembled by neutral mechanisms (Hubbell, 2001). Curtis and Sloan (2006) used the stochastic approach based on random-assembly to describe autotrophic populations of ammonia oxidizing bacteria (AOB) in wastewater systems. Ofiţeru et al. (2010) suggested that neutral community models should form the foundation of any description of open biological system. There are also room for a spectrum of theoretical mechanisms between niche-assembly and random-assembly approach if one considers that the niche exclusion principle does not absolutely limit community diversity and that even competition may be compatible with more neutral community assembly models. For example, competition dynamics can generate periodic or chaotic oscillations in the abundances of species that can generate niches with more species than | 1059 limiting resources (Huisman & Weissing, 1999). This last finding was also supported by a modeling approach based on game theory which showed that the diversity in a local niche is more dependent on the meta-community diversity for a given function than the number of limiting resources associated with the niche (Allesina & Levine, 2011).
These considerations suggest that much more data (i.e., limited number of samples in this experimental design along with missing data) and theoretical development will be necessary to understand community assemblies in biological wastewater treatment systems.
Finally, this study provides important validations of pilot-scale studies in understanding the microbial community dynamics. In this study, the community compositions between the full-scale and pilotscale reactors at LaPrairie-WWTP were highly similar compared to variations between plants (Fig. 1), and the diversities was essentially the same at both scales (Table 2)