**DOI:**10.1128/AEM.03286-14

## ABSTRACT

Biodiversities can differ substantially among different wastewater treatment plant (WWTP) communities. Whether differences in biodiversity translate into differences in the provision of particular ecosystem services, however, is under active debate. Theoretical considerations predict that WWTP communities with more biodiversity are more likely to contain strains that have positive effects on the rates of particular ecosystem functions, thus resulting in positive associations between those two variables. However, if WWTP communities were sufficiently biodiverse to nearly saturate the set of possible positive effects, then positive associations would not occur between biodiversity and the rates of particular ecosystem functions. To test these expectations, we measured the taxonomic biodiversity, functional biodiversity, and rates of 10 different micropollutant biotransformations for 10 full-scale WWTP communities. We have demonstrated that biodiversity is positively associated with the rates of specific, but not all, micropollutant biotransformations. Thus, one cannot assume whether or how biodiversity will associate with the rate of any particular micropollutant biotransformation. We have further demonstrated that the strongest positive association is between biodiversity and the collective rate of multiple micropollutant biotransformations. Thus, more biodiversity is likely required to maximize the collective rates of multiple micropollutant biotransformations than is required to maximize the rate of any individual micropollutant biotransformation. We finally provide evidence that the positive associations are stronger for rare micropollutant biotransformations than for common micropollutant biotransformations. Together, our results are consistent with the hypothesis that differences in biodiversity can indeed translate into differences in the provision of particular ecosystem services by full-scale WWTP communities.

## INTRODUCTION

The microbial communities residing within wastewater treatment plants (WWTPs) provide ecosystem services that are important for maintaining the quality of our environment. They consume carbon-, nitrogen-, and phosphorous-containing substrates and biotransform organic pollutants present within WWTP influent (1, 2), thus mitigating the potentially deleterious effects of these chemicals on receiving waters. The ability to perform these ecosystem functions is to some extent related to the biodiversity present within WWTP communities (3). WWTP communities are estimated to contain many thousands of different microbial strains and to express tens of thousands of different genes (e.g., 4–7), thus providing a supply of functional traits that could enable and facilitate particular ecosystem functions.

While biodiversity is thought to be important for enabling and facilitating particular ecosystem functions (3), biodiversity can also differ substantially both among different WWTP communities (4, 5, 7, 8) and over time for a single WWTP community (9, 10). In one recent investigation, the taxonomic richness measurements among 10 full-scale WWTP communities differed by 4- to 6-fold (4). These differences in biodiversity underscore one of the most general and debated questions in ecology: how do differences in biodiversity translate into differences in the rate of a particular ecosystem function (11, 12)? Stated alternatively, what is the shape of the association between biodiversity and the rate of a particular ecosystem function? Determining the shape could be important for understanding why some WWTP communities perform a particular ecosystem function at higher rates than others. It could also be important for predicting how differences or changes in biodiversity are likely to affect the rate of that particular ecosystem function.

The emerging view from experiments with other types of communities is that the association between biodiversity and the rate of an ecosystem function is positive (11, 12) (Fig. 1). The positive association is thought to result from the accumulation of strains that have unique niche partitioning or facilitative interaction effects on a particular ecosystem function (i.e., complementarity effects) (11, 13, 14). Niche partitioning effects occur when different ecological niches are inhabited by different specialist strains rather than by any single generalist strain, and the collection of specialist strains achieves a greater aggregate rate of a particular ecosystem function than any single generalist strain. Facilitative interaction effects occur when the presence of some strains stimulates the activities of other strains (e.g., by providing growth-limiting nutrients, detoxifying local environments, etc.), thus leading to higher rates of a particular ecosystem function. The positive association could also result from sampling effects (11, 13, 14). Briefly, communities with more strains are more likely to contain strains that provide niche partitioning or facilitative interaction effects or have large trait values for a particular ecosystem function. Sampling effects not only could modify the rates of particular ecosystem functions but could also result in qualitatively different capabilities among different communities (i.e., sampling effects could affect the presence or absence of a particular ecosystem function).

While the association between biodiversity and the rate of a particular ecosystem function is thought to be positive, the association could have one of two different shapes (11, 12) (Fig. 1). A nondecelerating shape is expected when different strains have complementarity and sampling effects that are unique rather than the same as those of other strains (Fig. 1, dashed line). This might occur if there were far less biodiversity present than would be needed to begin saturating the set of possible positive effects. Conversely, a decelerating shape is expected when different strains have to some extent the same complementarity and sampling effects as other strains (i.e., functionally redundant effects) (Fig. 1, solid line). This would occur if there were sufficient biodiversity present to begin saturating the set of possible positive effects. An important consequence of deceleration is that at some point, differences in biodiversity might no longer translate into differences in the rate of that particular ecosystem function (Fig. 1, far right region of the solid line).

