To encompass the full range of diversity of the nosZ gene, all gene features and CDS (coding sequences) of annotated nosZ nucleotide sequences between 500–3,500 bp were downloaded from the National Center for Biotechnology Information (NCBI). Genbank database on June 1^st, 2021 (28,866 sequences). These nosZ gene sequences were extracted from whole genomes as well as fragments from environmental clones and metagenomes. To extract these sequences Hidden Markov Models (HMM) for nosZ sequences were generated using a custom curated subset of nosZ sequences from NCBI and used to search assembled reads of metagenomes from the Arabian Sea (AS, during late Monsoon/intermonsoon transition in September–October 2007) and Eastern Tropical North (ETNP in 2016 and 2018) and South Pacific (ETSP in 2013) using HMMR (v3.3.2) (Eddy, 1998) (Supplementary Table 1, Supplementary Figure 1). HMMR recovered 255 new sequences from these sources with the inclusion threshold of an E value of 0.01 or less and a per-domain conditional E-value of 0.01 or less.

Additional nosZ sequences from the Chesapeake Bay in 2020 and ETNP in 2018 were obtained from clone libraries of PCR products amplified using clade specific primers. PCR products (~699 bp) for Clade I were amplified using NosZ1F (5′-WCS YTG TTC MTC GAC AGC CAG-3′) and NosZR (5′-CAT GTG CAG NGC RTG GCA GAA-3′) (Henry et al., 2006) and for Clade II, a ~745 bp product using nosZ-II-F: (5′-CTI GGI CCI YTK CAY AC-3′) and nosZ-II-R36 (5′-GCI GAR CAR AAI TCB GTR C-3′) (Jones et al., 2012). This resulted in 107 and 268 novel nosZ sequences found from the Chesapeake Bay in 2020 and ETNP in 2018, respectively. Lastly, sequences from Jayakumar et al. (2018) and Bertagnolli et al. (2020) were added to the compilation. Duplicate sequences were then removed from the compilation of all the sequences listed above using SeqKit (v.2.5.1; Shen et al., 2016) by keeping only the first instance of each sequence and its accession number, resulting in a total of 11,095 nosZ sequences.

The compilation of nosZ sequences was clustered into operational taxonomic units (OTUs) using CD-HIT (Li and Godzik, 2006; Fu et al., 2012) with a defined similarity threshold of 87% sequence identity (Taroncher-Oldenburg et al., 2003; Ward et al., 2007; Bulow et al., 2008), which generated 2,282 OTUs. This threshold was chosen in order to cluster homologous genes at the level found to differentiate functional genes at the approximate level of species distinction (Purkhold et al., 2000; Taroncher-Oldenburg et al., 2003). Six hundred seventy-four OTUs contained >2 sequences. The remaining OTUs containing one or two sequences were treated as individual sequences, which are identified by their accession numbers.

The environmental source of the sequences clustered within each OTU were determined by searching NCBI and published literature. Sequences were labeled with the following environmental sources: marine sediment, marine water column, marine OMZs, hydrothermal vent systems, terrestrial, animal (i.e., from the animal microbiome), salt marsh, freshwater and estuary systems, wastewater systems, and aquatic other. Other aquatic environments mainly included hypersaline and/or soda lakes and contaminated waters. These environments were classified as other because of the great variability between them. There were not enough sequences from each kind of these other environments to make clear conclusions. OTUs that did not contain any marine or aquatic sequences or sequences from Jayakumar et al. (2018) and Bertagnolli et al. (2020) were removed from further study. The remaining OTUs containing more than three sequences were aligned with MAFFT v7.407 (Katoh and Standley, 2013), and a consensus sequence for each OTU was generated with emboss v6.6.0 (Stamatakis, 2014).

