Microbial controls on net production of nitrous oxide in a denitrifying woodchip bioreactor

Denitrifying woodchip bioreactors (DWBs) are potential low-cost technologies for the removal of nitrate (NO3 - ) in water through denitrification. However, if environmental conditions do not support microbial communities performing complete denitrification, other N transformation processes will occur resulting in the export of nitrite (NO2 - ), nitrous oxide (N2 O), or ammonium (NH4 + ). In order to identify the factors controlling the relative accumulation of NO2 - , N2 O, and/or NH4 + in DWBs, porewater samples were collected over two operational years from a DWB designed for removing NO3 - from mine water. Woodchip samples were collected at the end of the operational period. Changes in the abundances of functional genes involved in denitrification, N2 O reduction, and dissimilatory nitrate reduction to ammonium were correlated with pore water chemistry and temperature. Temporal changes in the abundance of the denitrification gene nirS were significantly correlated with increases in porewater N2 O concentrations, and indicated the preferential selection of incomplete denitrifying pathways ending with N2 O. Temperature and the TOC/NO3 - ratio were strongly correlated with NH4 + concentrations and inversely correlated with the ratio between denitrification genes and the genes indicative of ammonification (∑nir/nrfA), suggesting an environmental control on NO3 - transformations. Overall, our results for a DWB operated at hydraulic residence times of 1.0 - 2.6 days demonstrate the temporal development in the microbial community and indicate an increased potential for N2 O emissions with time from the DWB. This article is protected by copyright. All rights reserved.

In DWBs, C availability is controlled by woodchip degradation, and the relative availability of different C substrates may change throughout DWB operations (Grießmeier & Gescher, 2018;Grießmeier et al., 2017;Nordström & Herbert, 2018). It can therefore be expected that differences in relative abundances of the functional groups involved in the different N transforming processes develop over time during DWB operation. The proportion between these functional groups ultimately controls the export of N species from DWBs, but little is known about the temporal development of the Ntransforming community and the associated temporal changes in the production of NO 2 − , N 2 O, and NH 4 + in DWBs.
In this study, temporal and spatial changes of the abundances of functional groups performing denitrification and DNRA in the porewater were studied with the objectives to relate these patterns with changes in the concentrations of NO 2 − , N 2 O, and NH 4 + in the porewater and overall reactor performance in a previously described DWB (Nordström & Herbert, 2018). We hypothesized that temporal and spatial changes in the genetic potential for denitrification and DNRA, determined as abundances of functional genes in den-

Core Ideas
• A high degree of spatial and temporal variability for functional gene abundances was noted in porewater. • Temperature dependence was exhibited especially for the gene nrfA. • N 2 O production correlated with nirS abundance and Σnir/ΣnosZ ratio. • TOC/NO 3 − ratio positively correlated with nrfA abundance and NH 4 + concentration. • A truncated denitrification pathway was promoted with time in the bioreactor.
itrifying and DNRA bacteria, control concentrations of NO 2 − , N 2 O, and NH 4 + and that their abundances are a consequence of changes in C/NO 3 − ratios. Because the porewater community is more easily sampled for temporal studies than the woodchip-associated community, we primarily monitored the development of the N-reducing community in the porewater. However, for the last sampling episode, we compared the spatial distribution of the abundance of denitrifying and DNRA bacteria in the woodchip matrix with their spatial distribution in porewater.

