Metabolic Diversity and Aero-Tolerance in Anammox Bacteria from Geochemically Distinct Aquifers

ABSTRACT Anaerobic ammonium oxidation (anammox) is important for converting bioavailable nitrogen into dinitrogen gas, particularly in carbon-poor environments. However, the diversity and prevalence of anammox bacteria in the terrestrial subsurface—a typically oligotrophic environment—are little understood. To determine the distribution and activity of anammox bacteria across a range of aquifer lithologies and physicochemistries, we analyzed 16S rRNA genes and quantified hydrazine synthase genes and transcripts sampled from 59 groundwater wells and metagenomes and metatranscriptomes from an oxic-to-dysoxic subset. Data indicate that anammox and anammox-associated bacteria (class “Candidatus Brocadiae”) are prevalent in the aquifers studied, and that anammox community composition is strongly differentiated by dissolved oxygen (DO), but not ammonia/nitrite. While “Candidatus Brocadiae” diversity decreased with increasing DO, “Candidatus Brocadiae” 16S rRNA genes and hydrazine synthase (hzsB) genes and transcripts were detected across a wide range of bulk groundwater DO concentrations (0 to 10 mg/L). Anammox genes and transcripts correlated significantly with those involved in aerobic ammonia oxidation (amoA), potentially representing a major source of nitrite for anammox. Eight “Candidatus Brocadiae” genomes (63 to 95% complete), representing 2 uncharacterized families and 6 novel species, were reconstructed. Six genomes have genes characteristic of anammox, all for chemolithoautotrophy. Anammox and aerotolerance genes of up to four “Candidatus Brocadiae” genomes were transcriptionally active under oxic and dysoxic conditions, although activity was highest in dysoxic groundwater. The coexpression of nrfAH nitrite reductase genes by “Candidatus Brocadiae” suggests active regeneration of ammonia for anammox. Our findings indicate that anammox bacteria contribute to loss of fixed N across diverse anoxic-to-oxic aquifer conditions, which is likely supported by nitrite from aerobic ammonia oxidation. IMPORTANCE Anammox is increasingly shown to play a major role in the aquatic nitrogen cycle and can outcompete heterotrophic denitrification in environments low in organic carbon. Given that aquifers are characteristically oligotrophic, anammox may represent a major route for the removal of fixed nitrogen in these environments, including agricultural nitrogen, a common groundwater contaminant. Our research confirms that anammox bacteria and the anammox process are prevalent in aquifers and occur across diverse lithologies (e.g., sandy gravel, sand-silt, and volcanic) and groundwater physicochemistries (e.g., various oxygen, carbon, nitrate, and ammonium concentrations). Results reveal niche differentiation among anammox bacteria largely driven by groundwater oxygen contents and provide evidence that anammox is supported by proximity to oxic niches and handoffs from aerobic ammonia oxidizers. We further show that this process, while anaerobic, is active in groundwater characterized as oxic, likely due to the availability of anoxic niches.

KEYWORDS aquifer, groundwater, anammox, "Candidatus Brocadiae", ammonia oxidizers, aero-tolerance F or decades, loss of fixed nitrogen (N) from freshwater environments, including aquifers was exclusively attributed to heterotrophic denitrification (1). Ammonium (NH 4 1 ) was deemed inert under anoxic conditions (2), and denitrification was considered the dominant anoxic sink for inorganic nitrogen (1). This view was altered when anaerobic ammonium oxidation (anammox), mediated by autotrophic bacteria from a deep-branching class within the phylum Planctomycetes, was discovered (3). Currently, only 10 anammox species have been enriched in culture (4); none are axenic. Their physiological characteristics and niches are defined using Monod model parameters (5). Anammox bacteria oxidize ammonium with nitrite (NO 2 2 ) to N 2 gas without nitrous oxide emissions (6) and are considered chemoautotrophic (7). Therefore, they tend to outcompete heterotrophic denitrifiers in low-carbon systems (8). Here, we explored the significance of anammox in aquifers. As typically oligotrophic environments (9), aquifers may represent ideal ecosystems for chemolithoautotrophic anammox bacteria.
A growing body of research now shows that anammox bacteria play an important role in the global aquatic N cycle, based on their discovery in diverse aquatic environments, including oceanic oxygen-minimum zones (10,11), meromictic lakes (12), marine subsurface sediments (13), uncontaminated suboxic groundwater (8), and contaminated groundwaters (14)(15)(16). Of these, aquifers harbor the largest freshwater store on Earth and can be characterized by long groundwater residence times (17), suitable for slow-growing anammox bacteria (with doubling times of 7 to 22 days) (18). Moreover, while groundwater often contains low levels of naturally occurring N species, such as nitrate (,0.25 mg/L) (19) and ammonia (,0.2 mg/L) (20), and frequent inputs of N from agricultural sources (21), aquifers tend to contain little organic carbon (9), potentially limiting heterotrophic denitrification (1) and favoring anammox bacteria. Anammox bacteria may also be able to compete with denitrifiers for organic carbon. Enrichment experiments have demonstrated that some anammox bacteria have the ability to oxidize small organic acids, such as propionate and acetate, with NO 3 2 /NO 2 2 (22), pointing to greater metabolic versatility than previously assumed for these chemolithoautotrophic organisms.
Anammox relies on the presence of both oxidized and reduced inorganic N compounds and, subsequently, is highly active in redox transition zones, where N species fluctuate between oxidation and reduction reactions (23). Anammox bacteria are therefore found naturally occurring at the interface of aerobic-anaerobic conditions, where NO 2 2 -containing water and anaerobic water with NH 4 1 mix (24). N fluctuation frequently occurs in groundwater due to interactions with surface water or to oxygen penetration from above in unconfined shallow aquifers (25). Elevated groundwater tables, which promote oxygenation through interactions with the overlying soil and groundwater recharge, stimulate a variety of anammox bacteria (16). Surface water contaminated with nitrogen species can be a source of NH 4 1 in aquifers. This ammonium can be partially absorbed by clay minerals (26) or oxidized to NO 2 2 or NO 3 2 by nitrification (1). NO 2 2 produced by aerobic ammonia oxidizers can serve as an electron acceptor in anammox (27). While anammox bacteria have been observed previously in groundwaters (8), their prevalence and genetic diversity are unclear. It is also not known how they are impacted by the heterogeneous chemical conditions found in aquifers.
To provide insights into anammox ecology and activity across a diverse range of groundwater ecosystems, we collected groundwater from 59 sites over 4 geographic regions across New Zealand. To the best of our knowledge, this is the largest survey of anammox communities in groundwater to date, encompassing a range of characteristics, including variations in nitrogen species and in organic carbon and dissolved oxygen concentrations. By reconstruction and metabolic characterization of diverse novel "Candidatus Brocadiae" genomes, and by quantifying the distribution, abundance, and activity of anammox bacteria across distinct groundwater conditions, we reveal the diverse ecological niches of these bacteria in groundwater. Results provide evidence in support of anammox as an important biological process in aquifers and a major N sink (16).

