AGS sequencing batch reactor (SBR) operation is detailed in Bodle et al., 2022. In brief, AGS was grown in two identical glass SBRs with a working volume of 3.4 L. Both SBRs were operated in repeating three-hour cycles: 72 minutes anaerobic feed, 100 minutes aeration at a gas flow rate of 5 L/min, three minutes settling, and five minutes effluent discharge. The hydraulic residence time in both reactors was approximately 6.4 hours, and the solids residence time was controlled at approximately 25 ± 5 days. SBR operating parameters are consistent with those in other lab-scale studies (e.g., Kong et al., 2015; Margot et al., 2015; Rodriguez-Sanchez et al., 2017). During aeration, pH and dissolved oxygen were controlled with LabVIEW software (National Instruments) at 7.0 ± 0.3 and 1.75 ± 0.25 mg/L, respectively. Influent media were identical to those described by de Kreuk and van Loosdrecht, 2004, except that influent sodium acetate was increased to 10.3 mM, resulting in an organic loading rate of 2.5 g C/L*d. Both SBRs were initially seeded with AGS from an AquaNereda^® treatment plant in Utrecht, The Netherlands and operated at steady state (i.e., complete nitrogen and phosphate removal) for over 300 days prior to beginning experimentation. Immediately before starting experimentation, both SBRs were emptied and AGS was combined, mixed, and redistributed so that granule qualities would be as similar as possible in both reactors.

For 80 days, the test reactor received 46 mL of pharmaceutical media with the influent medium, resulting in an influent concentration of approximately 150 µg/L of each pharmaceutical. Individual stock solutions of each pharmaceutical were prepared first in methanol at 10 g/L, then diluted into nanopure water to obtain two pharmaceutical media: one solution (“DG”) containing 17.86 mg/L each of diclofenac sodium (Acros Organics) and gemfibrozil (Acros Organics), and a second solution (“ERY”) containing 35.8 mg/L erythromycin (TCI Chemicals). Twenty-three mL of each solution (46 mL total) were delivered with the influent medium throughout the anaerobic feed period. Approximately 0.13 mg/L methanol was also present in the combined influent stream. The pharmaceutical media were protected from light to prevent photolytic degradation and prepared fresh every 8-10 days.

Results discussed in Bodle et al., 2022 showed that the pharmaceuticals tested herein sorbed to different lab materials with different affinities; therefore, to ensure the accuracy of pharmaceutical dosing, ERY solution was pumped into the reactor using silicone tubing (Masterflex). DG solution was pumped into the reactor using PharmaPure tubing (Masterflex). Despite different pharmaceutical stock solution concentrations, sorption to tubing resulted in fairly stable influent concentrations of approximately 150 µg/L of each pharmaceutical. Although each pharmaceutical is typically measured in wastewater treatment plant influents at approximately 1-10 µg/L (Fang et al., 2012; Schafhauser et al., 2018; Sathishkumar et al., 2020), pharmaceutical sorption to tubing drove usage of this elevated influent concentration: 150 µg/L was the minimum concentration of pharmaceuticals that could be consistently dosed to the test reactor without significant losses from sorption. Influent samples were taken from a sampling port in the tubing located at the base of the reactor (reactor schematic available in Supplementary Information Figure S1 ) and were extracted and quantified per methods detailed in Section 2.3.

Influent and effluent samples from both reactors were regularly taken and filtered through 0.45 µm regenerated cellulose filters prior to analyses for ammonia, nitrite, nitrate, phosphate, and dissolved organic carbon. Ammonia was quantified with Hach kit TNT 832 and a Hach DR 3900 spectrophotometer, equivalent to US EPA Method 350.1. Other anions listed were quantified with a Dionex ICS-1100 anion chromatography system equipped with an IonPac AS22 RFIC column (4 x 250 mm) and IonPac guard column (4 x 50 mm). Dissolved organic carbon (DOC), defined as that which could pass through a 0.45 µm filter, was measured with a SKALAR Formacs^HT Total Organic Carbon analyzer system.

