We sampled a pilot-scale advanced purification facility in El Paso, Texas that operated from June 8, 2015 to January 29, 2016 and treated 0.14 million gallons per day. The feed to the plant was secondary-treated wastewater, which was ozonated to a target concentration of 3.5 mg/L (∼5 min storage time, after which there was no detectable residual) and chloraminated to a target residual of 2–4 mg/L as Cl₂ to reduce fouling of the MF membranes. The flow was split and treated by parallel MF units (Pall module type UNA-620A and Evoqua Memcor CPII L40N), and effluents were recombined in a storage tank. Water leaving the storage tank was then split equally and fed into parallel NF (Dow, model NF90-400/34i) and RO (Hydranautics, model ESPA2-LD) processes. The recombined flow was treated by an ultraviolet/advanced oxidation process (UV/AOP). UV was supplied via Trojan UV Swift SC B08^TM at a dose of 840 mJ/cm², and oxidation was achieved by addition of hydrogen peroxide to a target concentration of 4 mg/L. UV/AOP treated water was then split for parallel filtration by three types of GAC (see Supplementary Methods).

On alternating days, effluents from each of the GAC filters were collected in sterile bottles, transported on ice to the laboratory, and chlorinated within 6 h in a 5 L glass reservoir using sodium hypochlorite (ACS grade, Spectrum Chemical MFG Corp.) to achieve a CT of 30 mg-min/L and a free chlorine residual of approximately 0.8–1.0 mg/L as Cl₂ at the end of the 30-min batch chlorination (unless otherwise noted). Chlorinated water was then transferred to three separate storage containers and continuously pumped to three annular reactors to simulate distribution (referred to throughout as SDS). Prior to operation, the SDS were sterilized and prepared as previous described (Wang et al., 2013). SDS were operated in the dark at a hydraulic residence time of 18 h (flow rate ∼0.92 mL ⋅ min⁻¹) and ambient temperature of ∼22°C. The inner cylinder rotation speed was set to 50 rpm, resulting in a shear stress of ∼0.15 N ⋅ m⁻² on the inner cylinder surface. All components of the simulated system (including the chlorination and storage reservoirs) were enclosed in aluminum foil to block light. Reactor operation is further described in Supplementary Methods. Several days before the end of SDS operation, SDS reactors were intentionally left stagnant for 24 h, after which feeding and mixing were resumed to allow stagnated reactor effluent to be sampled. Stagnation was induced only this once in the entire SDS operation.

Bulk water was sampled weekly from November 2015 through January 2016, during months 6–8 of pilot facility operation. Samples were collected in autoclave-sterilized 1-L Nalgene bottles from secondary wastewater feed to the pilot facility and from taps located after each treatment step: after chloramination, post-MF storage tank, NF permeate, combined storage after parallel NF/RO, after each of three parallel GAC column filters, and after each of three parallel SDS (Figure 1 and Table 1). Samples with residual chloramine disinfectant were quenched with excess sodium thiosulfate. Sampling of just the RO permeate did not yield sufficient DNA for sequencing, and therefore sampling was changed to the NF/RO combined tap in an attempt to concentrate more biomass onto a single filter. Sample taps were flushed for following periods: secondary wastewater and chloramine (>5 min), MF/UF filtrate (>15 min), NF and NF/RO combined permeate (>30 min), and GAC filtrate (>15 min). Water samples were filtered onto 0.22 μm mixed cellulose ester filters (Millipore) (see Supplementary Table S1 for filtration volumes). All filtration equipment (glass vacuum flask, glass funnel, glass filter supports) was autoclave-sterilized prior to use. Filters were handled using flame-sterilized tweezers. Filtered samples were stored in sterile 5-mL screw cap tubes (Axygen^TM, Thermo Fisher Scientific).

Total and intact cell counts in bulk water grab samples were determined by flow cytometry as described in Supplementary Methods. Biofilm samples were collected during plant shutdown by swabbing the internal surfaces of pipes downstream from each treatment process using sterile cotton swabs, but samples did not yield enough DNA for further analysis (data not shown). GAC media samples were collected by disassembling GAC columns immediately following the end of plant operation on January 28, 2016 (day 234 of the study). Media was poured into 50 mL tubes and stored at −80°C until DNA extraction. Media from columns 1 and 2 was from the top of the column, while column 3 was sampled from the middle. Sampling information is summarized in Table 1 and presented in detail in Supplementary Table S1.

