Fifty-one polar TOrCs which have been typically found in municipal wastewater were selected as target compounds in this study (Table 1), including pesticides, pharmaceuticals, and personal care products with different biodegradabilities in CAS treatment. For example, acyclovir, atenolol, caffeine, and ibuprofen can be relatively easily degraded (Prasse et al., 2011; Ferrando-Climent et al., 2012; Xu et al., 2017a; Chtourou et al., 2018), while carbamazepine and diclofenac are persistent and difficult to be biotransformed (Zhang et al., 2008; Kruglova et al., 2014), and some compounds such as emtricitabine and cetirizin have only been rarely studied so far. All TOrCs reference standards used were of analytical grade (Sigma-Aldrich) with a minimal purity of 98%.

Six different sand materials from rapid sand filters of municipal WWTPs and two materials from slow sand filters used for water treatment were sampled and then stored at 4°C for ~3 months. The batch experiments were performed in triplicates with all eight sand filters under the same conditions: 20 g of sand material was added to each 250 mL bottle and 80 mL of treated wastewater from the WWTP Koblenz. All batches were aerated with 7 mg·L⁻¹ oxygen and continuously shaken at 22°C and the speed of 100 rpm in the dark. To acclimate the microorganisms, all batches were equilibrated to the experimental conditions for ~24 h before the experiment started. The experiment was started by spiking the bottles with a mixture of 51 TOrCs at the initial concentration of 0.5 μg·L⁻¹ for each compound, which was within the range of TOrCs concentration in actual wastewater (Luo et al., 2014).

Over 72 h incubation, water samples were withdrawn from each bottle at regular intervals resulting in eight time points (0, 1, 4, 8, 12, 24, 48, 72 h) for k_biol calculation (see below). Samples were immediately filtered using 0.45 μm regenerated cellulose membrane (Macherey-Nagel, Germany) and stored at 4°C for maximum 2 weeks until liquid chromatography tandem mass spectrometry (LC-MS/MS) analysis. At the end of the experiment, the remaining water was removed and the sand materials were collected and stored at −80°C for DNA extraction. In addition, each sand material was fixed 1:1 with ethanol for fluorescence in situ hybridization (FISH) analysis.

In contrast to the DNA based analysis and the monitoring of the compound removal (see below), the biomass measurements and the cell counting was performed as one measurement per filter material (n = 8). Loss on ignition (LOI, a measure for total organic matter content) was determined by drying 15 g of sand filter material at 105°C overnight, and then placing in a muffle furnace at 440°C for 4 h, weighing the samples at each step. ATP from 1 g sand filter material was measured by the single tube luminometer Sirius FB12 (Titertek Berthold, Germany) using the ATP Biomass Kit HS (BioThema, Sweden) according to the manufacturer's instructions with eight measurement points prior and after ATP standard addition. The sand samples were sent to Vermicon AG (Munich, Germany) for cell count measurement and FISH analysis. FISH counting was performed with the sand samples using fluorescence microscopy for following groups: Alphaproteobacteria, Betaproteobacteria, Gammaproteobacteria, Deltaproteobacteria, Epsilonproteobacteria, Actinobacteria, Firmicutes, Cytophaga-Flexibacer-Subphylum, Planctomycetes, Chloroflexi, Nitrospirae, TM7, Archaea. Probe sequences and hybridization conditions were taken from Daims et al. (2006). FISH protocol was performed following Snaidr et al. (1997) with following modifications for sample preparation: 10 g of sand samples were thoroughly vortexed for 5 min. After settlement of the samples at room temperature 200 μL of supernatant were transferred into a 2 mL vial and mixed with 1 mL of PBS buffer. Samples were vortexed, followed by centrifugation at 5,000 rpm for 1 min. Supernatant was extracted and transferred into new 2 ml vial, followed by centrifugation at 5,000 rpm for 5 min. Supernatant was discarded and remaining pellet was solved in 200 μL of 50% EtOH. In situ hybridizations were performed on the EtOH-fixed samples using the gene probe technology VIT® (vermicon identification technology) by Vermicon AG (Munich, Germany). At least 10 fields of view were counted per measurement using a Zeiss Axioscop 2 epifluorescence microscope equipped with fluorescence filter sets for the dyes FAM, Cy3, and 4,6-Diamidino-2-phenylindol-dihydrochloride (DAPI). More information can be obtained by request from Vermicon AG (Munich, Germany).

