Samples were collected on the R/V Nathaniel B. Palmer during June to July 2013 (cruise NBP 1305) in the OMZ of the ETSP at the off-shore station (BB1; 14.0°S, 81.2°W) and the coastal station (BB2; 20.50°S, 70.70°W) (Supplementary Figure 1). Particulate material was collected in Niskin bottles mounted on the standard conductivity-temperature-depth (CTD) rosette system (Seabird Electronics, Seattle, WA, United States) at four depths at each station (BB1: 60, 130, 300, and 1000 m; BB2: 60, 115, 300, and 1000 m) and concentrated by filtration (up to 4 L) through Sterivex filters (0.22 μm). Filters were flash frozen in liquid nitrogen onboard and stored at -80°C until DNA and RNA extraction was performed.

Temperature, salinity, sigma theta, bottom depth and pressure at each station were measured on the SBE 911+ CTD system. Fluorescence, representing chlorophyll a, was measured using a single channel fluorometer (Wet labs, Philomath, OR, United States) mounted on the CTD. Oxygen distributions were determined using the STOX sensor (detection limit = 10 nM) mounted on the CTD rosette (Revsbech et al., 2009). Ammonium, nitrite and nitrate concentrations were measured using standard colorimetric protocols (UNESCO, 1994). N₂O concentration was determined using mass spectrometry (Ji et al., 2015). N^∗ is the deviation of measured dissolved inorganic nitrogen (DIN = nitrate + nitrite + ammonia) from predicted DIN by Redfield ratio and the world-ocean nitrogen to phosphate regression relationship (Deutsch et al., 2001). Environmental data were reported by Ji et al. (2015) and are provided in Supplementary Table 1.

Both DNA and RNA were extracted from eight Sterivex filters using the plant tissue protocol of the All Prep DNA/RNA Mini Kit (50) using a QIAcube (Qiagen). Reverse transcription from RNA to cDNA was performed using SuperScript^® III First-Strand Synthesis System for RT-PCR (Invitrogen^TM by Life Technologies^TM). Excess RNA was removed by RNase at the end of the synthesis.

The abundance of total and active nosZ assemblages were estimated by quantitative PCR (qPCR) using SYBR^® Green based assays using protocols described previously (Jayakumar et al., 2013). Primers nosZ1F and nosZ1R (Henry et al., 2006) were used to amplify a 259-bp conserved fragment of the nosZ gene. Known quantities (∼20–25 ng) of DNA and cDNA samples were assayed along with a minimum of five serial dilutions of plasmids containing nosZ gene, no template controls and no primer controls, all in triplicate on the same plate. To maintain continuity and consistency among qPCR assays, a subset of samples from the first qPCR assay was run with subsequent assays and fresh standard dilutions were prepared for each assay. DNA, cDNA and the concentrated standards were quantified prior to every assay using PicoGreen fluorescence (Molecular Probes, Eugene, OR, United States) calibrated with several dilutions of phage lambda standards, to account for DNA loss due to freeze thaw cycles. qPCR assays were run on a Stratagene MX3000P (Agilent Technologies, La Jolla, CA, United States). Automatic analysis settings were used to determine the threshold cycle (Ct) values. The copy numbers (number of copies of the gene sequence detected in the sample) were calculated according to: Copy number = (ng ^∗ number/mole)/(bp ^∗ ng/g ^∗ g/mole of bp) and then converted to copy number per ml seawater filtered, assuming 100% extraction efficiency.

DNA and cDNA qPCR products were used as targets for microarray experiments to characterize the community composition of total and active nosZ assemblages, respectively. Triplicate qPCR products from each depth were pooled. nosZ gene targets were purified and extracted from agarose gels using the QIAquick gel extraction kit (Qiagen). Purified DNA qPCR products from eight depths and cDNA qPCR products from seven depths were used to prepare targets for microarray analysis.

Microarray targets were prepared from the qPCR products following the protocol of Ward and Bouskill (2011). Briefly, dUaa was incorporated into purified DNA and cDNA during linear amplification using the BioPrime kit (Invitrogen^TM). The dUaa-Klenow product was labeled with Cy3 (dissolved in dimethyl sulfoxide), purified using QIAquick columns (Qiagen) and quantified by Nanodrop 2000 (Thermo Scientific). Duplicate Cy3 products for each sample were hybridized at 65°C overnight (16 h) onto replicate microarrays under ozone free conditions. Hybridized microarrays were washed and scanned with an Axon 4300 laser scanner.

