Metagenomic analysis of an ecological wastewater treatment plant’s microbial communities and their potential to metabolize pharmaceuticals

Pharmaceuticals and other micropollutants have been detected in drinking water, groundwater, surface water, and soil around the world. Even in locations where wastewater treatment is required, they can be found in drinking water wells, municipal water supplies, and agricultural soils. It is clear conventional wastewater treatment technologies are not meeting the challenge of the mounting pressures on global freshwater supplies. Cost-effective ecological wastewater treatment technologies have been developed in response. To determine whether the removal of micropollutants in ecological wastewater treatment plants (WWTPs) is promoted by the plant-microbe interactions, as has been reported for other recalcitrant xenobiotics, biofilm microbial communities growing on the surfaces of plant roots were profiled by whole metagenome sequencing and compared to the microbial communities residing in the wastewater. In this study, the concentrations of pharmaceuticals and personal care products (PPCPs) were quantified in each treatment tank of the ecological WWTP treating human wastewater at a highway rest stop and visitor center in Vermont. The concentrations of detected PPCPs were substantially greater than values reported for conventional WWTPs likely due to onsite recirculation of wastewater. The greatest reductions in PPCPs concentrations were observed in the anoxic treatment tank where Bacilli dominated the biofilm community. Benzoate degradation was the most abundant xenobiotic metabolic category identified throughout the system. Collectively, the microbial communities residing in the wastewater were taxonomically and metabolically more diverse than the immersed plant root biofilm. However, greater heterogeneity and higher relative abundances of xenobiotic metabolism genes was observed for the root biofilm.

The treatment of human wastewater by ecological systems predates the advent of engineered wastewater treatment plants (WWTPs). While exposure to human pathogens is greatly reduced in communities where modern wastewater treatment technologies have been implemented 1 , widespread detection of micropollutants in the environment 2 raises serious concerns about the efficacy of modern WWTPs to treat this class of contaminants. Moreover, with about two-fifths of the world's population experience health effects due to poor sanitary conditions 3 . The emerging field of ecological engineering has provided a variety of viable, costeffective wastewater treatment designs 4 . The organizing principle of ecological wastewater treatment is the construction of "task oriented mesocosms" 5,6 of eutrophic ecosystems that, like conventional systems, primarily rely on microbial metabolic processes to achieve water quality goals. Where ecological WWTPs and conventionally engineered WWTPs differ significantly is in their reliance on ecological processes to assimilate nitrogen, phosphorous, and carbon from the wastewater into biomass. Whereas, conventional WWTPs utilize mechanically assisted, microbial processes to evolve gaseous C, N, and frequently chemical precipitation of phosphorous. A number of ecological systems have been operating around the world for decades including constructed wetlands 7,8 , Eco-Machines TM4 , and biofilters 9,10 for residential, industrial, and municipal wastewater. These systems perform reliably based on tertiary wastewater standards 4 , while reducing operational costs 11 and environmental and human health impacts of wastewater 12 . The unique potential provided by ecologically engineered waste management is direct conversion of a liability (i.e. wastewater) to an asset (sequestered carbon, biomass, products, biodiversity, etc.) 13 .
At the core of wastewater treatment is the biodegradation, oxidation, and reduction of organic macromolecules and inorganic chemical species primarily by resident microbial communities. While the microbial communities of conventional WWTPs have been thoroughly studied 14 , very little is known about microbial communities in existing ecological WWTPs despite the fact that they are central to the functions these systems provide 6 . The introduction of activated sludge and environmental media from diverse sources is thought to provide essential microbial functional groups 4,15 . It is not known whether these "seeding" events provide microbial functional groups with the capacity for biodegradation of micropollutants.
The promotion of microbial biodegradation of recalcitrant xenobiotic pollutants by plant roots has been well documented 16,17 . The "rhizosphere effect" 18-20 , driven by the release of plant metabolites from plant roots, accelerates microbial biodegradation of recalcitrant pollutants in soil and water 21,22 . In some cases, microbial biodegradation is promoted by a nonspecific increase in microbial metabolic activity in the area surrounding roots 23 , yet other studies have shown a relationship between specific plant metabolites and certain pollutant degrading organisms 24 . The interaction has been described as co-metabolic induction, or "co-metabolism", where metabolism for one compound is promoted in the presence of other compounds 21,22 . This phenomenon has been successfully employed to accelerate the removal of a variety of recalcitrant pollutants from the soil and water including polychlorinated biphenols (PCBs) 19,20,24 , polycyclic aromatic hydrocarbons 22 , and chlorinated solvents such as trichloroethylene 25 . However, little research has been done on whether co-metabolism is occurring in ecological WWTPs as a result of plant-microbe feedback processes. Using whole metagenome sequencing (WMS) we have examined whether the microbial populations residing on the plant roots immersed in wastewater of an ecological WWTP showed evidence of the capacity for micropollutantant biodegradation. These populations were compared to microbial communities free-floating in the wastewater, enrichment cultures growing on individual pharmaceutical compound carbon sources, as well as PPCP concentrations throughout the treatment system to determine whether plant-microbe feedback processes are supporting PPCP biodegradation in ecological WWTPs.
EcoMachine TM sampling An ecological WWTP (Eco-Machine TM ) in Sharon, Vermont at the Vietnam Memorial rest area on northbound Interstate 89 (43.727896, -72.425564) was sampled on June 30 and July 1, 2013. Wastewater from the toilets, urinals, and sinks, is collected in a holding tank and then treated in a series of tanks ( Figure 1). These consist of an anoxic tank (ANOX), a closed tank (CLO) and planted aerobic tanks (HR1, HR2 and HR3). These are followed by a clarifier and final treatment by a sand filter (SAND  40 was used to test statistical significance of differentially abundant taxonomic groups and functional categories for 1) WWTP sample groups (aqueous and biofilm phases) and 2) enrichment culture sample groups (carbamazepine (C), sulfamethoxazole (S), and trimethoprim (T)). LCA taxonomic profiles and KEGG including abundances, were imported to STAMP for each sample set. Two-sided Welch's t-test was used to compare aqueous and biofilm phases with a confidence interval of the effect size and multiple test correction using the Benjamini-Hochberg FDR method. One-way ANOVA was used to compare enrichment culture groups with an effect size (Eta-squared) and multiple test correction using the Benjamini-Hochberg FDR method. Tukey-Kramer post-hoc test (0.95) was used to determine which means were significantly different when an ANOVA produced a significant p-value.
To visualize the distribution of microbial taxa in WWTP samples the Circos software package v0.69 41 was used to depict the location and relative abundances of microbial taxa at the class rank identified in MEGAN using the LCA algorithm.