To determine the concentrations of the 51 spiked TOrCs, all filtered samples of the eight time points were used for LC-MS/MS analysis via an Agilent 1260-LC coupled to a SCIEX QTrap 5500-MS according to the method previously described by Falås et al. (2016). Along the sampled 72 h, biotransformation rate constants were estimated assuming first-order kinetics and negligible sorption according to Equation (1), whereby S is the TOrCs concentration (μg·L⁻¹), X_Biomass is the biomass concentration represented by the ATP (pmol·g⁻¹), and k_biol is the normalized pseudo-first-order rate constant in g·(pmol·d)⁻¹.

The DNA of the microorganisms was extracted from the sand filter material, and was used to analyze the present taxa by 16S rRNA amplicon sequencing, and by whole genome shotgun sequencing. Approximately 1.5 g of the frozen filter material from the end of the experiment (24 samples) was used for DNA extraction using the FastDNA SPIN Kit for soil (MP Biomedicals, Eschwege, Germany) according to the manufacturer's instructions. Concentrations and purity of the individual DNA extracts were measured by microspectrophotometry (NanoPhotometer P330, Implen). 16S rRNA gene amplification, library preparation and sequencing were performed by IMGM Laboratories GmbH (Planegg, Germany). In short, primer pair 341F (CCTACGGGNGGCWGCAG) and 805R (GACTACHVGGGTATCTAATCC) for V3-V4 hypervariable region was used for amplification of the 16S gene. Next, the purified PCR products were normalized to equimolar concentration with the SequalPrep Normalization Plate Kit (Thermo Fisher Scientific) and the DNA library was sequenced on the Illumina MiSeq next generation sequencing system (Illumina Inc.) in paired end mode with 2 × 250 bp. Whole metagenome sequencing was performed by IMGM Laboratories GmbH (Planegg, Germany). In brief, the extracted DNA samples were checked for quality and quantity by a 1% agarose gel (Midori Green-stained) electrophoresis and a Qubit dsDNA HS Assay Kit (Invitrogen). Prior to library preparation, the high molecular weight gDNA was sheared to target fragment size of 550 bp using a Covaris M220 Focused-ultrasonicator (Brighton, United Kingdom). In the following, DNA was prepared for sequencing using the NEBnext UltraTM II DNA library preparation kit for Illumina (New England Biolabs). Sequence data were generated on the Illumina HiSeq 2500 next generation sequencing system with the 2 × 250 bp paired end mode.