Aqueous influent and effluent samples were prepared for mass spectrometry (MS) analyses by solid phase extraction (SPE) as detailed in Bodle et al., 2022. It is important to note that the periodic flow conditions that are characteristic of SBR systems allow constant sorption and desorption of pharmaceuticals within influent tubing. For this reason, the most accurate method of quantifying influent concentrations would have required extraction of the entire influent volume, which would been too destructive to reactor operation. For that reason, 100 mL influent or effluent sample were filtered with a 1.5 µm glass fiber filter (Hach) to remove solids and loaded on preconditioned Waters Oasis HLB cartridges (30 mg, 20 mL) at 10 mL/min using a vacuum manifold system. Loaded cartridges were washed, dried, and frozen at -18°C until elution (no longer than 14 days). Pharmaceutical losses to the lab materials used during extraction (glassware, syringes, pipette tips, and glass fiber filters) were tested for and were found to be minimal (Bodle et al., 2022).

Influent and effluent pharmaceutical samples were periodically taken in duplicate to account for possible pharmaceutical losses during the extraction process, as detailed in Bodle et al., 2022. In brief, one sample from each duplicate set was pre-spiked with pharmaceutical stock solution to a final nominal concentration of 100 µg/L prior to extraction. After extraction, the unspiked sample was split in half, and one half was post-spiked to a final nominal concentration equal to 100 µg/L multiplied by each samples’ concentration factor. Recovery was then determined as follows:

Measured pharmaceutical concentrations in unspiked samples were thereby corrected for recovery of each analyte. Recovery was 97 ± 9% for DCF, 117 ± 17% for ERY, and 98 ± 11% for GEM ( Table 1 ). Triplicate influent and effluent samples were taken once per month to assess accuracy. Pharmaceutical quantification and detection were performed with an Acquity I Class Plus ultra-performance liquid chromatograph (UPLC) coupled to a Waters Synapt XS quadrupole time-of-flight mass spectrometer (QToF-MS) in positive ion mode. Chromatographic analysis methods were adapted from Sodré and Sampaio, 2020. In brief, pharmaceuticals were separated with an Agilent Eclipse Plus C18 column (2.1 x 100 mm) at 30°C under gradient elution. Formic acid-enriched methanol and ultrapure water (0.1% v/v) were used as mobile phases at 0.7 mL/min, and the concentration of organic solvent was increased from 1% initially to 95% at 10.1 minutes, then returned to initial conditions and re-equilibrated for 3 minutes (total run length of 13 minutes). Sample injection volume was 8 µL. Retention times are summarized in Table 1 .

A calibration curve relating each pharmaceuticals’ peak area with its nominal concentration in analytical standards was used to determine parent compounds’ concentrations. Standards were prepared in both methanol and water. The instrument limits of quantification were under 10 µg/L for all compounds. The open-source software MZmine and R were used to analyze and compile mass spectra data (RStudioTeam, 2009; Pluskal et al., 2010).

All samples were also screened for DCF, ERY, and GEM biodegradation products using a personal compound database and library (PCDL) developed from the literature ( SI Table S1 ). Degradation products are reported when mass errors were less than 5 ppm and the signal-to-noise ratio of peaks was greater than or equal to 10. The goal of this approach is not to quantify degradation products’ concentrations but simply to use the presence of degradation products as an indicator of biodegradation. The relative concentrations of aqueous degradation products were tracked over time by calculating corrected peak areas as follows:

where the “relevant standard” term refers to the peak area of the related parent compound at 100 µg/L. For example, the peak area of an aqueous ERY degradation product would be divided by the peak area of the 100 µg/L ERY standard measured during the same mass spectrometry run and multiplied by the sample’s concentration factor (e.g., a 100 mL sample extracted and concentrated to 1 mL had a concentration factor of 100).

Pharmaceuticals were also extracted from granules to quantify solid phase concentrations per methods adapted from Martin et al., 2010. Extraction methods are summarized in the Supplementary Information. All solid phase samples were screened for degradation products as described above. Solid phase samples were extracted in triplicate once per month and periodically evaluated for recovery using spike-recovery testing, as described above. Average recoveries are summarized in Table 1 . Peak areas of degradation products in the solid phase were corrected as follows:

Aqueous and solid samples were periodically analyzed from the control reactor for pharmaceutical parents and degradation products. Except when noted otherwise, pharmaceuticals were not detected in any phase in the control reactor. Likewise, degradation products were generally not found in the influent to the test reactor.

Lastly, the toxicities of detected degradation products were estimated using the US Environmental Protection Agency’s Toxicity Estimation Software Tool (TEST). TEST estimates chemical toxicity using quantitative structure activity relationships (Martin et al., 2020).