PPCPs concentrations
There were 11,568 visitors during the week in which sampling was conducted (June 25-July 1, 2013). Each visitor used an average of 2.27 liters of water contributing 26,452 liters of water to the wastewater treatment system 42 . The wastewater used to isolate microbial DNA samples contained detectable concentrations of caffeine, carbamazepine, DEET, gemfibrozil, ibuprofen, naproxen, sulfamethoxazole, thiabendazole, and trimethoprim. Of the compounds detected in influent water (from facility toilets, urinals, and sink drains), the cholesterol medication gemfibrozil was detected at the highest concentration (1.5 × 10 5 ng L -1 ), followed by caffeine, ibuprofen, and naproxen (9.5, 6.6, and 5.5 × 10 4 ng L -1 , respectively) ( Table 1). The concentrations of PPCPs in the wastewater samples generally decreased the further through the treatment process ( Figure 1) the sample was obtained. However, gemfibrozil, caffeine, and ibuprofen were detected at higher concentrations in the sand filter or effluent water samples than the preceding tank.
The concentrations of carbamazepine, DEET, and trimethoprim did not change substantially over the entire treatment process.

Metagenome sequencing and sequence processing
Whole metagenome shotgun sequencing of 12 WWTP samples generated more than 388 million paired-end reads, 101 bp in length, with an average depth of 32.4 million reads per sample (range: 2.35-53.3 million) (Supplementary Material ST13). Eighty-eight percent of raw reads (343,355,560) were retained after qualitytrimming and were aligned to the NCBI-NR protein database. Of the 175,238,945 reads with at least one hit to NR proteins, approximately 84% were assigned taxonomy by the lowest common ancestor (LCA) algorithm in MEGAN. Over half of quality-trimmed reads in our samples (51%) had no protein hits in NCBI-NR and 19.3% of reads with protein hits could not be classified by the LCA algorithm. As a result, the latter were designated "Not Assigned" reads in MEGAN.
Whole metagenome shotgun sequencing of the 8 enrichment culture samples generated more than 177 million paired-end reads, 101 bp in length, with an average depth of 22 million reads per sample (range: 17-29 million) (Supplementary Material ST14). Eighty-nine percent of raw reads (157,598,266) were retained after quality-trimming and were aligned to the NCBI-NR protein database. Of the 96,824,600 reads with at least one hit to NR proteins, over 99% were assigned taxonomy by the LCA algorithm in MEGAN. Nearly 39% of quality-trimmed reads in our samples (60,773,666) had no protein hits in NCBI-NR and 0.4% of reads with protein hits could not be classified by the LCA algorithm and were designated "Not Assigned" reads in MEGAN. The minimum-support percent threshold in MEGAN for both WWTP and enrichment culture analyses was set to 1.0% based on our bioinformatics workflow results from the HMP and Delaware River control samples.