The microarray (BC016) contains 114 nosZ archetype probes. Each probe is a 90-bp sequence comprised of a 70-bp nosZ gene fragment and a 20-bp control region. Each archetype probe represents, and hybridizes with, all nosZ sequences with >85% identity, based on published sequences available in 2013. There are 71 NosZ archetypes, which represent typical or Clade I nosZ genes, and 43 WNZ archetypes, which represent the atypical or Clade II nosZ genes. The development of the microarray is described in Jayakumar et al. (in preparation) and the sequences are shown in Supplementary Table 2.

Fluorescence signal intensities for nosZ probes hybridized to the microarrays were obtained using GenePix Pro 7 software. The fluorescence ratio (FR) of each feature is defined as the ratio Cy3/Cy5 (70-mer probe/20-mer standard for each feature). The FR for each nosZ archetype was calculated as the average of probe signal intensities for duplicate features on the same microarray. Normalized fluorescence ratio (FRn) was calculated by dividing the FR of each nosZ probe by the maximum nosZ FR on the same microarray. FRn is the proxy of the relative abundance of each archetype and was used for further analyses.

Detrended correspondence analysis (DCA) was performed to analyze the overall microbial community composition. A dissimilarity test was performed using Permutational Multivariate Analysis of Variance (adonis). α-diversities (Shannon diversity indices) of total and active nosZ assemblages were calculated. β-diversities (Bray–Curtis dissimilarities) of total and active nosZ assemblages between different sites (i.e., depths) were calculated to perform a Mantel test. The Mantel test was used to determine significant environmental variables correlated with microbial community composition. These analyses were carried out using the vegan package in R (version 3.3.1). A maximum likelihood phylogenetic tree was built from aligned archetype sequences with MEGA 7 software. FRn values for each archetype at different depths from both stations were visualized on the phylogenetic tree by iTOL¹. The copy number of nosZ genes at each depth is given as mean (± standard error) of the qPCR triplicates.

At stations BB1 and BB2, the continuously undetectable O₂ concentration, the local nitrite maximum and the nitrate deficit at intermediate depths (130–370 m at BB1; 75–400 m at BB2) all indicated the presence of ODZs (gray areas in Figure 1). Sampling depths were chosen to represent water column features defined by oxygen concentration, as measured with the in situ STOX sensor: oxygenated surface water, upper oxycline [characterized by sharp O₂ concentration gradient ranging from saturation to below detection limit (<10 nM)], top of the ODZ (O₂ concentration <10 nM), core of the ODZ (O₂ concentration <10 nM) and lower oxycline (O₂ concentration >10 nM). The abundance of the total nosZ genes ranged from 24.1 (±1.4) copies mL^-1 in a sample from the lower oxycline to 636.4 (±28.3) copies mL^-1 in a sample from the ODZ (Figure 1). As for the active nosZ assemblage, the lowest abundance of active nosZ genes was 5.1 (±0.5) copies mL^-1 in a sample from the lower oxycline and the highest abundance was 604.6 (±103.7) copies mL^-1 in a sample from the surface water.

The abundance of total and active nosZ assemblages showed different distribution patterns at the two stations (Figure 1). At station BB1, the abundance of both total and active nosZ genes was highest in a sample from the surface water, and decreased with depth. At station BB2, the abundance of both total and active nosZ genes peaked in samples from the ODZ and was lowest in samples from the lower oxycline. The active nosZ genes were most abundant in the sample from 300 m. However, the total nosZ genes were most abundant in the sample from 115 m, where the abundance of the active nosZ genes was only 3% of the total.

The distribution of FRn of the total or active archetypes was similar across all depths within the same station (Figure 2). The average α-diversity was not significantly different (student’s t-test, P = 0.102) between the total assemblages (3.21) and the active assemblages (2.60) (Table 1). The least diverse total assemblage was from the lower oxycline (1000 m of station BB2), but the two least diverse active assemblages were from the ODZs (130 m of BB1 and 300 m of BB2).

The FRn distribution of nosZ archetypes showed that a very limited number of archetypes dominated the total or the active nosZ assemblages (Figure 2). Dominant archetypes were affiliated with bacteria from various environments, including salt marsh, soil, marine sediment, marine hot spring and activated sludge of a wastewater treatment plant (Supplementary Tables 3, 4). The FRn of the top five most abundant archetypes accounted for 48.9 to 83.3% of the total nosZ hybridization signal (Figure 3A and Supplementary Table 3). Notably, the highest percentage (83.3%) was from the sample from the lower oxycline (1000 m) at station BB2 and the most abundant archetype (NosZ42, an uncultured clone of nosZ gene derived from salt marsh sediments; Kearns et al., 2015) accounted for 31.6% of total FRn. However, this archetype was not among the top five archetypes of the active nosZ assemblage in the same sample (Figure 3B and Supplementary Table 4). The two most dominant typical nosZ archetypes in the total assemblage were NosZ6, derived from an uncultured clone from salt marsh sediments (Kearns et al., 2015), which is closely related to Marinobacter, and NosZ65, derived from Marinobacter sp. BSs20148 from marine sediment (Song et al., 2013) (Figure 3A and Supplementary Table 3). In contrast, NosZ6 and NosZ65 were not dominant in the active assemblage (Figure 3B and Supplementary Table 4).