All 16S rRNA gene amplicons data were processed in R (v 3.6.0) (R Core Team, 2016) using the DADA2 (v 1.9.0) pipeline (Callahan et al., 2016). Quality trimming, denoising, error-correction, paired-end read merging, chimera removal, and dereplication steps were performed according to the default setting (Supplementary Table 1). The amplicon sequence variants (ASVs) were taxonomically classified with a naïve Bayesian classifier using the SILVA training dataset (v 137) (Quast et al., 2012). The ASV and taxonomy tables, along with associated sample metadata were imported into phyloseq (v 1.22.3) (McMurdie and Holmes, 2013) for community analysis. Whole metagenome sequencing data was processed using Trimmomatic (v 0.39) (Bolger et al., 2014), and the quality of the reads was checked with FastQC (v 0.11.8) (Andrews, 2010). The metagenome data was used for (a) validating the amplicon results with METAXA2 software (v 2.1.1) (Bengtsson-Palme et al., 2015); (b) calculating the metagenomic distances between the samples by making a gene-pool comparison/relatedness determination with the k-mer based weighted inner product, kWIP (Murray et al., 2017); (c) metagenome assembly and functional analysis. For these latter two analysis, we used MEGAHIT (v 1.2.9) (Li et al., 2015) to assemble the data, and QUAST (v 5.0.2) (Gurevich et al., 2013) to evaluate the assembly quality. For the subsequent functional analysis, SUPER-FOCUS (v 0.31) (Silva et al., 2016) was applied using the aligner DIAMOND (v 0.9.34) to obtain the enzyme commission (EC) catogories (Buchfink et al., 2015). In addition, taxonomic informed pathways classifications by sample groups were determined by HUMAnN2 (v 2.8.1) (Franzosa et al., 2018). Metagenome-assembled genomes (MAGs) were obtained by MetaBAT2 (v 2.15) (Kang et al., 2019) and manually refined after assessing the completeness (>70%) and contamination (<10%) by CheckM (v 1.0.13) (Parks et al., 2015). Taxonomic classification was conducted by GTDB-Tk (Chaumeil et al., 2020). Prodigal (v 2.6.3) (Hyatt et al., 2010) was used for open reading frames (ORFs) prediction. KofamKOALA (https://www.genome.jp/tools/kofamkoala/) (Aramaki et al., 2020) was used to obtain KO annotations for genes predicted by Prodigal. Xenobiotics metabolism was performed by the “Reconstruct Pathway” tool in KEGG mapper (https://www.genome.jp/kegg/mapper.html, accessed February 2021).

All statistical analysis was performed in R environment (v 3.6.0). cor.test function was used to test the correlation coefficient r between the TOrCs biotransformation performance and different biomass indicators based on the Pearson method. Linear regression was used to model the global compound removal by ATP concentration. Residuals were tested for normality (Shapiro-Wilk normality test), and the distribution was inspected through QQ (quantile-quantile) plot (Supplementary Figure 1). A paired two-sample Mann-Whitney-Wilcoxon test (non-parametric) was used to identify significant differences between k_biol values of two sand groups. Kruskal-Wallis test was used to evaluate the differences of individual compound k_biol among all sand samples and between categories. The relationships between DADA2 and METAXA2 data were examined by Mantel tests. The microbial diversity indices were analyzed using the vegan package (v 2.5-6) (Oksanen et al., 2019). The species richness was determined by rarefying the amplicon dataset to the smallest sample (5,799 reads) through the “rrarefy” function (Supplementary Figure 2). Chao1 and Inverse Simpson index (performed on full dataset prior to rarefaction) were used to present community richness and alpha diversity, respectively. Community compositions were compared using Bray-Curtis dissimilarities on ASV abundances and presented using NMDS ordinations. Ordinations and heatmaps were done in the R package “ampvis2” (v 2.4.6) (Andersen et al., 2018). The analysis of correlations between TOrCs biotransformation rates and microorganisms or functional genes was based on Pearson correlations outlined in the Rhea script collection (Lagkouvardos et al., 2017). Before doing the correlation analysis, we tested the differences between each compound k_biol and zero (criterion 1), and set the minimum removal percentage to 10% (criterion 2). As a consequence, the k_biol of carbamazepine, fluconazole and xipamide were excluded from further analysis, since either criterion 1 or 2 was not met in any of the batches. Differential functional genes and pathway analyses were conducted by DESeq2 (v 1.29.5) (Love et al., 2014).

Physical, chemical and biological characteristics of eight sand filters are shown in Table 2. Materials of these filters included quartz gravel, quartz sand, anthracite and pumice. Loss on ignition, ATP, DNA concentration and cell counts served as biomass indicators. Biomass and organic matter content was high in Friedrichshafen, Eriskirchen and Wangen (>25 pmol/g ATP, >40 mg/g Loss on ignition), whereas the remaining five biofilters were low in biomass and organic matter, with the lowest measurements in the filter materials from drinking water systems (BWA and IFW; <1.5 pmol/g ATP, <9 mg/g Loss on ignition) (Table 2, Figure 1A). The biotransformation rate constants of TOrCs were calculated excluding six compounds (i.e., acridone, caffeine, ibuprofen, levetiracetam acid, ramiprilat and valsartan) which did not fit the first-order kinetics. Moreover, carbamazepine, fluconazole, and xipamide had non-significant k_biol and showed <10% removal in all sand samples. Wangen achieved the best overall TOrC parent compound removal with a global mean removal percentage of 58 ± 0.64%, followed by Friedrichshafen (49 ± 3.4%) and Eriskirchen (45 ± 1.2%) (Figure 1B). Hungen and Stuttgart showed similar global removal percentages of 28 and 26%, followed by Moos with 14%, respectively. In the drinking water filters IFW and BWA, there were almost no TOrCs degraded (removal below 0.5%). These differences in global transformation potential were highly correlated with biomass expressed as ATP concentration (Pearson r = 0.92, p < 0.001), followed by living cell counts (Pearson r = 0.76, p < 0.001), while other parameters showed no correlation (Table 2). The mean global k_biol could be predicted by ATP concentration [logarithmic: ln(ATP), LM, t = 8.3, adjusted R² = 0.76, p < 0.001; or linear: LM, t = 10.9, adjusted R² = 0.84, p < 0.001, Figure 1C]. The linear model showed a better fit, but the distribution of residuals was non-normal (Supplementary Figure 1C). In addition to k_biol, the overall average removal percentage showed a clear relationship with ATP concentration, this time on a natural logarithmic scale (LM, t = 13.9, adjusted R² = 0.90, p < 0.001, Figure 1D). Therefore, we used ATP for XBiomass normalization (Equation 1) of k_biol for the clustering analysis and the following correlation analysis. According to the ATP biomass, the samples with high biomass (Eriskirchen, Friedrichshafen and Wangen, named as HBG) had a high general biotransformation ability, while the biotransformation potential of the samples with low biomass (Hungen, Moos and Stuttgart, named as LBG) was relatively low, which can be also seen across individual compounds (Figure 2A). After ATP normalization, however, the individual differences between filter materials and compounds became less pronounced (Figure 2B, except IFW). When comparing drinking water sand group (IFW and BWA, named as DW) with the other two groups, their k_biol (mean value of individual substances in each group) remained different (Mann-Whitney-Wilcoxon test, W = 1305, p = 0.02 for HBG; W = 1318, p = 0.01 for LBG). However, when testing for differences in biotransformation of LBG and HBG, there was no significant difference in the mean global k_biol between the two groups (Mann-Whitney-Wilcoxon test, W = 1305, p = 0.33). This indicated their mean biotransformation performance were comparable per biomass unit, independent of other parameters such as microbial community composition. Differences between the k_biol of individual compounds persisted when comparing seven sand samples (IFW was excluded: the values close or below zero caused superimposed signals when normalized, Figure 2B) or HBG and LBG. Overall, 19 compounds showed significant differences in seven sand samples, and 15 substances showed significantly differential k_biol between HBG and LBG (Kruskal-Wallis test, adjusted p < 0.05) (Table 1). For comparison, we found 35 and 38 substances to differ in their k_biol before normalization for the sand filter materials and HBG vs. LBG, respectively.

We investigated the microbial community composition of all incubations using 16S rRNA amplicon sequencing and metagenomics. For the amplicon data, all sequences were clustered into a total of 23147 bacterial ASVs after the filtering step, showing obvious differences among the sand filters at the phylum level (Figure 3A). Apart from a high abundance of Proteobacteria in all filter materials, Nitrospirae dominated in the HBG with mean read abundances of 27.6%, followed by Bacteroidetes (11.5%). In the LBG, Nitrospirae only accounted for 0.7%. Microbial community composition was further cross-validated by microscopy (FISH) and metagenome classifications (METAXA2) (Figure 3A). In particular, the community composition matrix that resulted from the metagenome (METAXA2) was highly correlated with the composition matrix generated (at much deeper resolution) with the ASV amplicon data (Mantel tests, r = 0.96, p = 0.001). The FISH results confirmed that the relative proportion of cells and their respective biomass (assuming the same cell size) follow the results of the DNA based data (with exceptions; e.g., for the Nitrospirae community in Stuttgart).

At the genus level, a high degree of variation in the taxonomic composition between the sand samples was observed (Figure 3B, Supplementary Figure 3). In the HBG, there was a high amount of Nitrospira in Wangen and Eriskirchen with mean relative abundances of 45.2 and 39.0%, respectively. In Friedrichshafen (which was operated for post-denitrification), Nitrospira only accounted for 3.6% and other lineages became more abundant (e.g., Propionivibrio, Geobacter, Romboutsia, Legionella). In LBG, the three sand samples exhibited largely different microbial composition without sharing abundant genera. Stenotrophobacter (24.4%) was prominent in Hungen and Dokdonella (17.8%) dominated the Stuttgart filter materials. Pseudomonas (20.1%) and Janthinobacterium (17.3%) were dominant in Moos. In contrast to HBG and LBG, DW showed a stark dominance of Pseudomonas with a mean relative abundance of 60.1%, which also resulted in the lowest biodiversity indices (Table 2). In a multivariate analysis of the community composition (presented as ASV or as k-mers derived from the metagenome) we could recover a separation into HBG and LBG (Supplementary Figure 4, adonis: for ASV, R² = 0.21, p = 0.001; for metagenomics k-mers (kWIP), R² = 0.19, p = 0.001).

Finally, we confirmed that the community composition may be linked to the biotransformations observed in the systems using a Mantel test for the ASV matrix and the (normalized) k_biol values of the incubations (r = 0.50, p = 0.001). Subsequently, we were interested if there are any specific linkages between microbes and biotransformations. For this, Pearson's coefficients r were used to find hypothetical linkages between the k_biol and relative abundances of (a) microbial genera (Figure 4A) or (b) functional genes (Figure 4B, see paragraph below). In total, there were 62 genera (as a data reduction step, all ASV were collapsed to genera for this purpose) that showed either a significantly positive or negative correlation with the biotransformation rates of TOrCs (with a cutoff of abs(r) > 0.7, adjusted p < 0.05, observation > 9) (Figure 4A). All of these (except Pseudomonas) ranked below 2% relative abundance and can be regarded as rare taxa. Thirty genera showed highly positive correlations with more than one TOrC, such as Denitratisoma, Hydrogenophaga and Ideonella. Four genera correlated with single TOrCs, i.e., Litorilinea, Novosphingobium, Paludibaculum, Phaeodactylibacter were only positively associated with the biotransformation of diclofenac, acyclovir, sulfamethoxazole, and bezafibrate, respectively. Nitrospira, the dominant genus in the HBG, had no significant correlation with any compound (only when non-normalized k_biol values were considered we found multiple negative correlations). Examples for the negative correlation were Conexibacter, Ferruginibacter, Intestinibacter, Methylibium, Pseudomonas, and Terrimonas.

By comparison with the integrated enzyme database of KEGG, we identified 20409 EC numbers from the metagenome sequences. Principle component analysis (PCA) based on the relative abundance of annotated functions demonstrated, similar to the taxonomic patterns described above, distinct clustering of different sand samples (Figure 5). We can observe a separation between the HBG and LBG on axis 2 (PC2, 12%) with the exception of Friedrichshafen that shifted toward LBG, potentially caused by the low abundance of Nitrospira. PC1 (78%) distinguished two clusters of HBG, LBG and DW, which may due to the microbial structure differences between drinking water filters and wastewater filters. When comparing HBG with LBG, 1017 functions were over-represented in HBG and 1702 were over-represented in LBG, respectively (Supplementary Figure 5). From these, we selected those enzyme commission categories responsible for biocatalysis/biodegradation according to EAWAG-BBD database and analyzed their associations with TOrCs biotransformation. This resulted in 30 functions that showed a correlation with TOrCs biotransformation [with a cutoff of abs(r) > 0.7, adjusted p < 0.05, observation > 9] (Figure 4B). Notably, enoyl-CoA hydratase (EC 4.2.1.17)/3-hydroxyacyl-CoA dehydrogenase (EC 1.1.1.35) encoding genes were positively correlated with eight compounds (i.e., benzotriazole, carboxy acyclovir, ceterizin, DEET, lamotrigine, mecoprop, terbuthylazine, and torasemid). Both, dissimilatory sulfite reductase (EC 1.8.99.3) and nitrite reductase (EC 1.7.2.1) genes, which encode for universal and essential enzymes in the sulfur and nitrogen cycle, were also correlated to the biotransformation of carboxy acyclovir, clarithromycin, furosemide and mecoprop. Further, a few functions showed a correlation to only a single compound removal rate. For example, amidase (EC 3.5.1.4) encoding genes were found to be only correlated to the k_biol of carboxy acyclovir; threonine dehydratase (EC 4.3.1.19); acyl-CoA dehydrogenase, benzoate degradation ring-cleavage hydrolase, nitrate reductase (EC 1.7.2.1), S-formylglutathione hydrolase (EC 3.1.2.12) encoding genes were only correlated with the removal rate of atenolol.

Overall, there were 163 significantly differential abundant pathways annotated by HUMAnN2 comparing the HBG with LBG (see Supplementary Figure 6 for an overview excluding biosynthesis pathways). In the HBG, sulfate reduction, aromatic biogenic amine degradation and pathways regarding carbohydrate degradation (i.e. starch, glycogen, stachyose, D–galactose, galactose degradation) were overrepresented, the involved microorganisms were identified to be mainly Nitrospirae. In the LBG, energy metabolism pathways, such as the TCA cycle, Calvin-Benson-Bassham cycle, NAD/NADP–NADH/NADPH cytosolic interconversion, or the octane oxidation pathway were overrepresented.

We further inspected MAGs that taxonomically matched the genera for which we found correlations with individual k_biol (see section above, Figure 4A: genera with green bars) in order to see if these statistically as potential relevant genera contain biotransformation pathways. We annotated 37 MAGs in the KEGG Mapper for pathway reconstruction and we estimated the completeness of the pathways individually (Figure 6). In the category of xenobiotics metabolism, MAGs classified to Hydrogenophaga showed the most complete pathways, especially in benzoate and furfural degradation. Steroid degradation pathway was prominent in Sphingorhabdus and Pseudomonas.

We observed largely differing taxonomic and genetic profiles within each sand filter, accompanied by stark differences in the overall biotransformation performance (−1.4 vs. 58% removal) across all investigated TOrCs. The sand filter materials are influencing the biofilm formation and thus ATP and potentially microbial community composition, but they were not considered to have a direct association with (abiotical) TOrC transformations, for which their sorption effect can be ignored due to a) the polar substances and the rather low sorption affinity and b) the adsorption equilibrium from previous long-term exposure in the environment. This may point to the conclusion that functional microbial communities are the main predictor of TOrC transformations, however, most of these differences were eliminated after accounting for biomass, measured as ATP. Our results imply that the effect of biomass on the global TOrC transformation potential is foremost independent of the microbial community composition. We therefore had to partially reject our hypothesis that the microbial community composition is the major driver between global biotransformation potential in our experimental setup. However, on the other hand, the clear correlation of the compound matrix with the taxonomic matrix found by the Mantel test and later with single genera after normalization is evidence that although the global transformation potential is determined by biomass, the transformation of individual compounds is related to the taxonomic composition of the biological active sand filter system. After normalization, we could still identify significant differences for more than 15 different compounds (Table 1, e.g., clarithromycin, mecoprop, metoprolol, trimethoprim, furosemide, atenolol), which may be good candidates for a rather taxa or community specific degradation. Other studies already found indications that biomass is important for overall system performance, e.g., Liang et al. (2021) observed that the biomass increased with the running time of the reactors when also the TOrC removal increased, and Torresi et al. (2016) observed an enhanced TOrC removal with an increased biofilm thickness, although this feature was simultaneously attributed to increased diversity. In our study, however, alpha-diversity showed no such correlation (Table 2).

