Microbiome Structure and Function in Woodchip Bioreactors for Nitrate Removal in Agricultural Drainage Water

Woodchip bioreactors are increasingly used to remove nitrate (NO3–) from agricultural drainage water in order to protect aquatic ecosystems from excess nitrogen. Nitrate removal in woodchip bioreactors is based on microbial processes, but the microbiomes and their role in bioreactor efficiency are generally poorly characterized. Using metagenomic analyses, we characterized the microbiomes from 3 full-scale bioreactors in Denmark, which had been operating for 4–7 years. The microbiomes were dominated by Proteobacteria and especially the genus Pseudomonas, which is consistent with heterotrophic denitrification as the main pathway of NO3– reduction. This was supported by functional gene analyses, showing the presence of the full suite of denitrification genes from NO3– reductases to nitrous oxide reductases. Genes encoding for dissimilatory NO3– reduction to ammonium were found only in minor proportions. In addition to NO3– reducers, the bioreactors harbored distinct functional groups, such as lignocellulose degrading fungi and bacteria, dissimilatory sulfate reducers and methanogens. Further, all bioreactors harbored genera of heterotrophic iron reducers and anaerobic iron oxidizers (Acidovorax) indicating a potential for iron-mediated denitrification. Ecological indices of species diversity showed high similarity between the bioreactors and between the different positions along the flow path, indicating that the woodchip resource niche was important in shaping the microbiome. This trait may be favorable for the development of common microbiological strategies to increase the NO3– removal from agricultural drainage water.


INTRODUCTION
Leaching of nitrate (NO 3 − ) from agricultural soils to aquatic ecosystems is a growing environmental concern due to the globally increasing use of nitrogen (N) fertilizers in agriculture (Howarth et al., 2002). An effective way to mitigate NO 3 − leaching is to treat agricultural drainage water before it reaches recipient waters, e.g., by use of woodchip bioreactors (Schipper et al., 2010b). In such facilities, woodchips provide the organic carbon (C) substrates and electron donors for anaerobic microorganisms that convert NO 3 − to gaseous N species, thereby mitigating the N input to aqueous ecosystems. However, the NO 3 − removal efficiency varies among bioreactors and may not always reach the environmental goals (Carstensen et al., 2020). Optimization of woodchip bioreactors has focused mainly on abiotic factors influencing NO 3 − removal, such as dimensioning, hydraulic residence time (HRT), and physical matrix composition of the bioreactors (Cameron and Schipper, 2010;Chun et al., 2009). Yet, even under optimized conditions, NO 3 − removal varies between woodchip bioreactors, indicating additional controls, most likely linked to the microbiology of these constructed ecosystems. Therefore, management of the bioreactor microbiomes to increase the NO 3 − removal efficiency is an option that needs to be further explored and implemented (Grießmeier et al., 2017;Jang et al., 2019;Anderson et al., 2020). This necessitates a better understanding of the complex microbial populations and functional capacity of operating woodchip bioreactors that treat agricultural drainage water. So far, a number of studies have addressed the microbial composition of woodchip bioreactors under various conditions and at various scales using 16s rRNA and microbial cultivation methods (Grießmeier et al., 2017;Jang et al., 2019;Abdi et al., 2020;Hellman et al., 2021;Nordström et al., 2021). These studies allowed to link a number of factors, such as wood type, NO 3 − concentration, and pesticide contamination with the presence and activity of specific microbes. The previous studies generally highlighted the importance of the phylum Proteobacteria in the bioreactor ecosystems, but provided limited insight on the entire woodchip microbiome.
Here we used metagenomics to obtain a more complete overview of the taxonomic composition and functional potential of woodchip bioreactor microbiomes. We report the microbiome composition at a given time point in three stable operating fullscale woodchip bioreactors (BR1-BR3), where NO 3 − removal efficiencies were previously monitored (Bruun et al., 2016;Carstensen et al., 2019;Audet et al., 2021) and where concurrent metadata on water chemistry were available. The aim was to determine if the microbiomes differed among the bioreactors due to different design parameters and location in different agricultural settings. Because of a close physical location and similarity in design, the microbiomes of BR1 and BR2 were hypothesized to be more similar to each other than to BR3, which had a smaller ratio of bioreactor area to drained area and achieved greater NO 3 − removal efficiency (Audet et al., 2021). Our second aim was to address how the microbiome and functional potential differed along the flow path within the bioreactors due to gradual NO 3 − depletion, which may regulate the potential of anaerobic respiration based on other terminal electron acceptors, such as sulfate (SO 4 2− ) and carbon dioxide (CO 2 ) with resulting environmental losses of hydrogen sulfide (H 2 S) and methane (CH 4 ), respectively.