The sample from the ODZ (130 m) at station BB1 illustrates the contrasts observed between total and active nosZ assemblages. The FRns of the top three dominant archetypes (WNZ21, WNZ16 and NosZ65) in the total nosZ assemblage were comparable to each other, and they constituted 56.1% of the total community (Supplementary Table 3). However, NosZ65 was nearly undetectable in the active assemblage and the top two dominant archetypes (WNZ21 and WNZ16) accounted for 87.8% of the active assemblage in the same sample (Supplementary Table 4). The total nosZ assemblage was much more diverse than the active assemblage at this depth (Table 1). The representative sequences of WNZ21 and WNZ16 archetypes were derived from nosZ gene sequences of Anaeromyxobacter dehalogenans strain DCP18 (Chee-Sanford et al., unpublished) and an uncultured bacterium clone obtained from agricultural soils (Sanford et al., 2012), respectively. WNZ21 and WNZ16 archetypes were not only dominant in the active assemblage in one sample from the ODZ, but were among the top five abundant archetypes of both total and active assemblages in almost all samples (Figure 3).

Functional gene microarrays were used to describe the community composition of nosZ assemblages. FRn values from duplicate microarrays replicated well (r² = 0.802–0.997) (Supplementary Figure 2) and each pair of duplicates clustered together in the DCA plots (Figure 4).

The two-dimensional DCA model including both DNA and RNA microarray results explained 43.7% of the community composition of nosZ assemblages with 31.6% explained by the first axis and 12.1% explained by the second axis (Figure 4A). The clearest pattern was the clear separation of total (filled symbols) and active (open symbols) nosZ assemblages, indicating that they were different from each other. The significance (P < 0.001) of the difference between total and active assemblages was confirmed by the dissimilarity test (Table 2). Therefore, the community composition of total and active nosZ assemblages was further analyzed by two DCA models separately to better examine other patterns.

The DCA model of DNA microarray results explained 63.9% of the composition of the total nosZ assemblage (Figure 4B). The total nosZ community composition revealed site difference and O₂ dependence. Samples from the same station clustered together in the DCA model, indicating community composition was different between the off-shore station (BB1) and the coastal station (BB2). The site difference was statistically significant (P < 0.001) based on the dissimilarity test (Table 2). Besides the geographical pattern, composition of the total nosZ assemblage was also affected by O₂ concentration (Figure 4B). Samples from the ODZs clustered together, while samples with higher O₂ concentrations were distinct from the ODZ samples and different from each other. O₂ concentration of the seawater might not be the most important driver of microbial community composition, however, since the O₂ pattern was not captured by either axis of the DCA model.

The DCA model of RNA microarray results explained 40.2% of the community composition of the active nosZ assemblage (Figure 4C). Significant site difference (P = 0.039) of the active community composition was also revealed by the DCA model (Figure 4C) and the dissimilarity test (Table 2). However, the clustering based on O₂ concentration that was observed in the total nosZ assemblage was not observed for the active nosZ assemblage.

The composition of total and active nosZ assemblages was correlated with different environmental variables based on a Mantel test (Table 3). Relative depth, nitrate concentration, temperature, density (sigma theta) and pressure were significantly related to the β-diversity of the total nosZ assemblage. However, the β-diversity of the active nosZ assemblage was significantly related to N₂O concentration, nitrite concentration and fluorescence. Bottom depth, which was dramatically different between two stations and is a proxy for important ecological differences between the two sites, was a significant factor for both total and active nosZ assemblages.

Oxygen minimum zones are sites of high N₂O flux to the atmosphere (Law and Owens, 1990; Arévalo-Martínez et al., 2015). N₂O consuming organisms are the only biological sink for N₂O. Hence their abundance and community composition in the OMZ may be important in understanding the N₂O flux. The abundance of total and active N₂O consuming bacteria in the OMZ of the ETSP was estimated by measuring nosZ gene copy number (Figure 1). The relationship between abundance of N₂O consuming bacteria and depth in this study differed from that of denitrifiers indicated by nirS gene copy number at the same stations (Ji et al., 2015): the abundances of the total and active N₂O consuming bacteria in the surface water were similar to or higher than those in the ODZs (Figure 1), but the abundance of denitrifiers in the surface water was two orders of magnitude smaller than that in the ODZs. This difference suggests that the two genes represent functionally different groups.