Our study indicates that we can predict the global transformation potential of (ex situ) materials from biofilters largely (explaining 76–84% of the k_biol variation, or 90% of the variation for compound removal) by measuring the living biomass as ATP concentration. Further studies will be needed to assess if the relationship between ATP and k_biol is linear or logarithmic [the two models were not significantly different (ANOVA, p = 0.24)]. Biomass is used as a linear parameter for TOrCs biotransformation kinetics (see Equation 1), and originally, the k_biol was determined by the absolute abundance of the functional degraders of the respective TOrC (Bekins et al., 1998). In mixed microbial communities, however, this relationship may have to be re-evaluated in the future.

When comparing frequently used biomass estimators for k_biol, ATP, a biomass parameter that is frequently used in microbiology, has not been considered for k_biol calculation of TOrC removal rates before. The biotransformation rate constants k_biol were previously normalized by attached or suspended biomass as dry weight (Casas et al., 2015; Mazioti et al., 2015; Torresi et al., 2016), total suspended solids (Achermann et al., 2018), or DNA concentration (Liang et al., 2021). However, most methods are accompanied by certain biases. For example, dead and dormant cells will strongly influence DNA measurements and cell counts; or the biofilm EPS will bias DNA measurement (through external DNA), dry weight, and total organic carbon measurement (LOI). In our study, living and total cells (biased only by the different cell volumes) were still the best predictors after ATP (Table 2). Our findings suggest that at least for biofiltration systems the global potential for the biotransformation of TOrCs is more dependent on ATP than on other biomass indicators. Hence, ATP could have a profound impact when comparing biotransformation results across studies and treatments.