Bioreactors and Woodchip Sampling
Woodchips were sampled from three woodchip bioreactors in Denmark (2019-11-20), designed conceptually as shown in Figure 1. The annual mean NO 3 − concentrations at the inlet of BR1, BR2 and BR3 in 2019 were 49, 52 and 61 mg NO 3 − L −1 , respectively, and annual NO 3 − removal efficiencies were typically 26-56% at BR1 and BR2 and 40-70% at BR3 (Carstensen et al., 2019;Audet et al., 2021). The bioreactors had horizontal flow design (Figure 1) and received subsurface drainage water from agricultural fields mainly with annual cropping systems. BR1 and BR2 (56.214 • N,9.743 • E) had dimensions of 10 × 10 × 1.2 m (W × L × D) and were filled with a matrix of willow woodchips (2-32 mm diameter; Ny Vraa Bioenergi I/S, Denmark) mixed with seashells (Danshells A/S, Denmark). The woodchip/seashell ratio (v/v) was 50:50 for BR1 and 75:25 for BR2. Seashells were used as a material to improve the physical matrix stability and enhance water conductivity (Bruun et al., 2016). Dimensions of BR3 (55.988 • N, 10.081 • E) were 8 × 31 × 1.2 m and it was filled with willow woodchips. At all sites, the woodchip matrix was water saturated, except for the top 0.2 m, which remained unsaturated and was exposed to free air (i.e., the bioreactors were not covered by a layer of soil). BR1 and BR2 had been operating since 2012 (Carstensen et al., 2019;Hoffmann et al., 2019) and BR3 had been operating since 2015. The ratio between the bioreactor area and the drained agricultural area was 1:1300 for BR1 and BR2 and 1:830 for BR3.
Woodchip materials for metagenomic analyses were collected at three positions along the flow path in each bioreactor, i.e., near the inlet, middle and outlet (Figure 1). The positions were about 2.5, 5, and 7.5 m from the water inlet at BR1 and BR2 and about 8.5, 14.5, and 22.5 m from the water inlet at BR3. At each position, a composite sample was obtained by pooling three subsamples (taken across the bioreactor) from the upper 0.3 m of the watersaturated woodchip layer (i.e., sampled at the depth interval from ca. 0.2-0.5 cm below the bioreactor surface). The samples were handled using sterile gloves and immediately transferred to ziplock plastic bags, transported in a cooling box, and stored at −20 • C (2 days) before analysis.