An alignment including all of the consensus sequences and the individual sequences could not be made using the nucleotide sequences. Therefore, DNA sequences were converted to protein sequences using gene prediction program GeneMarkS (v.4.28; Besemer, 2001). Protein sequences were then aligned using Clustal Omega (v.1.2.4; Sievers et al., 2011). Because most of the sequences are nosZ fragments (usually cloned PCR products), they fell on either the 5′ or 3′ end of a full nosZ gene with only a small (~42 bp) region of overlap. Consequently, sequences were further separated into 5′ and 3′ groups for further analysis (Supplementary Figure 2). In the absence of complete continuous sequences the phylogenetic coherence between the 5′ and 3′ regions cannot be ascertained; therefore, the two regions were analyzed separately. Twenty-four clusters and 127 individual sequences, which were mainly extracted from whole genomes, spanned both the 5′ and 3′ regions and were included in both groups. The longest possible aligned 5′ and 3′ region within each alignment was extracted and used for further analyses.

The extracted aligned protein regions were converted back to nucleotide sequences and analyzed with the probe finding algorithm to form archetypes (Bulow et al., 2008). These aligned nucleotide sequences for each region were analyzed using a pairwise distance matrix. The 70 bp region which had the greatest average distance between sequences was selected for further analyses. Therefore, this algorithm determines the most distinct 70 bp region within the aligned region and was used to group the sequences into archetypes, each of which is represented by a unique 70-mer region. The 70-mer archetype region with an identity threshold of 87% provides the optimal discrimination between related sequences at approximately the species level for previously studied functional genes (Ward et al., 2007; Bulow et al., 2008). A 70-mer region is long enough to represent and distinguish amongst related types, but also allowed us to incorporate shorter sequences into the analyses, therefore, allowing us to obtain a high resolution of divergence and diversity. Compared to the reference sequence of the Pseudomonas fluorescens's nosZ gene consisting of 1,920 bp (Philippot et al., 2001; AF197468.1), the 5′ 70-mer region is located around ~792–862 bp, and the 3′ 70-mer region is ~1,322–1,392 bp (Supplementary Figure 2).

The archetypes were identified as Clade I, II, or unknown based on phylogenetic inference by using NCBI's blastn and blastp (Altschul et al., 1990). Additional metadata were collected by searching with the archetype using blastn optimized for highly similar sequences on NCBI (Camacho et al., 2009). Archetype sequence information and metadata are listed in Supplementary Tables 2, 3. Phylogenetic trees were generated with the aligned 5′ and 3′ 70-mer archetype regions using RAxML v8.2.12 (Stamatakis, 2014) and viewed and edited in iTOL v6 (Letunic and Bork, 2006).

A total of 11,095 nosZ sequences were collected as described above from the databases and from other published and newly obtained clone library analysis and metagenome searches. In addition, we searched assembled reads of metagenomes from the Arabian Sea (in 2007), ETSP (in 2013), and ETNP (in 2016 and 2018), and found 146, 16, 52 and 41 new nosZ sequences respectively. Furthermore, small clone libraries from the Chesapeake Bay and ETNP (2018) added 107 and 268 nosZ sequences respectively (Supplementary Table 1). Database compilation resulted in 3,942 total unique marine or other aquatic nosZ sequences that were further used for analyses. Most of the sequences were derived from PCR amplification but several new sequences were discovered in metagenomes.

All the nosZ sequences (i.e., prior to separation into 5′ and 3′ groups) were clustered into 2,282 OTUs using CD-HIT (Supplementary Figure 3). However, only 592 clusters contained sequences from marine or other aquatic environmental sources and were included in further analyses. Of those, 179 OTUs contained more than 2 sequences. These 179 OTUs represent only ~30% of the clusters, but account for ~87% of the sequences (3,442 sequences). The largest cluster contains 643 sequences (~16% of the sequences). The remaining OTUs, the unique individual sequences, represent ~70% of the clusters and ~13% of the sequences (500 sequences). This distribution of sequences shows that while the majority of the clusters are unique, most of the sequences group into larger OTUs (Figure 1).

