Filter material and water samples were collected from Islevbro waterworks (Copenhagen, Denmark), a groundwater-fed rapid sand filter system described previously (Lee et al., 2014). Briefly, raw water abstracted from a deep chalk aquifer undergoes aeration, iron oxidation in a retention tank (~20 min contact time), and dual-stage filtration (pre-filter for Fe-hydroxide retention; after-filter for biological ammonium removal). After-filters (~0.7 m depth, 1 mm grain diameter, ~30 years operation, backwashed every 14 days) receive influent with ~0.13 mg NH₄⁺-N/L, 9.3 mg/L dissolved oxygen (DO), pH 7.3, and 9–11 °C temperature, achieving ammonium removal to <0.01 mg/L. Filter material from one after-filter was core-sampled midway between backwash events using a 60 cm plexiglass cylinder, extruded, and aseptically sliced into depth sections: 0–10 cm (top), 10–40 cm (middle), and 40–50 cm (bottom). Top and bottom materials were used for column experiments and stored wet at 4 °C during transportation. Subsamples for DNA extraction were stored at −80 °C.

A lab-scale column assay, detailed elsewhere (Tatari et al., 2013), was used to investigate ammonium removal under varied conditions. Briefly, plexiglass columns (5 cm bed height, 2.6 cm inner diameter) were packed with sand and continuously fed with water. Two independent experiments, each comprising eight parallel columns, were operated for 30 days (Table 1 and Supplementary Figure S1). Column systems were placed in a temperature-controlled room set at either 10 °C (Experiment 2) or kept at room temperature (mean ≈ 20 °C; Experiment 1). For Experiment 1 (20 °C), columns were packed with either top- or bottom-layer material and fed with after-filter effluent water spiked with ammonium (as NH₄Cl) at a reference loading rate of 35 g NH₄⁺-N/m³ filter material/d (1.46 g NH₄⁺-N/m³/h; 1 mg/L NH₄⁺-N at 0.96 L/d flowrate), equivalent to full-scale conditions. Additional loading rates included 0.1 × reference (0.1 mg/L NH₄⁺-N) and 5 × reference, achieved either by increasing influent NH₄⁺-N concentration to 5 mg/L (oxygen-limited condition) or by increasing flowrate fivefold (non-limiting oxygen condition). The high-concentration treatment was designated oxygen-limited because complete nitrification of 5 mg/L NH₄⁺-N requires 22.9 mg O₂/L (stoichiometric demand of 4.57 mg O₂ per mg NH₄⁺-N oxidized), exceeding the measured influent dissolved oxygen of 10.6 ± 0.8 mg/L. For Experiment 2 (10 °C), columns packed with top-layer material were fed with three different waters to evaluate organic carbon effects on nitrifier-heterotroph interactions: after-filter effluent spiked with NH₄Cl (minimal organic carbon control), after-filter effluent spiked with NH₄Cl + 1.5 mg/L sodium acetate (labile carbon amendment), or pre-filter effluent spiked with NH₄Cl (containing natural groundwater dissolved organic matter (DOM)). In this experiment, 5× loading rates were set by increasing NH₄⁺-N concentration to 5 mg/L (oxygen-limited condition as described in Experiment 1). Four times during each experiment, 2 g of filter material was sampled from each column for microbial characterization (qPCR, 16S rRNA gene amplicon sequencing), increasing the volumetric loading rate by 52 ± 13% over the experimental period due to material removal.

Duplicate water samples were collected along the experiments every 2–4 days, filtered (Sartorius Minisart 0.20 μm), and stored at −20 °C. In Experiment 1, NH₄⁺, NO₂⁻, and NO₃⁻ were analyzed using an autoanalyzer (Bran+Luebbe Analytics, 2012). In Experiment 2, NH₄⁺ was measured using a salicylate-hypochlorite method (Bower and Holm-Hansen, 1980), NO₂⁻ via an adapted Grasshoff et al. (1983) method, and NO₃⁻ using a Merck Spectroquant test kit 109,713. Different methods were used due to equipment availability, with cross-calibration ensuring consistency. DO (influent and effluent) was periodically measured with a handheld meter (WTW, Multi 3,430, with FDO® 925) during Experiment 2.