While this conclusion will hopefully change the way we analyze TOrCs biotransformations, ATP itself has been frequently used as an important biomass indicator in drinking water biofiltration systems for natural organic matter (NOM) transformation (Pharand et al., 2014; Chen et al., 2016; Kirisits et al., 2019). This is not surprising since NOM is used as a substrate by the microbial communities in various environments such as soil, sediment, marine, and freshwater (Kolehmainen et al., 2007; Huang et al., 2011; Diem et al., 2013; Simon et al., 2013). However, TOrCs are not considered as substrates per se and their concentrations are by definition generally considered as too low for biomass maintenance and growth and therefore efficient degradation. Co-metabolism through promiscuity or mixed substrate use can be a mechanism that lead to TOrC removal (Rauch-Williams et al., 2010; Hellauer et al., 2019).

In general, we divided the filters into three groups according to the clustering results, microbial community composition, and biomass estimates. The taxonomic composition between the HBG, LBG, DW, and also between single sand filters was different. HBG, though, shared Nitrospirae as highly abundant lineage (Figure 3), a diverse and widespread group of often autotrophic, nitrite-oxidizing bacteria (Koch et al., 2015). Previous studies found that biotransformation of certain TOrCs can be related to ammonia oxidation activity of nitrifying activated sludge and biofilms in WWTPs (Helbling et al., 2012; Rattier et al., 2014; Men et al., 2017; Xu et al., 2017b). Moreover, asulam, carbendazim, fenhexamid, mianserin, and ranitidine showed biotransformation (16–85%) by the isolate Nitrospira inopinata (Han et al., 2019). Metagenome data of rapid gravity sand filter microorganisms also suggested that Nitrospira may serve as keystone species that drives the microbial ecosystems by providing organic carbon compounds and enable heterotrophic ammonium and carbon cycling (Palomo et al., 2016). In our system, we found no direct statistical linkage between Nitrospira and the biotransformation rate constants, however, it is possible that at least in the investigated biofilters Nitrospirae lineages as autotrophic species enables the establishment of a high microbial biomass, and thus enables the TOrC transformations through other microbial members. This also matches the observations by Liang et al. (2021), which indicates that the biotransformation may rely on rare community members of TOrCs-specific degraders, while the community of their moving bed reactor followed a progressive succession toward a Nitrospirae based climax community. The hypothesized positive and negative correlations of (biomass normalized) TOrC biotransformation with microbial genera pointed to rare biosphere microorganisms (<2%), which exhibited most of the correlations (>98%) with individual TOrCs removal rates. Although these correlations require experimental verification to test for causal relationships, these correlation-based hypotheses are in line with previous reports indicating that a small fraction of highly-specialized microorganisms, accounting for <0.1% of the microbial communities in biofilms, may be responsible for the TOrC transformation (Falås et al., 2018). One of the highly correlating genera was Hydrogenophaga, which showed a putative linkage to the removal of amisulprid, clarithromycin, climbazole, flecainide, oxypurinol, sitagliptine, sulpirid, terbutryn and venlafaxine (the top left cluster in Figure 4A). This genus has also been previously reported to efficiently remove diclofenac, metoprolol, clarithromycin, erythromycin, atenolol and codeine (Kanaujiya et al., 2019).