A total of 866 sequences aligned at the 5′ region of the nosZ gene and were grouped into 250 archetypes, which were further examined. One hundred eighty-one archetypes represented 252 unique individual sequences only and therefore, did not contain clusters. However, the remaining 69 archetypes contained the majority of the sequences (614 sequences). The largest archetype represented 75 sequences. Archetypes were classified into Clade I (169 archetypes, 67.6%) and Clade II (78, 31.2%). Clade I archetypes represent 681 sequences, while Clade II represent 182 sequences. The clade affiliation of three archetypes could not be determined so they were classified as unknown. The unknown sequences were all single sequences obtained from marine OMZs. The clades therefore show very clear phylogenetic separation (Figure 2).

The majority of archetypes contain sequences obtained from only one of the environmental sources (205 archetypes, 82%); no archetype contains sequences representing more than three different sources (36 represent two environmental sources and 8 represent three different environmental sources). Figure 2 depicts the phylogeny of the nosZ sequences of the 5′ region of the gene and shows the spread of environmental sources across archetypes and clades. The environmental sources for each archetype (based on the 70-mer archetype analysis) were identified, separated by clade, and the total number of archetypes representing an environmental source within each clade was determined (Figure 3A). Archetypes within Clade I were predominately obtained from marine sediment (62.7%) while archetypes of Clade II sequences derive from a much more diverse range of environmental sources, although the marine water column represents the largest representative environmental source (27.6%). Environmental sources of terrestrial, animal, wastewater systems, marine water column, hydrothermal vent systems, freshwater & estuary systems are found more often in Clade II nosZ genes.

Archetypes containing sequences retrieved from marine OMZs were examined further to resolve if there was smaller scale regional biogeography by comparing their occurrence within and among the three major OMZs (Table 1). The majority (76%) of archetypes containing OMZ sequences contained sequences from only one of the OMZ regions (i.e., either AS, ETNP, or ETSP). Twenty-four percent of the archetypes contained sequences from both the AS and ETNP, and no archetypes contained sequences from all three major OMZs. Therefore, most of the OMZ sequences appear to be unique to a specific region.

The taxonomy of each archetype's representative sequence was determined through phylogenetic inference (Figure 4A). The phylum/class of both Clade I and II could be determined for a majority of the sequences, 70.4%/66.9% and 89.6%/83.1% respectively. However, only 20.7% of Clade I sequences could be assigned to the genus/family level and 64.9% for Clade II. Clade I primarily consisted of sequences from the Pseudomonadota phylum (65.7%), with most classified as Alphaproteobacteria (36.1%). On the other hand, Clade II was dominated by Bacteroidota sequences (41.6%). In addition, Clade II was more diverse, containing sequences from 16 different phyla, while Clade I only had five phyla, three of which were found in both Clades.

The two dominant phyla, Bacteroidata and Pseudomonadota, were examined further by determining what percent of archetypes came from each environmental source within each phylum (Figure 5A). Bacteroidota, which was represented in of 33 archetypes (59 sequences), was primarily obtained from the marine water column (~48.5%), followed by freshwater and estuarine systems (21.2%). On the other hand, Pseudomonoadota (125 archetypes representing 447 sequences) was mostly retrieved from the marine sediment (53.6%). The Pseudomonoadota classes were further analyzed. Alpha-, Gamma-, and Deltaproteobacteria were all predominately obtained from marine sediment (71.9%, 31.9%, and 71.4% respectively). Betaproteobacteria's main environmental source was terrestrial; however, it only comprises 7 archetypes and 18 sequences.

The majority of the nosZ sequences aligned on the 3′ region of the gene. A total of 3,190 nosZ sequences grouped into 280 archetypes, which were further analyzed. Two hundred and twenty-four individual sequences made up 189 archetypes while most of the sequences (2,966 sequences) were clustered into 171 archetypes. The largest archetype represented 849 sequences. Archetypes divided almost equally between Clade I and II, 47.9% (134) and 51.8% (145) respectively, with only one archetype in which the clade affiliation could not be determined. Despite a nearly equal divide between Clade I and II archetypes, there are significantly more sequences represented within the Clade I archetypes vs. Clade II, 2,758 and 431 sequences respectively. Similar to the 5′ region of the nosZ gene, the clades separate out phylogenetically (Figure 6).

