Metagenomes from microbial communities were obtained from publicly available datasets across four station locations in the Pacific ( Figure 1 ). Detailed depth profiles of cellular metagenomes from St ALOHA were obtained as part of the Hawaii Ocean Time-series (HOT) (22°45°N and 158°W) for May (HOT 272), August (HOT 275), and November 2015 (HOT 278) and were downloaded from Bioproject PRJNA352737 (Luo et al., 2020). Nutrient and CTD measurements for these cruises can be downloaded with the Hawaii Ocean Time Series Data Organization and Graphical System (HOT-DOGS) application at the University of Hawai’i at Mānoa (https://hahana.soest.hawaii.edu/hot/hot-dogs/) and in the supplement to the original paper (Luo et al., 2020). For Eastern Tropical South Pacific (ETSP) Station 9 (13°S and 82.2°W) and Station 17 (16.7°S and 79°W) sampled in July 2013, metagenomes were obtained from Bioproject PRJNA704804 (Fuchsman et al., 2022). Hydrographic and nutrient data from this ETSP cruise have previously been published (Peters et al., 2018; Fuchsman et al., 2022). Cyanate data from this ETSP cruise were published, but urea was not measured (Widner et al., 2018b) ( Supplementary Figure S3 ). Metagenomes from ETNP Station 136 (17.04°N, 106.54°W) in April 2012 can be found in Bioproject PRJNA350692 (Fuchsman et al., 2017), and nutrient data for ETNP ST136 can be seen in Fuchsman et al. (2017). Cyanate and urea concentration data from this ETNP cruise are published (Widner et al., 2018a), and data for ST136 is included here ( Supplementary Figure S4 ). Temperature, salinity, and nutrient metadata for all samples can be found in Supplementary Table S1 .

In this paper, we add five metagenomes to existing ETNP ST136 (600 m [deep ODZ], 800 m [ODZ boundary], 1,000 m [hypoxic]) and ETSP ST17 (80 m [oxic], 125 m [hypoxic]) depth profiles. As described in Fuchsman et al. (2017, 2022), 4 L of water from the CTD was filtered through 0.2 µm SUPOR filters on board the ship and immediately frozen at -80°C. These filters were extracted in June 2023, following the protocol in Fuchsman et al. (2017); DNA was extracted from filters using freeze–thaw followed by incubation with lysozyme and proteinase K and phenol/chloroform extraction. Libraries were created and run on a NovaSeq 6000 at the Northwest Genomics Center (Seattle, WA). The two new ETSP metagenomes can be found with the rest of the 2013 ETSP metagenomes at Bioproject PRJNA704804, but the three new 2012 ETNP metagenomes are in Bioproject PRJNA1083228.

All ETSP samples, new and old, were assembled separately, and assemblies are available at Bioproject PRJNA704804. Samples were quality-screened using Trimmomatic v0.39 (Bolger et al., 2014) and assembled using Megahit v1.2.8 (Li et al., 2015). To find the genes of interest in the assembled contigs, genes were called and annotated with Prokka 1.14.6 (Seemann, 2014). Sequence information about the new metagenomes can be found in Supplementary Table S2 .

Transcripts of sequences corresponding to station P2 (16°55 N, 107°09 W) in the ETNP in May 2018 were obtained from the oxycline and top of the ODZ (76–150 m) (Mattes et al., 2022), downloaded from BioProject PRJNA727903, and processed following the same methodology described below. These transcript samples were sampled over multiple days at the station, and thus, due to variations such as internal waves, oxygen concentrations from the samples do not change linearly with depth (Mattes et al., 2022). STP2 is geographically close to ETNP ST136, examined for metagenomes.

Phylogenetic trees for urease (ureC) ( Supplementary Figure S1 ) and cyanase (cynS) ( Supplementary Figure S2 ) from Widner et al. (2018a) and nitrite oxidoreductase (nxrB) from Fuchsman et al. (2017) were updated using assembled sequences from the ETSP, marine single-cell genomes (Berube et al., 2018), and eukaryotic assemblies from the MMETSP (Keeling et al., 2014). Hydrazine oxidoreductase (hzo) data have already been published (Fuchsman et al., 2022). Nitrate reductase gene (narG) data for Thioglobaceae was obtained using a phylogenetic tree modified by Fuchsman et al. (2017).

The RNA polymerase (rpoB) phylogenetic tree has already been published (Fuchsman and Hays, 2023). A set of reference amino acid sequences (encoding for the genes nxrB, ureC, and cynS separately) were used as a query to be blasted against custom databases for the assembled protein ETSP 2013 dataset and for single-celled genomes and eukaryotic assemblies using blastp (Altschul et al., 1997). References and assembled sequences were aligned using MUSCLE v3.8.1551 (Edgar, 2004). The alignment was used to build a phylogenetic tree with bootstrapping using RAxML-ng (Kozlov et al., 2019) and later visualized using FigTree v1.4.4 (https://github.com/rambaut/figtree/releases). The tree constructed using RAxML-ng was used as the reference frame to place the metagenomic reads from the ETSP, ETNP, and HOT datasets.

As described previously (Fuchsman et al., 2023), read placement was done by recruiting short metagenomic reads via a tblastn search of the metagenomes (using an e-value of ≤ 5). The reads were trimmed to remove Ns and converted to amino acid sequences. After the quality trimming, only sequences longer than 100 bp (33 amino acids) were kept. The amino acid-translated reads were then aligned against the reference sequences using PaPaRa 2.0 (Berger and Stamatakis, 2011). The nonoverlapping paired read ends were then combined into one sequence in the same alignment using a Python script, and they were later placed in the tree using EPA-ng v0.3.6 with filter-max as 1 (Barbera et al., 2019). Placed reads have a pendant length, indicating the similarity between a query read and the location it places on the tree. Reads that were placed with a pendant length greater than 2 were removed; 1% or less of reads were removed at this step. To sort the reads into taxonomic groups, the reads in each group were enumerated using the assigned subcommand of Gappa v0.6.1 and a taxonomy file listing the taxonomy of the tree reference sequences (Czech et al., 2020). Taxonomic read counts were normalized using the method previously described by Fuchsman et al. (2019) where normalization factors for each sample were determined by dividing 48,556,135 (the number of reads in the 100-m ETNP sample) by the number of good-quality reads in the sample. The read counts were multiplied by the sample normalization factor, divided by the gene length (1,704 base pairs (bp) for ureC, 468 bp for cynS, 3847 bp for rpoB, and 1,275 bp for nxrB), and then multiplied by 100 to make visualization easier. We refer to these numbers as normalized reads. We then calculated the percentage of the total microbial community containing cynS and ureC by dividing the normalized reads for cynS or ureC by the combined Bacteria- and Archaea-normalized reads for single-copy core gene RNA polymerase (rpoB). Similarly, we calculated the proportions of individual taxa containing cynS or ureC by dividing by rpoB normalized reads for that taxa. The numerical results from these analyses can be found in Supplementary Table 3 .

The errors in these analyses can be broken into three types of error. For HOT 272 (May 2015), 700 m and 1,000 m samples were sequenced twice. Our analysis of these duplicate samples produced almost identical results; the average standard deviation between ureC groups was 0.05 normalized reads for 700 m and 0.15 normalized reads for 1,000 m. These data indicated that we have good reproducibility. However, there are two systemic caveats to this analysis: (1) some organisms may be missing from the phylogenetic trees, and (2) gene lengths vary between organisms. In the case of incomplete phylogenetic trees, it is less accurate to place sequences when there are no near relatives to that sequence. We have minimized this source of error by building trees both with known genomes and with environmentally assembled contigs. In terms of length variability, any differences between organisms’ genes and our estimates of gene length will result in misestimated gene ratios. However, the magnitude of such an error is small: the average length of reference genes composing the trees here were cynS = 460 bp ± 40 bp, ureC = 1,706 bp ± 23 bp, and rpoB = 3,815 bp ± 333 bp. RNA polymerase (rpoB) was the longest gene with the most variability. This variability in rpoB would change a proportion of ureC/rpoB of 100% by ± 10%.

In this section, we will describe the ratios of taxon-specific functional genes ureC and cynS to taxon-specific rpoB reads for a variety of groups important to biogeochemical cycling. We also explore the proportions of the taxon-specific functional genes nitrite oxidoreductase (nxrB; Nitrospina), nitrate reductase (narG; Thioglobaceae) and hydrazine oxidoreductase (hzo; Anammox) to support our taxon-specific rpoB data.

Although the transcript data and the metagenomic reads data from the ETNP shown here do not correspond to the same sample collection date (April 2018 versus April 2012, respectively), we observed similar trends between datasets. When ureC transcripts from STP2 in the ETNP were also analyzed, we observed several groups: Picocyanobacteria had the highest abundance of ureC transcripts, with 16.4 normalized ureC transcripts at depths of 106 m and 150 m, and the general alphaproteobacteria and gammaproteobacteria groups also had significant numbers of transcripts. Key N cycling microbes had low abundances of ureC transcripts: Thaumarchaeota had 1 to 2.5 normalized ureC transcripts at oxic depths; Anammox bacteria (Cand. Scalindua) had 2.4 normalized ureC transcripts at 112 m depth and 2.7 at 150 m depth; and Nitrospina ureC transcript abundances were low but present at all the sampled depths, with abundances between 0.2 and 1.8 normalized ureC transcripts ( Figure 11A ). Verrucomicrobia had low numbers of transcripts throughout, with three normalized reads at 150 m. The cynS transcripts in station P2 were dominated by transcripts of Anammox bacteria (Cand. Scalindua) ( Figure 11B ), with 104.4 normalized transcripts at 112 m and 83.2 normalized transcripts at 150 m depth. After Anammox bacteria (Cand. Scalindua), the second-highest number of cynS transcripts belonged to "other bacteria". Nitrospina cynS transcripts were found at 106 m and 112 m with a maximum of 2.9 normalized cynS transcripts at 112 m depth ( Figure 11B ).

Unlike Mao et al. (2022), who found that cynS in larger-size fractionated samples (>0.8 µm) was dominated by eukaryotic phytoplankton, we only observed small numbers of algae-derived cynS at all stations. The highest normalized reads were found at the DCM for each station and corresponded to Pelagomonales. The lack of algae-derived cynS in the ODZs was not surprising, as eukaryotic algae do not thrive in the ODZs (Fuchsman et al., 2022; Wong et al., 2023), but algae cynS reads were also not particularly abundant at HOT ( Supplementary Figure S5 ). This finding could be because HOT is dominated by Picocyanobacteria rather than eukaryotic algae, though eukaryotic algae are definitely present (Rii et al., 2016). Part of the difference between our results and those of Mao et al. (2022) could also be because our data were obtained from bulk water samples dominated by bacteria, while Mao et al. examined > 0.8 µm fractions enriched in eukaryotic algae. In either case, our results imply that eukaryotic algae are present but not dominant consumers of cyanate in the systems studied here.

Both Nitrospina and Cand. Scalindua Anammox bacteria have representatives with cyanase and urease, but the two bacterial groups appear to have different preferences for urea and cyanase ( Figures 3 , 4 ). The Anammox (Cand. Scalindua) bacterial community, which can use these compounds for dissimilatory energy production, had a higher proportion of Anammox bacteria containing the cynS gene compared to ureC, indicating a preference for cyanate over urea ( Figure 3 ). These results are consistent with rate data from our 2013 ETSP cruise, which showed that Anammox N₂ production rates could be supported by the N in cyanate while urea could only stimulate Anammox N₂ production rates after a 1.5-day lag time (Babbin et al., 2017), indicating that the Anammox bacteria were not utilizing the urea in situ at the time of sampling. Similar results were observed with the ETNP transcripts ( Figure 11 ), where the highest abundance of the cynS transcript in the ODZ core corresponded to Cand. Scalindua and the abundance of ureC in the transcript data was low in numbers compared to cynS. These transcripts imply that Anammox bacteria (Cand. Scalindua) were actively transcribing the cynS gene to metabolize cyanate in the ETNP, and assimilation rate data from the ETNP ODZ indicate that cyanate was used both for C assimilation and N dissimilation (Widner et al., 2018a). Contrastingly, the proportion of Nitrospina with ureC was higher than for cynS both at HOT and in ODZ waters, indicating a metabolic preference for urea as a reduced nitrogen form over cyanate. Additionally, expression of Nitrospina ureC transcripts was similar in numbers and present at all depths, but Nitrospina cynS transcripts showed lower expression (< 1 cynS transcript) in oxic water but increased to 2.2 and 2.9 at 106 and 112 m depth in the ODZ ( Figure 11 ). These observations suggest that Nitrospina could be actively transcribing both genes to possibly metabolize both cyanate and urea in the upper ETNP ODZ. The transcript numbers were low, but they contrasted with the distribution of ureC and cynS metagenomic reads ( Figure 4 ). A possible explanation for the presence of cynS transcripts from Nitrospina in the upper ODZ could be due to the high cyanate concentrations (40 nM) at the top of the ETNP ODZ ( Supplementary Figure S3 ; Widner et al., 2018a). The differences in the organic N preferences of Anammox and Nitrospina bacteria may prevent them from competing where their depth ranges overlap, except when an abundantly reduced N resource is available.

Although Thaumarchaeota does not possess the cynS gene, members of the group can use cyanate and urea both for nitrification (dissimilatory) and for assimilation (Kitzinger et al., 2019). Here we only examine urease in Thaumarchaeota since the cyanase gene in these archaea is unknown. At HOT, Thaumarchaeota abundance increased from depths of 100 to 500 m but then decreased again in deeper waters ( Figure 9A ). The proportion of Thaumarchaeota with ureC reached 100% at ~ 175 m but then decreased to 60%–70% at 500 to 1000 m ( Figure 9B ). However, at 4,000 m, 100% of Thaumarchaeota again had the ureC gene ( Supplementary Figure S10 ). In the ODZ stations, the abundances of Thaumarchaeota and its ureC gene were high in hypoxic areas above the ODZ ( Figure 9 ). ureC transcripts from the ETNP mirrored these findings, with abundant ureC transcripts for Thaumarchaeota in the upper oxycline but fewer in the ODZ ( Figure 11 ). However, similar to at HOT, in the 500- to 1,000-m range below the ODZs, 60%–70% of Thaumarchaeota had ureC. There is little previous data about Thaumarchaeota ureC in deep waters. Thaumarchaeota subpopulations (ecotypes) vary with depth and with coastal to offshore regions (Luo et al., 2014; Santoro et al., 2017, 2015). In the Arctic Ocean, the abundance of Thaumarchaeota ureC increased with depth from the surface to 100 m, but there was no data for deep waters (Shiozaki et al., 2021). In the more coastal Gulf of Mexico, only between 10% and 15% of Thaumarchaeota cells contain an ureC gene (Kitzinger et al., 2019), while 60%–100% of Thaumarchaeota contain ureC at SPOT (Ahlgren et al., 2017). It seems likely that only some ecotypes of Thaumarchaeota have ureC. In the open ocean datasets examined here, Thaumarchaeota is the dominant microbe with ureC at depth, but a smaller proportion of the population contains ureC in the mesopelagic.

Contrastingly, SAR11 bacteria were abundant (estimated by rpoB) throughout the water column, but the abundance of ureC from SAR11 was highest in surface waters and decreased with depth ( Figure 7 ). Urea is likely used for assimilation in SAR11. The abundance of the ureC reads gives an insight into the potential this microbial group has to use urea in oxic waters with low concentrations of other N-containing nutrients. At HOT, nitrate, an oxidized form of N, was not detected in surface waters ( Supplementary Figure S11 ; Supplementary Table S1 ). In fact, the proportion of SAR11 containing ureC negatively correlated with nitrate concentrations (log (nitrate) = − 0.053 * proportion of SAR11 with urease + 1.1487; R ² = 0.89, p-value = 1E ⁻¹⁷) ( Supplementary Figure S12 ). In contrast, in the ETNP and ETSP, the SAR11 ureC read abundances were low but tended to increase slightly (5%–10%) in the ODZ regions. The phylotypes of SAR11 ureC were completely different in oxic waters and in the ODZ, with an ODZ-specific ureC phylotype. SAR11 ureC was identified in ETNP ureC transcripts, but only at 150 m, the deepest depth sampled for transcripts. Ammonium is often undetectable (< 10 nM) in ODZ regions, but nitrate concentrations are high (~ 20 µM), so oxidized sources of N are present (Fuchsman et al., 2018; Widner et al., 2018b, 2018a). Urease appears to be less advantageous to SAR11 when oxidized sources of inorganic N are present, even though energy is needed to convert oxidized N to reduced N.

In Picocyanobacteria, we observed different proportions of cynS and ureC from different Picocyanobacteria groups occurring at different depths. For example, at HOT 275 (August 2015), on the surface, the majority of ureC reads corresponded to HL Prochlorococcus ( Figure 10C ), followed by a peak of the LLI Prochlorococcus at 100 m depth; the uncharacterized Prochlorococcus peaked at 150 m depth, whereas in the ODZs, LLIV Prochlorococcus ureC had the highest abundance. We hypothesize that uncharacterized Prochlorococcus corresponds to the uncultured NC1/LLVII ecotype of Prochlorococcus found in ITS data from these stations (Fuchsman et al., 2023). Our ureC data corresponds to the typical cascade of Prochlorococcus ecotypes with depth (Fuchsman et al., 2023).

The distribution of Picocyanobacteria cynS varied between stations and between depths ( Figure 10D ). At HOT, the abundances of picocyanobacteria cynS were lower in the top 50 m (10%–30%) but increased, reaching 70% at 150 m at HOT 272 and 100 m at HOT 275 and 278 ( Figure 10D ). The majority of these reads correspond to members of the Prochlorococcus group, but it is not possible to assign cynS to specific ecotypes ( Supplementary Figure S2 ). Some Prochlorococcus can utilize nitrate under N-limited conditions (Martiny et al., 2009; Berube et al., 2016). However, the distribution of cyanase in Prochlorococcus is not consistent with a negative correlation with nitrate. Rather, cyanase becomes the most abundant at the depth where nitrate first becomes measurable ( Supplementary Figure S11 ; Supplementary Table S1 ), which is also in the region of the DCM ( Supplementary Figure S5 ), where LLI Prochlorococcus is dominant (Fuchsman et al., 2023). While we do not have cyanate concentrations from HOT, in the North Atlantic, there are often cyanate maxima at the primary nitrite maxima, right below the DCM (Widner et al., 2016). Thus, more Prochlorococcus may be able to use cyanate at depths where it is available. The DCM region is also where eukaryotic algae have the cyanase gene ( Supplementary Figure S5 ). For the ODZ stations, a larger proportion of picocyanobacteria had cynS in the ETSP compared to the ETNP; Synechococcus was the main picocyanobacteria containing cynS in the ETSP ( Figure 10E ). This is consistent with published ITS data, which indicated that Synechococcus was more abundant in our ETSP dataset than in the ETNP (Fuchsman et al., 2023). Indeed, transcripts from the ETNP indicated that Prochlorococcus was using urea, but not cyanate, in and above the ODZ ( Figure 11 ). Thus, the ODZ Prochlorococcus are not competing with Anammox for cyanate in the ODZs. Stable isotopes of ODZ Prochlorococcus cells indicate that they primarily utilize nitrite, a partially oxidized form of inorganic N, which reaches 2–3 µM in ODZs (Aldunate et al., 2020). Thus, Picocyanobacteria and SAR11, which use urea or cyanate for assimilation rather than dissimilation, can use oxidized forms of nitrogen instead of urea and cyanate and appear to prefer these more abundant oxidized N compounds.

We identified two new groups of bacteria that use urea in low oxygen conditions: the Gammaproteobacteria Thioglobaceae and an unknown Verrucomicrobia. The Thioglobaceae clade has two distinct published subclades. A cultured strain EF1 of the SUP05 subgroup is capable of respiration with either oxygen or nitrate to oxidize sulfide to fuel carbon fixation (Shah et al., 2017). The EF1 culture can perform nitrate reduction and some steps of the denitrification pathway but cannot produce nitrogen gas (Shah et al., 2017). A cultured strain of the Arctic96BD-19 subgroup oxidized reduced sulfur, particularly thiosulfate, with oxygen but was a mixotroph, increasing growth in the presence of glucose (Marshall and Morris, 2013). In the ocean, the Arctic96BD-19 subgroup primarily lives where oxygen is present, and SUP05 primarily lives where sulfide is available (Walsh et al., 2009). Due to its ability to store elemental S, SUP05 can metabolize for some time after being taken out of more coastal sulfidic regions (Callbeck et al., 2018; Shah et al., 2019). In offshore ODZs, sulfide is not present. Instead of two clades of Thioglobacaea, by looking at the rpoB phylogenetic tree ( Supplementary Figure S7 ), we observed three subclades: Arctic96-BD19, SUP05, and a new ODZ clade of Thioglobaceae. The ODZ clade of Thioglobaceae dominated in ODZs ( Figure 6 ). While we have defined this clade using rpoB data, not genomic data, we can see that all the Thioglobaceae in the ODZ have a narG-type of nitrate reductase ( Figure 5 ). Additionally, based on the fact that the cultured Thioglobaceae clades are both autotrophic (Marshall and Morris, 2013; Shah et al., 2017), we assume that the ODZ clade is also autotrophic.

The members of the ODZ-specific subclade of Thioglobaceae appear to all be capable of using urea in the ODZ (100%), while the Arctic96BD-19 subclade only had 10%–20% of its members having ureC in oxic waters ( Figures 5 , 6 ). Thioglobaceae abundance (from rpoB) was only high in deep water (> 250 m) at HOT ( Figure 5A ). In contrast, at the ODZ stations, the abundance of Thioglobaceae was higher in the hypoxic and anoxic water columns ( Figure 5 ). The ODZ clade was the clade with the highest abundances in the ODZ core, so this clade likely contributed to the high number of ureC reads there. Thus, the ODZ Thioglobaceae may have a greater preference for using urea in the ODZ than does the oxic Thioglobaceae Arctic96-BD19 at HOT. We did not see many ureC transcripts for Thioglobaceae in the ETNP ODZ ( Figure 11 ), but the depth range of the transcripts only reached 150 m, and the maximal proportions of Thioglobaceae with ureC in the metagenomes were deeper than that (300 m) in the ODZ ( Figure 6 ).

Verrucomicrobia has the potential to use urea. At HOT, the abundance of Verrucomicrobia was low (< 2%) ( Figure 8A ), but the proportion of Verrucomicrobia with ureC was high (29.8%–150%) ( Figure 8B ). Contrastingly in the ODZ, the abundance of Verrucomicrobia examined with the rpoB gene is slightly higher (2%–4%), but the proportion of ureC is low (2%–25%) in the ODZ core ( Figure 8 ). In the hypoxic waters below the ODZs, Verrucomicrobia abundance increased to 5%–7% of the community, and all of the Verrucomicrobia contained the ureC gene ( Figure 8 ). Verrucomicrobia had low numbers of transcripts throughout the ETNP oxycline and upper ODZ, with three normalized reads at 150 m. However, activity below the ODZ, where ureC in Verrucomicrobia was the most abundant, is still unknown ( Figure 11 ). Verrucomicrobia is a phylum of bacteria whose members have many biogeochemical functionalities. Some Verrucomicrobia are heterotrophs that degrade polysaccharides and cellulose compounds and are often nitrogen fixers (Nixon et al., 2019; Delmont et al., 2022). Other members of the Verrucomicrobia phylum are methanotrophs (Schmitz et al., 2021; Howe et al., 2023). MAGs from the ETNP ODZ indicate that Verrucomicrobia had the napA type of nitrate reductase but did not have the rest of the denitrification pathway (Zhang et al., 2023). Verrucomicrobia MAGs from the ETNP were related to methane-oxidizing Verrucomicrobia of the Pedosphaerales group (Howe et al., 2023; Zhang et al., 2023), but methane monooxygenase is not present in these incomplete MAGs. Methane oxidizers are heterotrophs, and their carbon source is methane, so they need an external source of N (Sieburth et al., 1987). We cannot be sure of the metabolism of the Verrucomicrobia in our systems, but we can see that under low but not zero oxygen conditions, Verrucomicrobia were present and could utilize urea.

Though the majority of the unknown contigs on the urease (ureC) tree have now been identified ( Supplementary Figure S1 ), and we have a good understanding of which microbes are using urea in ODZs, the cyanase tree (cynS) still has an Other Bacteria group with unclear taxonomy ( Supplementary Figure S2 ). This group had significant numbers of metagenomic reads in the oxyclines above and below the ODZs and was fairly abundant in the ETNP transcripts ( Figures 2 , 11 ). Thus, there appears to be another key player still to be identified using cyanase in the oxyclines around ODZs.