While both nondecelerating and decelerating associations have been observed for some types of ecosystem functions and microbial communities (e.g., see references 3 and 15 to 18), associations have not been observed for other types (e.g., see references 19 to 21). Moreover, the associations among full-scale WWTP communities remain largely unexplored. It is plausible that full-scale WWTP communities are sufficiently biodiverse to nearly saturate the set of possible complementarity and sampling effects for a particular ecosystem function. If this were the case, then a decelerating association would be expected, and differences in biodiversity might not translate into differences in the rate of that particular ecosystem function (Fig. 1, far right region of the solid line). On the other hand, full-scale WWTPs contain considerable spatial, temporal, and resource heterogeneity, thus providing a large set of possible complementarity and sampling effects. It is therefore plausible that full-scale WWTPs are insufficiently biodiverse to substantially saturate the set of possible complementarity and sampling effects. If this were the case, then a nondecelerating association might be expected, and differences in biodiversity would translate into differences in the rate of that particular ecosystem function (Fig. 1, dashed line). Indeed, recent studies demonstrated that nondecelerating associations are more likely to occur as spatial and temporal heterogeneity increase (e.g., see references 16, 17, and 22).

Our main objective was to test for and measure the shapes of the associations between biodiversity and the rates of 10 different micropollutant biotransformations among 10 independent and full-scale WWTP communities, thus addressing an important gap in our knowledge about their ecology and functioning. We use the term micropollutant to refer to a synthetic organic chemical that is present at trace concentrations in WWTP influent (23). We previously measured the rate at which each WWTP community biotransforms each micropollutant (24). In parallel, we measured the taxonomic and functional biodiversity of each WWTP community from bacterial 16S rRNA and mRNA sequence reads, respectively (4). In this study, we tested for associations between different aspects of biodiversity and the rates of each individual micropollutant biotransformation. We then asked whether a decelerating or nondecelerating shape more accurately describes each of the observed positive associations. We finally proposed a hypothesis and obtained evidence for why the positive associations are stronger for some micropollutant biotransformations and weaker for others.

## MATERIALS AND METHODS

WWTP communities and biotransformation rate constants.We previously collected WWTP communities from 10 independent and full-scale WWTPs located across north-central Switzerland. We sampled each WWTP using a standardized protocol (24). Briefly, we collected 1 liter of activated sludge directly from the biological aeration basin of the WWTP, transported the activated sludge to the laboratory in a loosely capped 2-liter amber glass bottle, and added a magnetic stir bar to the 2-liter bottle immediately upon arrival at the laboratory.

We initiated a biotransformation experiment with each WWTP community within 3 h of collecting the activated sludge samples from the WWTP (24). We dissolved a defined mixture of 10 micropollutants in methanol to obtain a stock concentration of 100 mg per liter of each micropollutant, added 70 μl of the stock mixture into empty triplicate-stirred batch reactors, and allowed the methanol to completely evaporate. We next seeded the batch reactors with 70 ml of activated sludge to obtain an initial concentration of 100 μg per liter of each micropollutant (24). We then periodically removed 2-ml liquid samples from the stirred batch reactors over 4 days, filtered the liquid samples using a glass-fiber filter, stored the liquid samples at 4°C, and analyzed the samples using high-performance liquid chromatography-mass spectrometry (HPLC-MS) chemical analyses (24). We performed abiotic controls in parallel using twice-autoclaved activated sludge (24). We quantified the biologically mediated rate of disappearance of each micropollutant and determined the observed first-order rate constant using Bayesian parameter inference and Monte Carlo sampling as previously described (24). We normalized the biotransformation rate constants to empirical estimates of active biomass as described by Majewsky and colleagues (25) (see Table S1 in the supplemental material). We further normalized the biotransformation rate constants to total suspended solids (TSS) and obtained qualitatively similar outcomes. Because active biomass and TSS normalization generated qualitatively similar results, we only present data using active biomass normalization for the remainder of the article.