nirS and nosZ also differed in their absolute abundance. The highest abundance of N₂O consuming bacteria was only a few hundred copies mL^-1, which was three orders of magnitude smaller than the highest abundance of denitrifiers measured at the same stations (Ji et al., 2015). It is assumed that both genes are present in the genome as single copy genes, although there are exceptions for nosZ (Sanford et al., 2012). One possible explanation for the differences in both distribution and abundance of the nosZ assemblage and the nirS assemblage is that not all N₂O consuming bacteria contain the complete denitrification gene sequence (Sanford et al., 2012). The atypical nosZ genes are associated with bacteria that lack the other steps in the conventional denitrification pathway. Notably, bacteria with only nosZ genes but no other denitrification genes were overrepresented in the genomes of marine bacteria compared to other ecosystems (Graf et al., 2014). nirS, however, was preferentially associated with bacteria that contained a complete denitrification pathway (Graf et al., 2014).

Another contributing factor may be the specificity or bias of the PCR primers. The nosZ primers used in this study were optimized to amplify all known nosZ sequences as of 2006, and should therefore represent the large database of both terrestrial and marine sequences available at the time. However, it is clear that they might underrepresent the atypical N₂O consuming bacteria, which were not known at the time. The nirS primers used in the previous analysis of these samples (Ji et al., 2015) are potentially biased toward marine sequences (Braker et al., 1998) and may underrepresent more diverse sequences now available from other environments. One way to improve the nosZ coverage is to use multiple primer sets targeting different groups of nosZ archetypes.

The N₂O consuming bacteria are a small component of the total microbial assemblage, but are still quite diverse (Jones et al., 2013), so they are difficult to characterize by pure culture or metagenomics. The microarray, which was designed specifically to target N₂O consuming bacteria using more than 100 nosZ gene probes (Figure 2), may be a better tool to capture these underrepresented organisms without cultivation or detection of rare sequences in complex metagenomic datasets. The high reproducibility of microarrays reported previously (Bulow et al., 2008) was confirmed in this study in that duplicates for each sample run on two different microarrays clustered together in DCA (Figure 4) and had high r² of linear regressions (Supplementary Figure 2).

Based on the FRn values of diverse nosZ archetypes determined by microarray hybridization (Figure 2), a very limited number of archetypes dominated the total or the active assemblages. Moreover, the top five active archetypes accounted for larger percentage of the assemblage than that of the total archetypes (Figure 3), consistent with the less diverse active assemblage compared to the total assemblage (Table 1). These findings imply that although the total nosZ assemblage is very diverse, a few nosZ archetypes might be the major contributors of N₂O consumption at the study sites. The relative abundance of active archetypes, however, might uncouple that of enzymes of different archetypes and/or the contribution of different archetypes to the N₂O consumption rate due to different stabilities of enzymes from different archetypes.

The active nosZ assemblage was different from the total assemblage in both abundance profiles and community composition as detected by qPCR and nosZ microarray hybridization analysis, respectively. The highest abundance of the total N₂O consuming bacteria indicated by nosZ DNA copy number was 636.4 (±28.3) copies mL^-1 in the sample from the ODZ (115 m) at station BB2, but the abundance of the active bacteria in the same sample was only 21.1 (±5.2) copies mL^-1 (Figure 1). In the sample from oxygenated surface seawater (60 m) at station BB1, the abundance of the active N₂O consuming bacteria was 604.6 (±103.7) copies mL^-1, but the abundance of the total bacteria was only 357.2 (±12.5) copies mL^-1 (Figure 1).

Significant differences between the community composition of total and active N₂O consuming bacteria were indicated by the results of DCA and dissimilarity test (Figure 4A; Table 2). The different community composition was attributed to the differences in the distribution of FRn of more than 100 archetypes, especially the dominant ones. NosZ6 and NosZ65 were dominant in the total nosZ assemblage but were minor components in the active assemblage. NosZ42 was the most abundant archetype in the total nosZ assemblage from the lower oxycline (1000 m) at station BB2, but was not among the top five archetypes of the active assemblage. These differences between active and total nosZ communities are consistent with observations from soil and salt marsh sediments in which the active component of the microbial assemblage was apparently more responsive to environmental conditions (Barnard et al., 2013; Kearns et al., 2016).