The environmental source for each archetype was determined. Like the 5′ region of the nosZ gene, most of the archetypes from the 3′ region (210 archetypes, 75%) represent only one environmental source. However, there are two archetypes represented by sequences from five different environmental sources (Labrenzia sp. and uncultured gammaproteobacterial), and one archetype representing six different environmental sources (uncultured alphaproteobacterial). The remainder of the archetypes represent three or fewer different environmental sources. The phylogenetic tree of the representative nosZ archetype sequences in the 3′ region of the gene in Figure 6 portrays the spread of the environmental sources across the clades. Figure 3B depicts spread of environmental sources for each archetype's representative sequence. Similar to the 5′ region of the nosZ gene, archetypes within Clade I were predominately obtained from marine sediment (38.8%), followed by salt marshes (21.6%). On the other hand, Clade II archetypes have a wider distribution across environmental sources. Among the Clade II archetypes the terrestrial and marine OMZ environmental sources were equally represented (21.38%) followed by freshwater and estuarine systems (17.2%). Environmental sources of terrestrial, marine OMZs, animal, wastewater systems, marine water column, hydrothermal vent systems, wastewater systems, and freshwater & estuary systems were found more often in Clade II nosZ genes.

Biogeography of the nosZ sequence was examined at a smaller scale by focusing on the archetypes that contained sequences retrieved from marine OMZs. Within each of these archetypes, the source of the sequences was determined, i.e., which of the major OMZs (AS, ETNP, or ETSP) (Table 1). Similar to the 5′ region, archetypes containing OMZ sequences from the 3′ region were predominately obtained (59.2%) from only one of the major OMZs. Fewer archetypes were retrieved from two distinct OMZs (36.7%). Only 4.1% archetypes contained sequences from all three of the major OMZs. This suggests the majority of nosZ sequences are unique to their environment.

Like the 5′ region of the nosZ gene, the majority of the sequences could be classified into phylum/class for Clade I (60.4%/58.2%) and Clade II (63.2%/54.2%) (Figure 4B). Additionally, only 23.1% of Clade I and 32.6% of Clade II sequences could be assigned to the genus/family. Like the 5′ region, most of the Clade I sequences are within the Pseudomonadota phylum (59.7%) with Alpha- and Gamma-proteobacteria dominating (26.1% and 23.9%). Clade II was again predominately Bacteroidota sequences (27.1%). Clade II had greater diversity consisting of sequences from 18 different phyla compared to sequences from Clade I which were found within only 2 different phyla.

As noted above, Bacteroidota and Pseudomonadota were the largest phyla containing nosZ sequences, as was also observed from the 5′ region. Within each phylum, the percent of archetypes derived from each environmental source was determined (Figure 5B). Similar to the 5′ region, Bacteroidota's (39 archetypes representing 61 sequences) main environmental source was the marine water column (30.8%); however, this is followed by marine sediment (25.6%) instead of freshwater and estuarine systems. As with the 5′ region, archetypes within the phylum Pseudomonadota (92 archetypes representing 2,371 sequences) were mostly obtained from the marine sediment (30.8%). Alpha- and Gammaproteobacteria were both predominately retrieved from marine sediment (30.7% and 32.4% respectively). Betaproteobacteria was mainly obtained from freshwater and estuarine systems; however, it was the smallest Pseudomonadota class in terms of both number of archetypes and sequences (13 and 78 respectively).

Although the nosZ gene has been identified and studied for decades, knowledge of the function, distribution, and taxonomy of nosZ containing organisms in the environment is still fragmented (Hallin et al., 2018; Bertagnolli et al., 2020). This is in part due to the high taxonomic diversity of nosZ-containing organisms revealed through phylogenetic- and high throughput sequencing analysis of nosZ genes in natural communities (Jones et al., 2012; Sanford et al., 2012; Orellana et al., 2014; Hallin et al., 2018). The nosZ sequences analyzed in this study were so distinct at the nucleotide level, they could not be well aligned, further highlighting the nosZ diversity. However, the sequences were similar enough at the protein level such that they could be aligned, suggesting despite the great DNA level diversity, the genes encode fundamentally the same enzyme. The high level of similarity at the protein level implies strong functional constraints on the enzymes, which probably extends to ecologically significant features such as oxygen sensitivity and the requirement for copper at the active sites. While the gene sequence has evolved consistently with the host genome, the functionality of the protein has been conserved.