RNA sample collection and processing.We collected 1.5-ml liquid samples for RNA analyses from each of the WWTP communities at 2 h after mixing the WWTP communities with the micropollutants (4). We analyzed RNA rather than DNA because RNA analyses consider only metabolically active strains and transcriptionally expressed genes, thus excluding metabolically inactive strains and nonexpressed genes, which are unlikely to have contributed to the measured micropollutant biotransformations. We collected samples for RNA analyses at 2 h for three reasons. First, 2 h provided time for the WWTP communities to respond to the addition of the micropollutants. Second, because the micropollutants were typically biotransformed without an observable lag period and followed first-order kinetics, 2 h coincided with a rapid period of micropollutant biotransformation. Third, the micropollutants were not completely biotransformed after 2 h (24). We provide full descriptions of the RNA isolation and purification methods elsewhere (4).

Preparation, sequencing, and analysis of bacterial 16S rRNAs.We reverse transcribed bacterial 16S rRNAs, amplified the B27F-B534R region of the reverse transcription products, and sequenced the reverse transcription products using the 454-GS FLX platform (Roche 454 Life Sciences, Branford, CT, USA) as described elsewhere (4). All of the resulting 16S rRNA sequence reads are publically available in the NCBI Sequence Read Archive (see below). We analyzed the 16S rRNA sequence reads using mothur software (26) and assigned those reads that passed all of the quality control steps into operational taxonomic units (OTUs) using a sequence similarity threshold of 97%. We chose 97% because this threshold has been widely used in comparable studies and reduces the probability that sequencing errors could lead to artifactual estimates of biodiversity. We provide a full description of the 16S rRNA sequencing methods and quality control parameters elsewhere (4).

Preparation, sequencing, and analysis of mRNAs.We enriched mRNAs and sequenced the enriched mRNAs using the HiSeq2000 platform (Illumina, San Diego, CA, USA) as described elsewhere (4). All of the resulting mRNA sequence reads are publically available in the MG-RAST database (see below). We analyzed the mRNA sequence reads using the MG-RAST metagenomics analysis server (version 3.2) (27) and assigned those reads that passed all of the quality control steps into level 3 SEED subsystems functional categories (28). We provide a full description of the mRNA sequencing methods, quality control parameters, and annotation parameters elsewhere (4).

Biodiversity measurements.Prior to measuring biodiversity, we corrected for differences in sequencing depth across the 10 WWTP communities by rarefying the 16S rRNA sequence reads that were assigned to OTUs to 2,500 and the mRNA sequence reads that were assigned to level 3 functional categories to 300,000. We measured the observed richness of each WWTP community as the observed number of unique OTUs or unique level 3 functional categories. We measured the Chao1 extrapolated richness (29), the abundance-based coverage estimation (ACE) extrapolated richness (29, 30), the Shannon diversity, and the Shannon evenness of each WWTP community from the same sets of OTUs or level 3 functional categories. We reported the Shannon diversity and Shannon evenness measurements in terms of the Hill numbers *H*_{1} and *H*_{0}, which provide more intuitive interpretations of the measurements (31). All of the reported biodiversity measurements are the average values from 1,000 independently rarefied 16S rRNA or mRNA sequence read data sets. We further rarefied the 16S rRNA sequence reads that were assigned to OTUs to 1,000, 1,500, and 2,000 and found that the rank orderings of the biodiversity indices were identical to those when we rarefied the16S rRNA sequence reads to 2,500 (data not shown). We therefore present only the data obtained when we rarefied the 16S rRNA sequence reads to 2,500. We performed rarefaction and calculated all biodiversity measurements in the R software environment (32) using functions from the vegan package (33).

Micropollutant multifunctionality measurements.We measured the collective rates of multiple micropollutant biotransformations using the multifunctionality measure described by Zavaleta and colleagues (34). Briefly, we scaled the rate constants for each individual micropollutant biotransformation to an average value of 0 and a standard deviation of 1. We then calculated micropollutant multifunctionality as the average value of the scaled rate constants among the 10 individual micropollutant biotransformations. The advantages and limitations of this measure of multifunctionality have been discussed in detail elsewhere (34, 35).

Linear and logarithmic models.We fit linear (*y* = [*m* × *x*] + *b*) and logarithmic (*y* = [*m* × LN(*x*)] + *b*) models to the biodiversity measurements in order to gain insight into the underlying relationships between biodiversity and micropollutant biotransformations, where the dependent variable *y* is the rate constant for an individual micropollutant or micropollutant multifunctionality, the independent variable *x* is a measure of biodiversity, *m* is the slope, and *b* is the intercept. We fit each model, estimated the model parameters, and calculated the coefficient of determination (*R*^{2}) and the root mean squared error (RMSE) in the R software environment using functions from the stats package (32).