An optimized nucleic acid extraction protocol was used on the concentrated bulk water samples because very low biomass was expected in some pilot facility and annular reactor bulk water samples. Briefly, filter membranes were torn and placed into Lysing Matrix A tubes (MP Biomedicals) with flame-sterilized tweezers. A modified phenol-chloroform-isoamyl alcohol extraction was used as previously described, with the exception of the 70°C water bath pre-heat of the filter (Pinto et al., 2012; Bautista-de Los Santos et al., 2016). Field blanks were created by filtering 1 L of autoclaved de-ionized water onto 0.22 μm mixed cellulose ester filters (Millipore) and then processed identically to field samples.

Genomic DNA (gDNA) was extracted from the GAC media samples using the MoBio DNeasy PowerSoil Kit (Qiagen) according to manufacturer’s instructions, with the following modifications: samples were incubated in lysis buffer at 65°C for 15 min prior to proceeding with the extraction protocol, bead-beating was performed at maximum speed for 10 min using a vortex adapter, the two inhibitor removal steps were combined into a single step to limit loss of gDNA, and the spin filter was incubated at 55°C for 10 min with nuclease-free water prior to elution to increase gDNA eluted. Extraction blanks were created by processing PCR-grade water alongside with and identically to field samples. Media samples yielded DNA for both amplicon and metagenomic sequencing. The amplicon libraries suffered from contamination and were excluded from analysis.

Extraction yields were assessed via Qubit dsDNA HS assay kit (Thermo Fisher Scientific). When necessary, gDNA was concentrated using a SpeedVac, and inhibitors were reduced with OneStep PCR Inhibitor Removal Kit (Zymo Research, Irvine, CA, United States).

Library preparation for amplicon sequencing followed the Schloss Lab MiSeq wet-lab protocol for amplification of the V4 region of the 16S rRNA gene¹ (Kozich et al., 2013) with slight modification (Supplementary Table S2). Briefly, the V4 region was amplified using uniquely barcoded 515F and 806R primers with Phusion Hot Start II polymerase, HF buffer, dimethyl sulfoxide (3.75% final concentration) and 1–2 μl of genomic DNA. Triplicate 25 μl reactions were combined and concentrated by SpeedVac. For difficult reactions, up to five replicates were pooled. The dual-barcoded libraries (including sampling, extraction, and PCR negative control libraries) were normalized using the SequalPrep Normalization Plate kit (Invitrogen), pooled, and sequenced on a MiSeq with V3 chemistry for 600 cycles, yielding paired-end 250 bp reads. Mock-community control libraries were constructed using the ZymoBIOMICS Microbial Community DNA Standard (Zymo Research, Irvine, CA, United States) and included in amplicon and metagenomic sequencing (Supplementary Figure S1).

Library preparation for metagenomic sequencing was performed at the Functional Genomics Laboratory at UC Berkeley. Briefly, 500 bp insert libraries were constructed using the KAPA Hyper prep kit with a PCR amplification step. Sequencing was performed on an Illumina HiSeq 4000 at the Vincent J. Coates Genomics Sequencing Laboratory at UC Berkeley, yielding 150 bp paired-end reads.

Reads were demultiplexed and mapped to PhiX (21.7%). Forward reads were processed using the DADA2 pipeline v1.4.0 (Callahan et al., 2016) according to http://benjjneb.github.io/dada2/bigdata.html to produce ASVs. Briefly, based on evaluation of read quality using FastQC², reads were truncated to 230 nucleotides (nt) and left-trimmed to remove 15 nt from the 5′ end, resulting in 215 nt sequences. Reads were also truncated where quality dipped below Q = 11, and any reads <200 nt were removed. On average 70% of reads remained after trimming. Error modeling used the first 2 million filtered reads, and samples were denoised individually using this error model. The denoised sequences were screened for chimeras with removeBimeraDenovo, which removed 3153 sequences accounting for 4.3% of all reads. Taxonomy was assigned using the Ribosomal Database Project’s naive Bayesian classifier (Wang et al., 2007) within DADA2. The classifier was trained on Silva v128 (Quast et al., 2013) and set to use a minimum bootstrap confidence of 50 to assign a given rank (default minBoot = 50). Error-corrected, taxonomically classified data were analyzed with the Phyloseq (v1.20.0) package in R (McMurdie and Holmes, 2013) (see Supplementary Data Sheet S2 and https://github.com/rosekantor/16S_metagenomics_notebooks).