In the drinking water sand filters DW, the biotransformation performance significantly differed from the wastewater filters. Here, Pseudomonas was found dominating DW, a common inhabitant of aquatic environments, including oligotrophic drinking water, lake water and surface water (Lopez et al., 2005; Huang et al., 2015; Nasreen et al., 2015). Although Pseudomonas has been reported to be able to degrade certain micropollutants (Li et al., 2010; Tezel et al., 2012; Devi et al., 2019), in our case, it also clustered in proximity to Hydrogenophaga, however, with more pronounced negative correlations with TOrCs.

The relative abundance of functional genes and encoded biotransformation pathways also correlated with individual compounds and resulted in a clearer clustering than the microbial lineages. Enoyl-CoA hydratase/3-hydroxyacyl-CoA dehydrogenase was overrepresented in the HBG and was positively correlated with the k_biol of eight TOrCs in our study. The enoyl-CoA hydratase is known to catalyze a β-oxidation substrate by adding hydroxyl groups and a proton to an unsaturated β-carbon of the molecule (Salgado et al., 2020). For instance, the degradation of ibuprofen can be initiated by enoyl-CoA hydratase with the introduction of hydroxyl groups, this enzyme was found to be upregulated during the biodegradation process (Almeida et al., 2013). Enoyl-CoA hydratase was also identified in steroid estrogen degradation by bacterium Serratia nematodiphila (Zhao et al., 2020). Furthermore, Cameron et al. (2019) reported enoyl-CoA hydratase contributed to biofilm formation and the antibiotic tolerance, which also supported our findings that the high biomass promoted the TOrCs removal. Therefore, we suggest that a high abundance of enoyl-CoA hydratase/3-hydroxyacyl-CoA dehydrogenase might indicate a favorable functional potential for TOrCs biotransformation. Dissimilatory sulfite reductase (EC 1.8.99.3) and nitrite reductase (EC 1.7.2.1) were found to be correlated to the biotransformation of carboxy acyclovir, clarithromycin, furosemide and mecoprop in our study. In previous studies, sulfite reductase was reported to catalyze the cleavage of isoxazole and piperazinyl rings (Jia et al., 2019), and nitrite reductase was correlated with the biotransformation rate constant of some compounds like sulfamethoxazole, erythromycin and trimethoprim (Torresi et al., 2018).