Bioinformatic analyses.We used the UM-BBD database (36) to predict sets of bacterial enzymes that could hypothetically catalyze each of the experimentally measured micropollutant biotransformations. We used the NCBI protein search tool (www.ncbi.nlm.nih.gov/protein/) to retrieve all taxonomically identified protein sequences for each predicted enzyme. We performed the protein search across eight publically available sequence databases (RefSeq, GenBank, EMBL, DDBJ, PDB, UniProtKB/Swiss-Prot, PIR, and PRF).

Statistical methods.Because the underlying distributions of our measurements are unknown, we used nonparametric methods for all of our statistical hypothesis tests. We used the two-sided Spearman rank correlation test to test whether the biodiversity measurements associate with the rate constants of each individual micropollutant biotransformation or with micropollutant multifunctionality. We used the Benjamini-Hochberg method to account for multiple comparisons (37). We used the one-sample two-sided Wilcoxon test to test whether the Spearman rank correlation coefficients of all 10 of the individual micropollutant biotransformations significantly deviate from zero and whether the values of *R*^{2}_{linear} − *R*^{2}_{logarithmic} significantly deviate from zero. For the last two tests, we assumed that the individual micropollutant biotransformations are independent of each other (see Results). We performed all of the statistical tests in the R environment using functions from the stats package (32).