Immediately after pharmaceutical dosing to the test reactor began, nitrogen removal dropped to approximately 70% ( Figure 1 ). There was a brief recovery in nitrogen removal from day 6 to 20, during which removal peaked at 93%, though removal declined thereafter and stabilized at approximately 70%. Poor nitrogen removal was due mainly to incomplete ammonia oxidation ( SI Figure S2 ). Nitrite and nitrate concentrations in both control and test reactors were approximately equal for the entire experimental duration ( Figure S3 ), indicating that nitrite oxidation and denitrification were minimally impacted by pharmaceutical exposure. Nitrogen removal data for the control reactor are not plotted past day 40 because the software controlling both reactors experienced an error on day 40, stopping dissolved oxygen control in the test reactor and resulting in an acid overdose to the control reactor. DO control to the test reactor was interrupted for just 24 hours; however, the acid overdose to the control reactor caused shutdown of that reactor for 10 days. It is reasonable to expect that complete nitrogen removal would have continued in the control reactor if not for this disruption.

Phosphate removal also declined sharply in the test reactor and remained noisy for the entire experimental duration ( Figure 1 ). Despite this, DOC consumption was over 90% in both reactors ( Figure S4 ), and most carbon continued to be consumed anaerobically in the test reactor ( Figure S5 ). Near-complete anaerobic carbon consumption suggests the potential activity of glycogen accumulating organisms, discussed more in section 4.1.

All pharmaceuticals were partially removed in the first 10-12 days of dosing, evidenced by lower effluent than influent concentrations ( Figure 2 , Figure S6 ). Removal was calculated by performing a mass balance on influent and effluent concentrations at each time point (SI Equations S1 and S2). Degradation products were also detected during the first 12 days ( Figure 3 ). However, from day 12-23, effluent concentrations of all pharmaceuticals spiked to approximately twice influent concentrations, indicating a release of retained pharmaceuticals from the reactor. Solid phase data collected over this time frame also appear to indicate desorption of pharmaceuticals from AGS ( Figure S6 ). Pharmaceutical degradation is discussed in more detail in the following sections.

Metagenomic sequencing generated approximately 25 million high quality reads, resulting in a total of 19 bins with greater than 90% completeness and less than 5% contamination, and an additional 38 with completeness greater than 50% with less than 10% contamination. All MAGs were classified as members of the bacteria, including known members of glycogen accumulating organisms (GAOs) such as species in the Competibacter and Contendibacter. Phosphorus accumulating organisms (PAOs) were also present, such as species within the Accumulibacter. A MAG identified as a species in the Accumulibacter was the most abundant organism detected, at 19% of the population based on the number of raw reads that mapped to the MAG.

Despite the presence of multiple MAGs identified as known ammonia oxidizing bacteria (AOB), including Nitrosospira and Nitrosomonas, this pathway was not uncovered in the KEGG analysis of constructed MAGs, though the complete ammonia oxidation pathway was found in the total shotgun community ( Figure 4 ). The scaffolds with identified ammonia monooxygenase (amo) genes (amoA, amoB, and amoC) were relatively short (< 6000 bp). One scaffold contained all three amo genes, and another scaffold contained just amoA and amoC. Based on BLAST (Altschul et al., 1990), both were highly similar to sequences identified from Nitrosomonas species; one scaffold was similar to a MAG identified as Nitrosomonas from a metagenomic analysis of biofilms in wastewater treatment plants (Suarez et al., 2022). The absence of this pathway in putative AOB in the assembled MAGs could be due to the relatively low genome completeness calculated for these taxa, as the Nitrosomonas MAG (Granule57, Figure 4 ) was 87% complete; however, the relative abundance of this organism was 0.66%, and was not included in the overall KEGG analysis. In addition, this MAG had the pathway for nitrite reduction and nitric oxide reduction associated with NOB. We also identified the complete pathway for hydroxylamine oxidation, the second step in nitrification, in the total community, but not in the assembled MAGs.

Nitrate reductase genes were widespread in the granule community, as seen in previous metagenomic studies of wastewater (Singleton et al., 2021). Twelve MAGs contained genes for the complete dissimilatory nitrate reduction pathway and of those, ten also contained the full nitrite oxidation pathway ( Figure 4 ) including members of the Rhodocyclaceae, a Competibacter, and a MAG classified as a member of the Thermoanaerobaculia. Only a single MAG (Granule56) had the pathway for dissimilatory nitrate reduction to ammonia (DNRA), identified as a species in the Thermoanaerobaculia in the phylum Acidobacteriota, which also contained the dissimilatory nitrate reduction pathway. An organism related to the Thermoanaerobaculia was also detected in the active community, but at very low levels (< 0.1% relative abundance).