Despite the differences between the total and the active N₂O consuming assemblages, their distributions both depended on geographic location. The composition of total and active nosZ assemblages was significantly different between the coastal and the off-shore stations as indicated by the results of DCA and dissimilarity test. Since the two stations shared similar dominant nosZ archetypes, the geographical divergence mainly reflected differences among the large number of rare archetypes between the two stations. The more negative N^∗ at the coastal station BB2 (Supplementary Table 1) indicates more intense nitrogen-loss fueled by more organic matter. The different amount of organic matter, which supports the metabolism of heterotrophic nosZ bacteria, might partially contribute to the geographical differences of nosZ assemblages. Geographical differences might also result from different nutrient sources at the two stations, since their distance to the sediment (bottom depth) and to the shore were dramatically different. The dependence of geographic location was also observed for ammonia oxidizing archaea (Peng et al., 2013) and for nirS denitrifiers (Jayakumar et al., 2013) in the ETSP and Arabian Sea OMZs.

In addition to geographical patterns, the total and active assemblages were correlated with different environmental variables (Table 3). Depth and environmental parameters that co-varied with depth (including temperature, density and pressure) were major drivers of the β-diversity of the total N₂O consuming assemblages, implying different organisms coexist in the water column by occupying different ecological niches. However, N₂O concentration difference was a major driver of the β-diversity of the active nosZ assemblages, implying the active nosZ community was a better indicator for N₂O consumption potential.

The role of the nosZ assemblage in oxygenated seawater has been ignored because N₂O consumption is considered the least oxygen tolerant anaerobic step in the conventional denitrification pathway (Zumft, 1997). However, nosZ genes were abundant in oxygenated surface water in the Southern Indian Ocean (Raes et al., 2016) and nosZ mRNAs were detected in the oxic regions in the Arabian Sea (Wyman et al., 2013). Our study confirmed that a nosZ assemblage was not only present but also active in oxygenated surface water in the OMZ of the ETSP. In particular, atypical nosZ archetypes, usually associated with N₂O consuming bacteria lacking a complete denitrification pathway, were present and active in surface waters. In addition, the most abundant archetypes of total and active nosZ communities were both atypical nosZ archetypes (WNZ21 and WNZ16), implying the significant contribution of atypical archetypes to the nosZ communities and the necessity to consider atypical archetypes while analyzing the potential of N₂O consumption.

N₂O reductase enzymes from denitrifiers had very low O₂ tolerance (Bonin et al., 1989; Körner and Zumft, 1989); on the contrary, nosZ assemblages were detected in the oxygenated surface waters and O₂ concentration was not significantly correlated with the active microbial community, as indicated by DCA and Mantel test. The survival of N₂O consuming bacteria in oxic layers and their O₂-independence might be attributed to anoxic micro-environments created by phytoplankton microaggregates or particles. Free-living and particle-associated microbes from the same seawater sample can have different community compositions (Delong et al., 1993). More specifically, a recent study in the OMZ of the ETSP showed that nosZ mRNAs were 28-fold more abundant on particles (>1.6 μm) compared to free-living microbes (0.2–1.6 μm) (Ganesh et al., 2015). Additionally, nosZ mRNA co-occurred with the cyanobacterium Trichodesmium in oxic water in the Arabian Sea (Wyman et al., 2013). Consistently, fluorescence, a proxy for chlorophyll a, was significantly correlated with the β-diversity of the active nosZ assemblages in this study (Table 3).

The active nosZ community in oxygenated surface water might capture N₂O produced in deeper seawater and thus reduce the flux into the atmosphere. Thus, evaluating the nosZ community is essential to the prediction of the oceanic N₂O emissions. Moreover, the oceanic N₂O emissions represent net fluxes, which are controlled by both N₂O production and N₂O consumption. Some N₂O flux models (Suntharalingam and Sarmiento, 2000; Martinez-Rey et al., 2015; Trimmer et al., 2016) do not parameterize N₂O consumption, and other models either consider N₂O consumption only in suboxic or anoxic waters (Cornejo and Farías, 2012; Babbin et al., 2015) or estimate N₂O consumption assuming it is constrained by O₂ concentration (Zamora et al., 2012). Failing to consider the O₂-independent, non-denitrification N₂O consumption potential in these O₂ forcing models might contribute to their uncertainty and the variation among different models. Additionally, the N₂O consuming organisms have not been fully investigated. Besides denitrifiers and atypical N₂O consuming bacteria analyzed in this study, other organisms (Trichodesmium and Crocosphaera) also exhibited N₂O consuming capacity under laboratory conditions (Farías et al., 2013), suggesting that their significance in the environment warrants further investigation.