Metagenomic Analyses and Microbial Diversity
Metagenomic analyses were performed after mixing the composite samples to obtain a representative woodchip subsample of 15-20 g. The woodchip material was placed in 50-mL test tubes with 35 mL of milliQ water and shaken overnight (150 rev min −1 , 20 • C). After settling, the supernatant was recovered and centrifuged to harvest the microorganisms (14,000 g, 10 min). The pellets were recovered and treated with lysozyme (Sigma Aldrich, United States) to facilitate cell lysis. DNA was extracted as previously described using QIAamp kit for QIACUBE (QIAGEN, Germany, cat no. 51126) and sequenced using an Illumina MiSeq system (Petersen et al., 2020). The sequencing depth was between 18 and 1611 mega base pairs (Mbp). The fastq sequences of the metagenomes (Supplementary Table 1) were uploaded on the metagenomics RAST (MG-RAST) server (Meyer et al., 2008) and analyzed using default parameters with a maximum E value of 10 −5 , minimum identity cut-off value of 60%, and minimum alignment length of 15 bp (Randle-Boggis et al., 2016). The sequences were compared with the Reference Sequence database (O'Leary et al., 2016) for taxonomical identification and with the Kyoto Encyclopedia of Genes and Genomes Orthology database (Kanehisa et al., 2004) for identification of functional genes associated with nitrate reduction (Zumft, 1997), including dissimilatory nitrate reduction to ammonium (DNRA) and denitrification (Supplementary Table 2). The genes chosen to assess denitrifying potential in the samples were markers of larger complexes as described in Kraft et al. (2011). More specifically, narG was chosen as marker of the abundance of nitrate reductase related to the membrane bound complex NarGHI; napA was marker for the periplasmic nitrate reductase encoded by the gene cluster napABCDE; nir, representing the sum of nirS and nirK, was marker for nitrite (NO 2 − ) reductase; norB was marker for the membrane bound nitric oxide (NO) reductase gene; and nosZ was marker for the periplasmic nitrous oxide (N 2 O) reductase.
Unconstrained non-metric multidimensional scaling (NMDS) and constrained canonical analysis of principal coordinates (CAP), both using Bray-Curtis distance, were used to analyze the dissimilarity in taxonomic composition and abundance of the identified microbial genera between the bioreactors and between the different positions along the flow path of the bioreactors (Anderson and Willis, 2003;Oksanen et al., 2020). The data was rarefied using the function vegdist from the R package vegan (Oksanen et al., 2020) before conducting CAP (function capscale) using the same R package and NMDS with the function nmds from the R package ecodist (Goslee and Urban, 2007). These two analyses focused on variations along the flow path (from inlet to outlet) and between different woodchip bioreactors. Also, total species richness (S), Shannon diversity (H ) and Pielou evenness (J) were used as ecological indices of similarity (Pielou, 1975) and tested for differences between the bioreactors and between the different positions along the flow path of the bioreactors using the procedure for two-way ANOVA analysis without replication (Zar, 2010) and reporting of p-values (Hurlbert et al., 2019). All analyses were performed using R version 3.6.1 (R Development Core Team, 2020) including tests of normal distribution and homoscedasticity.

Environmental Metadata
Sampling and analyses of environmental metadata and water chemistry at the bioreactor inlets and outlets were conducted as described in Audet et al. (2021). The bioreactors were continuously monitored and the environmental metadata retrieved for this study represented a period from within 10 days before to 1 day after the time of woodchip sampling for metagenomic analyses.

Microbiome Composition and Functional Groups
The microbiomes of the woodchip bioreactors were dominated by the phylum Proteobacteria (70-91%) with Pseudomonas as the most abundant genus, ranging from a relative abundance of 26% at the middle of BR1 to 61% at the outlet of BR3 (Supplementary Figure 1). Two other core taxa in the bioreactors were Actinobacteria and Firmicutes (Supplementary Figure 1). The high abundance of Actinobacteria was driven by multiple genera, but mainly by the genus Cellulomonas. The phylum Firmicutes was mainly represented by the genera Bacillus and Exiguobacterium.
The functional metabolic groups in the bioreactors were dominated by heterotrophic NO 3 − reducers, which in addition to Pseudomonas were represented by, e.g., Burkholderia and Enterobacter (Figure 2). The second most abundant metabolic groups were heterotrophic ferric iron (Fe 3+ ) reducers (Geobacter and Shewanella) and ferrous iron (Fe 2+ ) oxidizers with Acidovorax as the dominating genus.
FIGURE 2 | Heatmap of the relative abundance of microbial genera and functional metabolic groups in the woodchip bioreactors BR1, BR2, and BR3 close to the zones of inlet (in), middle (mid) and outlet (out) of the agricultural drainage water flow. The results are presented as the relative abundance of the hits attributed to each genus in percent of the total amount of hits obtained from the sample.
Dissimilatory sulfate-reducers were represented mainly by Desulfovibrio and endospore-forming Desulfotomaculum (Figure 2). Archaea mediating both hydrogenotrophic and acetoclastic methanogenesis were found in all bioreactors, with Methanosarcina as the dominating genus. Aerobic methanotrophs and fungi were also found in all samples, with fungi dominated by the phylum Ascomycota, especially the genus Aspergillus in both asexual and sexual (Neosartorya)states.

Functional Genes
The relative abundance of functional nitrate reduction and denitrification genes was rather similar across the woodchip samples (Figure 4), and always with negligible contribution from nrfA, which is involved in DNRA (10 times lower relative abundance than nir). Nitrate reductase genes (narG and napA) associated with reduction of NO 3 − to NO 2 − were found in highest relative abundance (Figure 4), and exceeded the abundance of NO 2 − reductases (nir) by an order of magnitude. The nosZ genes for N 2 O reductase were generally found in lowest relative abundance, except when compared to nrfA. Ratios of nir to nosZ (nir/nosZ) ranged from 1.0 to 2.7 with a mean of 2.0 ± 0.6 FIGURE 3 | Plots of non-metric multidimensional scaling (NMDS) and canonical analysis of principal coordinates (CAP) using Bray-Curtis distance based on taxonomic diversity and abundance for the three bioreactors (A,C) and for the three zones of woodchip sampling (B,D).
(mean ± standard deviation, n = 9), indicating comparable numbers of genes coding for enzymes related to production and consumption of N 2 O.

Environmental Metadata
Based on the environmental metadata (Supplementary Table 4) there were clear differences between the bioreactors. First, BR3 had a much higher HRT than BR1 and BR2, which was influenced by the hydrological settings and the bioreactor volumes. Linked to the higher HRT for BR3, this bioreactor showed lower total N and NO 3 − concentrations at the outlet, as well as higher CH 4 production as compared with BR1 and BR2. In addition, BR3 reduced the concentrations of N 2 O in the inlet water from 17.8 to 2.3 µg N 2 O-N L −1 , whereas N 2 O concentrations were lower at the inlet of BR1 and BR2 (4.0 µg N 2 O-N L −1 ), but increased toward the outlet (5.3-12.2 µg N 2 O-N L −1 ).

Functional Groups
The prevalence of the genus Pseudomonas, combined with low abundance of DNRA genes, supported previous studies showing heterotrophic denitrification as the main NO 3 − converting process in woodchip bioreactors (Schipper et al., 2010a). Thus, pseudomonads are known to thrive in soils and freshwater environments and contribute to denitrification as facultative anaerobes, since many species possess all the genes associated with complete denitrification (Lalucat et al., 2006). The relative abundance of functional denitrification genes showed a similar pattern across and within the bioreactors (Figure 4), and indicated relatively low nir/nosZ ratios, consistent with a low risk of N 2 O emissions (Conthe et al., 2019). This was in agreement with field measurements of dissolved and gaseous N 2 O fluxes from the three mature woodchip bioreactors, which showed a mean N 2 O-N emission factor of 0.6% of NO 3 -N removal (Audet et al., 2021). In studies of lab-scale woodchip bioreactors, ratios of nir/nosZ were initially high (> 10), but decreased after 6 months , which supported that the genetic potential for N 2 O reduction develops toward generally low emissions from denitrifying bioreactors (Aalto et al., 2020).
The genera Geobacter and Shewanella, found in all bioreactors, are versatile heterotrophs that may oxidize a range of organic substrates under anoxic conditions by dissimilatory reduction of Fe 3+ to Fe 2+ (Bird et al., 2011). In the bioreactors, the iron reducers co-existed with Acidovorax, known to comprise Fe 2+ oxidizing anaerobes that may utilize NO 3 − or NO 2 − as electron acceptor (Liu et al., 2019), thereby potentially contributing to iron-mediated denitrification (Bryce et al., 2018). Thus, microbial groups were identified, which could drive a dynamic cycling between Fe 3+ and Fe 2+ in the bioreactors, fueled by electrons from woodchip C and leading to denitrification. Such FIGURE 4 | Functional genes of denitrification and dissimilatory nitrate reduction to ammonium (nrfA) in samples from bioreactors BR1, BR2, and BR3 at the zones close to the inlet (in), middle (mid) and outlet (out) of the agricultural drainage water flow. The results are presented as the relative abundance of the hits attributed to each gene in percent of the total amount of hits obtained from the samples with the KO database. Encoding genes: Nar and Nap, nitrate reductase; nir, nitrite reductase; nor, nitric oxide reductase; nosZ, nitrous oxide reductase; nrfA, nitrite reductase by formate (see Supplementary Table 2). denitrification mediated by iron cycling has been reported for nitrate-rich groundwater settings (Hill, 2019) and wetland soil where the process was suggested to be of ecological relevance (Petersen et al., 2020). However, it is uncertain to what extent the entire process would be microbially mediated, since abiotic denitrification may occur by chemical Fe 2+ reactions with N intermediates, such as NO 2 − (Klueglein and Kappler, 2013;Jones et al., 2015). Such chemodenitrification has been suggested to increase the risk of N 2 O emissions from other C-rich ecosystems (Zhu-Barker et al., 2015), but so far, both the pathways, Fe availability, and ecological importance of iron-mediated denitrification in woodchip bioreactors remain to be documented.
The presence of dissimilatory sulfate-reducers in the woodchip bioreactors substantiated the risk of H 2 S formation, which has previously been indicated at high HRT when NO 3 − is fully depleted (Carstensen et al., 2019). Likewise, the presence of methanogens supported previous reports of CH 4 emissions from all three bioreactors (Bruun et al., 2017;Carstensen et al., 2019), with a higher production from BR3. Hence, the metabolic diversity of the bioreactors included anaerobes using SO 4 2− or CO 2 as electron acceptors, which thermodynamically would be favorable only after depletion of NO 3 − (Thauer et al., 1977).
Yet, there was no clear stratification of these metabolic types along the flow path of the bioreactors, indicating a latent potential of H 2 S and CH 4 formation at times of complete NO 3 − removal. Moreover, it was surprising to find similar abundances of methanogenic populations between the bioreactors since a higher abundance could have been expected in BR3, because of the higher HRT and CH 4 production. The omnipresence of fungal genera suggested a potential catabolic role related to degradation of wood-derived lignocellulose to simple C substrates for denitrification. The activity of fungi in denitrifying bioreactors may be restricted by low oxygen (O 2 ) availability, but due to the inflow of oxygenated drain water, oxic conditions may prevail at least temporarily in the woodchip matrix at concentration gradients, which decrease along the flow path and depth in the bioreactors (Jéglot, unpublished results). So far, however, little is known about the life and ecological role of fungi in denitrifying woodchip bioreactors. Remarkably, and in addition a role in lignocellulose degradation, Aldossari and Ishii (2020) newly reported the occurrence of denitrifying fungi from a woodchip bioreactor in Minnesota (United States), where these organisms seemingly contributed directly to anaerobic N transformations, thus indicating a novel physiological and ecological role to be further examined.
Although fungi have been commonly associated with degradation of lignocellulose, an increasing number of bacterial enzymes have been identified, which are also able to degrade these compounds. Proteobacteria, Firmicutes, Actinobacteria, Verrucomicrobia, and Bacteroidetes have all been found to contain genera that may attack lignin structures (López-Mondéjar et al., 2019). Therefore, the relative importance of fungi and bacteria for lignocellulose degradation in woodchip bioreactors may depend on competition and environmental constraints that regulates the activity of the different groups. In a metagenomic and metatranscriptomic study of a denitrifying beech woodchip bioreactor in Germany, Grießmeier et al. (2021) found that bacteria contributed to the wood decomposition process to similar extent as fungi. Nevertheless, a general scheme was suggested where fungi in cellulose hydrolyzing biofilms on woodchip surfaces provided soluble labile C compounds as electron donors for denitrifying bacteria predominantly occurring in the planktonic phase of the bioreactor (Grießmeier et al., 2021). In the present study, the genus Cellulomomas was prevalent, which may produce starch-, xylan-, and cellulosedegrading enzymes under microaerobic or even anaerobic conditions (Stackebrandt and Schumann, 2014) and thereby contribute to initial woodchip degradation. Based on genome analyses of Cellulomonas sp. strain WB94, isolated from a woodchip bioreactor in Minnesota, Jang et al. (2019) recently reported the occurrence of denitrification genes in Cellulomonas, which indicated that this genus might contribute to both initial degradation of complex C compounds and complete C mineralization by denitrification in woodchip bioreactors.

Microbiome Diversity
The bioreactors BR1 and BR2 deviated from BR3 in design, environmental location and characteristics of NO 3 − removal efficiency (Carstensen et al., 2019;Audet et al., 2021). Yet, a high similarity of the microbial communities between the bioreactors and between the different positions within the bioreactors was found. This indicated (i) that difference in design parameters, and especially HRT, did not significantly influence the microbiome of the woodchip bioreactors, and (ii) that microorganisms transported in the drain water from the surrounding agricultural fields were either rather similar or did not have a strong impact on the microbiomes established in the bioreactors. A metagenomic study with more than 700 semi-aquatic bacterial communities, sampled over five orders of spatial distance, emphasized the general impact of resource niches over more stochastic effects in shaping the communitylevel signatures of bacterial communities (Pascual-García and Bell, 2020). Also, a strong role of substrate type in shaping the microbial community composition was indicated in studies of lab-scale denitrifying bioreactors based on sedge, barley straw and pine woodchips . Therefore, fullscale willow woodchip bioreactors receiving nitrate-rich drainage water from agricultural fields could provide consistent resource niches for congruent microbiomes that develop across different bioreactors. If so, one could also expect that common strategies of microbial management (e.g., bioaugmentation or biostimulation) might be exploited across such woodchip bioreactors to increase the efficiency of NO 3 − removal from agricultural drainage water. Only few studies have tried to reveal how microbial communities fluctuate seasonally within individual field-scale denitrifying bioreactors. Grießmeier et al. (2021) indicated a stable microbial community composition (dominated by phylum Proteobacteria) over different seasons in a beech woodchip bioreactor, i.e., with minor influence from HRT and temperature changes. Yet, in another study, Jéglot et al. (2021) indicated divergent temperature responses of denitrifiers from a willow woodchip bioreactor, thus suggesting a potential adaptation of the woodchip microbiome to seasonal temperature changes. Thus, further studies with sampling at multiple time points for metagenomic and metatranscriptomic analyses are required to better characterize the microbial functioning of woodchip bioreactors in response to environmental temperature fluctuations.

CONCLUSION
The present study showed high similarity in the microbiomes between bioreactors and between different positions within three operating full-scale woodchip bioreactors at a given time point, indicating minor importance of design parameters and location in different agricultural settings. Major microbial groups and genes in NO 3 − reduction pathways were associated with denitrification with only minor contribution from DNRA. Functional diversity related to denitrification further included microbial lignocellulose degradation and bacterial Fe cycling, which are processes that should be further examined for their ecological role in efficient NO 3 − removal from agricultural drainage water.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/ repositories and accession number(s) can be found in the article/Supplementary Material.