RESULTS AND DISCUSSION
Prevalence and activity of anammox in chemically diverse groundwaters. (i) Anammox bacteria are widespread in groundwater and phylogenetically diverse. All known anammox bacteria associate exclusively with order "Candidatus Brocadiales" (28) in the class "Ca. Brocadiae." Based on 16S rRNA gene amplicon analysis of 80 groundwater samples collected across 10 aquifers, "Ca. Brocadiae" was present in 60 samples (from 47 wells) across all four regions assessed in this study. Results demonstrate "Ca. Brocadiae" prevalence within aquifers comprising unconsolidated materials (sandy gravel = 55/71; sand/silt = 1/1), which are a common aquifer type globally (29), and consolidated volcanic rock (basalt/ignimbrite = 4/7) ( Fig. 1a; Table S1). Furthermore, "Ca. Brocadiae" was present in groundwaters encompassing a range of nitrate N (0-22 g/m 3   g/m 3 ), and dissolved organic carbon (DOC, 0 to 26 mg/L) concentrations (based on bulk groundwater measurements) (Table S1). The wide range of groundwater DO concentrations across which members of the class "Ca. Brocadiae" were detected is consistent with findings from two pristine carbonate-rock aquifers in Germany (8). There, anammox was presumed to outcompete denitrification due to low organic carbon and contributed to 83% of total N loss. Overall, "Ca. Brocadiae" represented 0.37% of groundwater prokaryotic communities (3,940/1,071,440 sequences; ;49.3 sequences 6 202 [standard deviation {SD}] per sample) and ranked 36th in abundance out of all 205 classes of bacteria and archaea. Research has shown that anammox is globally common in soils overlying aquifers, periodically saturated by high groundwater tables (16). Results here demonstrate that anammox-related bacteria are also common in the sediments and rock constituting aquifer saturated zones.
(ii) "Ca. Brocadiae" distribution or relative abundances were greater in suspended than attached aquifer fractions. Microbial community compositions can vary considerably between suspended or planktonic biomass in groundwater and biomass attached as biofilms on aquifer sediments and rocks (32). When comparing groundwater (the planktonic fraction) and biomass-enriched groundwater postsonication (the planktonic 1 biofilm fraction), we found that class "Ca. Brocadiae" 16S rRNA gene amplicon sequences were recovered from almost all groundwater samples (seven of eight wells sampled for both fractions) but only half postsonication ( Fig. 1c; Table S1). This difference is mostly explained by higher "Ca. Brocadiae" relative abundances in groundwater communities in the planktonic versus biofilm fractions at higher DO concentrations. Thus, "Ca. Brocadiae" represented a greater fraction of subsurface communities in groundwater than attached to aquifer materials. However, we identified no significant difference between these sample types (planktonic and biofilm) in terms of "Ca. Brocadiae" composition based on linear discriminant analysis effect size (LEfSe) analysis (Kruskal-Wallis test, P . 0.05). Instead, oxygen was associated with a greater effect overall, as sites with lower DO concentrations contained greater "Ca. Brocadiae" diversity in both sample types, along with a substantially greater number of "Ca. Brocadiae" amplicon sequences (Fig. 1b and c and 2a).
(iii) "Ca. Brocadiae" relative abundance and community composition vary strongly with dissolved oxygen. DO concentrations ranged from oxic to anoxic and are characterized here as anoxic (0 mg/L), suboxic (,0.3 mg/L), dysoxic (0.3 to 3 mg/L), and oxic (.3 mg/L) (33). They, along with other geochemical parameters, were found to be generally stable through time in a subset of eight groundwater wells that were sampled twice, 3 to 4 months apart (Table S1). DO was positively correlated with oxidation-reduction potential (ORP), sulfate, and nitrate N and negatively correlated with ammoniacal N (detected in only 24/80 samples), phosphate, dissolved reactive phosphorus (DRP), pH, temperature, and conductivity (Spearman's rank correlations) (Table S2). While DO was not correlated with borehole depth in this study, borehole depth was positively correlated with ammoniacal N, pH, and total dissolved solids (TDS) and negatively correlated with nitrate N and ORP, reflecting an increase in reducing conditions with well depth. DOC is generally associated with lower oxygen availability in groundwater (34); however, concentrations were near or at those consistent with oligotrophic groundwater conditions and did not exhibit a significant relationship with any of the parameters, possibly reflecting rapid utilization in groundwater (Table S1) (9).
The overall proportional abundance of class "Ca. Brocadiae" 16S rRNA gene sequences was significantly and negatively correlated with DO concentrations (r = 20.28) (Table S2), a trend also demonstrated in batch culture experiments (35) and observed in carbonate-rock aquifers (8). Significant positive correlations were instead found with DOC (r = 0.29), which was uncorrelated with DO, and also conductivity (r = 0.27), phosphate (r = 0.31), DRP (r = 0.31), and NO 2 2 (r = 0.31, although only 14% of samples had concentrations above the detection limit for NO 2 2 ). Of these, as indicated above, conductivity and phosphate/DRP concentrations negatively correlated with DO overall (P , 0.05), and NO 2 2 was also weakly negatively correlated (P = 0.05). Correlations between "Ca. Brocadiae" and other substrates, including those critical for ammonia oxidation (NH 4 1 and NO 3 2 ), were not significant, indicating that parameters such as DO, nitrite, and DOC are more important factors for determining the presence of anammox bacteria. In addition to these parameters, multivariate regression tree analysis showed that redox potential (ORP threshold, 175.9 mV) was the main factor discriminating "Ca. Brocadiae" communities ( Fig. S2).
As indicated by comparisons of planktonic and attached "Ca. Brocadiae" fractions ( Fig. 1c), analysis of anammox community diversity by oxygenation regime (anoxic, suboxic, dysoxic, and oxic) across all samples revealed that OTU richness and alpha diversity (inverse Simpson) were significantly higher at the suboxic and dysoxic sites than oxic sites (Fig. 2a), suggesting that groundwater with low oxygen supports more anammox species. The spatial distribution of different "Ca. Brocadiae" taxa was significantly associated with changes in DO and borehole depth (distance-based redundancy analysis [dbRDA]) (Fig. 2b), along with aquifer location and lithology, phosphate, and pH, which collectively explained 16% (R 2 adjusted) of the variation (permutation test; P , 0.05, permutations = 999). When only DO was considered, dramatic differences in composition were evident (Fig. 1b). Notably, the relative abundance of the genus "Ca. Brocadia" fraction was positively correlated with DO concentrations (r = 0.46) (Fig. 1b). Three "Ca. Brocadia" and "Ca. Brocadiaceae" OTUs made up 30.4 to 89.0% of the anammox-associated community in groundwater with .4 mg/L DO ( Fig. S3-S4). Oxygen tolerance in anammox bacteria has been demonstrated previously (36,37), particularly for "Ca. Brocadia." For example, "Ca. Brocadia sinica" exhibited tolerance in the presence of .1 mg/L DO in vitro (38), and anammox bacteria closely related to "Ca. Brocadia fulgida" have been reported in groundwater with up to 3 mg/L oxygen (8) (.3 mg/L DO is generally classified as oxic [33]). A recent study analyzing the metabolism of "Ca. Brocadia" species identified genes encoding superoxide dismutase and cytochrome c peroxidase genes that are involved in oxygen detoxification (39), as might be expected for organisms living at reduction-oxidation interfaces. "Ca. Brocadiae" overall were detected in groundwater with much higher DO concentrations in this study (up to 10.6 mg/L), and results indicate that "Ca. Brocadia" taxa are better adapted to higher-oxygen conditions. However, as aquifers are spatially heterogenous and anammox bacteria require access to reduced forms of nitrogen, it is likely that "Ca. Brocadiae" inhabit relatively lower-oxygen niches within the sampled aquifers than is reflected by bulk groundwater measurements. As such, the upper limits of DO tolerance for anammox bacteria require further validation. In comparison, a further 25.9% of "Ca. Brocadiae" comprised Planctomycetes RIFCSPHIGHO2_02_FULL_50_42 (mostly OTU 74). The overall relative abundance was significantly and negatively correlated with DO (r = 20.38) and ORP (r = 20.45), indicating a preference for reducing conditions (P , 0.05) (Table S2). It was originally recovered from an aquifer adjacent to the Colorado River (Rifle, CO, USA) (40), with conditions similar to site Wel13 here, where the genus was most abundant, including identical dysoxic conditions (0.78 mg/L DO) with low NO 2 2 and NO 3 2 .
(iv) DOC did not negatively impact "Ca. Brocadiae" community abundance. Anammox generally occurs at low organic carbon levels (41). As reported above, DOC and "Ca. Brocadiae" community abundances were instead positively correlated. The DOC composition was not analyzed in this study, which makes further comparisons between sites difficult. However, given the importance of redox conditions (and hence DO) for the anammox process, the observed association may be best explained by the negative relationship typically shared by DO and DOC in groundwater, where organic carbon availability regulates aerobic metabolism (34) and is associated with anoxia (although DO and DOC were not significantly correlated in this study [r = 20.07, P = 0.66]) (Table S2). Nevertheless, some anammox bacteria (e.g., "Candidatus Anammoxoglobus propionicus") have a higher affinity for NO 2 2 in the presence of organic acids, such as propionate, or can assimilate formate via the Wood-Ljungdahl pathway (as recently shown for "Candidatus Kuenenia stuttgartiensis") (42), implying that groundwater DOC could support the growth of some anammox bacteria.
(v) Hydrazine synthase capacity, but not activity, was correlated with DO. To confirm anammox potential and activity in groundwater, we quantified genes and transcripts encoding hydrazine synthase subunit hzsB. Together, genes hzsA, hzsB, and hzsC encode hydrazine synthase, which is a key enzyme in anammox metabolism, converting nitric oxide and ammonium into hydrazine (6). The relative abundance of "Ca. Brocadiae" 16S rRNA gene sequences correlated strongly with hzsB gene copies (r = 0.62, P , 0.05) (Fig. 3a), consistent with the expectation that "Ca. Brocadiae" undertake anammox (43). The concentration of hzsB genes ranged from 3.64 Â 10 1 to 1.23 Â 10 7 genes/L of groundwater (excluding 11.1% of samples below detection) with an average of 4.9 Â 10 5 6 2.1 Â 10 6 (SD) genes/L. This is consistent with concentrations of ;1.5 Â 10 4 to ;1.0 Â 10 8 genes/L found elsewhere in the subsurface, including in soil horizons interacting with the water table (8,44). The average number of hzsB transcript copies was significantly lower than that of hzsB gene copies (Wilcoxon test, P , 0.01; 4.5 Â 10 5 copies/L fewer, ;13.2 times lower) (Fig. 3b), suggesting a latent capacity of 92.4%. Overall, hzsB gene copy numbers were negatively correlated with DO (r = 20.42, P , 0.05), reflecting the trend seen with "Ca. Brocadiae" 16S rRNA gene relative abundance and DO, and as expected for an anaerobic process. While DO is expected to decrease with borehole depth, favoring anaerobic processes such as anammox, gene copy numbers were negatively correlated with borehole depth (r = 20.34). This may be explained by the weak relationship observed between DO and borehole depth in this study: DO tended to decrease with depth, but 22% of wells #10 m deep had DO concentrations of ,3 mg/L (Table S1). No significant relationship was evident between hzsB gene transcripts or anammox activity and oxygen availability (r = 20.16, P = 0.65) or bore depth (r = 0.14, P = 0.70).
Co-occurrence of aerobic and anaerobic ammonia oxidizers. The nitrificationanammox process allows for elevated N removal efficiencies from NH 4 1 in naturally ammonium-rich or contaminated water (45), by ensuring that a portion of ammonium is converted to nitrite. Potential combined nitrification-anammox was determined by quantification of archaeal and bacterial ammonia monooxygenase genes (amoA). Results show that hzsB and archaeal and bacterial amoA genes (r = 0.5) and transcripts (r = 0.71) were strongly

FIG 3
Correlations between 16S rRNA gene amplicon data, aerobic ammonia oxidizers (archaea and bacteria), and abundance of hydrazine synthase genes and transcripts from ddPCR data. (a) Correlation of hzsB gene copies per liter of groundwater (log 10 ) and number of "Ca.
Brocadiae" (log 10 ) gene amplicon sequences from the rarefied OTU table with samples containing the class "Ca. Brocadiae," showing strong agreement between quantitative data for the hzsB anammox functional marker gene and "Ca. Brocadiae" relative abundances. (b) Box plot representing log 10 abundance of hzsB genes for genes and transcripts. Black open circles represent means and horizontal lines represent the medians of gene copies per liter. Purple gradients for panels a and b correspond to DO concentration. (c) Correlation of amoA genes (log 10 ) and hzsB genes (log 10 ) copies per liter of groundwater using Spearman's correlation for genetic potential (DNA). (d) Correlation of amoA genes (log 10 ) and hzsB genes (log 10 ) copies per liter of groundwater using Spearman's correlation in transcripts (RNA). The gray dashed line shows the theoretical 1:1 ratio.
correlated, suggesting co-occurrence of aerobic and anaerobic ammonium oxidation (P , 0.05) ( Fig. 3c and d). Co-occurrence has been observed previously within laboratorybased models (46) and an aquifer system, where anammox was the dominant process (8).
Together, these results suggest that anammox in aquifers relies, at least partially, on nitrite derived from aerobic ammonia oxidation. Results also indicate that anammox is more frequently the dominant process, in terms of gene expression, when a wide range of aquifers is considered (Table S1). Although aerobic ammonia monooxygenase gene copies were more abundant in 83% of samples than hydrazine synthase genes (detected across 10 aquifers), hydrazine synthase transcripts were more abundant than ammonia monooxygenase transcripts in 69% of samples (detected across nine aquifers and higher in eight).
Phylogenomic diversity of anammox bacteria in aquifers. Genomes were reconstructed from 16 samples (gwj01 to gwj16) across eight wells (groundwater with or without biomass enrichment using sonication to detach biofilms from the surrounding aquifer materials) spanning a range of oxygen (0.37-7.5 mg/L) and nitrate (0.47-12.6 g/m 3 ) concentrations (Table S1). Of 541 recovered metagenome-assembled genomes (MAGs) that were .50% complete with #5% contamination, eight were identified as "Ca. Brocadiae" (63.4 to 95.6% complete) (Table S3). Of these MAGs that contained partial or (near) full-length (445 to 1,496 bp) 16S rRNA gene sequences (nzgw513, nzgw516, and nzgw517), those from nzgw516 and nzgw517 aligned fully with OTU729 (272 bp long), which was present in 21/80 samples (total of 92 sequences). Both the nzgw516 and nzgw517 and OTU729 sequences had best hits to Planctomycetes spp., including Planctomycetes strain Pla86 and "Ca. Kuenenia stuttgartiensis," when searched against the NCBI NR database (BLASTN; 79.54% to 82.85% identity), consistent with phylogenomic analyses showing that these genomes comprise a novel clade (clade II) (Fig. 4a). OTU729 was classified only to the domain Bacteria using USEARCH with the SILVA SSU Ref NR99 database and RDP Classifier, suggesting that these novel "Ca. Brocadiae" could be underestimated in the environment when a 16S rRNA gene-based approach is used.
A pairwise comparison with 16 other "Ca. Brocadiae" genomes from various environments yielded average nucleotide identity (ANI) values below the proposed species cutoff of 95% (47) (Fig. S5), implying that the newly recovered MAGs represent novel species. The phylogenomic tree (Fig. 4) revealed that six NZ aquifer-derived MAGs form three (sub)clades with recently recovered genomes from the aquifer in Rifle, CO (40), suggesting that these are groundwater-adapted species. One (nzgw512) clustered with the well-characterized anammox bacterium "Candidatus Jettenia caeni" (48). An eighth MAG, nzgw510, was phylogenetically distinct from other "Ca. Brocadiae." MAGs clustering with the Rifle MAG MHYJ01 (from aquifer sediment [40]) constitute clade II, which is distinct from well-known anammox bacteria (clade I) (Fig. 4a). They contained a much higher GC content (.67%) than others from this study, excluding the outlier nzgw510 (average 6 SD, 42.5% 6 3%) (Table S3). A higher GC content may indicate a higher evolutionary rate for these organisms (49). All lower-GC "Ca. Brocadiae" MAGs from this study, including two aquifer subclades, coclustered with well-known anammox Candidatus genera from enrichment studies (clade I). The estimated genome sizes were significantly smaller than the high GC group (clade I, 3.22 6 0.92 Mbp; clade II, 5.16 6 1.24 Mbp), further suggesting a divergence in evolutionary strategy employed by the two groups.
As for 16S rRNA gene analyses (Fig. 1b), "Ca. Brocadiae" MAG relative abundances varied across sites, corresponding to distinct geochemical conditions ( Fig. 4b; Table S3) but not strictly phylogenetic relatedness (Fig. 4a). Three of four MAGs associated with well-characterized clade I anammox bacteria (nzgw511 to nzgw513) were most abundant at sites with low DO (,1.1 mg/L) and high DOC (3 to 26 g/m 3 ). MAG nzgw514 is one of two closely related to "Ca. Brocadia." Consistent with 16S rRNA gene data (Fig. 1b), it demonstrated a preference for oxic groundwater (Fig. 4b). Together with the high-GC MAGs (nzgw510 and clade II), it was more abundant at sites with high DO (5.6 to 7.5 mg/L) and low DOC (0 to 0.8 g/m 3 ). Factors identified that influence ecological niche differentiation among anammox species include microbial growth kinetics, organic matter, oxygen tolerance, aggregation ability, and interspecific competition (5). Collectively, results here suggest that DO availability may be a major driver of anammox bacterial diversification in groundwater.
Reinforcing genomic evidence, metatranscriptomics data for six samples (gwj09, gwj11, and gwj13 to gwj16) across two sites revealed the transcriptional activity of hzoA by MAGs nzgw511 to nzgw514. This activity was 356-fold higher at the dysoxic site (wells E1 and N3) than the oxic site (Fig. 5b) and likely contributed in part to excess N 2 measured at those sites (Fig. 6). Despite the compositional bias in anammox bacteria associated with differences in DO ( Fig. 1b and 4), results confirm that in groundwater with low oxygen concentrations, anammox bacteria are more active (8). Contemporaneous measurement of excess N 2 indicated active N 2 generation in dysoxic groundwater (from wells E1 and N3) due to denitrification or anammox. In contrast, oxic groundwater was either at the bounds of uncertainty (wells BW8 and RF3) or devoid of measurable excess N 2 (wells SR1-2, BW19, and RF2) (Fig. 6), Arrows indicate competing processes that can alter gas concentrations, and gray dashed lines indicate excess air in groundwater relative to atmosphere (with upper and lower lines representing addition of unfractionated excess air relative to equilibrium concentrations at 10 and 15°C). The black horizontal arrow depicts additional excess N 2 inferred to be from biological processes (denitrification or anammox). Reconstructed N 2 data (in equilibrium with inert atmospheric gases), based on groundwater recharge temperatures and excess air concentrations derived from dissolved Ne and Ar data (shown as black circles). Recharge temperature is the temperature of recharging water at the time it enters the groundwater system. Excess air is dissolved air in excess of the equilibrium soluble amount at the given recharge temperature (thought to originate from processes such as bubble entrapment occurring during recharge and subsequent dissolution under increased hydrostatic pressure). The difference between these and the measured N 2 data (numbered blue circles and their shift along the x axis relative to the corresponding numbered black circles) indicates the amount of N 2 in excess, formed via denitrification and/or anammox at each site. Error bars show the combined statistical standard uncertainty from all processes and calculations contributing to the measurement uncertainty, expressed as 1 standard deviation. Groundwater from wells SR1-2, BW8, BW19, and RF2-3 is characterized as oxic, while groundwater from E1 and N3 is dysoxic-suboxic (Table S1). which, accordingly, included groundwater with relatively little observed hydrazine synthase/ dehydrogenase activity (wells SR1 and SR2). Additional quantification via ddPCR of hydrazine synthase transcripts from each site, sampled at an earlier date, further demonstrated that anammox was active in both dysoxic and oxic groundwater (Fig. 7) associated with all "Ca. Brocadiae" MAGs, which were distributed across eight wells over four sites ( Fig. 4b; Table S3). However, quantitative results were consistent with relative metatranscriptomics data, indicating that transcriptional activity was 100-fold higher in dysoxic than oxic groundwater in the four samples analyzed.
It is suggested that anammox bacteria adapt a dissimilatory nitrate reduction to ammonia (DNRA)-like metabolism in the absence of NH 4 1 , using nitrate reductase, followed by NrfAH cytochrome c nitrite reductase (28). The respective products, NO 2 2 and NH 4 1 , can then drive anammox. We found nrfAH genes in phylogenetically diverse genomes (nzgw514, nzgw516, and nzgw517). Further, genomes nzgw512 to nzgw514 contain a tetraheme cytochrome c nrfH gene, with CxxCH repeated sequence motifs similar to "Ca. Brocadia" (62) and "Ca. Jettenia" (60). Both nrfA and nrfH genes were observed to be expressed by these taxa in dysoxic groundwater (Fig. 5b). Together, these genes and transcripts indicate the active reduction of NO 2 2 to NH 4 1 in a NAD(P) H-dependent, six-electron transfer reaction, via a mechanism resembling DNRA (50), potentially fueling active anammox in the same sites. NAD(P)H nitrite reductase genes (nirD and/or nirB), which also reduce NO 2 2 to NH 4 1 , were found in clade II genomes, nzgw515 to nzgw517. The large subunit encodes NirB, an iron-sulfur protein. All protein-coding sequences had the highest identity (79 to 84%) to NirB from a Planctomycetes genome (GCA_016200125.1), recently recovered from a pristine aquifer in the United States (64). This had a geochemical profile similar to that of the site where nzgw516 and nzgw517 were most abundant (low NH 4 1 , ,0.05 mg/L, and NO 3 2, 0.23 mg/L). The nirB gene was recently also identified in a marine-derived anammox bacterium, related to "Ca. Scalindua" (65). We also identified a small subunit, NirD, in nzgw516 and nzgw517. The sequence contains a Rieske nonheme iron oxygenase family domain with closest similarity to Verrucomicrobia species (63.5%). NH 4 1 is often limiting in groundwater (20). In this study, NH 4 1 was undetectable in groundwater from which MAGs were recovered and was detected at only 31% of sites overall. The presence of these Nrf and Nir proteins could therefore mediate an important alternative process providing NH 4 1 for anammox (50).
(iii) Mechanisms of oxygen tolerance in groundwater anammox bacteria. All "Ca. Brocadiae" genomes encode proteins for oxygen protection, such as BatD and BatA, which are involved in oxygen tolerance in Bacteroides species (66), alkyl hydroperoxide reductases, and thioredoxin (67). Thioredoxin and peroxiredoxin genes were expressed in MAG nzgw510 at the oxic site 0.06 to 0.15 transcripts per kilobase per million reads mapped ([TPM]), and BatD, BatA, catalase, peroxiredoxin and thioredoxin were expressed by MAGs nzgw511 to nzgw514 at the dysoxic site (0.112 to 2.81 TPM). Genes encoding respiratorychain cytochrome c oxidase, cbb 3 type (ccoN, ccoO, and ccoP), were all or partially present in four MAGs (nzgw511 to nzgw514). The cbb 3 -type oxidase is a multichain transmembrane protein located near the anammoxosome membrane (68). This enzyme is believed to have evolved to perform a specialized function in microaerobic energy metabolism; however, it could also merely act in oxygen detoxification (28). Additionally, three clade II genomes encode a caa 3 -type cytochrome oxidase, containing unique CoxAB subunits. In Acidithiobacillus ferrooxidans, coxAB genes encode the terminal part of the ferrous iron "downhill" pathway, which shuttles electrons from extracellular iron to oxygen (69), and thus represents another potential detoxification strategy.
(iv) Alternative metabolic pathways in groundwater anammox bacteria. Anammox bacteria grow chemolithoautotrophically by performing CO 2 fixation via the Wood-Ljungdahl pathway (70), and genes encoding the pathway were found in all of our groundwater genomes, including the complete pathway in nzgw511 and nzgw512 ( Fig. 5; Table S4). The products of this reaction, acetyl coenzyme A (acetyl-CoA) and pyruvate, enter the oxidative tricarboxylic acid (TCA) cycle and gluconeogenesis ( Fig. 5a; Table S4). However, in "Ca. Kuenenia stuttgartiensis," as in four of our genomes, the oxidative branch of the TCA cycle is likely mediated by a Re-citrate synthase rather than a citrate synthase. MAGs nzgw511 to nzgw514 encode a Re-citrate synthase that was highly similar to that identified in Clostridium kluyveri (55 to 58%) (71). It operates incompletely to synthesize alpha-ketoglutarate, similar to other anaerobic bacteria (72).
In addition, several ABC transport systems were present within the genomes, revealing variations in substrate importation such as phosphate, cobalt, nickel, iron(III), zinc, sulfate, molybdate, lipoproteins, ribose, rhamnose, polysaccharides, and oligopeptides (Table S4). Oligopeptide transport systems were present in three genomes (nzgw515 to nzgw517), as seen in "Ca. Scalindua profunda," and suggest these genomes are capable of oxidizing decaying organic matter (53), in addition to carbon fixation via the Wood-Ljungdahl pathway (Fig. 5a). Additionally, iron(III) transporters were present in three genomes, which may enable the coupling of formate oxidation to iron(III) reduction in the absence of NO 2 2 or NO, as previously identified in other species (73). Hydrogenase genes were also present in three of the genomes. The first encodes a [Ni-Fe] group 3b hydrogenase, found in nzgw511, which is proposed to harbor (sulf)hydrogenase activity (74). This genome also encoded a sulfate adenylyltransferase, which may play a role in sulfide oxidation (75). The second is a group 4 hydrogenase present in nzgw513. These hydrogenases carry out coupled oxidation of NADPH to the evolution of H 2 (76), indicating hydrogen turnover and/or an alternative energy source. Our results therefore suggest these groundwater anammox bacteria are metabolically versatile, containing various hydrogenases and ABC transporters for organic compounds and genes encoding DNRA, which would allow for growth under substrate-limiting conditions (such as low ammonium concentrations).
Implications for groundwater. Findings presented here indicate that natural and agriculturally derived nitrate and ammonium (77) may be removed from aquifers with a range of bulk oxygen concentrations by endogenous anammox bacteria. In this study, gene transcripts indicative of anammox were positively correlated with those for aerobic ammonia oxidation, which could fuel the anammox process by providing nitrite, across anoxic-to-oxic conditions. These two processes are considered to have opposing oxygen requirements; however, some oxygen in groundwater may facilitate loss of fixed N by stimulating aerobic ammonia oxidation. Additionally, although anammox bacteria have been shown to outcompete denitrifiers in low-carbon groundwaters (8), we found no evidence that increased DOC availability negatively impacted anammox bacteria (overall) or anammox-related activity. While DOC is required by heterotrophic denitrifying bacteria (1), its influence on denitrifiers (autotrophic or heterotrophic) and anammox bacteria is likely complex. For instance, DOC plays a role in the regeneration of ammonium (78), essential for anammox, and oxygen depletion stimulated by DOC creates conditions favorable for both competing pathways leading to N 2 . Moreover, genomic data indicate that anammox bacteria in groundwater may be capable of supplementing their metabolism through the acquisition of exogenous organic carbon. Further studies are needed to validate the oxygen tolerances and metabolic versatility of anammox bacteria in the terrestrial subsurface.
Conclusions. This study shows that anammox bacteria are prevalent and active across wide-ranging aquifer chemistries and lithologies, including in oxic groundwater, although this is predicted to be unfavorable. While we found significantly more "Ca. Brocadiae" diversity in anoxic groundwater, some taxa were positively associated with DO concentrations. Of the eight novel "Ca. Brocadiae" genomes reconstructed, three belong to undercharacterized clades previously recovered only from aquifers. Most possess genes involved in signature pathways of anammox, including novel hydrazine synthase genes. A co-occurrence of anaerobic and aerobic ammonia oxidizers at many sites suggests metabolic handoffs (such as nitrite) between these processes. Furthermore, genomic characterization of "Ca. Brocadiae" identified a range of potential aero-tolerance mechanisms, explaining our finding of anammox in oxic groundwaters. Results indicate that niche differentiation occurs among anammox bacteria based on oxygen concentrations and that anammox is a common mechanism for nitrogen removal in aquifers.

MATERIALS AND METHODS
Study sites and sample collection. Groundwater was collected from 59 wells spanning 10 aquifers in the Waikato, Wellington, and Canterbury regions of New Zealand (Table S1), from several aquifer lithologies: alluvial sandy gravel (71 samples), sand-silt (1 sample), basalt (1 sample), shell bed (1 sample), peat (1 sample) and ignimbrite (6 samples). Wells were purged (;3 to 5 bore volumes); then, 3 to 90 L of groundwater or 0.5 to 15 L of biomass-enriched groundwater was filtered on site. Biomass-enriched groundwater was collected from 8 sites in Canterbury (Table S1) directly following standard groundwater collection, and low-frequency sonication (2.43 kW) for 2 min to detach biofilms and aquifer particles (also covered in biofilms) (79). Biomass was captured onto mixed cellulose ester membrane filters (1.2-mm-pore-size prefilter over a 0.22-mm filter), using a 14-mm stainless steel filter holder (Merck Millipore Ltd., Cork, Ireland). Both filters were immediately submerged in RNA Later (Thermo Fisher Scientific, Waltham, MA, USA), transported on dry ice, and stored at 280°C.
Chemical analysis of water samples. Dissolved oxygen (DO), water temperature, pH, oxidationreduction potential (ORP), and specific conductance (SPC or conductivity) were measured on-site using a flowthrough cell and field probes (YSI EXO sonde 2, YSI PRO1, and YSI ProDSS; YSI, Yellow Springs, OH, USA). Samples were categorized by DO concentration as follows: anoxic (0 mg/L), suboxic (,0.3 mg/L), dysoxic (0.3 to 3 mg/L), and oxic (.3 mg/L) (33). Unfiltered groundwater was analyzed for P, N, C, S, Mn, and alkalinity at Hill Laboratories (Hamilton, New Zealand). Total phosphorus and phosphate were determined according to American Public Health Association method 4500-P (APHA 4500-P), parts B and E, (modified to include an acidic ammonium persulfate to convert organophosphates and polyphosphates to orthophosphate) (80). Dissolved reactive phosphorus (DRP) was determined according to APHA 4500-P, part G (sample was reacted with ammonium molybdate and ascorbic acid to form molybdenum blue and then detected at 880 nm). Total ammoniacal N was determined according to APHA 4500-NH3, part H (using a phenol/hypochlorite reaction forming a complex that was detected at 630 nm), and calculated as NH 4 -N = NH 4 1 -N 1 NH 3 -N. Nitrite N and nitrate N 1 nitrite N were determined according to APHA 4500-NO 3 , part I (NO x N is measured via automated cadmium reduction and Griess reaction; NO 2 -N is calculated by automated Griess reaction [sulfanilamide] detected at 540 nm). Nitrate N was calculated as (nitrate N 1 nitrite N) 2 nitrite N. Total organic carbon (TOC) and dissolved organic carbon (DOC) were measured according to APHA 5310, part C (analyzed using supercritical persulfate oxidation with phosphoric acid and sodium persulfate), where TOC is calculated as total carbon 2 total inorganic carbon. For DOC, groundwater was filtered first using a 0.45-mm polypropylene filter. Total suspended solids were measured by first evaporating groundwater samples in an oven at 105°C until dry. Dried solids were then weighed and normalized to the total water volume analyzed. Sulfate was measured according to APHA 4110 B. Alkalinity and dissolved manganese were analyzed according to APHA 2320 B and APHA 3125 B, respectively.
To assess denitrification/anammox, excess N 2 gas was quantified using methods by Martindale et al. (81). Briefly, dissolved Ar, Ne, and N 2 were measured by a standard curve via a system comprising two detectors, a pulsed-discharge helium ionization detector (Valco Instruments D-4-I-SH14-R), a thermal conductivity detector (Shimadzu TCD-2014), and ultrahigh-purity helium (He) gas. The measurement of a sample uses the principles of headspace analysis and Boyle's law (81). Excess N 2 , attributable to denitrification/anammox, was determined by comparison to inert Ar and Ne gas concentrations. The uncertainty was reported for each measurement of the original sample concentration as the standard measurement error (combined standard uncertainty).
Nucleic acid extraction and sequencing. (i) Nucleic acid extraction. RNA and DNA were extracted using RNeasy PowerSoil Total RNA and DNeasy PowerSoil Pro kits (Qiagen, Valencia, CA, USA), respectively, with nuclease-free glycogen added to aid RNA precipitation (0.1-mg/mL final concentration; Roche Diagnostics, Basel, Switzerland). DNA extractions used 0.14 to 0.89 g of crushed filter (1 to 47 extractions per sample). For RNA, 2.12 to 3.90 g was used per extraction (1 extraction per sample). RNA was DNase treated using the Turbo DNA-free kit ("rigorous" protocol; Invitrogen, Carlsbad, CA, USA). DNA removal was verified via 16S rRNA gene amplification (as described below, but over 55 cycles) and gel electrophoresis. The products of replicate DNA extractions were pooled and concentrated using sodium acetate (0.3 M final concentration) and ethanol (2Â volume) precipitation with 0.1 mg/mL glycogen (Roche Diagnostics) via ethanol precipitation. RNA extractions were concentrated using the Zymo RNA Clean and Concentrator-5 kit (Zymo Research, Irvine, CA, USA).
High-molecular-weight DNA for whole-genome shotgun (WGS) sequencing was verified via 1% agarose gel electrophoresis. Nucleic acids were quantified with a Qubit 3.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) using double-strand DNA (dsDNA) HS and RNA HS assay kits and quality checked using a Nanophotometer (Implen, Munich, Germany). RNA was checked using an Agilent Bioanalyzer with RNA 6000 Nano and Pico chips (Integrated Sciences, NSW, Australia). Samples with an RNA integrity number of 6 or more or a fragment distribution value (DV200) of .30% (percentage of RNA fragments above 200 nucleotides) were used to quantify transcripts and for transcriptomics. Between 24 pg and 1.1 mg of total RNA was converted to cDNA using Superscript III Supermix (Invitrogen, Carlsbad, CA, USA).
(ii) Metagenome, metatranscriptome, and amplicon sequencing. DNA libraries for 15 samples (gwj01 to gwj16) were prepared using the TruSeq Nano DNA kit (Illumina, San Diego, CA, USA), with a targeted insert size of 550 bp, at the Otago Genomics Facility (University of Otago, New Zealand), except low-yield sample gwj02, which was prepared with the ThruPLEX DNA-seq kit (TaKaRa Bio USA, Inc., Mountain View, CA, USA). Then, 2 Â 250 bp sequencing was performed using the Illumina HiSeq 2500 V4 platform. RNA libraries were prepared using the Ovation SoLo RNA-Seq system (NuGEN, Redwood City, CA, USA) using custom probes for rRNA depletion at the Otago Genomics Facility. Custom rRNA probes were designed by the manufacturer using small and large ribosomal subunit sequences reconstructed from the 16 metagenomes with EMIRGE (82) over 40 iterations with clustering at 97% identity and using the SILVA 132 database (83). Ribosomal sequences generated were used as target sequences in the design of custom AnyDeplete probes using NuGEN's proprietary algorithm. rRNA and genomic DNA were removed as a part of the library preparation kit, using the custom probes and DNase treatment, respectively. Paired-end 2 Â 125 bp reads were generated from RNA libraries using the Illumina HiSeq 2500 V4 platform. PCR amplification of 16S rRNA genes used modified Earth Microbiome Project primers EMP-16S-5159F and EMP-16S-8069R primers (84,85) with Illumina Nextera adapters and MyTaq HS Red mix (Bioline, London, UK). PCR cycling conditions were as follows: 95°C for 5 min, 35 cycles consisting of 45 s at 95°C, 60 s at 50°C, and 90 s for 72°C, followed by 72°C for 10 min. Amplicons were purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA). Barcoded libraries, prepared according to Illumina's 16S metagenomic sequencing library preparation manual, were loaded with 10% PhiX for 2 Â 250 bp sequencing via Illumina MiSeq with V2 chemistry (Auckland Genomics, University of Auckland, Auckland, New Zealand).
Phylogenetic and protein sequence trees. Bacterial core gene alignments (75 to 114 genes per genome) generated from GTDB-Tk were used to construct a maximum-likelihood phylogenomic tree in IQ-TREE (v1.6.9) (115) using ModelFinder (116) best-fit model LG1F1R5, and annotated with iTOL (117). EMIRGE, Metaxa2, and amplicon sequences were aligned using MUSCLE (118) with default parameters and trimmed to 295 bp using Geneious 11.1.2. A maximum-likelihood tree was constructed as with IQ-TREE using ModelFinder best-fit model TIM31F1I1G4. Protein sequences of hydrazine synthase subunits A, B, and C were aligned with MUSCLE and trimmed to remove columns with .50% gaps using trimAl (119). A maximum-likelihood tree was constructed as with IQ-TREE using ModelFinder best-fit model WAG1G4.
Metatranscriptome processing. Adapters were removed from transcriptomic reads and then quality trimmed and checked as described above for metagenomic reads. Residual rRNA sequences were removed using SortMeRNA v2.1 (120). Additional checks were carried out to ensure that filtered reads were still paired, and read pairs were ordered using a repair script from BBMap v38.81 (121). Filtered transcriptomic reads were mapped to contigs from the set of dereplicated MAGs using Bowtie2 (103) (v2.3.5, -end-to-end -very_sensitive). Read counts were determined using featureCounts (122) (v1.6.3, -F SAF). Singleton reads per gene were removed, and the remaining read counts were normalized to transcripts per kilobase per million reads mapped (TPM) using the following equation: (number of reads mapped to gene) Â (1,000/gene length) Â (1,000,000/library size) (123).

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.