Data from bulk water samples and corresponding field, extraction, and amplification negative controls were prefiltered for successful sequencing (>10000 reads), and low-biomass samples were removed prior to analyses due to obvious contamination with secondary wastewater (two samples) and MF filtrate (one sample). This prefiltering yielded 37 samples and 3 negative controls. DESeq2 analysis was performed with all ASVs found in any prefiltered negative control or sample. Of 108 ASVs that were shared between any sample and any negative control, 7 were significantly enriched in samples over negative controls (DESeq2 p-adjust <0.01). These may represent cross-contamination of negative controls with samples or different but closely related organisms present in samples and negative controls. We removed from samples all ASVs found in negative controls except for these 7. After filtering contaminant reads, one sample and one negative control with fewer than 300 reads were removed, and the two remaining amplification and extraction controls contained 2 and 3 ASVs, respectively.

Rarefaction curves were plotted for bulk water samples using rarecurve in the vegan R package (v2.4.4). The DADA2 pipeline removes singletons because they cannot be reliably distinguished from sequencing errors, precluding the use of alpha diversity metrics that rely on singletons. Given that the rarefaction curves plateau for all samples (Supplementary Figure S2), total observed ASVs were used as a metric of alpha diversity. Beta diversity was analyzed via principle coordinates analysis (PCoA) using Bray-Curtis dissimilarity and variables contributing to differences between communities were identified with PERMANOVA using adonis in the vegan package for R. Additional PERMANOVA testing using Aitchison distance found the same level of significance as was observed when using Bray-Curtis dissimilarity of percent relative abundance data. Calculating and testing for homogeneity of group dispersions used the betadisper and permutest functions within the vegan package in R (groups were sample locations).

To estimate abundances of genera containing opportunistic pathogens, the relative abundances of all ASVs (>0.05% relative abundance in any sample) belonging to the genera Legionella and Mycobacterium were summed for each sample at each location in the treatment train. The minimum, median, and maximum relative abundances were then determined for each location. To generate the range of absolute abundances, we multiplied minimum, median, and maximum relative abundances by the minimum, median, and maximum total cell counts ever observed at a given location across the duration of the study.