Study site and DWB system
The subsurface DWB described by Nordström and Herbert (2018) was constructed at the Kiruna iron ore mine, northern Sweden (67˚51′ N, 20˚13 ′E), with the purpose to reduce NO 3 − concentrations in mine and process water originating from the use of ammonium nitrate-based explosives. For this study, samples from the inlet, outlet, and five porewater sampling points along the bottom centerline of the DWB were used (Supplemental Figure S1). The DWB was filled with decorticated pine woodchips. To increase the initial abundance of denitrifying microorganisms, digested sewage sludge mixed with water was added while the DWB was filled with woodchips. The center areas of the DWB were capped by glacial till (Supplemental Figure S1) with the intention to restrict oxygen diffusion into the DWB.
Mine drainage from the clarification pond at the mine site was pumped to the DWB, where it entered through a perforated drainage pipe near the surface of the DWB and extending across its width. The hydraulic residence time (HRT) of the DWB was adjusted several times during the operational periods by changing the pump discharge (Supplemental Table S1), with 2.6 d being the most common HRT (varying between 1 and 2.6 d). The choice of HRT was based on a previous laboratory-scale experiment where relatively long HRTs at low temperature provided nearly complete removal of NO 3 − (Nordström & Herbert, 2017 Porewater, inlet water, and outlet water were sampled once a month (referred to as "profile sampling") to analyze water chemistry (Nordström & Herbert, 2018). Briefly, a peristaltic pump was used for porewater sampling, and the first ∼2 L of water was discarded prior to sample collection to ensure a representative sample. Inlet and outlet water was grab-sampled. For analyses of microbial communities involved in N transformation processes, water sampled on Days 57,85,113,365,400,428,456, and 477 was filtered using Sterivex filter units, with 0.22-μm Millipore Express polyethersulfone membranes, attached to 60-ml syringes. Inlet water samples from Days 57 and 85 were discarded because of technical problems during sampling. Between 540 and 1,680 ml of water (average, 990 ml) was required to saturate the Sterivex filter units. The syringes were rinsed in sample water three times prior to filtration, and new polyvinyl chloride plastic tubing for the peristaltic pump was used for each sample. The Sterivex filter units were stored on ice for ∼6 h following collection and frozen at −20˚C until DNA extraction. Denitrifying woodchip bioreactor porewater temperature was obtained from thermistors attached to porewater sampling points at the base of the DWB.

Woodchips and sewage sludge media
Sewage sludge samples used as inoculum were collected in sterile 50-ml plastic tubes at the time of DWB construction and frozen at −20˚C until analysis. Woodchips were sampled following the termination of DWB operations (Day 490). The sampling focused on the deepest regions of the DWB because a previous DWB study indicated a significantly greater abundance of 16S rRNA, nirS, nirK, and nosZI genes at the greatest depths in a DWB (Herbert, Winbjörk, Hellman, & Hallin, 2014). Water was drained from the DWB, and traverse trenches were excavated at 2.7, 11.2, 19.7, 28.2, and 33.9 m from the inlet (Supplemental Figure S1). Woodchip samples were collected along a center line at bottom depth (corresponding to the five porewater sampling points) and at 0.4 m above the bottom. Samples were also collected at depths of 11.2, 19.7, and 28.2 m at ±1.35 m from the center line and 0.4 m above bottom (i.e., the deepest possible depth due to the trapezoidal shape of the reactor). At the inlet and outlet, woodchip samples were collected from the surface of the bioreactor. All samples were collected in triplicate. The samples were stored in sterile 50-ml plastic tubes and placed on ice for ∼6 h following collection, frozen at −10˚C for 4 d, and stored at −20˚C until DNA extraction.

Chemical analyses of water samples
We used porewater chemistry data from Nordström and Herbert (2018)

DNA extraction and quantitative real-time polymerase chain reaction
DNA from the Sterivex filter units was extracted using the MoBio PowerWater Sterivex DNA kit following the manufacturer's instructions, including incubation at 90˚C (Mobio Laboratories Inc., 2018). For DNA extraction from the woodchips (4 g) and sewage sludge (0.2 g), the DNeasy PowerMax Soil Kit was used according to the manufacturer's instructions (Qiagen GmbH). Both the woodchips and digested sewage sludge were freeze-dried prior DNA extraction.
The 16S rRNA gene (Muyzer, Dewaal, & Uitterlinden, 1993) was used as a proxy for the abundance of the total bacterial community. Each 15-μl qPCR reaction contained 2-5 ng (water samples) or 0.1-0.3 ng (woodchip samples) of DNA, 0.25-2 μM of each primer, 15 μg bovine serum albumin, and 1x iQ SYBR Green Supermix (BioRad Laboratories). Two independent quantifications per gene were performed using the BioRad CFX Real-Time System (BioRad Laboratories). Potential PCR inhibition in the samples was tested by spiking each sample with the pGEM-T plasmid (Promega Co.) and amplifying it using plasmid-specific primers. Amplification was compared between samples and nonsample (water control) reactions, and no inhibition was present in the samples with the amounts of DNA extract used. Cycling protocols and primer concentrations are described in the supporting information (Supplemental Table S2).

Statistical analysis and parameter estimation
Differences in the abundances of all functional genes and the 16S rRNA genes between porewater (including outlet samples), inlet water, woodchips, and sewage sludge were tested using Dunn's test, which is appropriate for groups with unequal numbers of observations (Zar, 2010), in R package 'FSA' (Dinno, 2017), which performs a Kruskal-Wallis test (normality could not be assumed based on a Shapiro-Wilks test; Supplemental Table S3) followed by pairwise comparisons. Corrections for multiple comparisons were done by false discovery rate (Benjamini & Hochberg, 1995). For comparisons between two groups, Wilcoxon rank sum tests were used.
Nonmetric multidimensional scaling (NMDS; using the R package 'vegan' [Oksanen et al., 2018]) was used to illustrate the structural differences in the concatenated N reducing communities (i.e., based on all functional genes) between samples from the porewater (including outlet samples), inlet water, woodchips, and the inoculum (digested sewage sludge). The abundances of the studied marker genes from each individual sample in the DWB were assumed to represent a specific N-reducing community. Gene abundances were square-root transformed and submitted to Wisconsin double standardization. Then, a community matrix with Bray-Curtis dissimilar-ities was created and the NMDS was run with a maximum of 250 iterations.
We compared the similarity between the N reducing community in the porewater samples with those in the three potential sources: the inlet water, the woodchip media, and the inoculant sewage sludge. Community matrices based on Bray-Curtis dissimilarities were generated for each sample source using the same procedure as described above and compared using a permuted (n = 999) analysis of similarity (anosim; R-package 'vegan') (Oksanen et al., 2018). To

Porewater chemistry
The NO

Gene abundances and N-reducing community structure
The abundances of the total bacterial community and genes coding for the different N reducing pathways differed among sample types, but the hdh gene coding for anammox was not detected in any of the samples (Figure 2a). An assessment of the difference between woodchip samples collected at different depths indicated that there was not a significant difference (p < .05, one-sided Wilcoxon test) in the abundances of the 16S rRNA genes and the functional genes between the two sampling depths, with the exception of nirS. For nirS, there were significantly lower abundances at the bottom compared with 0.4 m above (p < .01) (data not shown). Regardless, woodchip samples from these two sample depths are considered together for the rest of this study.
There were significant differences in the abundances of the 16S rRNA and functional genes between the porewater, inlet water, woodchips, and digested sewage sludge (Figure 2a; Supplemental Table S4). The sample source with the highest abundance of functional genes differed depending on the gene, although gene abundances in the inlet water were most often the lowest among the different sample sources (Figure 2a). For the different genes and sample sources, nirS was most abundant in the woodchips and sludge, whereas nrfA in sludge had the greatest abundance. Among the woodchip samples, the 16S rRNA and nirS abundance varied the most along the flowpath from inlet to outlet and the nrfA abundance did not vary (Figure 2b). nirK and nosZI displayed the same pattern with the highest abundances in the middle section of the reactor, whereas nirS and nosZII had contrasting patterns, with nirS lower at the inlet and nosZII lower at the outlet.
The NMDS analysis (Figure 2c) showed that nirS-type denitrifiers were indicative of the N-reducing community structure in the woodchip media. By contrast, nosZI and nosZII involved in N 2 O reduction were relatively more enriched in the porewater and inlet water samples in comparison to the woodchip media (Figure 2c), and nrfA (DNRA) characterized the N-reducing community in the digested sewage sludge.
Despite differences in absolute abundances, the distribution of gene abundances relative to each other in the porewater was comparable with that in the woodchip media (nirS > nrfA ≥ nirK > nosZII > nosZI) (Figure 2a). However, when comparing the N-cycling community in porewater at Day 477 and woodchip samples from Day 490, there was a significant difference between these two sample pools (R = .87; p = .001), demonstrating that these two matrices represent different sampling environments.
increasing symbol sizes imply later sampling dates. Gene names in black boxes denote species scores

Temporal changes in the N-reducing community and associated biogeochemical changes in the porewater
We observed differences in gene abundances between the operational years and along the flowpath of the DWB (Figure 3). In general, there was a significant increase in gene abundance between the corresponding periods in the first and second operational year with respect to nirS, nirK, nosZII, and nrfA, whereas the total bacterial community remained approximately the same and nosZI decreased significantly (Supplemental Figures S2 and S3). However, the opposite changes in nosZI and nosZII abundances more or less canceled out the total change in potential N 2 O reduction capacity between years. The structure of the N reducing community in the porewater also changed over time and differed between the two operational years, as shown with NMDS (Figure 4), and the change was significantly associated with changes in porewater concentrations of NO 3 − -N (R 2 = .34; p = .001), NO 2 − -N (R 2 = .30; p = .003), N 2 O-N (R 2 = .17; p = .022), NH 4 + -N (R 2 = .14; p = .040), TOC (R 2 = .28; p = .003), DWB temperature (R 2 = 0.46; p = .001; see also Figure 2), and TOC/NO 3 − ratio (R 2 = .3627; p = .001) but not with porewater pH (R 2 = .01; p = .901). The differences in the structure of the N reducing community between years were mainly split in relation to axis 1 (Figure 4), corresponding to changes in the porewater concentrations of NO 2 − -N and N 2 O-N. Although there was no significant difference in the TOC/NO 3 − ratio between the two operational years (Supplemental Figure S2), the NMDS demonstrates that the TOC/NO 3 − ratio is likely controlled by DWB temperature (Figure 4). Because the gene abundance data from the first operational year were limited to three sampling dates, it was not possible to make conclusive interpretations of temporal trends during the first year. The investigation of temporal changes over an operational year was hence restricted to data from the second year. During this period, gene copies representing the total bacterial community (Figure 3a) decreased in abundance by approximately 50%. Maximum gene abundances were observed during the summer months (Days 400 and 428) for nirS and nrfA and to a certain extent nirK, whereas nosZII abundance decreased and nosZI was relatively constant over this period (see Figure 3a-f). Indeed, abundances of nrfA closely followed porewater temperature variations during the second year, with the greatest covariation existing for sampling points towards the outlet (Figure 3f). The abundances of nirS also followed porewater temperature variations, but the peak in abundance appeared to lag behind the temperature maximum (Figure 3b). However, contrary to nrfA, nirS abundance peaked at locations in the first half of the DWB and increased nearest the outlet toward the end of the sampling period.
Among the investigated genes coding for NO 2 − reductases in denitrifiers (i.e., nirS, nirK), the abundance of nirS consistently exceeded that of nirK in the porewater (Figure 3b,c). The nirS/nirK ratio attained a maximum value at locations close to the inlet and during the summer months (data not shown). Relative to nrfA, the abundance of nirS + nirK was consistently greater during the two operational years (Figure 3g), with the exception of Day 400 near the outlet where nrfA > nirS + nirK. The abundance of nirS and nirK genes was also greatly in excess of the sum of the genes coding for N 2 O reductase (i.e., nosZI, nosZII; Figure 3h). Correlation analyses using data from the second year indicated that significant correlations (r s > |0.6|; p < .05) existed among geochemical parameters and between geochemical parameters and gene abundances (Table 1)

DISCUSSION
The higher abundance of functional genes involved in the denitrification pathway compared with DNRA and the absence of anammox in both the porewater and woodchip media agrees with denitrification being the major pathway for NO 3 − reduction and N removal in the DWB studied here (Nordström & Herbert, 2018) as well as in DWBs in general . Based on gene abundance data, the structures of the N-reducing communities in the woodchips, sludge, inlet water, and the porewater were significantly different from one another in many instances. The original intention of inoculating with digested sewage sludge was to promote the rapid development of a denitrifying community in the woodchip material. However, the sludge also provided a high-abundance source of organisms with the nrfA gene (i.e., genetic capacity for DNRA, an undesired side-reaction in a denitrifying bioreactor). The gene abundances determined in the woodchip samples on Day 490 likely reflect the accumulated contribution from the sludge inoculant, the original community in the woodchips (not analyzed), the inlet water, and the temporal changes of the community in response to selective pressures. However, development of the N-reducing community is mainly determined by the type of substrate used in reactors (Hellman et al., 2020) and, as shown in Figure 4, the porewater TOC concentration, the TOC/NO 3 − ratio, and DWB temperature.

F I G U R E 3
Gene abundances (a-f) and gene abundance ratios (g and h) as a function of time since start of denitrifying woodchip bioreactor operations and as a function of distance from the inlet. Abundances of (a) 16S rRNA genes and functional genes (b) nirS, (c) nirK, (d) nosZI, (e) nosZII, and (f) nrfA. Ratios (g) (nirS + nirK)/nrfA (h) and (nirS + nirK)/(nosZI+nosZII). Bioreactor temperature is plotted on secondary y axis. The x axis corresponds to day of porewater sampling (Days 57,85,113,330,365,400,428,456,and 477) T A B L E 1 Spearman's rank correlation between chemistry, temperature, gene abundances, and ratios of different gene abundances in the porewater (bioreactor centerline)

F I G U R E 4
Nonmetric multidimensional scaling ordination based on Bray-Curtis dissimilarities for the abundance of the functional genes nirS, nirK, nosZI, nosZII, and nrfA in porewater samples. Arrows show significant correlations (p < .05) between the N cycling community structure across samples and porewater chemistry. Gene names in black boxes denote species scores. TOC, total organic C The surface-bound N reducing community associated with the woodchips and containing nirS, nirK, and nosZI increased in abundance from the bioreactor inlet to the middle of the DWB (Figure 2b), suggesting that environmental conditions promoting the growth of organisms with nirS, nirK, and nosZI were better in the central region of the DWB compared with conditions closer to the inlet. A similar spatial development in gene abundance was observed in the porewater on Day 477, just prior to woodchip sampling (i.e., nirS and nirK gene abundances increased through bioreactor). Because temperature did not vary greatly within the bioreactor for any given date (data not shown), the increase in nirS and nirK gene abundances along the bioreactor flowpath was likely controlled by the relatively low concentration of organic C, elevated concentrations of NO 3 − and NO 2 − , and consequently the low TOC/NO 3 − ratio (Figure 1f). These results are contrary to our previous study (Herbert et al., 2014) that indicated a decrease in nirS abundance in a sawdust DWB with distance from the bioreactor inlet; however, this system (Herbert et al., 2014) was NO 3 − limited in regions further from the inlet but was not C limited (i.e., high TOC/NO 3 − ratio).
Based on the porewater samples (Figures 2 and 3), we clearly detected a shift in the structure of the N reducing community over time, as characterized by the change in the relative abundances of functional groups (Figure 3; Supplemental Figure S2). The increased abundance of the capacity for NO 2 − reduction relative to N 2 O reduction during the second operational year suggests a preferential selection of nirS-type denitrifiers with a truncated denitrifying pathway terminating with N 2 O. Dissolved N 2 O porewater concentrations were most strongly correlated (Table 1) with pH and the type of NO 2 − reductase (positively with nirS and nirS/nirK ratio and negatively with nirK; c.f. Barrett et al., 2016;Jones et al., 2014) and also with the relative abundance of NO 2 − reducers to N 2 O reducers (Σnir/ΣnosZ). Furthermore, the abundances of nirS and nirK (and also nrfA) during the second operational year ( Figure 3) demonstrated a clear temperature dependence, whereas abundances of nosZI and nosZII appeared to be independent of temperature, which led to a correlation between temperature and Σnir/ΣnosZ (Table 1). Warneke, Schipper, Bruesewitz, McDonald, and Cameron (2011a) showed that an increased abundance of NO 2 − reducers over N 2 O reducers at higher temperatures, including a relative enrichment in nirS-type to nirK-type denitrifiers, was associated with increased N 2 O emissions in denitrifying bioreactors. Our observations imply an increased genetic potential for the production of N 2 O from DWBs over time and with increasing temperature at HRTs of 1-2.6 d. In addition, the genetic potential for N 2 O production was spatially variable within the bioreactor, implying a geochemical control (e.g., TOC/NO 3 − , see below) as well. These implications are important contributions to our general understanding of greenhouse gas emissions from DWBs over longer time scales because they indicate the importance of understanding the temporal development of these systems. The preferential promotion of denitrifiers with a denitrification pathway terminating with N 2 O may have been an effect of the increased competition for electron donors. The simulations by Nordström and Herbert (2018) suggested a decrease in the export of acetate from the DWB studied here with time from start-up, which is consistent with the observations of other studies (Grießmeier & Gescher, 2018;Grießmeier et al., 2017). When there is competition for available organic C by denitrifiers capable of complete denitrification, electrons are preferentially directed to the NO 2 − reductase rather than the N 2 O reductase (Pan et al., 2013). This suggests that low TOC/NO 3 − ratios could lead to the preferential promotion of microorganisms with denitrification pathways truncated to N 2 O, in this case nirS denitrifiers. This was supported in this study by the negative correlation between TOC/NO 3 − ratio and the nirS/nirK ratio and its positive correlation with nirK abundances. Promotion of organisms lacking the nosZ gene is also supported by the spatially dependent relationship between N 2 O concentrations and Σnir/ΣnosZ (Figures 1d  and 3h). Further, clade II nosZ is frequently associated with nondenitrifying bacteria that do not have the genetic makeup needed for N 2 O production (Graf et al., 2014;Jones et al., 2014). The disproportional abundance of nos and nir genes suggests that the denitrification pathway is split among community members. This decoupling of the intermediate steps of denitrification onto several populations would reduce intraorganism competition (Lilja & Johnson, 2016;Pan et al., 2013). Such decoupling and metabolic specialization among truncated denitrifiers and N 2 O reducers has been inferred from metagenomes in tidal sediments and grassland soils (Diamond et al., 2019;Marchant et al., 2018). The abundance of nrfA genes detected in the porewater suggests DNRA bacteria may be competing for NO 3 − with the denitrifiers (Figures 2a and 3f). The genetic potential for NO 2 − reduction relative to DNRA (i.e., Σnir/nrfA) was negatively correlated with the TOC/NO 3 − ratio (Table 1), in agreement with studies showing that DNRA is favored by high C/N ratios (e.g., Kraft et al., 2014). Both temperature and the TOC/NO 3 − ratio had an important control on nrfA abundance ( Figure 3f) and hence NH 4 + -N concentrations, as indicated by significant correlations between these parameters (Table 1). However, the NH 4 + concentrations were overall low, suggesting limited importance of this unwanted process. Interestingly, many DNRA bacteria are also fermenting (Muyzer & Stams, 2008;Van Den Berg, Elisário, Kuenen, Kleerebezem, & van Loosdrecht, 2017), and hence the DNRA bacteria may also contribute to fermentation and thereby support denitrification despite being competitors for NO 3 − .

CONCLUSIONS
Denitrification was the major pathway for NO 3 − reduction and N removal in the DWB over two operational years at HRTs of 1-2.6 d but was associated with an increased genetic potential for N 2 O production with time. We conclude that pH and temporal and spatial changes in the relative abundance of different denitrifier genotypes, indicated by abundances of genes involved in different steps in the denitrification pathway, controlled porewater concentrations of N 2 O. Bioreactor temperature and TOC/NO 3 − ratio had a strong control on the occurrence of bacteria capable of NO 2 − reduction and those capable of DNRA (preferring high TOC/NO 3 − ). A spatially variable community, likely dominated by nirS-type denitrifiers with a truncated pathway terminating with N 2 O, developed over time in the DWB, where the supply of electron donors from substrate decomposition may be a controlling parameter on community development. Considering the significant differences in DWB chemistry and functional gene abundance between the two operational years, this study highlights the importance of distinguishing between initial variations during a start-up period, which may be extensive in length, and the long-term performance of a DWB.

A C K N O W L E D G M E N T S
This work was supported by the Luossavaara-Kiirunavaara Aktiebolag (LKAB) and VINNOVA, The Swedish Innovation Agency, under Grant 2014-011334. The authors thank Christopher Jones, Carlos Palacin-Lizarbe, and Robert Almstrand for help related to the use of Sterivex filter units. The authors thank the anonymous reviewers for all their comments and advice.

C O N F L I C T O F I N T E R E S T
There are no conflicts of interest.