Taxonomic classification of sequences
The LCA algorithm provided a microbial taxonomic profile of the 12 WWTP samples and the 8 enrichment culture samples ( Figure 2 and Figure 3, respectively). The read counts were normalized to  the sample with the least number of total reads to allow relative abundance to be depicted and shown as bar-charts at the leaves of each phylogram. The terminal taxa from various ranks identified by the LCA algorithm in the 12 WWTP samples ranged from a low of 8 taxa (HR1_B sample) to 17 taxa (ANOX_W sample) ( Figure 2; Supplementary Material ST2). Members of Mycobacterium, Pseudomonas, and Verrucomicrobia were identified in all of the aqueous phase (_W) samples as well as the biofilm in the anoxic tank (ANOX_B) and the sand filter (SAND_B). Additionally, the families Rhodobacteraceae, Burkholderiaceae, Comamonadaceae and Xanthomonadaceae were identified in all of the aqueous phase samples as well as in at least one, but not all immersed biofilm samples. The family Xanthomonadaceae was identified in 10 of the 12 WWTP samples, followed by Rhodobacteraceae (9), Comamonadaceae and Rhizobiaceae (7), and Burkholderiaceae (6).
According to the LCA algorithm taxonomic assignments by MEGAN, the enrichment cultures originating from the wastewater effluent inoculant produced mixed cultures ranging from 2 (T3B and T3C) to 8 (S3D) taxa identified in each culture (Figure 3; Supplementary Material ST4). Bacteria were identified in all enrichment cultures and ascomycete fungi were in all carbamazepine cultures and one trimethoprim culture (T3D A graphical representation of the pairwise distance matrix 39 generated using normalized Goodall's similarity index 38 of the 12 WWTP samples is shown as an unrooted phylogenetic neighbor network in Figure 4. The LCA taxonomic assignments of the microbial populations residing in the aqueous phase samples cluster near one another in the neighbor network, while the immersed biofilm samples showed greater dissimilarity. In STAMP, the aqueous and biofilm samples were compared at different taxonomic ranks to highlight differences between the WWTP sample sets. Figure 5 depicts the percent relative abundances of microbial classes identified in the 12 WWTP samples. In general, Actinobacteria, Alphaproteobacteria, Betaproteobacteria, and Gammaproteobacteria had higher relative abundances in the aqueous samples, whereas Bacilli were more abundant in the biofilm samples. The relative abundances of Bacilli in the biofilm samples were highest in the first three tanks (ANOX, CLO, HR1) and sharply reduced thereafter (Figure 5a, Figure 9). The relative abundances of Deltaproteobacteria, Betaproteobacteria, Gammaproteobacteria, and Bacilli were significantly different between the two physical phases (aqueous and biofilm) sample sets (Figure 5b).
To assess similarities in taxonomic composition between the pharmaceutical compound enrichment cultures and the WWTP samples, we examined the co-occurrence of LCA taxa among datasets ( Figure 6). The taxonomic similarities of the five major sample types were examined by combining all aqueous samples, all immersed biofilm samples, and the replicates of each pharmaceutical carbon source enrichment cultures into five datasets. As is typical of enrichment cultures originating from complex environmental samples, taxonomic similarities were limited to a few shared taxa. When combined, twelve microbial taxa were identified in the carbamazepine enrichment cultures. However, only three taxa were identified in at least two out of three replicate carbamazepine cultures. Of the twelve taxa identified in the combined carbamazepine cultures, only one each was identified in the biofilm and aqueous samples (Bacillus cereus, Mycobacterium spp., respectively) obtained from the wastewater treatment plant. Bacillus cereus was identified in the biofilm samples obtained from the first three tanks of the wastewater treatment system (ANOX_B, CLO_B, and HR1_B). It has been associated with human feces 43   To assess the influence of the carbon sources on enrichment culture taxonomic composition, LCA-based taxonomic profiles of the enrichment culture replicates were compared in STAMP. The most numerous significant differences with strong effect sizes were found at the family level (see supplementary material ST35). The relative abundance of reads from Bradyrhizobiaceae, Hyphomicrobiaceae, Microbacteriaceae, and Rhizobiaceae were significantly greater in sulfamethoxazole cultures than in trimethoprim or carbamazepine cultures (with FDR corrected p-values of 2.98 × 10 -8 , 0.0274, 0.0199, 0.381, respectively).

Microbial metabolism of PPCPs
To investigate the role of the WWTP's microbial communities in the metabolism of PPCPs, KEGG was queried using MEGAN for sequences identified as involved in xenobiotic metabolism. Figure 7 depicts the percent relative abundances of sequences identified in the metagenomes obtained from each sample as involved in xenobiotic biodegradation and metabolism KEGG pathways for WWTP samples. The most abundant xenobiotic biodegradation and metabolism subcategory in all cultures and WWTP samples, with a few exceptions (HR1_B, C3B and C3D), was benzoate degradation comprising 14.5%-17.7% of reads in this subcategory for WWTP samples and 11.1%-23.6% in the enrichment culture samples. The biofilm sample obtained from the sand filter (Sand_B) at the end of the WWTP treatment train produced the greatest number of sequences (18,766) aligning to the "xenobiotic biodegradation and metabolism" category in the KEGG database followed by the aqueous samples of the sand filter (SAND_W) (18,318), and the aqueous phase of the third planted aerobic tank (HR3_W) (18,002). The sample with the greatest number of reads associated with KEGG pathway "drug metabolism -cytochrome P450" was the aqueous phase of the third planted aerobic tank (HR3_W) (1,384) followed by aqueous samples of the sand filter (SAND_W) (1,372) and the immersed biofilm sample obtained from the first planted aerobic tank (HR1_W) (1,330). The sample with the greatest number of reads associated with the KEGG category "drug metabolism -other enzymes" was the biofilm in the second planted aerobic tank (HR2_B) (2,181) followed by aqueous phase of the anoxic tank (ANOX_W) (1,904), and the biofilm sampled from the sand filter (SAND_B) (1,827).
The percent relative abundances of xenobiotic metabolism-associated sequences identified in enrichment cultures' metagenomes are illustrated as a heat map in Figure 8 (see Supplementary materials ST6). Reads associated with benzoate degradation were most abundant for trimethoprim cultures, whereas chloroalkane and chloroalkene degradation reads were most abundant  for sulfamethoxazole cultures. The most abundant xenobiotic metabolism category for the carbamazepine cultures varied from culture to culture with aminobenzoate degradation, benzoate degradation, and chloroalkane/chloroalkene degradation categories for each of the three replicates. The samples C3A and T3D produced the greatest number of reads associated collectively with xenobiotic metabolism genes at 232, 201 and 218, 591, respectively. Of these, C3A and T3D metagenomes contained 10,151 and 11,889 reads associated with "drug metabolism -other enzymes" in KEGG, respectively. In this category, the sulfamethoxazole culture (S3B and S3D) had the most assigned reads at 16,102 and 13,730, respectively. For the category "Drug metabolism -cytochrome P450" the carbamazepine culture sample C3A contained the greatest relative number of sequences (16,199) followed by the sulfamethoxazole culture S3D (11,847).

Discussion
To our knowledge, this is the first study to use WMS to characterize the microbial communities of an ecological WWTP. While a taxonomic description of the microbial communities is provided, here we have focused on microbial metabolism of PPCPs in the wastewater. Using WMS we have identified representation and relative abundance of microorganisms in the six major treatment tanks. To examine the role plants have on the structure of the microbial communities, we compared the communities found in the wastewater and those attached to plant roots immersed in the wastewater. Microbial metabolic pathways for most emerging pollutants, including micropollutants such as PPCPs, have not been characterized. Therefore, we have focused on the xenobiotic metabolic capacity as represented by the location and abundance of known genes represented in the KEGG database.

PPCPs in the wastewater
The concentrations of PPCPs in the samples obtained from the Sharon, VT WWTP indicate that some of the PPCPs are being effectively removed by the system (naproxen and thiabendazole), while others are accumulating in the recirculating system (caffeine, carbamazepine, DEET, trimethoprim, and sulfamethoxazole) ( Table 1). The increasing concentrations of caffeine, gemfribozil and ibuprofen in or after the sand filter suggest partitioning to the aqueous phase from attenuated organic matter may be occurring in the sand filter. However, based on these data, biotic (biodegradation) and abiotic (partitioning to the primary sewage sludge) processes occurring in the influent holding tank and first treatment tank (ANOX) account for significant reductions in concentrations for the PPCPs gemfibrozil, naproxen, thiabendazole, and to a lesser degree caffeine, ibuprofen, sulfamethoxazole, whereas carbamazepine, DEET, and trimethoprim were not removed.
Partitioning to the solid phase (sewage sludge) is an important aqueous phase removal processes that is driven by a compound's hydrophobicity 44,45 . For certain chlorinated organic compounds, the octanol-water partitioning coefficient (KOW) correlates positively with sorption to biosolids when log KOW values range from 1.26 to 5.48 46 . Carbamazepine (log K OW 2.3-2.5), DEET (log KOW 2.18-2.44), and trimethoprim (log KOW 0.9-1.4) show limited removal in the Sharon, VT WWTP (Table 1). These findings are consistent with similar studies on partitioning of PPCP in conventional WWTP [47][48][49] . Partitioning rates from aqueous influent to biosolids (sewage sludge) is variable with log K d (solid-water distribution) values ranging from <0.7 to 4.2 50 depending on the compound. The removal of the parasiticide, fungicide thiabendazole (log K OW 2.47) from the wastewater between the influent sample location and the anoxic treatment tank was likely due to strong partitioning to the primary sludge in the holding tank. With this exception, all other compounds were detected in more than one tank of the system. The aqueous concentrations of the cholesterol-lowering drug gemfribozil (log K OW 4.77) also indicated strong partitioning to the primary sludge showing reduction in concentration by three orders of magnitude between the influent and the anoxic treatment tank.
Pharmaceuticals and PPCPs in wastewater can undergo a number of processes that contribute to their complete or partial aqueous phase removal in wastewater treatment systems 45 . These include chemical and or physical processes such as sorption to organic matter 46,48,51 , photolysis 52,53 , volatilization 50 , and biological transformation 54 . Biological transformation is unique among these as it provides WWTP operators the potential to increase the removal of PPCPs from wastewater while partially or completely mineralizing the compounds thereby eliminating any risks associated with their release to the environment. In contrast, sorption to biomass (primarily sewage sludge) results in decreased mineralization 55 , which when applied to land (dominant disposal method of the processed sewage sludge) is likely a significant source of PPCPs in the environment 56 . It is unclear whether thermal treatment and dewatering, as is commonly done to biosolids prior to land application, alters the mass of PPCPs in this media.
The influent concentrations of the PPCPs were in many cases an order of magnitude or greater than the concentrations reported in conventional WWTP 49,57 . The Sharon, VT ecological WWTP recirculates the effluent onsite as flush water (sterilized and dyed blue prior to being used as toilet and urinal flush-water). The recirculation and reuse of the effluent is likely to result in an additive or concentrating effect for compounds with low removal and/or partitioning rates. Additionally, the concentration of PPCPs in municipal wastewater is likely diluted by the mixing of non-human wastewater such as wash water, storm water (in combined sewer systems), industrial process water and a variety of other sources. While the higher concentration of PPCPs in this recirculating ecological WWTP may present elevated exposure risks to operators and the environment if materials are discharged, conventional WWTP, which do not recirculate wastewater, are likely to discharge greater mass of PPCPs per liter wastewater treated. Additionally, retaining the PPCPs in the WWTP through recirculation is preferential to releasing them into receiving water bodies.
The results reported here are initial findings on the removal of PPCPs from the wastewater processed by this system as the species and concentrations of detected PPCPs are likely to change with time, fluctuating with the changing population of visitors. Significant variability in PPCPs species and mass loading into the system is likely responsible for the non-detects (Table 1 ND's) of individual compounds detected later in the treatment train. As we staggered our sampling of the initial three and latter three treatment tanks by 48 hours to accommodate the residence time of the wastewater in the treatment system, it is reasonable to assume that the concentrations quantified here represent the flux of PPCPs through the WWTP. Therefore, changes in the concentration of individual PPCPs throughout the system are the result of biotic and abiotic aqueous phase removal processes. The trends observed here could be influenced by fluctuations in PPCPs inputs. These fluctuations are likely responsible for the non-detects observed in the holding tank (INF) for carbamazepine, DEET and trimethoprim, while these compounds were detected in the next treatment tank (ANOX).

Microbial communities of the ecological WWTP
The immersed biofilm and aqueous phase microbial communities exhibited two distinct taxonomic structures. According to the LCA algorithm, this difference was most evident at the class level (Figure 5a). The relative abundances of Deltaproteobacteria, Betaproteobacteria, and Gammaproteobacteria were significantly higher in the microbial metagenomes of the aqueous phase, while Bacilli were observed in greater abundance in the immersed biofilm microbial communities (Figure 5b). The relative abundance of Bacilli in the immersed biofilm communities was highly variable, with their dominance diminishing significantly after the first three tanks. Figure 9 illustrates the microbial taxonomic spatial variability of the ecological WWTP. The sample order in the diagram reflects wastewater movement through the treatment system starting with the anoxic tank (ANOX) moving clockwise to the sand filter (SAND). The width of each ribbon extending from the taxon to the sample represents relative abundance based on sequence counts. Bacilli dominated the immersed biofilm in the first three tanks (ANOX, CLO, and HR1), with representatives of other taxa present to a much lesser degree. The immersed biofilm samples also show greater microbial taxonomic richness in the latter phases of the treatment system (HR2, HR3, and SAND). This pattern is likely due to retention of fecal taxa by the earlier tanks as Bacilli are known to be abundant in human feces 43 . Bacteroidia, which is also abundant in human feces 58 , was not detected colonizing the immersed plant root surfaces. Only a small population was identified in the wastewater obtained from the anoxic tank.
Members of the Firmicutes (Bacilli) and Bacterioidia, which are common in human feces 58 , were identified in the ecological WWTP microbial communities. However, these were the only organisms identified in the samples that are associated with human feces. Of the organisms used in water quality criteria to indicate contamination by feces (others include Clostridium perfringens, enterococci, Escherichia coli, and fecal coliforms), Clostridium perfringens was detected in only one sample, HR1_B (Figure 2). This would indicate that the ecological WWTP is performing well with regard to attenuating fecal coliforms. Bacilli were dominant in the samples of immersed biofilm collected from the first three tanks in contrast to the taxonomic composition of the latter three treatment tanks (Figure 2, Figure 9). This could indicate colonization of the surfaces in these initial treatment tanks by organisms of fecal origin. The wastewater samples from these tanks as well as biofilm samples collected from tanks later in the treatment train did not show this pattern (Figure 9).
The microbial taxonomic composition of the immersed biofilm in the downstream tanks increased in richness with Actinobacteria, Alphaproteobacteria, Bacilli, Gammaproteobacteria, Sphingobacteria and to a lesser degree Betaproteobacteria, and Planctomycetia identified by LCA. This is in contrast to the aqueous phase microbial community, which showed taxonomically rich populations throughout the system with members of the Actinobacteria, Betaproteobacteria, Alphaproteobacteria, and Gammaproteobacteria in greatest abundance (Figure 9).
The taxonomic composition of the microbial communities forming biofilm on the plant root surfaces immersed in the wastewater changed significantly over the course of the wastewater treatment process. The taxonomic dissimilarity (Figure 4) observed among the immersed plant root biofilm samples showed three distinct communities: The anoxic (ANOX_B) and the closed aerobic (CLO_B); the first planted aerobic (HR1_B); the second and third planted aerobic tanks (HR2_B and HR3_B); and the sand filter (SAND_B). The changing characteristics of the first tanks of the treatment system (anoxic to aerobic conditions) is likely driving a transition from anaerobes such as Bacteriodes, and facultative anaerobes or microaerophiles such as Bacillus, Clostridium, Nocardia, Mycobacterium, to aerobic communities. The abundance of the facultative anaerobe group, Bacilli is likely promoted by the oxygen limiting environment in the first two tanks 43 . The immersed biofilm microbial community of the first planted aerobic tank (HR1) is taxonomically distinct from HR2 and HR3 despite identical physiochemical conditions. This pattern may be influenced by the composition of plant species in the individual tanks, which varies from tank to tank, or perhaps the diminishing influence of organisms contributed by feces.
There is limited ability to relate the results of our analysis of the microbial communities of this ecological WWTP to that of other systems serving different populations and geographic locations as this appears to be the first published findings on the topic. However, metagenomic analyses of the microbial ecology of conventional (activated sludge) wastewater systems have been described elsewhere [59][60][61][62] .
For example, Lee et al., 2014 59 employed 16S rRNA gene microarrays (PhyloChip) to establish a baseline microbial community structure of the municipal WWTP aeration basin. The microbial taxonomic composition of the aeration basin showed some similarities with that of the entire ecological WWTP sampled here. Specifically, Proteobacteria and Firmicutes were abundant in both the conventional system and ecological WWTP sampled here. To assess the microbial community seasonal variation of activated sludge over a four-year period Ju et al., 2014 63 employed WMS, as was used here. They showed variation in microbial taxonomic composition between the summer and winter samples. They also found variation in the microbial community composition over the four years sampled, irrespective of season. The metabolic structure of activated sludge according to SEED subgroups appeared to remain stable, in spite of variation in taxonomic composition, which suggests microbial community functional redundancy may be present in these systems. The ecological WWTP sampled for this work is housed in a climate-controlled glass house, which raises questions as to whether the microbial community varies from year to year and season to season. The microbial communities are not likely to exhibit temperature-dependent seasonal variation. Taxonomic variation may occur as a result of changes in the microbial communities contributed by the visitors.
The role plants have on influencing the structure of the rootcolonizing microbial communities appeared to increase after the first three treatment tanks. The attenuation of taxa contributed by feces (Bacillus) after the first three treatment tanks is reflected in the increased microbial taxonomic heterogeneity in the latter two planted aerobic tanks (HR2 and HR3), which is reflected in the branch lengths for these samples in the neighbor joining network shown in Figure 4. These results indicate the need to design ecological WWTPs with sufficient retention time to allow for the attenuation of stool microbial communities and the development of diverse microbial biofilm communities.

Removal of PPCPs by the ecological WWTP
The dominant PPCP removal process from the wastewater appears to be partitioning to sludge (biosolids) and biodegradation under nitrifying conditions, which are both reflected in the reductions in aqueous PPCP concentrations that occurred early in the treatment process. Primary sludge settles out of the wastewater in the holding tank and is periodically removed for off-site disposal. The concentration of some of the detected PPCPs continued to decline as the wastewater continued through the system (Table 1) indicating some continued removal beyond the first two or three tanks. For example, while the concentration of caffeine declined between the holding tank (9.5 × 10 4 ng L -1 ) and the anoxic tank (1.9 × 10 4 ng L -1 ), further reduction from the aqueous phase was observed in the subsequent three aerobic treatment tanks (CLO, HR1, HR2). Given that caffeine is a hydrophilic organic base (low KOW) only moderate partitioning to sludge is expected (86) and microbial biodegradation is likely to be responsible for the reduction in caffeine concentrations observed from the aerobic treatment tanks. Ibuprofen concentrations followed a pattern similar to caffeine's in that significant reductions were seen in the first three aerobic treatment tanks (1.1 × 10 4 ng L -1 to 5.6 × 10 2 ng L -1 ). Carbamazepine, DEET, and trimethoprim concentrations remained stable throughout the treatment process. The combination of low partition to primary sludge expected and metabolic recalcitrance accounts for their stability in system 47,50,55 .
Microbial biodegradation pathways for most PPCPs have not been characterized, which makes it difficult to directly detect responsible genes. Nevertheless, ammonia oxidizing bacteria have been associated with the biodegradation of some PPCPs, while others have been shown to be degraded by nitrite oxidizing bacteria 64,65 . The relative abundances of sequences that MEGAN associated with ammonia monooxygenases were very low throughout the system (ranging from 0 to 53 reads), but were highest in the biofilm samples obtained from the HR3 and HR2 tanks. Due to lower dissolved oxygen levels, genes involved in denitrification (nitrite reductases) were found to be in highest relative abundance in the anoxic (ANOX) and closed (CLO) tanks (644 and 688, respectively).
Functional attributes of detected taxa reported in the literature can be used to identify metabolic potential pertinent to uncharacterized xenobiotic metabolic pathways. For example, in the first three tanks, the Firmicutes colonizing plant root surfaces have been reported to metabolize xenobiotics. Bacillus cereus, B. megaterium and B. amyloliquefaciens have been reported to metabolize phenol 66 , crude oil 67 , textile dyes 68 , and other xenobiotics through the induction of cytochrome P450s 69,70 . Dehalobacter sp. FTH1, identified in the plant root biofilm sample obtained from the second planted aerobic tank (HR2), has been reported to dechlorinate a number of organohalide xenobiotics 71,72 . Clostridium, identified in the root biofilm sample obtained from the first aerobic tank (HR1), has been reported to be involved in metabolism of bromophenols as a member of a consortium including Delhalobater 73 . Entercoccus spp., identified in the plant root biofilm of HR2, has been reported to degrade azo dyes 74 .
Of the Actinomycetales identified, Mycobacterium spp., which have been reported to metabolize a variety of xenobiotics including polycyclic aromatic hydrocarbons 75 , biphenyls 76 , as well as various pharmaceuticals 77 were identified in low relative abundance in samples obtained from the biofilm growing on plant roots in the anoxic tank and in high relative abundance in the biofilm sampled from sand filter. Mycobacterium spp. was also identified in the aqueous wastewater throughout the system (Figure 2 &  Figure 9). While this metabolically plastic genus has been reported to be capable of metabolizing a wide variety of xenobiotics it should be noted that there are numerous pathogenic taxa including M. tuberculosis, M. bovis, and M. avium among others. The biofilm sample obtained from the sand filter contained 10,282 reads associated with human diseases and 18,766 reads associated with xenobiotic metabolism by the LCA algorithm in MEGAN.
Rhizobiales, which were identified in the aqueous phase throughout the system as well as in the biofilm sampled in the latter three treatment tanks (HR2, HR3, and SAND) have been reported as abundant in biofilm reactors treating sulfamethoxazole containing wastewater 78 . Also present throughout the system's aqueous phase were members of the Rhodobacteraceae, Burkholderiaceae, Comamonadaceae, Pseudomonas, and Xanthomonadaceae, which have been reported to metabolize aromatic hydrocarbons 79,80 . Among these, the genus Pseudomonas has been identified as capable of biodegradation a variety of xenobiotic including some pharmaceuticals in a number of other settings including membrane bioreactors 81 , cultures originating from pharmaceutical wastewaters 82 , and environmental samples 83 .
Greater xenobiotic metabolic heterogeneity was observed in the samples obtained from the plant root-associated biofilm as compared to the free-floating aqueous microbial community. The aqueous microbial metagenome, collectively, contained a greater total for xenobiotic metabolism gene copies (1.2 × 10 5 compared to 8.6 × 10 4 for the plant root biofilm) (see Supplementary material ST5). When comparing the proportion of sequences identified in the aqueous and biofilm phases of the system, which represent the non-root-associated and root-associated microbial populations, respectively, the xenobiotic metabolism gene categories nitrotoluene, benzoate, flurobenzoate and steroid degradation were found to be significantly higher in the aqueous phase samples (Welch's t-test p-values < 0.05) (see Supplementary material ST49) (Figure 7). However, when comparing the type and abundance of reads associated with xenobiotic metabolism (KEGG level 3) the aqueous phase samples all resembled one another whereas, the biofilm samples were heterogeneous (see Figure 7 and Supplementary material ST31).
For aqueous phase microbial communities, growth of microbial consortia with the capacity to metabolize a given PPCP is driven by the presence of the compound in the wastewater, which will fluctuate with time. Stationary biofilm communities are likely to be more stable populations. These communities can accumulate with time and potentially acquire metabolic genes by horizontal gene transfer 84 . In contrast to the aqueous phase xenobiotic metabolism, which was dominated by benzoate degradation genes, the plant root microbial biofilm metagenomes contained higher relative abundances of other categories including aminobenzoate, chloroalkane/chloroalkene, and to a lesser degree ethylbenzene categories (Figure 7). However, of these categories, only chloroalakane/ chloroalkene was significantly higher (Welches t-test p-value 0.031) in the biofilm samples collectively (see Supplementary material ST49).
While others have identified ammonia oxygenases and nitrite reductases as being involved in microbial PPCP biodegradation 64 , the relative abundances of these gene categories were very low throughout the system. The relative abundance of xenobiotic metabolism gene copies was highest for the sand filter samples at 18,318 for SAND_W and 18,766 for SAND_B (see Supplementary material ST30). The sand filter's ability to attenuate and accumulate sloughed off microbial cells as the wastewater passes through may be driving an accumulation of microbial biomass. If this is the case, sand filters are likely to have populations of the microbial communities found throughout the aqueous phase of the system, yet may not serve as a location of high metabolic activity, thereby contributing little to the metabolism of xenobiotics. The increase in concentrations observed for some of the PPCPs (caffeine, carbamazepine, DEET, gemfribozil, and ibuprofen) after the wastewater passed through the sand filter, if a real trend, could support this perspective.
Culture bias was reflected in the taxonomic composition of the enrichment cultures growing on the sole pharmaceutical carbon sources examined here. While culture bias is well known 85 and was expected, the enrichment of organisms capable of metabolizing individual pharmaceutical compounds from the WWTP effluent water reflects the ability for the ecological WWTP to support this metabolic capability. Only one taxon (Bacillus cereus) was identified in at least one enrichment culture (carbamazepine) and the biofilm sampled from the WWTP and one taxon (Mycobacterium spp.) was identified in the aqueous, biofilm and carbamazepine enrichment cultures (Figure 6). Given the selective pressure supplied by culturing, it is unlikely that these two taxa are solely responsible for the biodegradation of carbamazepine in the ecological WWTP. However, having isolated pharmaceutical metabolizing consortia from effluent water suggests the ecological WWTP supports microbial populations with the capacity remove recalcitrant micropollutants from wastewater.
Secondary plant metabolites contributed to the wastewater by the tropical species cultivated in planted tanks may help support populations of microbial communities with xenobiotic degradation capabilities by providing a more consistent supply of structurally diverse carbon sources. The most abundant category of xenobiotic metabolism genes was associated with benzoate degradation in nearly all the samples. Benzoate degradation is known to play a role in the degradation of a variety of aromatic compounds 23,86,87 . The dominance of this functional category suggests that microbial communities in the aqueous phase were metabolizing a variety of benzoate-containing compounds, which were likely to include metabolites of plant and xenobiotic origin. However, further work is needed to determine whether plant-microbe feedback processes promote PPCP biodegradation in ecological WWTPs.
Author contributions IB conceived of the study, designed the experiments and was awarded the funding that supported the work. HD contributed to the statistical and bioinformatics analysis. JV developed and refined the bioinformatics workflow as well as contributed bioinformatics analysis. ML created the enrichment cultures as a part of her senior thesis project. IB prepared the first full draft of the manuscript. IB and HD prepared the manuscript revisions. JV was involved in the revisions of the draft manuscript. All authors have contributed to and agreed to the final content of the manuscript.