Shifts in the active microbial community were evaluated using 16S sequencing of rRNA genes and transcripts throughout pharmaceutical exposure. As described in Section 2.4.4, the active microbial community was defined as ZOTUs with a ratio of 16S rRNA transcripts (i.e. cDNA) to rRNA genes greater than or equal to one. The 16S rRNA gene read numbers from ZOTUs identified as “active” were then used to calculate relative abundances within the total active community. Multidimensional ordination ( SI Figure S7 ) showed separation of the three pools (DNA, RNA, and active), but also showed separation of the treatment and control SBRs.

At day zero, the distributions of active families in both the control and test reactors were similar: Competibacteraceae, Rhodocyclaceae, and Chitinophagaceae were most abundant ( Figure 5 ). Notably, the Rhodocyclaceae family includes PAOs (Candidatus Accumulibacter) and a potential denitrifying genus (Zoogloea). The role of bacteria in the Chitinophagaceae family is unclear—multiple OTUs matching Niabella, Terrimonas, and Flavipsychrobacter genera were detected, but to date, the function of these genera in AGS is unknown.

After five days, active families in the test and control SBRs diverged greatly. In the test reactor, active Rhodocyclaceae spiked in abundance at day five, from 14 to 46% ( Figure 5 ), which may be linked with the brief recovery in phosphate removal between days four and 10 ( Figure 1 ). Active Rhodocyclaceae then declined sharply in abundance in the test reactor for the remainder of the dosing period, stabilizing at approximately 1.5 ± 0.7% abundance for the second half of the test. In contrast, active Rhodocyclaceae in the control reactor were present at 10.2 ± 5% for the entire test duration.

Despite large shifts in the test reactor community, significant correlations between microbial community composition and effluent pharmaceutical concentrations were generally not found. A significant, though small, correlation between effluent GEM concentrations and the microbial community was observed (R^{2 =} 0.28, p = 0.048).

Despite the lack of overall correlations, when times with the largest differences in pharmaceutical removal (specifically days 5 through 17, Figure 2 ) were focused on, the strongest declines in relative abundance were concentrated within the Gammaproteobacteria and primarily in Rhodocyclaceae ( Figure 5 , Figure S8 ). The most striking shift for any single ZOTU was observed for a sequence classified as a member of Azonexus (family Rhodocyclaceae), which made up 40% of the active community at day 5 but was absent by day 17. Large changes in relative abundance were also observed for two members of Accumulibacter, which were present at over 1% of the active community on day 5 but were not detected by day 17.

In contrast, positive responses were distributed across the phylogeny, with the largest shift occurring in a ZOTU identified as a Methylotenera species: relative abundance increased from under 0.1% at day 5 to 22% at day 17. Likewise, a ZOTU identified as a Bacteroidia species (UBA2475_sp013816615) increased from under 0.1% to 6% over the same time period. Positive shifts were also observed for members of the Chloroflexota, Gammaproteobacteria, and Acidobacteriota ( Figure S8 ).

We had limited success linking the ZOTUs to the organisms identified in the metagenome assemblies. Of the top ten most abundant ZOTUs in the active community, three had no close relatives in the metagenome. The most abundant active ZOTU in the test reactor (genus Methylotenera, order Burkholderiales), was not detected in the assembled MAGs; however, this is likely due to the length of time elapsed between metagenome sampling and the onset of experimentation as well as the presence of pharmaceuticals. Four ZOTUs were classified as Competibacter A denitrificans with varying levels of confidence. The remaining ZOTUs with a close relative in the metagenome included a Xanthomonadaceae (genus SCMT01), a species of Accumulibacter, and an Azonexus species.

Methylophilaceae became the dominant active family in the test reactor after day 17 and included Methylotenera and Methylophilus genera. Members of this family may have proliferated due to the presence of trace methanol in the influent from pharmaceutical dosing (0.13 mg/L influent concentration); however, some Methylophilaceae species can also consume acetate (Jenkins et al., 1987; Dedysh et al., 2005), and therefore it is more likely that Methylophilaceae proliferated due to increased acetate availability in the aerobic phase of SBR operation ( Figure S5 ). Some Methylotenera species are also capable of aerobic denitrification (Mustakhimov et al., 2013), which may explain why nitrate did not accumulate in the test reactor. It is worth noting that this family did not appear to be affected by pharmaceuticals’ presence, but also did not proliferate until after pharmaceutical removal stopped—Methylophilaceae were therefore likely not responsible for the pharmaceutical degradation observed in the first 12 days of the test.