Filter material from the top of the columns was sampled four times during the 30-day period, with the final sample at day 30. Samples were collected in cryotubes, flash-frozen in liquid nitrogen, and stored at −80 °C. DNA was extracted from the 0.5 g sand collected just before and at the end of the experiments using the FastDNA® Spin Kit for Soil (MP Biomedicals, Santa Ana, CA, USA) per the manufacturer’s instructions. DNA concentration and purity were assessed using a NanoDropTM 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA). Quantitative PCR (qPCR) assays were performed in duplicate using a Chromo 4 Thermal Cycler (Bio-Rad Laboratories, Hercules, CA, USA). Each reaction contained 12.5 μL 2 × iQ SYBR Green Supermix, 20 μM forward and reverse primers, 10 ng DNA template, and PCR-grade water. Targeted groups included total bacteria (16S rRNA gene), Nitrospira (specific 16S rRNA region), AOB (specific 16S rRNA region), and AOA (amoA gene) (Supplementary Table S1). Cell density calculations assumed one copy per cell, except for the 16S rRNA of the total community, where copy numbers were estimated using CaRcone (R script¹). DNA was PCR-amplified using primers PRK341F (5′-CCTAYGGGRBG CASCAG-3′) and PRK806R (5′-GGACTACNNGGGTATCTAAT-3′) for 35 cycles to target the V3-V4 hypervariable region (Yu et al., 2005). PCR products were purified and sequenced on the Illumina MiSeq platform at the DTU MultiAssay Core Centre (Lyngby, Denmark).

To evaluate physiological adaptation versus biomass growth, apparent cell-specific ammonium oxidation rates (SAOR; fmol NH₄⁺-N/cell/h) were estimated for the initial and final experimental phases. SAOR was calculated by dividing the volumetric ammonium removal rate (g NH₄⁺-N/m³/h) by the total nitrifier cell density (sum of AOB, AOA, and Nitrospira) converted to a volumetric basis using layer-specific filter material bulk densities of 1.1 g/cm³ for the top layer and 1.7 g/cm³ for the bottom layer, as previously reported for this system (Lee et al., 2014). Initial SAOR was calculated using average removal rates from the first 48 h and the nitrifier density of the corresponding inoculum (top or bottom material). Final SAOR used average removal rates from days 28–30 and the column-specific cell densities measured at day 30.

16S rRNA gene amplicon libraries were processed using the DADA2 pipeline (Callahan et al., 2016), which outputs the abundance of error-corrected amplicon sequence variants (ASVs). ASVs were classified with the SILVA prokaryotic reference database v132. After filtering, 1.2 million sequences were retained, averaging 65,000 sequences per sample. Further analysis was carried out in R packages phyloseq (McMurdie and Holmes, 2013) and ampvis2 (Andersen et al., 2018). Statistically significant differences in sequences before and after 30-day incubations were identified using DESeq2 (Wald significance test, parametric fit type, padj < 0.05) (Love et al., 2014).