Metagenomic read processing used FastQC (Babraham Bioinformatics, see text footnote 2) to inspect quality, bbmap (Bushnell³) to remove PhiX and adapter sequences, and sickle (Joshi and Fass, 2011) for quality trimming. Read-based phylogenetic characterization was performed with MetaPhlAn2 using default parameters (Truong et al., 2015). Assembly was performed independently for each sample using idba_ud (Peng et al., 2012) with the “–pre_correction” option. Annotation of scaffolds ( 1 kbp was performed with Prokka v1.12 (Seemann, 2014) and RNAmmer (Lagesen et al., 2007).

Metagenomic analysis was conducted in Anvi’o (Eren et al., 2015) according to the metagenomic workflow for v3 and v4. For each assembly (scaffolds ≥2.5 kbp), an Anvi’o contig database was generated and a profile database was created that included read-mapping information from each sample. Taxonomy was added using Centrifuge (v1.0.2-beta). Binning was performed manually based on hierarchical clustering by sequence composition and differential coverage in Anvi’o. Genome bins >70% complete with <10% contamination according to CheckM (Parks et al., 2015) were de-replicated to generate a non-redundant set of MAGs using dRep (Olm et al., 2017) with primary clustering at 95% ANI and secondary clustering at 99% ANI, requiring 60% coverage of the larger genome in each pairwise comparison. MAGs were taxonomically identified based on their position in a concatenated gene tree (see Supplementary Methods), Centrifuge results, and CheckM taxonomy placement (Supplementary Tables S3, S4). Non-redundant MAGs were curated to remove misassemblies using ra2.py (Brown et al., 2015)⁴ and open reading frames were predicted using Prodigal in single genome mode (Hyatt et al., 2010). Annotation of predicted proteins was conducted using USEARCH (Edgar, 2010) against the KEGG database (Kanehisa et al., 2016). C1-carbon utilization was determined by manual searches of KEGG-based annotations against select proteins representative of the pathways of interest (Supplementary Table S5).

To investigate shared microbial community members across samples, reads from each sample were stringently mapped to (1) assemblies of each sample (scaffolds ≥1 kbp), (2) scaffolds ≥1 kbp containing unique RpsC sequences, and (3) the non-redundant MAGs (see Supplementary Methods). Stringent mappings were also used to calculate indices of replication (iRep values) (Brown et al., 2016) for each MAG in each sample using iRep⁵ with default settings.

Antibiotic resistance genes were detected via HMMsearch of the predicted proteins against the ResFams core database of hidden Markov models (HMMs) (Gibson et al., 2015). Within each resistance gene family, ResFam hits were sorted by length and clustered at 99% identity using USEARCH -cluster_fast (Edgar, 2010). Tabulated results were searched for any clusters containing proteins from both secondary wastewater and post-NF/RO samples. Ammonium monooxygenase genes were identified using TIGRFAM HMMs, and custom HMMs were used to search for Adenovirus and JC Polyomavirus sequences (see Supplementary Methods).

A phylogenomic tree for the Mycobacteriaceae was constructed with the Anvi’o phylogenomics pipeline, using thirty reference genomes downloaded from NCBI representing each of the major clades described in Gupta et al. (2018). A subset of 108 of the single copy genes/domains from Campbell et al. (2013) were identified via HMMsearch (Eddy, 2009), and individual alignments were constructed with MUSCLE (Edgar, 2004). A concatenated alignment was manually curated to remove C- and N-termini overhangs and regions of poor alignment. The resulting alignment was used to construct a phylogenetic tree with RAxML (Stamatakis, 2014) using 100 bootstraps.

Absolute abundances and alpha diversity of the bacterial community followed similar trajectories throughout the treatment train. Flow-cytometric cell counts decreased from secondary wastewater through NF/RO and then increased after GAC filtration, indicating growth on filter media (Figure 2A). Cell counts in GAC filtrate and SDS were 3 to 4 log₁₀ lower than the secondary wastewater feed. Despite chlorination following GAC filtration, the median intact cell count in the three SDS was 2.7 × 10³ cells ⋅ mL⁻¹. Intact cell counts were often lower than total cell counts, indicating that some sequence information was derived from non-viable cells. Median richness decreased across the treatment train with a slight increase after GAC (Figure 2B). Observed ASVs were used as a proxy for richness within each sample because rarefaction curves plateaued for nearly all bulk water samples (Supplementary Figure S2) and singletons had been removed during data processing. Between secondary wastewater and SDS, median observed ASVs decreased from 536 to 34 (a 1.09 log₁₀ reduction) (Figure 2B).

There were significant differences in the taxonomic composition of bulk water bacterial communities from different sample locations (Figure 3A,B). By principal coordinate analysis (PCoA), bacterial communities clustered by location within the treatment train (Figure 4). This clustering was supported by PERMANOVA, which suggested differences between communities were partially explained by sample location (PERMANOVA R² = 0.57; p < 0.001; same level of significance using proportion-normalized or rarefied data), but not sample date (R² = 0.15; p < 0.08). Some of the significance of sample location in this test may be due to heterogeneous group dispersions of samples within different locations (permutest p = 0.053 with 10000 permutations) (Supplementary Figure S3). Notably, MF and GAC had the lowest distance to the centroid of their groups, suggesting the most consistent communities. NF/RO permeate and negative controls did not include enough successfully amplified and decontaminated samples to provide information about within-location variability (n = 2 for each).

After each treatment step, the taxonomy of the most prevalent ASVs changed, and overall membership also changed (Figure 3B). One exception was that secondary wastewater communities appeared similar to some samples taken after chloramine treatment. This effect might be expected as ozonation and chloramination inactivated cells (Figure 2A) but may not have significantly damaged their DNA. Relative abundances and phylogenetic classifications based on metagenomic data yielded a somewhat similar picture of the microbial community, except for among low-biomass samples. We used a read-based method (MetaPhlAn2) (Truong et al., 2015; Supplementary Figure S4), unique RpsC sequences from assembled data (Supplementary Figure S5), and non-redundant MAGs from across all samples (Figure 6A) to assess the community composition and structure within metagenomes compared with amplicon data.

In the secondary wastewater before and after chloramine treatment, the dominant ASVs (as high as 57% relative abundance), most abundant RpsC sequence, and most abundant MAGs corresponded to the family Neisseriaceae. Meanwhile, MetaPhlAn2 identified Rhodocyclaceae as the most abundant Family in secondary wastewater (Supplementary Figure S4). One post-chloramine sample and all MF filtrate samples were overwhelmingly dominated by a single ASV classified as Mycobacteriaceae. In the NF permeate sample, the most prevalent ASV was classified as Pseudonocardiaceae (61.0% relative abundance), while in the NF/RO permeate storage tank, an ASV classified as Methylobacteriaceae was highly abundant (92.1% relative abundance) (Figure 3A). Abundant ASVs, metagenomic reads, RpsC sequences, and MAGs in the GAC filtrates were classified as Rhizobiales (including Bradyrhizobiaceae, Hyphomicrobiaceae, Methylocystaceae, and Methylobacteriaceae) (Figure 3A, 6 and Supplementary Figures S4, S5). The SDS communities were more variable over time than GAC communities, likely due to variation of the chlorine residual and a 24-h stagnation event just prior to the final sampling (day 227). In the two sample days prior to stagnation, SDS 1 and SDS 3 communities contained abundant ASVs, RpsCs, and MAGs classified as Burkholderiales/Comamonadaceae, Methylobacteriaceae, and Bradyrhizobiaceae. By amplicon sequencing, the community in SDS 2 prior to stagnation (SDS 2, day 205) appeared similar to the GAC filtrate (Figure 3), but by all metagenomic analyses, the same sample appeared more similar to the other SDS samples prior to stagnation (Supplementary Figures S4, S5 and Figure 6A). Based on our amplicon data with limited sampling, all three SDS communities were different before and after stagnation (Figure 3A,B).

Seven ASVs found in negative controls and in samples from all locations were not removed during data decontamination because they were significantly enriched (adjusted p < 0.005) in samples relative to controls. Two of these ASVs, classified as Bradyrhizobium and Hyphomicrobiaceae, were abundant in low biomass samples, especially GAC filtrates (33 and 22% median relative abundances) and some SDS samples. In metagenomic data, these two taxa were also highly abundant in GAC media and filtrate but not SDS (Supplementary Figure S5 and Figure 6A).

In addition to bacterial genomes, several eukaryotic nuclear and plastid sequences were detected with metagenomics. Most notably, all three GAC media metagenomes contained identical sequences from a single Chlorophyta bin. Based on best BLAST hits of chloroplast 16S, 23S, and co-binned nuclear 28S rRNA genes to NCBI-nr, the genome was likely related to the Sphaeropleales. The chloroplast RpsC sequence from this bin was detected by read-mapping in all three GAC media metagenomes (Supplementary Figure S5), and the presence of algae is consistent with observations of green growth on walls of the media columns. An additional Chlorophyta species with 100% identity to Parachlorella kessleri across both the chloroplast 16S rRNA gene and RpsC sequence (NCBI accessions FJ968741.1 and YP_003058311.1) was present in the metagenome of only GAC media from column 3 (Supplementary Figure S5), and this column was designed to filter water 69% more slowly than the others (see Supplementary Methods).

To examine potential persistence of bacteria from secondary wastewater and early treatment steps into post-NF/RO samples, we investigated the reoccurrence of specific ASVs, and used metagenomic read-mapping to assemblies, to unique RpsC-containing scaffolds, and to non-redundant MAGs. ASVs were defined as recurring if they were detected with >50 reads in at least two samples from a given location. The 50-read cut-off was used as minimal protection against barcode bleed (also known as sample cross-talk) from Illumina sequencing (Mitra et al., 2015; Edgar, 2016). The criterion for recurrence in two samples was used to protect against one-off cross-contamination of samples or barcodes. There was no intersection between the recurring ASVs in the secondary wastewater or MF samples and ASVs found in the SDS (all three SDS treated as one location) (Figure 5). A total of four ASVs were shared between pre-NF/RO samples (secondary wastewater and MF) and post-NF/RO samples (GAC, all three treated as one location) (Figure 5 and Supplementary Figure S7). One of these four ASVs was the Neisseriaceae sequence highly abundant in secondary wastewater. This ASV was detected but not prevalent in samples after NF/RO, having a maximum of 160 reads or 1.3% relative abundance in any post-NF/RO sample (GAC column 3, day 184). Two other recurring ASVs, classified as Methylocystaceae and Bradyrhizobiaceae were most prevalent in GAC and SDS samples and were widespread in these samples. These ASVs were also detected before NF/RO with as many as 309 reads (1.5%) and 2168 reads (8.1%), respectively, in pre-NF/RO samples. The Bradyrhizobiaceae and Neisseriaceae ASVs were also present in negative controls, but were not removed by data decontamination (see section “Materials and Methods”). Lastly, the fourth ASV shared between pre- and post-NF/RO samples was classified as Comamonadaceae, and was relatively rare in all samples. This ASV constituted less than 0.44% relative abundance in any pre-NF/RO and less than 2.4% in any post-NF/RO sample (Supplementary Figure S7).

Metagenomic analyses with stringent read-mapping also revealed little to no overlap between secondary wastewater and post-NF/RO samples (GAC media, GAC filtrate, and SDS bulk water) (Supplementary Figure S8). At most, 0.011% of metagenomic reads from any given secondary wastewater sample were stringently mapped to the assemblies of any post-NF/RO samples (=1 mismatch per read across scaffolds ≥1 kbp). Likewise, fewer than 0.10% of reads from any post-NF/RO sample mapped to a secondary wastewater assembly (Supplementary Figure S8). Based on read-mapping to scaffolds containing the marker gene RpsC, no overlap was detected between the secondary wastewater and the post-NF/RO (Supplementary Figure S5). Mapping reads to 38 non-redundant MAGs gave the same result (Figure 6A).

As a group, GAC media samples shared the majority of unique RpsC sequences and MAGs with GAC filtrates. Individually, between 42 and 59% of GAC filtrate reads could be mapped to their respective GAC media assemblies (Supplementary Figure S8). Between GAC filtrates and SDS samples, overall community structure differed; while some MAGs were shared, others were detected in only GAC filtrate or SDS (Figure 5). Between these sample points, water underwent batch chlorination, storage, and simulated distribution, and there was a 10-fold drop in median intact cell counts (Figure 2A). Thus, some of the ASVs shared between GAC filtrate and SDS may represent chlorine-inactivated cells whose DNA remained intact (Figure 3B, 6A and Supplementary Figure S5).

To investigate whether the advanced treatment facility could be a source of opportunistic drinking water pathogens, we calculated the minimum, median, and maximum possible relative and absolute abundances of Legionella spp. and Mycobacterium spp. at each location in the treatment train (Table 2). Calculations were based on amplicon sequencing (including ASVs >0.05% relative abundance in at least one sample) combined with total cell counts quantified via flow cytometry. The estimated maximum Legionella spp. in secondary wastewater and SDS was 4.1 × 10⁵ cells ⋅ mL⁻¹ and 1.6 × 10³ cells ⋅ mL⁻¹, respectively (Table 2 and Supplementary Figure S9). Of the seven Legionella ASVs included in the calculations, four were from pre-NF/RO, two from post-NF/RO, and one shared between pre- and post-NF/RO samples (Supplementary Figure S9). In the genus Mycobacterium, maximum possible concentrations in secondary wastewater and SDS were 1.7 × 10⁶ cells ⋅ mL⁻¹ and 3.7 × 10⁴ cells ⋅ mL⁻¹, respectively (Table 2). Only one Mycobacterium sp. ASV was detected in both pre- and post-NF/RO samples. This was the same ASV that was highly abundant in MF filtrate samples and detected at very low abundances in every other location. The Legionella and Mycobacterium ASVs found in pre-and post-NF/RO were not identified in the analysis of persistent ASVs across the treatment train (see section “Persistence of Bacteria Through Advanced Treatment”) because they didn’t meet the criteria for consideration. Thus, while they were detected, their presence in the data may be due to Illumina barcode bleed or sample cross-contamination. The reported maximum opportunistic pathogen concentrations were likely overestimated because they were calculated using maximum relative abundances and maximum cell counts ever observed at a given location across the duration of the study, regardless of whether these occurred at the same time. Furthermore, abundance calculations used total cell counts and thus also included non-viable cells (Figure 2A) and assumed a single copy of the 16S rRNA gene per cell. Legionella spp. and Mycobacterium spp. have median 16S rRNA gene copy numbers of 3 and 1, respectively⁶.

Mycobacteriaceae were also identified in metagenomic analyses, where they could be more accurately classified. A recent phylogenomic study split the genus Mycobacterium into five new genera, of which one, the Mycobacterium spp., contains the majority of slow-growing pathogens (Gupta et al., 2018) (not reflected in SILVA v128 used for amplicon analysis). We reconstructed two unique, near-complete Mycobacteriaceae MAGs, each detected by read-mapping in at least three samples. Using Gupta et al. (2018) as a starting point, we conducted a phylogenomic analysis, placing the two MAGs within the Mycobacteriaceae. One MAG was nearest neighbor to M. gordonae (within the slow-growing, pathogen-containing Tuberculosis-Simiae clade, genus Mycobacterium), and the other closest to M. smegmatis (within the rapid-growing Fortuitum-Vaccae clade, genus Mycolicibacterium) (Supplementary Figure S10). Notably, the M. gordonae-related MAG was not detected by read-mapping in two of the SDS, and in the third SDS it was at low relative abundance (coverage of 0.24× per 10 million reads) (Figure 6A).

In addition to opportunistic pathogens, we searched metagenomic assemblies for protein sequences from human adenoviruses and polyomaviruses (both dsDNA enteric viruses) that could indicate passage through the treatment train. Using custom built HMMs (see Supplementary Methods), no hits were identified for any of these proteins within any metagenome, including secondary wastewater samples. This result was not surprising, given that pore size of membranes used for sample collection (0.22 μm) was much larger than these viruses and the sample concentration method was not optimized for retention of viruses.

We examined whether ARGs from secondary wastewater had passed through advanced treatment or were sourced from bacteria within the treatment train. Putative ARGs were identified in assembled metagenomic data by searching against the ResFams database (Gibson et al., 2015). Of the 123 ARG families in ResFams, only 34 were found, and these were mechanistically classified as efflux pumps, enzymes catalyzing structural modifications to antibiotics, beta-lactamases, and regulatory proteins (Figure 6C and Supplementary Figure S11). Of the 1609 protein hits across 12 metagenomes, no putative ARG sequences from secondary wastewater samples clustered with sequences in post-NF/RO samples at ≥99% amino acid identity. Putative ARGs were found at higher normalized coverages in GAC media, GAC filtrates, and SDS samples relative to wastewater, an effect primarily driven by ARG families associated with transport, and beta-lactamases of classes A and B (Supplementary Figure S11). Many ARGs co-occurred in the same MAGs (Figure 6C).

Given the increases in microbial cell counts after GAC filtration (Figure 2A), we were interested in microbial growth in oligotrophic conditions after NF/RO. We estimated TOC that could have passed through a single GAC column over 205 days prior to sampling. With a maximum flow rate of 4,400 L per day through a GAC column, and the typical 50–300 μg ⋅ L⁻¹ of TOC after RO in advanced treatment trains (Zeng et al., 2016), we estimated the maximum possible amount of TOC adsorbed in one column was 45–135 g. Chemical analyses have shown that aldehydes and other small carbon molecules are formed after ozonation of wastewater (Wert et al., 2007), and other studies have shown that small uncharged compounds may pass through RO membranes (Bellona et al., 2004). Other carbon sources for bacterial growth may have included carbon fixed during algal growth in the columns, and the GAC media itself.

Given these possible carbon sources, annotations of predicted proteins from the non-redundant MAGs were used to assess the capacity for C1-carbon compound utilization, including methanol and formaldehyde (Figure 6D). Seven MAGs recovered from post-NF/RO samples encoded putative methanol or methanol/ethanol dehydrogenases. Additional MAGs encoded the capacity to convert formaldehyde to formate, which could then be a substrate for energy generation via formate dehydrogenase. Lastly, incorporation of C1-carbon compounds into biomass may occur via the serine cycle, which was encoded in some MAGs (Figure 6D). It should be noted that the absence of specific pathways in MAGs does not necessarily indicate the inability of the corresponding organisms to grow on C1 compounds. These analyses were not exhaustive of all possible pathways, nor were all genomes complete.

Lastly, a proxy for bacterial growth rate, the index of replication (iRep), was calculated for each MAG in each sample. After filtering to remove MAG-by-sample combinations with low coverage, high genome fragmentation, and poorly fitted regressions (see Brown et al., 2016), iRep values were retained for 80 MAG-by-sample combinations (Figure 6B and Supplementary Table S5). Many genomes had low iReps (1.1–1.2) suggesting they were not actively replicating. The fastest replication rates (>1.5) were observed for four genomes in the secondary wastewater samples (Rhodocyclaceae_3, Flavobacteriaceae_1, Proteobacterium_1, and Acidimicrobidae_1), two genomes in GAC media (Hyphomicrobium_1 and Bradyrhizobiaceae_1), one genome in GAC filtrate (Bradyrhizobiaceae_1), and one genome in an SDS (Curvibacter_1).