Accession numbers.The 16S rRNA sequence reads obtained in this work are publically available in the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra) under BioProject number PRJNA232662. mRNA sequence reads are publically available in the MG-RAST database under project number 6015 (http://metagenomics.anl.gov/).

## RESULTS

Micropollutant biotransformation measurements.We previously measured the rates of 10 different micropollutant biotransformations by 10 independent and full-scale WWTP communities (24). We selected the WWTPs in order to obtain a set with substantial differences in their environmental and operational metrics, including total suspended solids, solid retention time, and influent source (24). We postulated that greater differences in those metrics would translate into greater differences in biodiversity and in turn into greater differences in the rates of micropollutant biotransformations. We summarized the environmental and operational metrics for each WWTP in Table S2 in the supplemental material and elsewhere (24). We selected the micropollutants to have diverse chemical structures and to putatively undergo different types of biotransformations (24). We summarized the experimentally observed biotransformation types and biotransformation products for each micropollutant in Table S3 and elsewhere (24). We measured the rate of each micropollutant biotransformation as the active biomass-normalized first-order rate constant (liters g^{−1} day^{−1}) (see Table S1). We used active biomass normalization because the biodiversity measurements were obtained for a fixed number of 16S rRNA or mRNA sequence reads, which are themselves surrogate measures of active biomass (38, 39).

Taxonomic biodiversity and individual micropollutant biotransformations.We first tested for associations between the taxonomic richness aspect of biodiversity and the rates of each individual micropollutant biotransformation. To accomplish this, we assigned the previously analyzed 16S rRNA sequence reads into operational taxonomic units (OTUs) using a sequence similarity threshold of 97% (4). We then measured the observed taxonomic richness of each WWTP community as the observed number of unique OTUs among 2,500 randomly sampled 16S rRNA sequence reads. We additionally measured the extrapolated taxonomic richness of each WWTP community using the Chao1 and abundance-based coverage estimation (ACE) indices (29, 30) from the same 2,500 randomly sampled 16S rRNA sequence reads. We summarized the taxonomic richness measurements in Table 1. We then used the Spearman rank correlation test to test whether the taxonomic richness measurements associate with the rate constants for each individual micropollutant biotransformation. This is a nonparametric test that does not make any assumptions about the underlying shape of the association; it is therefore valid for both nondecelerating and decelerating associations (Fig. 1). We plotted the taxonomic richness measurements against the rate constants for each individual micropollutant biotransformation in Fig. S1 to S3 in the supplemental material. We summarized the Spearman rank correlation coefficients in Fig. 2A.

We observed significant positive associations between the taxonomic richness aspect of biodiversity and the rates of specific individual micropollutant biotransformations, regardless of whether we used observed or extrapolated measures of taxonomic richness. This is based on two lines of evidence. First, the Spearman rank correlation coefficients are themselves significant and positive for 6 of the 10 individual micropollutant biotransformations (Fig. 2A) (Benjamini-Hochberg-adjusted two-sided *P* value of <0.05). Second, when all 10 of the individual micropollutant biotransformations are considered together, the Spearman rank correlation coefficients are more often positive than one would expect by chance (Fig. 2A) (one-sample Wilcoxon test, two-sided *P* value of 0.0039). These outcomes would only be expected if differences in taxonomic richness tend to translate into differences in the rates of specific individual micropollutant biotransformations. Taken together, our data indicate that WWTP communities with greater taxonomic richness are indeed more likely to biotransform specific individual micropollutants at higher rates.

We next tested for associations between the Shannon taxonomic diversity and Shannon taxonomic evenness aspects of biodiversity and the rate of each individual micropollutant biotransformation. As opposed to taxonomic richness, these aspects of biodiversity consider the proportional abundances of different taxa and place less weight on low-abundance taxa (40). We measured Shannon taxonomic diversity and Shannon taxonomic evenness using the same 2,500 randomly sampled 16S rRNA sequence reads that we used for measuring taxonomic richness. We summarized the Shannon taxonomic diversity and Shannon taxonomic evenness measurements in Table 1 and plotted the measurements against the rate constants for each individual micropollutant biotransformation in Fig. S4 and S5 in the supplemental material, respectively. We summarized the Spearman rank correlation coefficients in Fig. 2B and replotted the Spearman rank correlation coefficients for observed taxonomic richness to facilitate comparison.

We observed two important outcomes. First, the Spearman rank correlation coefficients for Shannon taxonomic diversity and Shannon taxonomic evenness are significant and positive for fewer of the individual micropollutant biotransformations than for observed taxonomic richness (Fig. 2) (Benjamini-Hochberg-adjusted two-sided *P* value of <0.05). Second, the median values of the Spearman rank correlation coefficients among all 10 of the individual micropollutant biotransformations are smaller for Shannon taxonomic diversity (0.63) and Shannon taxonomic evenness (0.61) than for observed taxonomic richness (0.70). Thus, the positive associations between biodiversity and the rates of individual micropollutant biotransformations appear to weaken when the proportional abundances of taxa are considered and less weight is placed on low-abundance taxa. Our data therefore lead to the hypothesis that low-abundance taxa are important for determining the rates of specific—but not necessarily all—individual micropollutant biotransformations.

Taxonomic biodiversity and multiple micropollutant biotransformations.While investigating individual micropollutant biotransformations is important for understanding the underlying relationships between biodiversity and the rates of ecosystem functions, the preferred outcome of wastewater treatment is typically not to remove any single micropollutant. Instead, the preferred outcome is typically to remove a large number and wide variety of different micropollutants. We therefore next tested for positive associations between the different aspects of biodiversity and the collective rates of multiple micropollutant biotransformations. To accomplish this, we measured the collective rate of multiple micropollutant biotransformations using the multifunctionality measure proposed by Zavaleta and colleagues (34). We summarized the micropollutant multifunctionality measurements in Table S4 in the supplemental material. We plotted the observed taxonomic richness, extrapolated taxonomic richness, Shannon taxonomic diversity, and Shannon taxonomic evenness measurements against the micropollutant multifunctionality measurements in Fig. 3. We summarized the Spearman rank correlation coefficients in Fig. 2.

We observed significant positive associations between all of the measured aspects of biodiversity and the collective rates of multiple micropollutant biotransformations (Fig. 3). The Spearman rank correlation coefficients are always significant and positive (Fig. 2) (two-sided *P* value of ≤0.016). Moreover, the Spearman rank correlation coefficients are typically larger than or equal in magnitude to the largest Spearman rank correlation coefficient among the individual micropollutant biotransformations (Fig. 2). Thus, the positive associations for the collective rates of multiple micropollutant biotransformations are typically stronger than or as strong as the positive associations for the rates of individual micropollutant biotransformations. Furthermore, the positive associations for the collective rates of multiple micropollutant biotransformations are weaker when the proportional abundances of taxa are considered and less weight is placed on low-abundance taxa (i.e., they are smaller for Shannon taxonomic diversity and Shannon taxonomic evenness than for taxonomic richness) (Fig. 2B). Thus, low-abundance taxa are also likely important for determining the collective rates of multiple micropollutant biotransformations.

Shapes of the observed associations.We next sought to gain insight into the underlying shapes of the observed positive associations. More specifically, we determined whether a nondecelerating or decelerating model explains more of the variance and enables more accurate predictions for the rates of individual or collective rates of multiple micropollutant biotransformations. To accomplish this, we used a linear model to describe a nondecelerating shape and a logarithmic model to describe a decelerating shape as proposed by Bell and colleagues (15). We fit the models to the micropollutant biotransformation measurements as a function of biodiversity and plotted the best-fit models in Fig. S1 to S5 in the supplemental material and Fig. 3. We summarized the best-fit models in Table S5. We then calculated the differences between the coefficients of determination for the best-fit linear and best-fit logarithmic models (*R*^{2}_{linear} *− R*^{2}_{logarithmic}) (Fig. 4). Thus, a positive difference indicates that the best-fit linear model explains more of the variance than the best-fit decelerating model. We further calculated the differences between the root mean squared errors for the best-fit linear and best-fit logarithmic models (RMSE_{linear} − RMSE_{logarithmic}) (see Fig. S6). Thus, a negative difference indicates that the best-fit linear model generates more accurate predictions than the best-fit logarithmic model. We excluded the biotransformation of diazinon from the analyses because the slopes of the best-fit linear and best-fit logarithmic models were always negative (see Table S5).

We were unable to detect a common shape for the associations between biodiversity and the rates of individual micropollutant biotransformations. This is based on two observations. First, among all of the individual micropollutant biotransformations, one type of model does not generally explain more of the variance (Fig. 4) or enable more accurate predictions (see Fig. S6 in the supplemental material). Linear models are marginally more effective for four of the individual micropollutant biotransformations (biotransformations of azoxystrobin, isoproturon, trinexepac ethyl, and venlafaxine), logarithmic models are marginally more effective for two of the individual micropollutant biotransformations (biotransformations of atenolol and valsartan), and neither model is clearly more effective than the other for three of the individual micropollutant biotransformations (biotransformations of ethofumesate, propachlor, and ranitidine) (Fig. 4; see also Fig. S6). Second, the values of *R*^{2}_{linear} − *R*^{2}_{logarithmic} among all of the individual micropollutant biotransformations do not significantly deviate from zero (one-sample Wilcoxon test, two-sided *P* value of >0.25). This is true regardless of the specific aspect of biodiversity that is considered. It is also true regardless of whether all of the individual micropollutant biotransformations are considered or only those individual micropollutant biotransformations that had both significant and positive Spearman correlation coefficients (Benjamini-Hochberg-adjusted two-sided *P* value of <0.05).

As opposed to the rates of individual micropollutant biotransformations, we detected clear nondecelerating shapes for the positive associations between biodiversity and the collective rates of multiple micropollutant biotransformations (Fig. 3). The best-fit linear models consistently explain more of the variance and enable more accurate predictions of micropollutant multifunctionality than the best-fit logarithmic models (Fig. 4; see also Fig. S6 in the supplemental material). Thus, at the levels of biodiversity that are present within full-scale WWTP communities, the positive associations between biodiversity and the collective rates of multiple micropollutant biotransformations do not appear to decelerate (Fig. 3).

Rarity of micropollutant biotransformations.We next asked why the positive associations between biodiversity and the rates of individual micropollutant biotransformations are stronger for some micropollutant biotransformations (e.g., venlafaxine) and weaker for others (e.g., trinexepac ethyl) (Fig. 2A). We hypothesized that the rarity of the micropollutant biotransformation (i.e., the number of different taxa that perform the biotransformation) is an important determinant. As a micropollutant biotransformation becomes increasingly rare, the probability that the set of possible complementarity and sampling effects is saturated should decrease, thus leading to a stronger positive association. To provide a first test of this hypothesis, we used the database UM-BBD (36) to identify bacterial enzymes that are predicted to perform each of the experimentally observed biotransformations (see Table S6 in the supplemental material). We next retrieved all publically available and taxonomically identified protein sequences for each predicted bacterial enzyme and quantified how many different bacterial taxa are predicted to perform each micropollutant biotransformation (Table 2). This quantity provides an estimate of how rare each micropollutant biotransformation is. We found that the positive associations for the observed taxonomic richness, extrapolated taxonomic richness, and Shannon taxonomic diversity measurements are indeed significantly stronger for rare micropollutant biotransformations than for common micropollutant biotransformations (Spearman rank correlation coefficients < −0.65; *P* < 0.05), thus providing evidence for this hypothesis.

Functional biodiversity and micropollutant biotransformations.The complementarity and sampling effects that give rise to positive associations are not directly caused by taxa but are instead caused by the phenotypes expressed by those taxa (41). These phenotypes could emerge via the expression of a wide range of different functions, including metabolic, motility, regulation, cell signaling, resistance, and stress response functions. Thus, we expect that the number of functions expressed by a WWTP community (i.e., functional richness) should also positively associate with the rates of individual and the collective multiple micropollutant biotransformations. Yet to our knowledge, this expectation has not been extensively tested. To test this, we measured the observed and extrapolated functional richness of each WWTP community from 300,000 randomly sampled and functionally annotated mRNA sequence reads that were assigned to level 3 functional categories. We summarized the functional richness measurements in Table 1. We then used the Spearman rank correlation test to test for associations between the functional richness measurements and the rate constants for each individual micropollutant biotransformation or for micropollutant multifunctionality. We summarized the Spearman rank correlation coefficients in Fig. 5.

We observed two important outcomes. First, as with taxonomic richness, the Spearman rank correlation coefficients for functional richness are themselves significant and positive for both the rate constants of specific individual micropollutant biotransformations (Fig. 5) (Benjamini-Hochberg-adjusted two-sided *P* value of <0.05) and for micropollutant multifunctionality (Fig. 5) (two-sided *P* value of <0.003). Second, when considering all 10 of the individual micropollutant biotransformations together, the Spearman rank correlation coefficients are more often positive than one would expect by chance (Fig. 5) (one-sample Wilcoxon test, two-sided *P* < 0.006). Thus, our data indicate that WWTP communities with greater functional richness are also more likely to biotransform specific individual and multiple micropollutants at higher rates.

Independence of micropollutant biotransformations.A key assumption of our experimental design is that each individual micropollutant biotransformation is independent of the others (i.e., each individual micropollutant biotransformation is performed by a different set of enzymes). This assumption is necessary, for example, to test whether the distributions of the Spearman rank correlation coefficients among all 10 of the individual micropollutant biotransformations significantly deviate from zero. To test their independence, we used a combination of bioinformatic and empirical approaches. We first analyzed the sets of bacterial enzymes that were predicted by the UM-BBD (36) to catalyze each of the experimentally measured micropollutant biotransformations (see Table S6 in the supplemental material). We found that unique sets of enzymes are predicted to biotransform six of the 10 individual micropollutants (see Table S6), thus providing bioinformatic support that these biotransformations are independent of each other. Of the remaining four, an identical set of enzymes is predicted to biotransform both azoxystrobin and trinexepac ethyl (see Table S6). However, the rate constants for these micropollutant biotransformations are not related to each other (Spearman rank correlation coefficient = 0.21; two-sided *P* > 0.50), thus providing empirical support that they are also independent of each other. Finally, an identical set of enzymes is predicted to biotransform both isoproturon and valsartan (see Table S6). In this case, we previously demonstrated that acetylene, which is an inhibitor of monooxygenase activity, completely eliminates the biotransformation of isoproturon but has no effect on the biotransformation of valsartan (24), thus providing empirical support that these micropollutant biotransformations are also independent of each other. In conclusion, our assumption that the different micropollutant biotransformations are independent of each other is supported by bioinformatic or empirical evidence.

## DISCUSSION

Our results revealed important insights into how differences in biodiversity translate into differences in the rates of micropollutant biotransformations among full-scale WWTP communities. While we observed significant and positive associations between biodiversity and the rates of some individual micropollutant biotransformations, the associations do not share a shape. Some associations are better described with a nondecelerating shape, while others are better described with a decelerating shape (Fig. 2 and 4). Thus, for any individual micropollutant biotransformation, one cannot naively assume whether and how differences in biodiversity will translate into differences in the rate of that ecosystem function. This outcome is consistent with those observed for comparable investigations of other types of microbial communities (e.g., see references 19 and 20).

While we were unable to identify a common shape for the positive associations between biodiversity and the rates of individual micropollutant biotransformations, we found that there is a clear nondecelerating shape for the positive association between biodiversity and the collective rates of multiple micropollutant biotransformations (Fig. 3 and 4). The absence of apparent deceleration suggests that greater levels of biodiversity are likely required to maximize the collective rate of multiple micropollutant biotransformations than are typically required to maximize the rate of any individual micropollutant biotransformation. In other words, WWTPs are insufficiently biodiverse to saturate the set of possible complementarity and sampling effects when multiple micropollutant biotransformations are considered. Thus, engineering strategies that maximize biodiversity could be particularly effective when the preferred outcome of wastewater treatment is the removal of a large number of different micropollutants. While there are only few comparable investigations of associations between biodiversity and multifunctionality among microbial communities (e.g., see references 42 and 43), this outcome is consistent with the outcomes observed for comparable investigations of plant communities (34, 35).

An important question that emerges from our investigation is why the observed positive associations are stronger for some micropollutant biotransformations than for others (Fig. 2). We tested the hypothesis that the strength of the positive association depends, at least in part, on the rarity of the particular micropollutant biotransformation. As a micropollutant biotransformation becomes increasingly rare, the probability that the set of possible complementarity and sampling effects is saturated should decrease, thus leading to a stronger positive association. Our estimates of the rarity of each micropollutant biotransformation are consistent with this hypothesis and therefore provide a first test of one of its main predictions. Such a hypothesis could reflect a general rule that has important applications beyond wastewater treatment. For example, it could be used to make specific predictions about how changes or differences in biodiversity are likely to affect the rates of particular types of ecosystem functions.

An important aspect of our experimental design is that all of the investigated full-scale WWTPs had differences in their operational and environmental metrics (see Table S2 in the supplemental material). One consequence of this aspect is that we cannot unequivocally test whether the observed positive associations resulted from causal linkages between biodiversity and the rates of particular micropollutant biotransformations. Nevertheless, we have two lines of indirect evidence that the observed positive associations do indeed result from causal linkages between these two variables. First, the positive associations are stronger for rare micropollutant biotransformations and weaker for common micropollutant biotransformations. If biodiversity were causally linked with the rates of particular micropollutant biotransformations, we would expect this relationship (see the discussion above). If biodiversity were not causally linked with the rates of particular micropollutant biotransformations, we would not expect any relationship between the strengths of the positive associations and the rarity of the measured micropollutant biotransformations, which is inconsistent with our observations. Second, we previously could not explain differences in the rates of particular micropollutant biotransformations by differences in any of the operational metrics of the WWTPs (24). Thus, a causal linkage between biodiversity and the rates of particular micropollutant biotransformations provides a plausible and parsimonious explanation for our results.

Our experimental design also allows us to address an important and largely unexplored gap in our knowledge about biodiversity: can we translate biodiversity-ecosystem functioning relationships derived from carefully controlled laboratory experiments into real-world applications (12, 22)? Consider that, in the natural environment, environmental conditions can differ dramatically across small spatial and temporal scales. The consequence is that for any set of real-world microbial communities (including full-scale WWTP communities), there is inherent uncertainty about differences in the environmental conditions from which those microbial communities were collected. If these differences in environmental conditions have opposite effects on the rates of micropollutant biotransformations from those of biodiversity, then the expected positive associations might be obscured and biodiversity would consequently be of limited predictive value. Our results suggest that this is not the case. Even though the 10 WWTPs had differences in their operational metrics, we were nevertheless able to observe the expected positive associations between biodiversity and the rates of individual and multiple micropollutant biotransformations (Fig. 2 and 3). Our results therefore suggest that biodiversity-ecosystem functioning relationships can indeed be of utility for predicting how differences in biodiversity translate into differences in the rates of micropollutant biotransformations among real-world WWTP communities.

In conclusion, we report three main outcomes that improve our general knowledge about the ecology and functioning of full-scale WWTP communities. First, differences in biodiversity translate into differences in the rates of specific individual micropollutant biotransformations. Second, there is a nondecelerating shape for the positive associations between biodiversity and the collective rate of multiple micropollutant biotransformations. Biodiversity is therefore likely to be particularly important for the removal of a large number and wide variety of different micropollutants. Finally, we present evidence for a potential general rule that determines how strongly biodiversity is likely to associate with the rates of different types of micropollutant biotransformations, which could have applications that extend well beyond wastewater treatment.

## ACKNOWLEDGMENTS

We thank Tom Curtis, Eberhard Morgenroth, Johanna Otto, Jan Roelof van der Meer, George Wells, and Paul Wilmes for useful discussions and for critically reading preliminary versions of the manuscript.

This work was supported by grants from Eawag (category: Seed) and the Korean-Swiss Science and Technology Cooperation.

## FOOTNOTES

- Received 7 October 2014.
- Accepted 11 November 2014.
- Accepted manuscript posted online 14 November 2014.
Supplemental material for this article may be found at http://dx.doi.org/10.1128/AEM.03286-14.

- Copyright © 2015, American Society for Microbiology. All Rights Reserved.