The 16S rRNA gene generally provides limited phylogenetic resolution for distinguishing comammox Nitrospira from canonical nitrite oxidizers (Daims et al., 2015; Lawson and Lücker, 2018). However, extensive prior characterization of the filter used for our column experiment, including metagenomic sequencing, genome binning, and recovery of Nitrospira metagenome-assembled genomes (MAGs) (Palomo et al., 2016, 2018, 2022a) and amoA-targeted qPCR and sequencing (Fowler et al., 2018), enabled us to develop a robust classification framework for Nitrospira ASVs based on phylogenetic placement and cross-validation with functional gene data. MAGs previously recovered from this system had been functionally annotated for comammox-specific genes, including ammonia monooxygenase (amoA) and hydroxylamine dehydrogenase (hao), enabling reliable distinction between comammox and canonical NOB, as well as between comammox clades A and B. For each MAG, we extracted the corresponding 16S rRNA gene sequence (when available) and aligned these sequences using MUSCLE v3.8.31 (Edgar, 2004) together with Nitrospira-affiliated ASVs from the present study and 16S rRNA gene sequences retrieved from publicly available Nitrospira genomes. These alignments were used to construct a maximum-likelihood phylogenetic tree using RAxML v8.2.11 (Stamatakis, 2014) with 550 rapid bootstraps (determined using the autoMRE option) and the GTRCAT substitution model (best model determined using jModelTest v.2.1.10; Posada, 2008). The tree was rooted using Leptospirillum species as outgroup and visualized using the Interactive Tree of Life (iTOL) web tool (Letunic and Bork, 2016). Each ASV was assigned to a category (comammox clade A, comammox clade B, or canonical NOB) based on its closest phylogenetic placement relative to MAG-derived 16S rRNA genes with known metabolic identity. Assignments were validated by BLASTn searches against the NCBI Whole Genome Shotgun contigs (WGS) database, confirming that top hits aligned with tree-based classifications (Supplementary Table S2). To further validate our classifications, we compared the relative abundances of comammox versus canonical Nitrospira, and comammox clade A versus clade B proportions, derived from 16S rRNA gene amplicon sequencing in this study with independent estimates from previous investigations of the same filter using: (1) metagenomic read recruitment, (2) qPCR targeting Nitrospira 16S rRNA genes and comammox-specific amoA genes to determine the proportion of comammox within total Nitrospira, and (3) comammox amoA amplicon sequencing to distinguish clade A versus clade B ratios. Amplicon-derived, qPCR-derived, and metagenome-derived abundance estimates showed strong concordance (Supplementary Table S3). Furthermore, direct comparison of individual MAG and ASV abundances revealed strong agreement (R² = 0.93, Supplementary Figure S2), including consistent rank ordering of dominant clade B, intermediate clade A, and low-abundance canonical NOB lineages. This multi-method validation supports the reliability of our 16S rRNA gene-based Nitrospira classifications in this well-characterized system.

Statistical analyses were performed in R version 4.5.1. Volumetric removal rates (g NH₄⁺-N/m³/h) were analyzed using linear mixed-effects models (lme4 package; Bates et al., 2015) with column identity as a random intercept to account for repeated measurements. A primary model included loading rate, temperature, layer origin, oxygen availability, and water source as fixed effects: Removal Rate ~ Loading_Rate + Temperature + Layer + Oxygen + Water_Source + (1|Column). Because the experimental design was partially confounded, focused subset analyses were conducted to obtain unbiased factor estimates, using Welch’s t-tests for pairwise contrasts. Fixed-effect significance was assessed by Type III ANOVA with Satterthwaite’s approximation (package lmerTest; Kuznetsova et al., 2017). Post hoc comparisons were performed using estimated marginal means with Tukey’s HSD adjustment (package emmeans; Lenth and Piaskowski, 2017). Model performance was evaluated based on marginal and conditional R² values. To analyze differences in microbial cell densities, data were log₁₀-transformed and analyzed using separate linear models for each microbial group. These models included temperature, layer origin, oxygen availability, loading rate, and water source as fixed effects. The significance of factors was assessed using a Type III ANOVA, with significant differences explored through pairwise comparisons of estimated marginal means adjusted with Tukey’s HSD. For all models, assumptions of normality (Shapiro–Wilk test) and homoscedasticity (Breusch-Pagan test) were met (all p > 0.17). Fold-changes in cell density relative to the source material (FC = Final/Source) were calculated, log₂-transformed, and analyzed using linear models. For all statistical tests, significance was set at α = 0.05.