This great diversity of nosZ at the DNA level introduced detection and analytical challenges, which suggests that the extent of the N2O-reducing microbial community is still underestimated (Chee-Sanford et al., 2020). Over 600 new and widely divergent nosZ sequences were found in our searches of marine metagenomes and even in small clone libraries from estuarine and oceanic environments (Supplementary Table 1). It is clear that new studies are continuously adding to the nosZ diversity, and likely much diversity is still undiscovered. This great diversity is evident in Figure 1, while the majority of sequences group together into larger clusters, most of the clusters, ~70%, are unique. Even after further grouping into archetypes, a similar pattern is evident: the majority of the archetypes represent individual sequences for both the 5′ and 3′ region of the nosZ gene (72% and 68% respectively), meaning most of the branches on the phylogenetic trees (Figures 2, 6) are unique. While this study uncovered an extensive diversity of nosZ sequences, this rate of discovery implies that additional diversity not covered in this study is yet to be revealed.

NosZ containing organisms are found inhabiting diverse environments and cover a wide range of genera. The prevalence of these organisms has the potential to regulate global N₂O emissions as nosZ is the only known biological sink of N₂O. Although nosZ was first characterized in 1997 (Zumft), much remains unknown about the distribution, taxonomy, and biogeochemical impacts in the environment. Great diversity of the nosZ gene was found at the nucleotide level (such that sequences could not be well aligned); however, nosZ was highly conserved at the protein level and functionally remains the same despite evolving in diverse environments. This apparent functional constraint on the enzyme implies further ecological constraints, such as oxygen sensitivity or copper requirements; however, this has yet to be proven and requires further investigation. In addition, the extensive nucleotide sequence diversity found introduces detection and analytical challenges. Universal PCR primers for nosZ are not possible and, in turn, it is likely that nosZ diversity and environmental abundance are greatly underestimated. The ecological significance of this underestimation should be further examined. Instead of PCR primers, applying HMM to metagenomic studies is likely a better approach to identifying and characterizing new nosZ sequences, as done in this study.

We found Clade II to have a wider taxonomic range and distribution across environmental sources compared to Clade I. Differences in community composition, environmental distribution, and/or kinetics between Clade I and II containing organisms are likely to impact environmental N₂O fluxes causing various ecological and environmental consequences. However, the extent of this is beyond the scope of this study, and field studies are required to understand the biogeochemical impacts of the two clades. The spread of the diverse taxa across environments has implications for N₂O consumption rates.

This study examined the environmental distribution of the nosZ gene with a focus on marine and other aquatic environments to explore the breadth of environmental sources across both clades. In both 5′ and 3′ regions the majority of archetypes contain sequences obtained from only one environmental source. This suggests nosZ sequences, regardless of clade, is typically unique to an environment. Environmental factors such as oxygen levels, redox, nutrient availability, temperature, pH, etc., are likely involved in creating this observed separation between environmental sources. Clade I was mostly retrieved from marine sediments while Clade II sequences were derived from a much more diverse range of environmental sources including terrestrial, animal, wastewater systems, marine water column, hydrothermal vent systems, freshwater, and estuary systems. We found Clade I sequences primarily consisted of sequences from the Pseudomonadota phylum (predominately alphaproteobacteria) and Clade II was dominated by Bacteroidota sequences. These Bacteroidota sequences were mainly obtained from the marine water column while the Pseudomonadota were predominately obtained from marine sediments. The mechanism for the environmental distribution of nosZ, such as N₂O affinity, nutrient availability, oxygen concentration, redox state, etc., is not well understood and should be examined further as it could explain variations in N₂O fluxes. NosZ has the potential to regulate N₂O emissions that contribute to ozone depletion and the greenhouse effect; this study characterized the diversity and environmental distribution of this gene to help further understand its biogeochemical impacts.