Competibacteraceae was the second most abundant active family in the test reactor, at 21.8 ± 10.5% abundance over the last 40 days of exposure. Bacteria in the Competibacteraceae family are glycogen accumulating organisms, or GAOs, that compete with PAOs for anaerobic carbon consumption but do not aerobically consume phosphate. For this reason, their proliferation is typically linked with reduced phosphate removal (McIlroy et al., 2014). It is likely that Competibacteraceae activity sustained anaerobic carbon consumption in the test reactor; however, it is surprising that the relative abundance of this family was lower in the test reactor than in the control reactor: active Competibacteraceae were present at 30.8 ± 18.6% abundance in the control reactor, despite complete phosphate removal ( Figure 1 ). It is possible that Competibacteraceae activity levels were similar in both reactors; the increased abundance of active Methylophilaceae in the test reactor may make the relative abundance of active Competibacteraceae appear smaller.

The relative abundance of active nitrifying families Nitrospiraceae and Nitrosomonadacea was higher in the test reactor than the control for the first 38 days of pharmaceutical exposure ( Figure 6 ), which may explain why nitrogen removal briefly recovered between days zero and 20 ( Figure 1 ). Nitrifier activity may also be responsible for pharmaceutical degradation in the first 12 days of the test: the ammonia monooxygenase enzyme used by both ammonia oxidizing bacteria and archaea (AOA) is known to react non-specifically with aromatic compounds, typically through an oxidation reaction (Keener and Arp, 1994; Schatteman et al., 2022). The primary degradation product of DCF, and all detected GEM degradation products, are formed through oxidation reactions. Aqueous GEM degradation products were not detected after day 34, approximately the same time that active nitrifier abundance dropped. Likewise, aqueous DCF1 levels fell and remained low after day 34. It is possible that decreased levels of DCF and GEM degradation products after day 34 are linked with decreased abundance of active nitrifiers over the same time period. It should be noted, however, that active Nitrosomonadaceae (the only AOB family detected) were not detected for most timepoints in both the control and test SBR. The lack of active AOB data, however, is contradicted by complete ammonia oxidation in the control reactor and partial ammonia oxidation in the test reactor, as well as active Nitrospiraceae abundance patterns over the dosing period. Several active AOA genera were also detected, which may further explain why ammonia oxidation continued throughout dosing. AOA abundance is discussed in more detail in section 4.2.

The addition of pharmaceutical compounds led to strong shifts in the microbial community and concomitant reductions in nitrogen and phosphorus removal. Based on the assessment of potential functional pathways, no high-quality MAGs were found with the ammonia oxidation pathway. Although AOB, such as Nitrosomonas, were detected in both the shotgun metagenomic and metabarcoding datasets, we were unable to construct a MAG with over 90% completeness, which limited our ability to determine the role of these organisms. However, despite their limited abundance overall, it is likely that Nitrosomonas is the main driver of ammonia oxidation in the reactors, given their presence in seed granules and the homology of the identified functional genes.

Several archaeal taxa were also found in the metabarcoding dataset though archaea were not detected in the granule metagenome. Archaea detected were primarily species within the Thermoproteota, including members of the order Nitrosophaerales. Nitrosophaerales contains numerous AOA, such as Nitrosotenuis and Nitrosocaldus genera. AOA have been posited to be important for the oxidation of ammonia to nitrite in commercial wastewater treatment plants (Limpiyakorn et al., 2013), though a study by Bai et al., 2012 found that AOA dominated in smaller scale treatment plants while AOB dominated in commercial ones. The authors hypothesized that AOA may be more sensitive to the toxic compounds more frequently measured in larger scale treatment plants.

Although archaea made up less than 1% of the sequences identified, 47 ZOTUs were identified as active taxa. We focused primarily on abundant taxa, but the ability of rare organisms to drive ecosystem function is well documented (Jousset et al., 2017); for example, rare taxa may be particularly important in pollutant degradation. The large shifts in relative abundance of wastewater-treating taxa, such as Rhodocyclaceae and Methylophilaceae ( Figure 5 ), suggest that additional studies in which microbe functionality is explored could provide important insights into community interactions that might otherwise be overlooked.