Rapid sand filters have been operated for groundwater treatment at Islevbro waterworks since 1923, achieving complete ammonium removal from influent concentrations of ~0.13 ± 0.05 mg NH₄⁺-N/L to < 0.01 mg NH₄⁺-N/L in effluents. Lab-scale columns packed with top-layer (0–10 cm depth) or bottom-layer (40–50 cm depth) filter material from Islevbro were operated under varying conditions to investigate RSF microbial community responses to short-term disturbances, including temperature (10 °C vs. 20 °C), oxygen availability (adequate vs. limiting), water source (after-filter effluent, acetate-amended after-filter effluent, or after-filter influent), and ammonium loading rates (0.1×, 1×, 5 × reference loading of 1.46 g NH₄⁺-N/m³/h, equivalent to full-scale conditions) (Table 1). Results showed that ammonium loading emerged as the primary determinant of volumetric removal capacity (p < 0.001), with rates scaling proportionally from 0.3 ± 0.1 g NH₄⁺-N/m³/h at 0.1 × loading to 7.1 ± 3.3 g NH₄⁺-N/m³/h at 5 × loading under non-oxygen-limiting conditions (Figure 1A). Columns packed with top-layer material achieved peak rates of 8.7 ± 2.9 g NH₄⁺-N/m³/h with 94 ± 3% removal efficiency under high-flow, 5 × loading conditions (Figure 1B). Under identical ammonium loading but oxygen-limiting operation, removal rates dropped to 3.7 ± 0.9 g NH₄⁺-N/m³/h (p < 0.05; Figure 1B). This oxygen limitation was stoichiometric rather than kinetic: complete oxidation of 5 mg NH₄⁺-N/L requires 22.9 mg O₂/L, whereas influent contained only 10.7 ± 0.2 mg/L of oxygen. Effluent oxygen concentrations confirmed near-complete oxygen depletion in O₂-limited columns (1.1 ± 0.7 mg/L) compared with oxic controls (4.8 ± 0.3 mg/L), suggesting operation at maximal O₂-limited capacity.

Filter layer origin also significantly affected nitrification performance (p < 0.05), with top-layer columns consistently outperforming those packed with bottom-layer material, particularly at the high loading (Figure 1). Notably, bottom-layer removal rates and efficiencies increased after the first few days (Supplementary Figure S3). By contrast, temperature and water source had no significant effects on ammonium removal (p > 0.6). Across all treatments, nitrification proceeded to completion with negligible nitrite accumulation, except during a transient lag phase in the first 24–48 h of high loading operation (Supplementary Figure S4).

The top and bottom sections of the full-scale RSF exhibited pronounced differences in community composition and microbial abundance prior to incubation (Figure 2A). Amplicon sequencing showed that the top-layer community was dominated by Proteobacteria (41%), Nitrospirae (29%), and Acidobacteria (13%). In contrast, the bottom layer community contained a higher proportion of Acidobacteria (22%) and a markedly lower abundance of Nitrospirae (4%) (Figure 2A). These differences were particularly evident among nitrifiers: Nitrospira spp. represented 28.8 ± 1.2% of the initial top-layer community but only 3.9% in the bottom layer (Figure 2B). Nitrosomonas spp. (AOB) were likewise in higher abundance in the top layer (0.9 ± 0.6%) relative to the bottom (0.1%), while ammonia-oxidizing archaea (AOA) constituted only a minor and similar fraction (0.03–0.1%) in both layers. Quantitative PCR confirmed that these compositional contrasts were mirrored by differences in absolute cell abundance. The top-layer material contained sixfold more total cells/g than the bottom layer (Figure 3). Among nitrifiers, Nitrospira abundance was 66-fold higher and AOB 11-fold higher in the top vs. bottom layer, whereas AOA differed only modestly (1.6-fold; Figure 3). These data establish that both the taxonomic composition and the biomass of nitrifying populations were stratified in the source filter material, forming the foundation for subsequent treatment-dependent dynamics described below.

Linear modeling of final cell densities showed that filter layer origin remained the dominant factor controlling total microbial abundance (p < 0.001), with top-layer columns maintaining significantly higher densities (average of 3.6 × 10⁹ cells/g) than bottom-layer columns (7.9 × 10⁸ cells/g) throughout the experiment (Figure 3). Temperature (p = 0.41), oxygen availability (p = 0.49), water source (p = 0.19), and loading rate (p = 0.68) had no significant effect on total microbial density (Supplementary Table S4). Fold-change analysis revealed distinct growth dynamics between layers. Bottom-layer columns at 20 °C showed a 49% increase in total cell density, coinciding with the progressive improvement in ammonium removal (Figure 3). In contrast, top-layer columns at 20 °C exhibited a slight decline, while those at 10 °C showed modest growth (Figure 3). Among nitrifiers, Nitrospira spp. were the most strongly influenced by layer origin (p < 0.001). Bottom-layer columns showed marked Nitrospira proliferation, with (3.4–4.4)-fold increases in cell density (Log₂FC = 1.7–2.1; p < 0.01; Figure 3). In contrast, top-layer columns at 20 °C exhibited slight declines, whereas those at 10 °C remained stable (Figure 3). Despite the substantial Nitrospira growth in bottom layers, cell densities at the end of experiments in top-layer columns (average of 7.9 × 10⁸ cells/g) were still 19-fold higher than those in bottom-layer columns (4.2 × 10⁷ cells/g). Water source showed a marginal effect on Nitrospira abundance (p = 0.13), with after-filter effluent-fed columns maintaining slightly higher densities compared to influent (pre-filter effluent) or acetate-amended after-filter effluent. Neither temperature (p = 0.76), oxygen availability (p = 0.67), nor loading rate (p = 0.32) significantly affected Nitrospira densities (Supplementary Table S4). Patterns of AOB abundance closely mirrored those of Nitrospira, with layer origin as the dominant factor (p < 0.001). Top-layer columns sustained 3.6-fold higher AOB densities (2.1 × 10⁷ cells/g) than bottom-layer columns (6.0 × 10⁶ cells/g; Figure 3). Similar to Nitrospira, bottom-layer communities showed substantial AOB proliferation, with (3.4–5.6)-fold increases (Log₂FC = 1.8–2.5; p = 0.017). Likewise, top-layer columns at 20 °C exhibited slight AOB declines, whereas those at 10 °C showed modest increases (Figure 3). Neither temperature, oxygen availability, water source, or loading rate significantly affected AOB densities (all p > 0.17) (Supplementary Table S4).

In contrast to bacterial nitrifier behavior, AOA exhibited a fundamentally different set of responses. Temperature was the dominant factor influencing AOA abundance (p = 0.002), with 10 °C columns supporting an average 3.5-fold higher cell densities (3.0 × 10⁶ cells/g) than 20 °C columns (8.3 × 10⁵ cells/g; Figure 3). Layer origin, which dominated bacterial nitrifier distributions, had no significant effect on AOA abundance (p = 0.70), but these nitrifiers responded significantly to oxygen availability (p < 0.05), with 1.7-fold higher abundances under oxygen-limited conditions. Additionally, water source also influenced AOA populations (p < 0.05), with acetate-amended after-filter effluent supporting slightly lower densities compared to unspiked after-filter effluent or after-filter influent water. Loading rate showed no significant effect on AOA abundance (p = 0.37) (Supplementary Table S4). Collectively, these data indicate concurrent growth of all three nitrifier guild members under most experimental conditions, with no evidence of short-term competitive exclusion. Bacterial nitrifiers were primarily structured by filter layer origin and varied little across temperatures, while AOA exhibited a layer-independent, temperature-driven dynamics and were favored under oxygen-limited conditions.

To determine whether enhanced ammonium removal was driven by increased nitrifier biomass or elevated cellular activity, we calculated cell-specific ammonium oxidation rates (SAOR). Initially, bottom-layer nitrifiers exhibited an order of magnitude higher per-cell activity than their top-layer counterparts across all loading conditions (Supplementary Figure S5). Despite this intrinsic physiological advantage, the substantially greater nitrifier biomass in the top layer drove superior bulk ammonium removal. Over the 30-day incubation, the response of SAOR was highly dependent on the ammonium loading. For bottom-layer communities, high loading triggered a significant physiological upregulation, elevating SAOR from 2.84 ± 0.34 to 4.27 ± 0.56 fmol NH₄⁺-N/cell/h. In contrast, at reference loading, SAOR remained statistically stable (from 1.10 ± 0.19 to 1.04 ± 0.14 fmol NH₄⁺-N/cell/h), indicating that improved ammonium removal rates (Supplementary Figure S3) were primarily driven by the observed nitrifier proliferation (Figure 3). For top-layer communities, where nitrifying biomass did not significantly change over the incubation period (Figure 3), the observed increase in per-cell activity (Supplementary Figure S5) appears to be the primary factor responsible for the improved ammonium removal rates (Supplementary Figure S3).

RSF community composition changed substantially over the 30-day period based on 16S rRNA gene amplicon sequencing, with the magnitude of restructuring mainly driven by temperature and layer origin (Figure 4). Columns operated at 20 °C experienced more pronounced compositional shifts from the initial communities than those at 10 °C, and bottom-layer columns underwent the most dramatic community reorganization (Figure 4). Sequencing results supported the qPCR-based quantitative trends and enabled resolution of clade-specific dynamics within Nitrospira, the dominant nitrifiers in all columns (Figure 5). The relative abundance of Nitrospira spp. increased significantly in bottom-layer columns (from 3.9% to 9.8 ± 2.2%), consistent with the substantial absolute growth detected by qPCR, while they remained similar in top-layer columns (from 28.8 ± 1.2% in the initial filter material to 26.0 ± 3.6% at the end of the experiment; Figure 2B). Comammox Nitrospira dominated the Nitrospira population according to 16S rRNA sequencing (Figure 5), consistent with previous metagenomic (92 ± 2% in top, and 95 ± 1% in bottom) and qPCR (70–99%) findings from the source filter (Fowler et al., 2018). Within this group, clade B was predominant (~85%), with clade A comprising ~15%, also in line with prior metagenomics and amoA amplicon data (Fowler et al., 2018). Clade-level dynamics, however, differed markedly between temperature treatments and filter-layer origin (Figure 5). In bottom-layer columns, differential abundance analysis identified two comammox clade B sequence variants (ASVs 16S_1076 and 16S_1077) that increased significantly during the incubation (padj < 0.05; Figure 5). In contrast, in top-layer columns we observed the opposite trend. At 10 °C, a consistent compositional adjustment occurred, characterized by a significant decline in the dominant clade B variant (16S_1076; from 14.1% to 9.6 ± 1.8%, padj < 0.0001) and a concurrent 2.4-fold enrichment of a clade A variant (16S_1079; from 1.4% to 3.4 ± 0.6%, padj < 0.0001), irrespective of ammonium loading, oxygen regime, or water source (Figure 5). Conversely, in top-layer columns operated at 20 °C, the same clade A variant slightly declined in relative abundance (from 2.5% to 1.4 ± 0.5%; Figure 5). No major compositional changes were detected among canonical Nitrospira ASVs.

Beyond nitrifiers, several taxa exhibited notable enrichment patterns. Hydrogenophaga, Sphingobium, and Ferribacterium showed increased relative abundances across all columns to varying degrees (Figure 2B). These genera, undetected in source filter material, showed substantial increases in columns operated at higher temperature, reaching highest relative abundances in the bottom-layer columns (Hydrogenophaga: 9.9 ± 4.5%; Sphingobium: 7.2 ± 2.5%; Ferribacterium: 4.5 ± 1.7%; Figure 2B). Other taxa exhibited condition-specific enrichment patterns. For instance, acetate amendment selectively stimulated Acinetobacter and Pseudomonas spp. in columns operated at reference loading rates (Figure 2B), while high ammonium loading conditions enriched Novosphingobium spp. (Figure 2B).

Together, these patterns reveal temperature- and substrate loading-dependent niche partitioning within comammox Nitrospira and accompanying shifts among heterotrophs, highlighting the rapid adaptive restructuring of RSF microbial communities under short-term environmental perturbations.