the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 25 Mar 2026
Global pattern of nitrogen metabolism in marine prokaryotes
Abstract. The ocean nitrogen cycle is driven by an ensemble of metabolic processes sustaining marine ecosystems and ocean productivity. However, the spatial distribution and environmental drivers of its major pathways, i.e., nitrogen fixation, denitrification, assimilatory and dissimilatory nitrate reduction to ammonium (ANRA, DNRA), and nitrification are not well known. Furthermore, the taxonomic composition of the prokaryotes supporting each pathway remain incompletely understood. Leveraging newly assembled global marine metagenomic datasets and a state-of-the-art machine learning framework, we inferred the global biogeography of the genomic potential for key metabolic pathways of the marine nitrogen cycle. This was achieved using a multi-output regression of gene read counts against environmental climatologies. Our results reveal distinct biogeographic patterns of genomic potential: anaerobic or light-inhibited pathways are enriched in high-latitude regions, eastern boundary upwelling systems, and deeper ocean layers, while nitrogen fixation and ANRA dominate in oligotrophic gyres. These patterns are consistent with known metabolic strategies, model-based estimates, and underlying taxonomy. Indeed, we identify that Cyanobacteria associate primarily with aerobic, biosynthetic pathways, while Gammaproteobacteria and Nitrososphaeria encode for nitrogen transformations related to energy requirements. By coupling microbial community composition with genome-level information, our approach advances understanding of the microbial foundations of nitrogen transformation pathways and offers new insights on underrepresented processes into biogeochemical models. We highlight the growing value of omic data to better understand marine ecosystem function in relation to environmental gradients and community composition, and their use as a potential observation-based alternative or complement to biogeochemical models.
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-1459', Anonymous Referee #1, 25 Apr 2026
-
AC1: 'Reply on RC1', Alexandre Schickele, 24 Jul 2026
Authors response to RC1
We thank the editor and the two reviewers for their feedback on the manuscript. Here we provide our answer to RC1. We numbered all comments and provided answers below and in the attached pdf file. Lines are numbered according to the new manuscript with track change disabled.
Reviewer 1
General comments:
The authors examine global patterns in key nitrogen cycling pathways in the ocean using genetic datasets. The tool used -- the statistical model CEPHALOPOD -- has been previously published and provides a thorough and statistically rigorous analysis of how the key genes vary with environmental factors. This is a step up from simple DIY correlations between datasets. Most of the results themselves are not novel -- we know that assimilation of N by cyanobacteria should dominate in oligotrophic gyres, and nitrification at depth, for example, but it is still useful to see how these known patterns emerge using this tool. It is very interesting to see the widespread but low abundance coverage of the denitrification genes. However, the conclusion of high relative abundance of DNRA across the ocean, especially at high latitudes in the surface, needs more work to be convincing. Not only is it very surprising (why would an anaerobic process show up in the most oxygenated waters at high latitudes?), but the genetic marker may not be appropriate here. I am not a geneticist or omics expert, but my understanding is that nir B and nir D can be involved in ANRA as well as DNRA (see for example this review: https://doi.org/10.1016/j.synbio.2025.12.017), and ANRA is actually what we expect in high latitude surface waters. I suggest the authors look for nrfA instead and reevaluate their conclusions about DNRA.
Thank you for this thoughtful overview of our work and constructive suggestions that have substantially improved the manuscript. All the specific points raised in your general comment are addressed in detail in the responses below, including the reevaluation of DNRA by considering the new marker genes you suggested.
Specific scientific comments:
1. Use of "aerobic" (ex: L. 17): I would not use the term aerobic for autotrophic processes/photosynthesis. I realize the ANRA may be associated with aerobic heterotrophy as well in the surface, which may explain why a small section of the pie in Fig. 6c for ANRA is not cyanobacteria. I and many others consider aerobic to mean O2-consuming, not O2-producing, so may be confusing.
Thank you for this comment. We agree that the term “aerobic” may be misleading in this context because it is often interpreted as referring to oxygen-consuming metabolisms. Our intention was instead to describe metabolisms associated with oxic environments, including processes that do not consume oxygen, such as photosynthesis. To avoid confusion, we have replaced “aerobic metabolism” with “metabolism associated with oxic environments” where appropriate.
L20-22: Indeed, we identify that Cyanobacteria primarily encodes for metabolic pathways associated with oxic environment and biosynthetic pathways, while Gammaproteobacteria and Nitrososphaeria encode for nitrogen transformations related to energy requirements.
L439-443: In the epipelagic layer, nitrogen cycling metabolic pathways associated with cellular biosynthesis and oxic conditions (i.e., ANRA, Nitrogen fixation) are primarily supported by Cyanobacteria (relative contribution between 0.6 and 0.9), in a lesser extent by non-cyanobacterial groups such as Gammaproteobacteria (relative contribution between 0.1 and 0.4; Fig. 6c).
2. Genes for nitrification: Again, I am not an expert, but why not look for amoA and nxr genes for ammonia and nitrite oxidation? This is what is used for this database on nitrifiers: Tang et al "Database of nitrification and nitrifiers in the global ocean," https://doi.org/10.5194/essd-15-5039-2023 -- which should be cited and with which results here could be compared.
Thank you for these suggestions regarding nitrification, following Tang et al. (2023). In our study, however, the amoA gene (K28504) was not available in the metagenomic annotation dataset used for our analyses (Paoli et al. 2022). In addition, the nxrA gene (K00371) is associated not only with nitrification but also with denitrification and DNRA (according to the KEGG database), making it unsuitable as a pathway-specific marker under our selection criteria. We therefore excluded nxrA to avoid ambiguity in pathway attribution. Consequently, nitrification was represented using the genes pmoB (K10945), pmoC (K10946), and hao (K10535; Fig. 2). We also removed pmoA (K10944) from the list of nitrification markers, as it is primarily associated with methane oxidation and was incorrectly assigned to nitrification in the previous version of the manuscript. We have now cited the nitrifier database of Tang et al. (2023) in the Methods and Discussion to place our results within the context of previous work on the global distribution of marine nitrifiers.
L200-204: For nitrification, although suggested by Tang et al. (2023), the amoA gene was not available in the metagenomic annotation dataset (Paoli et al. 2022), whereas nxrA was excluded because it is also associated with denitrification and DNRA (Kanehisa and Goto, 2000). Consequently, nitrification was represented by the pathway-specific marker genes pmoB (K10945), pmoC (K10946), and hao (K10535; Fig. 2b).
L488-493: The genomic potential, i.e., the relative abundance of those genes negatively correlates to PAR and positively correlates to nitrate and oxygen concentration in the epi- and mesopelagic layers (Fig. 6b), known drivers of nitrification rates (Tang et al., 2023; Ward, 2008; Zhang et al., 2020). Moreover, the estimated genomic potential distribution, showing elevated nitrification gene abundance in high-nutrient, low-chlorophyll (HNLC) regions such as the North Pacific and Atlantic, and Southern Ocean, aligns with higher measured nitrification rates (i.e., ammonia and nitrite oxidation; Tang et al., 2023).
L494-496:Therefore, the genomic potential of nitrification aligns with known nitrifiers distribution and nitrification rates (Pajares and Ramos, 2019; Tang et al., 2023).
3. Nitrification taxonomy: It is reassuring that the taxonomy associated with nitrification is Nitros*, but it is perhaps a bit circular to pitch it as a conclusion (considering the literature) since for nitrification, taxonomy and function are so related. It would make the paper stronger overall to pitch this more as a "validation" that the model can find known patterns such as this one.
Indeed, this result is not an original finding per se, but rather a validation that the taxonomy associated with our genomic potential matches current knowledge. We now state it as follow:
L486-488: Under oxic conditions, Nitrification encoding genes are predominantly associated to Nitrososphaeria archaea, which corresponds to established knowledge in the literature (Tang et al., 2023; Ward, 2008), thus validating taxonomic annotation of our selected marker genes.
4. Fig. 1: There are many wide black arrows, but then you say you did not consider them in this study. Perhaps there is another way to say this that is not so disserving? You are using methods that just identify correlations among datapoints, so probably OK to just not even mention that you cannot consider transport in the model. Same for anammox: you DID consider it, but it was just not sufficiently abundant to include in the statistical analysis.
Thank you for this comment. Our objective is to model the genomic potential associated with prokaryote genes and their environmental niches. Physical transport processes (Fig. 1, black arrows) are therefore outside the scope of the framework, as they are driven by ocean circulation and density gradients rather than by microbial genes. To clarify this point, we added the following statement to the Figure 1 caption:
L76-77: Note that physical transport processes are represented for completeness of the nitrogen cycle but were not modeled, as our study focuses on the environmental niches and the genomic potential of microbial communities.
Regarding anammox, we did evaluate this pathway but were unable to generate robust global projections because the available marker genes lacked pathway specificity and sufficient coverage in our dataset. To clarify this point, we added the following statement at the end of the Introduction, with a reference to the Methods section for further details:
L162-164: Although anammox was initially considered, it was not retained for global projections because the available marker genes lacked sufficient pathway specificity and coverage to support robust habitat modeling (see Sect. 2.1.2.).
We also clarified this point in the Methods section:
L204-208: Among the three enzymes associated with the Anammox pathway (Fig. 1, pathway 6) and detected in our metagenomic dataset (Paoli et al., 2022), the nitrite reductase NO-forming enzymes K00368 and K15864 did not meet our specificity criteria because they are also associated with denitrification, whereas hydrazine synthase (K20934) was detected only once in the eastern pacific tropical upwelling system (-84.61°E, -1.89°N). Consequently, the Anammox pathway (Fig. 1, pathway 6) was not retained as candidate for global-scale modeling.
5. L. 120: "Habitat modeling" seems like it should represent a much broader category of modeling, rather than just simply this statistical approach. I would assume that any ecosystem model that aims to identify the ecological niche of a population or functional type is a habitat model.
Thank you for this comment. We agree that habitat modeling encompasses a broad range of applications and is not limited to a specific statistical methodology. This is already stated in the manuscript as follows:
L135–136: Rooted in ecological niche theory, it is widely used to predict the spatial distribution of organisms based on the environmental conditions at observed locations (Peterson and Soberón, 2012).
Our approach falls within this broader framework because it uses environmental conditions to characterize the niches of microbial taxa and the genomic potential associated with the relative abundance of nitrogen-cycle marker genes. To clarify this point, we have added the following statement at the end of the Introduction in the revised manuscript:
L156-157: To address these questions, we apply the core principles of habitat modelling to model the environmental niches of nitrogen-cycle markers genes in prokaryotes, and to project the global distribution of their genomic potential.
6. Cutoff (50 observations; line 169): You say that a gene must be present in at least 50 observations "to ensure sufficient spatial coverage" -- this does seem incompatible with many of the anaerobic N cycling transformations in anoxic zones unless you have 50 samples IN anoxic zones. This meant that you couldn't show anammox patterns. Perhaps for these types that didn't make that cut (anammox, and perhaps nrf as a replacement for nirBD for DNRA), it's useful to explain that they do "show up where they should," i.e., the small volumes of anoxic zones, rather than saying "we didn't consider them." You did consider them, but they are just not candidates for global-scale representation which is reasonable!
Thank you for this comment, which is closely related to your comment #4. We agree that the previous wording may have suggested that anammox was not considered, whereas it was in fact evaluated during marker-gene selection. We have revised the text to clarify that the pathway was excluded from global-scale modeling because the only anammox-specific marker gene lacked sufficient spatial coverage, despite being detected in a location where anammox activity is expected.
L204-208: Among the three enzymes associated with the Anammox pathway (Fig. 1, pathway 6) and detected in our metagenomic dataset (Paoli et al., 2022), the nitrite reductase NO-forming enzymes K00368 and K15864 did not meet our specificity criteria because they are also associated with denitrification, whereas hydrazine synthase (K20934) was detected only once in the eastern pacific tropical upwelling system (-84.61°E, -1.89°N). Consequently, the Anammox pathway (Fig. 1, pathway 6) was not retained as candidate for global-scale modeling.
7. L. 225: "outliers in the biological inputs were removed." Is this reasonable, given that we do think there is much heterogeneity in microbial communities and BGC fluxes in the environment? What are the implications of this? Is it more like you are finding "average" patterns and ignoring the known deviations? (This seems reasonable to do.)
Thank you for this comment. We agree that microbial communities and biogeochemical fluxes exhibit substantial heterogeneity in space and time. However, the definition of outliers depends on the scale of analysis. In our case, the CEPHALOPOD framework operates at a 1° x 1° monthly climatological scale, with the objective of identifying large-scale climatological relationships rather than resolving sub-mesoscale or short-term variability. We have clarified this point in the manuscript as follows:
L253-255: This step also aims to improve the robustness of model fitting to large-scale patterns and to reduce the influence of extreme values that may reflect local or transient conditions not representative of the climatological scale at which CEPHALOPOD is applied.
8. DNRA results: As I state above, I'm deeply skeptical that DNRA is actually this abundant in oxygenated surface waters. It would be really amazing if these results WERE true, and so, you need to convince me and the other readers that I can believe this. See above major comment about nrf vs. nirBD.
Thank you for this major comment. Although nirB and nirD are annotated as DNRA-associated genes in KEGG, this is contested in the literature. In this context, and in agreement with your suggestion and reviewer #2 (RC2), we revised our definition of the DNRA pathway by replacing the previously used nirB and nirD markers with nrfA and nrfH. These are more widely recognized in the literature as specific markers of DNRA. All analyses, figures, and results were updated accordingly. With this revised and more specific definition, DNRA-related marker genes are predominantly detected in oxygen-depleted regions, particularly in the eastern tropical Pacific upwelling system, consistent with the current literature describing this pathway. The relevant changes are detailed below:
L328-331: We observe a clear dominance of nitrogen transformation pathways associated with bioavailable over bio-unavailable forms of nitrogen – namely ANRA, and nitrification – in both epi- and mesopelagic layers (Fig. 3). In the epipelagic layer, ANRA accounts for the largest share of the genomic potential (0.79), followed by nitrification (0.17). However, DNRA that is also associated with bioavailable forms of nitrogen only accounts for a minor share of the genomic potential (<0.01).
L356-360: DNRA, typically associated with heterotrophic activity in suboxic to anoxic environments, shows peak genomic potential in the tropical Pacific upwelling system, with values ranging from 0.015 in the epipelagic to 0.03 in the mesopelagic layers (Fig 4c). Additional hotspots are predicted in ODZ of the northern Indian Ocean and in the Benguela upwelling system, with values ranging from 0.01 in the epipalgic to 0.02 in the mesopelagic layers (Fig. 4c). These are associated with a high uncertainty however (NSD > 0.5 ; see Fig. S3 for detailed uncertainty pattern).
L364-373: see Figure in the supplement
Figure 4: Genomic potential for bioavailable nitrogen transformation pathways represented by their relative abundance of metagenomic reads. Nitrogen transformation pathways are displayed in rows with: (a) Assimilatory Nitrate Reduction to Ammonium (ANRA), (b) Nitrification and (c) Dissimilatory Nitrate Reduction to Ammonium (DNRA). Left panels correspond to the epipelagic projections and right panel correspond to the mesopelagic projections. Each panel represents the annual average spatial projection, with white stippling representing geographical areas associated with a high degree of uncertainty (i.e., standard deviation above 0.5 the global average). Latitudinal profiles are displayed in form of annual average (black), summer (red; northern hemisphere June to August) and winter (blue; northern hemisphere December to February). The projection uncertainty across bootstraps is presented as grey shading in the latitudinal profiles.
L446-448: Taxa that are obligatory or facultative anaerobes also show contributions to epipelagic DNRA (0.1 to 0.2 for each taxon) while mesopelagic DNRA is mainly associated with Marinisomatia in our dataset (0.8).
L597-502: The genomic potential for anaerobic nitrogen transformation pathways, such as denitrification and DNRA, positively correlates with low oxygen saturation and high apparent oxygen utilization (Fig. 6b). This is consistent with the dependence of Gammaproteobacteria on sub- or anoxic conditions to perform anaerobic respiration (Lam and Kuypers, 2011). The eastern Pacific upwelling, Arabian Sea, and Bay of Bengal ODZs show the highest genomic potential for DNRA and denitrification, reflecting known hotspots of anaerobic nitrogen transformation (Zhang et al., 2023).
9. L. 497: Levine et al 2025 is a review, not a model. Cite the models you are referring to directly.
Thank you for this comment. We clarified this by referring to the Régimbeau et al. (2022), Giordano et al. (2024) and Reynolds et al. (2024) studies directly. The statement is now as follows:
L526-530: Moreover, integrating metagenomics with genome-scale metabolic models (Giordano et al., 2024; Régimbeau et al., 2022; Reynolds et al., 2024), offers a promising avenue to link gene content with microbial metabolic fluxes (Levine et al., 2025). These metabolic models leverage genomic data to predict microbial nutrient transformations under varying environmental conditions, thereby enhancing the predictive capacity of biogeochemical cycles (Giordano et al., 2024; Régimbeau et al., 2022; Reynolds et al., 2024)
Smaller presentation comments:
10. L. 85: "ammonium that is being released" -- note that this may also be considered a key N transformation process! Perhaps also mention urea somewhere, as urea cycling is becoming more recognized currently.
Thank you for this comment. We agree that urea cycling is an important component of the marine nitrogen cycle and contributes to regenerated nitrogen, notably through the production of ammonium during urea degradation. This process may be included in the statement in L85 (“ammonium that is being released”). In the present study, however, we focus on redox-associated inorganic nitrogen transformation pathways that can be consistently resolved using available marker genes. In contrast, urea cycling belongs to the broader organic nitrogen cycle and is interconnected with multiple metabolic processes, including amino acid metabolism. For this reason, it was not explicitly included in the scope of our analysis. We clarified this statement and mention the urea cycle as follows:
L86-91: But nitrification thrives in the layers just below that, in the depth range between 100 and 500 m in the mesopelagic zone, where it oxidizes the ammonium released during the remineralization of the organic nitrogen that is sinking or transported downwards from the near-surface layers. These ammonium-releasing processes, that may include urea degradation among other pathways of organic nitrogen cycling (e.g., Wan et al., 2024), are not explicitly represented in this work, which focuses on redox-associated inorganic nitrogen transformations only.
11. L. 107: do you mean "such genes" instead of "such enzymes"?
We apologize for the confusion. While genes code for enzymes, we are referring to the genes. The statement has been clarified, also following reviewer #2 comments 4, 9 and 10 in our answer to RC2.
L115-118: The relative abundance of genes encoding a given metabolic function in prokaryotic metagenomes can be estimated from genome-wide metagenomic read coverage. Following Schickele et al. (2024), we refer to this relative abundance (e.g. of nifH,D,K) as the genomic potential of a prokaryotic community for that function (e.g., nitrogen fixation).
12. L. 113: "the one[s] associated" (not "the one")
Thank you for pointing this out. We corrected the text accordingly.
13. L. 145: "Doing so, [we] explore"
Thank you for pointing this out. We added the missing word [we]
14. L. 272: "The three features" -- could you name them here? I am curious at this stage!
Thank you for this suggestion. We added the mention of the features with reference to the corresponding result section:
L305-307: The three features explaining most of the variance in the observation presented a cumulative feature importance (> 50 %) of 75 % for the epipelagic model (POC, nitrate and oxygen concentrations) and 81 % for the mesopelagic model (POC, apparent oxygen utilization and PAR; see Sect 3.5.).
15. L. 302: Could you clarify/specify that "bio-unavailable forms" refers to N fixation and denitrification? This makes sense to look at these two separately from the others, but just took me a second to realize what was going on.
Thank you for your comment. We clarified this point as follows:
L337-340: This also means that the first order spatial gradient does not reflect whether nitrogen transformations are associated with cellular biosynthesis or energy metabolism, but rather whether they involve the bio-available nitrogen pool (ANRA, DNRA, and nitrification) or the exchange between the bio-unavailable nitrogen pool and atmospheric (N₂) nitrogen (nitrogen fixation and denitrification).
16. Figs. 4 and 5: Could you just label each of the maps with the process for much easier reference?
Thank you for this suggestion. We added labels referencing the process above each global projection for clarity.
17. L. 408: Again, I would not call cyanobacteria "aerobic."
Following our response to your comment #1, we do not refer to cyanobacteria as aerobic in the revised manuscript.
18. L. 410/11: Nitros* are named the way they are because they are associated with nitrification, so seems naive to state a result this way (see above comment).
Thank you for this comment. Following our response to your comment #3, we agree that this result is not an original finding per se, but rather a validation that the taxonomy associated with our genomic potential matches current knowledge. We state it as follow:
L486-488: Under oxic conditions, Nitrification encoding genes are predominantly associated to Nitrososphaeria archaea, which corresponds to established knowledge in the literature (Tang et al., 2023; Ward, 2008), thus validating taxonomic annotation of our selected marker genes.
19. L. 415: "In the lat[t]er"?
Thank you for pointing this out, however the corresponding statement has been removed in the present in the revised manuscript.
20. L. 416: Is it known that these are diazotrophs? Perhaps more fair to just say non-cyanobacterial groups?
Thank you for this comment. Non-cyanobacterial diazotrophs are reported in the literature, but we did not specifically test if the taxa reported in Fig. 6c are actively expressing diazotrophy. While the corresponding statement has been removed in the present in the revised manuscript, we applied your suggestion to a similar statement above:
L439-442: In the epipelagic layer, nitrogen cycling metabolic pathways associated with cellular biosynthesis and oxic conditions (i.e., ANRA, Nitrogen fixation) are primarily supported by Cyanobacteria (relative contribution between 0.6 and 0.9), in a lesser extent by non-cyanobacterial groups such as Gammaproteobacteria (relative contribution between 0.1 and 0.4; Fig. 6c).
-
AC1: 'Reply on RC1', Alexandre Schickele, 24 Jul 2026
-
RC2: 'Comment on egusphere-2026-1459', Anonymous Referee #2, 06 May 2026
General comments:
This work uses correlations between global genomic datasets and biogeochemical environmental data to validate and extrapolate on trends in prokaryotic nitrogen cycling predicted by the CEPHALOPOD statistical model. This is a proof of concept for using metagenomics to predict what the authors refer to as “genomic potential”, or the possibility that cells are expressing the proteins encoded by the genes sequenced in global metagenomic efforts. This work benefits from validation via comparison with known nitrogen cycling trends, and goes on to use the tool for prediction of unknown prokaryotic trends. This is most effective as a source of validation, where the CEPHALOPOD model successfully predicts known trends like increased assimilation of N in light, oxic oligotrophic equatorial waters. However, this doesn’t hold up when predicting N cycling in novel areas. There are base assumptions fed into the model that I think have produced spurious results. Like the first reviewer noted, I have doubts that DNRA and other anaerobic processes are comparatively more prevalent in well oxygenated upwelling systems and polar regions. These analyses hinge on the identification of specific genes to serve as a proxy for processes, which is a difficult thing to do even with metatranscriptomics. I think this paper would be served by a closer look at specific N gene function in eukaryotes, and comparing them to the known geography of prevalent N cycling processes, to use as further validation of the CEPHALOPOD model.
Specific scientific comments:
Light inhibition is not necessarily correlated to anaerobic conditions, especially in high-latitude regions. Those areas are often characterized as oxygen rich, high nutrient, low chlorophyll. This would make oxygen available as an electron receptor during photosynthesis, and eliminate the need for DNRA and other redox reactions. Deeper ocean layers are often oxygen rich while being nutrient poor, depending on depth below the euphotic zone, where oxygen tends to rebound as there are fewer organisms to consume it.
I would consider adjusting your depth-based strata based on region. In oligotrophic gyres, the chlorophyll maximum will often occur at 150m, meaning that the region above that will be well lit and nutrient depleted. Coastal, upwelling regions can have a chlorophyll maximum of 5-10m. Adjusting this will make your environmental data averages more accurate.
Has anyone shown that “genomic potential” has predictive power in prokaryotes? I don’t work with prokaryotes, so I’m not sure if there are studies tying metagenomics to metatranscriptomics on this scale in these organisms. If so, I would prominently include them to bolster your case for DNA as a proxy for likely expression.
I am slightly confused about the choice of nirB and nirD as a proxy for DNRA activity. As far as I know, nirB is also used for nitrate assimilation (ANRA), and some organisms use it concurrently with other forms of nitrite reductase in order to bolster N assimilation when “fresh” N is available. Its use alongside nirD, in the form sometimes referred to as nirBD, is documented as capable of DNRA alongside ANRA. The ratio of nirBD to nirA could be potentially used as evidence of DNRA activity, but I would include evidence that nirBD alone can be used as a DNRA proxy if you are convinced. In prokaryotes, nrfA is often used as the diagnostic gene for DNRA.
L58 -- generally, can you clarify your rates and put context to Tg N/year? It would be helpful to be able meaningfully compare these parts of the nitrogen cycle.
L94 -- denitrification has been shown to happen at some level in all oceans
L106 -- can you be more clear about what a “taxonomic unit” is? A whole genome? A gene? What is a metagenomic “in situ” observation?
125 -- it would be helpful for the reader to enumerate examples where this has been successful?
130 -- can you define the difference between species presence and genomic potential?
169 -- 50 individual observations per sample? Per station?
176 -- regarding enzyme selection: this would benefit from a negative control. Where you don’t see the enzyme, do you not see the activity? This would rule out muddled results caused by enzymes being used for concurrent N cycling activities. It would also help to cite papers where these genes are shown to be responsible for these processes in eukaryotes.
257 -- can you ground this .25 association in anything other than this past publication?
277 -- this is unclear, what makes high latitudes different?
305 -- as I said before, I suspect this DNRA pattern has to do with the boundaries of the euphotic zone being vastly different across these latitudes, where the defined “epipelagic” includes a a fair amount of darkness and anoxia in productive waters
309 -- enhanced nitrogen?
310 -- is the conclusion that more cells are in these regions expressing these genes, or are these genes being comparatively overexpressed? I know you can’t speak to expression data specifically, but it is difficult to grasp the meaning of “genomic potential” without acknowledging it.
464 -- the eastern Pacific upwelling zone, or the California Current Ecosystem, is a high productivity, highly oxygenated system. Maybe you meant the Eastern Tropical North Pacific ODZ?
470 -- how does figure 6B refer to diatoms?
471 -- These papers referring to the potential for diatoms (no capitalization) use of DNRA do not take place in high latitude conditions. Kamp and Stief have exclusively shown diatom DNRA to be associated with benthic conditions: complete anoxia, and induced nitrate storage. It is unlikely that these conditions are met in epipelagic, high latitude water.
480 -- this argument, that genes present don’t necessarily well represent expression, seems to undercut the base assumptions of this paper.
Summary: I like the idea of putting the ample amount of metagenomic data we have to use in models like CEPHALOPOD, and I’m glad to see efforts being put toward this. Novel results will only be credible if they are based on fine-tuned inputs developed from real world biogeochemistry and nitrogen cycling patterns, which I think are a bit lacking in this paper. If the paper pivoted toward using known patterns to validate the model input, I would support this. I would also like to note that some of the conclusions/introductions read as extremely text-generated, and I would encourage the authors to use the results to drive their own conclusions in correlation with existing literature.Citation: https://doi.org/10.5194/egusphere-2026-1459-RC2 -
AC2: 'Reply on RC2', Alexandre Schickele, 24 Jul 2026
Authors response to RC2
We thank the editor and the two reviewers for their feedback on the manuscript. Here we provide our answer to RC2. We numbered all comments and provided answers below and in the supplementary pdf. Lines are numbered according to the new manuscript with track change disabled.
Reviewer 2
General comments:
This work uses correlations between global genomic datasets and biogeochemical environmental data to validate and extrapolate on trends in prokaryotic nitrogen cycling predicted by the CEPHALOPOD statistical model. This is a proof of concept for using metagenomics to predict what the authors refer to as “genomic potential”, or the possibility that cells are expressing the proteins encoded by the genes sequenced in global metagenomic efforts. This work benefits from validation via comparison with known nitrogen cycling trends, and goes on to use the tool for prediction of unknown prokaryotic trends. This is most effective as a source of validation, where the CEPHALOPOD model successfully predicts known trends like increased assimilation of N in light, oxic oligotrophic equatorial waters. However, this doesn’t hold up when predicting N cycling in novel areas. There are base assumptions fed into the model that I think have produced spurious results. Like the first reviewer noted, I have doubts that DNRA and other anaerobic processes are comparatively more prevalent in well oxygenated upwelling systems and polar regions. These analyses hinge on the identification of specific genes to serve as a proxy for processes, which is a difficult thing to do even with metatranscriptomics. I think this paper would be served by a closer look at specific N gene function in eukaryotes, and comparing them to the known geography of prevalent N cycling processes, to use as further validation of the CEPHALOPOD model.
Thank you for this thoughtful overview of our work and constructive suggestions that have substantially improved the manuscript. Notable changes include a revision of the marker-gene selection, replacing nirB and nirD with the more pathway-specific nrfA and nrfH markers for DNRA, leading to a substantially revised and more realistic global DNRA distribution. We also redefined the epipelagic and mesopelagic layers using the climatological mixed layer depth, clarified the definition and interpretation of genomic potential, and better place our findings in the biogeochemical modelling context.
Specific scientific comments:
1. Light inhibition is not necessarily correlated to anaerobic conditions, especially in high-latitude regions. Those areas are often characterized as oxygen rich, high nutrient, low chlorophyll. This would make oxygen available as an electron receptor during photosynthesis, and eliminate the need for DNRA and other redox reactions. Deeper ocean layers are often oxygen rich while being nutrient poor, depending on depth below the euphotic zone, where oxygen tends to rebound as there are fewer organisms to consume it.
Thank you for this comment. We agree that low-light conditions do not necessarily imply low-oxygen or anaerobic conditions. Our intention was not to imply a direct relationship between light limitation and anaerobic metabolism, but rather to discuss light as one of several environmental factors associated with the distribution of the modeled pathways. In our results, nitrification marker genes are more abundant in low-irradiance but oxygenated waters, both in high-latitude epipelagic regions and throughout much of the mesopelagic ocean. This is consistent with the fact that nitrification is inhibited by high irradiance, while still requiring oxic conditions. The following statements support this:
L423-425: In the epipelagic layer, RDA1 and RDA2 positively associate with ANRA (0.08, 0.28), nitrogen fixation (0.39, 0.43), temperature (0.41, 0.10), and PAR (0.37, 0.08), and negatively with nitrification (–0.81, -0.01), nitrate (–0.32, -0.21), and oxygen (–0.-0.39, -0.10) concentrations.
L486-493: Under oxic conditions, Nitrification encoding genes are predominantly associated to Nitrososphaeria archaea, which corresponds to established knowledge in the literature (Tang et al., 2023; Ward, 2008), thus validating taxonomic annotation of our selected marker genes. The genomic potential, i.e., the relative abundance of those genes negatively correlates to PAR and positively correlates to nitrate and oxygen concentration in the epi- and mesopelagic layers (Fig. 6b), known drivers of nitrification rates (Tang et al., 2023; Ward, 2008; Zhang et al., 2020). Moreover, the estimated genomic potential distribution, showing elevated nitrification gene abundance in high-nutrient, low-chlorophyll (HNLC) regions such as the North Pacific and Atlantic, and Southern Ocean, aligns with higher measured nitrification rates (i.e., ammonia and nitrite oxidation; Tang et al., 2023).
L522-523: Our findings on nitrification, denitrification, and DNRA gene distributions correspond closely with known and decoupled oxygen and light controls (Lam and Kuypers, 2011)[…]
2. I would consider adjusting your depth-based strata based on region. In oligotrophic gyres, the chlorophyll maximum will often occur at 150m, meaning that the region above that will be well lit and nutrient depleted.
Thank you for this important suggestion. We agree that fixed depth strata may not consistently represent equivalent ecological environments across ocean regions. To account for this spatial variability, we redefined both the epipelagic and mesopelagic layers using a mixed-layer-depth climatology that was already included among the environmental predictors. This stratification was applied consistently to both the environmental and biological datasets. Consequently, we also revised the set of environmental predictors considered in the models and restricted surface-specific predictors to those appropriate for the newly defined layers (see Table S1). Relevant changes in the manuscript as stated below:
L194-196: Finally, we associated these observations with depth layers from 0 m to the Mixed Layer Depth (MLD; Boyer et al., 2018) and from the MLD to 1,000 m, here considered as the epipelagic (213 samples) and mesopelagic (180 samples) layers of the ocean, respectively.
L225-234: To capture the environmental drivers of gene distribution associated with nitrogen transformation pathways, we selected 12 monthly resolved, observation-based environmental climatologies (i.e., hereafter named features) with a spatial resolution of 1° x 1°, at the global scale (Table S1) (Schickele et al., 2025). These features describe the physical (e.g., sea water temperature), chemical (e.g., nutrient concentrations, oxygen chemistry), and biological (e.g., apparent oxygen utilization) properties of water masses, as well as their circulation and turbulence characteristics (e.g., eddy kinetic energy, EKE or Finite Size Lyapunov Exponent, FSLE). At climatological scales, EKE and FSLE do not aim to assess sub-mesoscale processes, but rather highlight large scale ocean currents (e.g., Gulf Stream, Alghulas current). When applicable (i.e., the 8 features retrieved from World Ocean Atlas)(Boyer et al., 2018), environmental fields were integrated over the depth ranges of biological sampling: from 0 m to the MLD (Boyer et al., 2018) for the epipelagic layer and from the MLD to 1000 m for the mesopelagic layer.
3. Coastal, upwelling regions can have a chlorophyll maximum of 5-10m. Adjusting this will make your environmental data averages more accurate.
Following your previous comment #2, the suggested modifications allow the modeled layers to better reflect region-specific environmental structure rather than fixed depth intervals.
4. Has anyone shown that “genomic potential” has predictive power in prokaryotes? I don’t work with prokaryotes, so I’m not sure if there are studies tying metagenomics to metatranscriptomics on this scale in these organisms. If so, I would prominently include them to bolster your case for DNA as a proxy for likely expression.
Thank you for this important comment. We acknowledge that the definition of genomic potential in the initial manuscript was not sufficiently clear (former L121: The relative abundance of such genes, referred to as the genomic potential for a given metabolic process […]). The genomic potential does not simply refer to the presence of a marker gene in a genome. Rather, it corresponds to the relative abundance of marker genes encoding a given metabolic function in prokaryotic metagenomes, estimated from genome-wide metagenomic read coverage and normalized for gene length and sequencing depth. Importantly, genomic potential should not be interpreted as a direct measure of metabolic rates or gene expression. Instead, it provides an indicator of where a given metabolic function is most encoded in the prokaryote genomes, thus likely to be expressed. We have clarified this point in the manuscript as follows:
L115-120: The relative abundance of genes encoding a given metabolic function in prokaryotic metagenomes can be estimated from genome-wide metagenomic read coverage. Following Schickele et al. (2024), we refer to this relative abundance (e.g. of nifH,D,K) as the genomic potential of a prokaryotic community for that function (e.g., nitrogen fixation). It does not simply describe the presence or absence of a gene in an individual genome, but the prevalence of genes encoding a given metabolic function within a sampled prokaryotic community, thereby linking prokaryotic identity and functional capacity (Sunagawa et al., 2015).
Several studies nonetheless support links between gene abundance, gene expression, and nitrogen-transformation rates in prokaryotes. In the marine environment, Salazar et al. (2019) showed that both nif transcript abundance and nitrogen fixation activity increase under low-nitrogen conditions in subtropical oceans, consistent with the environmental patterns of nif gene abundance estimated in our study. Similarly, Pogoreutz et al. (2017) reported a correspondence between nif gene abundance, nif expression, and nitrogen fixation activity in coral-associated microbial communities. Beyond marine systems, Morales et al. (2010) found relationships between denitrification gene abundance and denitrification activity in soils. We added these examples to the introduction to better position our approach relative to previous work:
L120-124: Although genomic potential should not be interpreted as a direct proxy for absolute metabolic rates, previous studies have reported similar environmental responses between nitrogen-cycle gene abundance, transcription, and process rates across marine (Pogoreutz et al., 2017; Salazar et al., 2019) and terrestrial (Morales et al., 2010) prokaryotic communities, supporting the use of genomic potential as an indicator of large-scale functional biogeography for the nitrogen cycle.
In addition, we performed an explicit comparison between the projected genomic potential hotspots for nitrogen fixation and denitrification and independent model-based rate estimates from both mechanistic and inverse biogeochemical approaches (Fig. S4)(Härri et al. 2026; Wang et al. 2019). We show that the spatial overlap between genomic potential and modeled rate hotspots is of the same order of magnitude as the overlap between the two modeling approaches themselves. While these results do not imply a direct one-to-one correspondence between gene abundance and process rates, they support the use of genomic potential as an informative indicator of the large-scale biogeography of nitrogen-cycling functions in prokaryotic communities. We therefore added the corresponding statements:
L398-405: The dominance of nitrogen fixation in epipelagic oligotrophic gyres and of denitrification in oxygen-deficient and upwelling regions is consistent with independent model-based estimates of these rates (Fig. S4). Using the 75th quantile as a threshold, the spatial overlap between genomic potential and model-based rate hotspots is of the same order of magnitude as the overlap between mechanistic and inverse model estimates themselves: 59 to 47% compared to 44% for nitrogen fixation, and 38 to 50% compared to 52% for denitrification (Fig. S4). In conclusion, metabolic requirements similarly shape the spatial distribution of the genomic potential for nitrogen transformation pathways for both bio-unavailable and bioavailable forms, that are consistent with the spatial distribution of independently modelled biogeochemical rates.
L462-467 : We showed elevated abundance of ANRA and nitrogen fixation encoding-genes in the epipelagic layer of tropical and subtropical regions, strongly associated with high temperature and low nitrate concentrations, consistent with known environmental controls on nitrate assimilation and nitrogen fixation (Luo et al., 2014; Zhang et al., 2020) and the identified nitrogen fixation hotspots from biogeochemical models (Fig. S4)(Härri et al., 2026; Wang et al., 2019). This agreement supports the interpretation of genomic potential as an indicator of the large-scale biogeography of nitrogen fixation capacity in prokaryotes.
(see Figure in the supplementary pdf attached)
Figure S4. Spatial overlap of hotspot regions defined by the 75th quantile of genomic potential and independently modelled biogeochemical rates from ROMS-BEC (Härri et al., 2026) and an inverse model approach (Wang et al., 2019). Overlap is shown for nitrogen fixation in the (a) epipelagic and (b) mesopelagic layers, and for denitrification in the (c) epipelagic and (d) mesopelagic layers.
5. I am slightly confused about the choice of nirB and nirD as a proxy for DNRA activity. As far as I know, nirB is also used for nitrate assimilation (ANRA), and some organisms use it concurrently with other forms of nitrite reductase in order to bolster N assimilation when “fresh” N is available. Its use alongside nirD, in the form sometimes referred to as nirBD, is documented as capable of DNRA alongside ANRA. The ratio of nirBD to nirA could be potentially used as evidence of DNRA activity, but I would include evidence that nirBD alone can be used as a DNRA proxy if you are convinced. In prokaryotes, nrfA is often used as the diagnostic gene for DNRA.
Thank you for this major comment, that has also been raised by reviewer #1 (RC1) in comment #8. Although nirB and nirD are annotated as DNRA-associated genes in KEGG, this is contested in the literature. In this context, and in agreement with your suggestion and reviewer #1, we revised our definition of the DNRA pathway by replacing the previously used nirB and nirD markers with nrfA and nrfH. These are more widely recognized in the literature as specific markers of DNRA. All analyses, figures, and results were updated accordingly. With this revised and more specific definition, DNRA-related marker genes are predominantly detected in oxygen-depleted regions, particularly in the eastern tropical Pacific upwelling system, consistent with the current literature describing this pathway. The relevant changes are detailed below:
L328-331: We observe a clear dominance of nitrogen transformation pathways associated with bioavailable over bio-unavailable forms of nitrogen – namely ANRA, and nitrification – in both epi- and mesopelagic layers (Fig. 3). In the epipelagic layer, ANRA accounts for the largest share of the genomic potential (0.79), followed by nitrification (0.17). However, DNRA that is also associated with bioavailable forms of nitrogen only accounts for a minor share of the genomic potential (<0.01).
L356-360: DNRA, typically associated with heterotrophic activity in suboxic to anoxic environments, shows peak genomic potential in the tropical Pacific upwelling system, with values ranging from 0.015 in the epipelagic to 0.03 in the mesopelagic layers (Fig 4c). Additional hotspots are predicted in ODZ of the northern Indian Ocean and in the Benguela upwelling system, with values ranging from 0.01 in the epipalgic to 0.02 in the mesopelagic layers (Fig. 4c). These are associated with a high uncertainty however (NSD > 0.5 ; see Fig. S3 for detailed uncertainty pattern).
L364-373: see Figure in the supplementary pdf
Figure 4: Genomic potential for bioavailable nitrogen transformation pathways represented by their relative abundance of metagenomic reads. Nitrogen transformation pathways are displayed in rows with: (a) Assimilatory Nitrate Reduction to Ammonium (ANRA), (b) Nitrification and (c) Dissimilatory Nitrate Reduction to Ammonium (DNRA). Left panels correspond to the epipelagic projections and right panel correspond to the mesopelagic projections. Each panel represents the annual average spatial projection, with white stippling representing geographical areas associated with a high degree of uncertainty (i.e., standard deviation above 0.5 the global average). Latitudinal profiles are displayed in form of annual average (black), summer (red; northern hemisphere June to August) and winter (blue; northern hemisphere December to February). The projection uncertainty across bootstraps is presented as grey shading in the latitudinal profiles.
L446-448: Taxa that are obligatory or facultative anaerobes also show contributions to epipelagic DNRA (0.1 to 0.2 for each taxon) while mesopelagic DNRA is mainly associated with Marinisomata in our dataset (0.8).
L597-502: The genomic potential for anaerobic nitrogen transformation pathways, such as denitrification and DNRA, positively correlates with low oxygen saturation and high apparent oxygen utilization (Fig. 6b). This is consistent with the dependence of Gammaproteobacteria on sub- or anoxic conditions to perform anaerobic respiration (Lam and Kuypers, 2011). The eastern Pacific upwelling, Arabian Sea, and Bay of Bengal ODZs show the highest genomic potential for DNRA and denitrification, reflecting known hotspots of anaerobic nitrogen transformation (Zhang et al., 2023).
6. L58 -- generally, can you clarify your rates and put context to Tg N/year? It would be helpful to be able meaningfully compare these parts of the nitrogen cycle.
Thank you for this comment. We agree that rates expressed in Tg N yr⁻¹ can be difficult to interpret without a broader biogeochemical context. To facilitate comparisons among nitrogen transformation pathways, we have added contextual information describing the relative magnitude of ANRA compared with other major marine nitrogen fluxes discussed in the introduction. We now write:
L58-59: This corresponds to an annual flux of approximately 850 to1700 Tg N yr⁻¹, making ANRA the largest nitrogen transformation processes in the global ocean.
7. L94 -- denitrification has been shown to happen at some level in all oceans
Thank you for this comment. We agree that our original wording was too restrictive and have revised the statement accordingly as following:
L98-102: Being an obligatory anaerobic process, open-water denitrification dominates in ODZs of the Arabian Sea, the eastern tropical Pacific, and parts of the eastern (sub)tropical Atlantic. Nevertheless, low-oxygen microenvironments within sinking organic particles can support denitrification in otherwise oxygenated waters, thereby extending the potential niche of denitrification across all oceans, beyond the ODZs (Bianchi et al., 2018; Deng et al., 2024; Wang et al., 2019)
8. L106 -- can you be more clear about what a “taxonomic unit” is? A whole genome? A gene? What is a metagenomic “in situ” observation?
Thank you for your comment. The taxonomic unit refers to the reconstructed genomes, each being assigned to a taxon. We clarified the statement as follows:
L113-115: Each of these reconstructed genomes is associated to a single taxon, thus corresponds to a taxonomic unit. They are also functionally annotated based on the enzymes they encode and the metabolic processes these enzymes are known to support (Kanehisa and Goto, 2000; Paoli et al., 2022).
9. 125 -- it would be helpful for the reader to enumerate examples where this has been successful?
Thank you for this comment. As stated in the response to your comment #4, we have added examples that support links between gene abundance, gene expression, and nitrogen-transformation rates in prokaryotes. In the marine environment, Salazar et al. (2019) showed that both nif transcript abundance and nitrogen fixation activity increase under low-nitrogen conditions in subtropical oceans, consistent with the environmental patterns of nif gene abundance estimated in our study. Similarly, Pogoreutz et al. (2017) reported a correspondence between nif gene abundance, nif expression, and nitrogen fixation activity in coral-associated microbial communities. Beyond marine systems, Morales et al. (2010) found relationships between denitrification gene abundance and denitrification activity in soils. We added these examples to the introduction to better position our approach relative to previous work:
L120-124: Although genomic potential should not be interpreted as a direct proxy for absolute metabolic rates, previous studies have reported similar environmental responses between nitrogen-cycle gene abundance, transcription, and process rates across marine (Pogoreutz et al., 2010; Salazar et al., 2019) and terrestrial (Morales et al., 2017) prokaryotic systems, supporting the use of genomic potential as an indicator of large-scale functional biogeography for the nitrogen cycle.
10. L130 -- can you define the difference between species presence and genomic potential?
Thank you for this comment. As stated in the response to your comment #4, we acknowledge that the definition of genomic potential in the initial manuscript was not sufficiently clear (former L121: The relative abundance of such genes, referred to as the genomic potential for a given metabolic process […]). The genomic potential does not simply refer to the presence of a marker gene in a genome. Rather, it corresponds to the relative abundance of marker genes encoding a given metabolic function in prokaryotic metagenomes, estimated from genome-wide metagenomic read coverage and normalized for gene length and sequencing depth. Importantly, genomic potential should not be interpreted as a direct measure of metabolic rates or gene expression. Instead, it provides an indicator of where a given metabolic function is most encoded in the prokaryote genomes, thus likely to be expressed. We have clarified this point in the manuscript as follows:
L115-120: The relative abundance of genes encoding a given metabolic function in prokaryotic metagenomes can be estimated from genome-wide metagenomic read coverage. Following Schickele et al. (2024), we refer to this relative abundance (e.g. of nifH,D,K) as the genomic potential of a prokaryotic community for that function (e.g., nitrogen fixation). It does not simply describe the presence or absence of a gene in an individual genome, but the prevalence of genes encoding a given metabolic function within a sampled prokaryotic community, thereby linking prokaryotic identity and functional capacity (Sunagawa et al., 2015).
11. L169 -- 50 individual observations per sample? Per station?
The station is used to define a location in space and time (longitude, latitude, climatological month). The sample correspond to the genomic material taken from the water at a certain depth, for a given station. We now clarify this as following:
L181-183: Hereafter, we consider a station as a location in space and time (longitude, latitude, climatological month), while the sample integrates the depth information as well (i.e., a station can have an epipelagic and a mesopelagic sample).
12. L176 -- regarding enzyme selection: this would benefit from a negative control. Where you don’t see the enzyme, do you not see the activity? This would rule out muddled results caused by enzymes being used for concurrent N cycling activities. It would also help to cite papers where these genes are shown to be responsible for these processes in eukaryotes.
Thank you for this suggestion. We agree that the ideal validation of enzyme selection would be to compare marker-gene abundance with transcript abundance and rate measurements from the same samples. However, such paired datasets are not currently available at the global scale for all nitrogen transformations considered here.
We therefore minimized ambiguity through a conservative marker-selection strategy. We retained only pathway-specific genes identified by KEGG annotations (Kanehisa and Goto, 2000) and previous marine prokaryotic metagenomic and metatranscriptomic studies (e.g., Salazar et al., 2019; Tang et al., 2023), and excluded markers known to be shared among multiple pathways (e.g., nxrA for nitrification and anammox-associated nitrite reductases). Following the reviewer #1 comments, we further refined marker selection by replacing nirB/nirD with nrfA/nrfH for DNRA and removing pmoA from the nitrification marker set. We have clarified this rationale in the manuscript as follow:
L189-191: (1) the gene must encode for an enzyme specific to one of the targeted nitrogen transformation pathways based on KEGG annotations and established literature (Kanehisa and Goto, 2000; Salazar et al., 2019; Tang et al., 2023).
L196-200: This filtering resulted in a set of 13 genes encoding enzymes specific to nitrogen fixation, ANRA, nitrification, DNRA, and denitrification, based on the KEGG classification (Fig. 2b)(Kanehisa and Goto, 2000). This conservative approach avoids selecting marker genes associated with concurrent nitrogen transformations, and is supported by marine prokaryotic metagenomic and metatranscriptomic studies, which have used these genes to characterize nitrogen transformation pathways (Salazar et al., 2019; Tang et al., 2023).
13. L257 -- can you ground this .25 association in anything other than this past publication?
We thank the reviewer for this comment. Similar thresholds have been used in other metagenomic-based studies (e.g. Faure et al., 2020), although these typically refer to univariate R² values, whereas our framework is multivariate. However, the 0.25 threshold is not based on this literature, but is directly derived from the validation and sensitivity analyses conducted in Schickele et al. (2025), where the CEPHALOPOD model was originally developed and benchmarked. In that study, this value was identified as a conservative cut-off, based on systematic evaluation of predictive skill across simulated datasets with varying levels of sampling effort and observational bias. Because CEPHALOPOD is applied here in the same model configuration and calibration framework, we retained this threshold to ensure consistency and comparability with the original benchmarked implementation. We have clarified this point in the manuscript as follows:
L286-289: A minimum threshold of 0.25 was required to ensure sufficient predictive power, following the benchmark evaluation of the CEPHALOPOD framework in Schickele et al. (2025), where this value was retained as a conservative cut-off for model performance. Comparable thresholds have also been adopted in other metagenomic-based modeling studies (e.g., Faure et al., 2021; Schickele et al., 2024).
14. L277 -- this is unclear, what makes high latitudes different?
Thank you for this comment. Here, our intention was to refer specifically to the extrapolation diagnosed by the MESS analysis. We have clarified that high-latitude and coastal regions are characterized by combinations of environmental values that fall outside the environmental range represented in the observations used to train CEPHALOPOD. Projections in these regions therefore correspond to model extrapolation and should be interpreted with caution. We revised the text as follows:
L310-314: The MESS analysis (see Sect. 2.2.1.) highlights however that environmental conditions characterizing high-latitude and coastal regions fall outside the range of environmental conditions represented in the observations used to train CEPHALOPOD. This means that projections in these regions represent model extrapolation, were not directly evaluated during model calibration, and should therefore be considered with caution (Fig. S2).
15. L305 -- as I said before, I suspect this DNRA pattern has to do with the boundaries of the euphotic zone being vastly different across these latitudes, where the defined “epipelagic” includes a a fair amount of darkness and anoxia in productive waters
Thank you for pointing this out, however the corresponding statement has been removed in the present in the revised manuscript, following our answer to reviewer #1, comment #8, and your comment #5.
In addition, as a response to your comment #2, we redefined the epipelagic and mesopelagic layers using the climatological mixed layer depth (MLD) rather than fixed depth intervals, to better account for regional variations in upper-ocean structure. Finally, following your suggestion and that of reviewer #1, we revised the DNRA marker-gene selection by replacing nirB and nirD with the more pathway-specific nrfA and nrfH gene markers, which substantially changed the projected DNRA distribution, consistent with current literature of marine DNRA biogeography.
16. L309 -- enhanced nitrogen?
Thank you for your comment. The statement was removed from the manuscript as DNRA is no longer the dominant pathway in high latitudes when considering the marker genes, you, and reviewer #1 suggested to use.
17. L310 -- is the conclusion that more cells are in these regions expressing these genes, or are these genes being comparatively overexpressed? I know you can’t speak to expression data specifically, but it is difficult to grasp the meaning of “genomic potential” without acknowledging it.
Thank you for your comment. This statement has been removed now that DNRA is not the dominant nitrogen pathway anymore in the high latitudes, when considering the new marker genes, you, and reviewer #1 proposed.
18. L464 -- the eastern Pacific upwelling zone, or the California Current Ecosystem, is a high productivity, highly oxygenated system. Maybe you meant the Eastern Tropical North Pacific ODZ?
Thank you for your comment. We do not see any mention of the “California current ecosystem” in the initial manuscript but agree that the mention of “eastern Pacific upwelling” is non-specific and confusing. We meant Eastern Tropical Pacific ODZs. This is now clarified as following:
L500-501: The eastern tropical Pacific, Arabian Sea, and Bay of Bengal ODZs show the highest genomic potential for DNRA and denitrification […].
19. L470 -- how does figure 6B refer to diatoms?
Thank you for pointing this out. The figure reference was indeed not justified. However, the updated model using the suggested DNRA marker genes does not present a hotspot in this area anymore, we therefore removed the statement.
20. L471 -- These papers referring to the potential for diatoms (no capitalization) use of DNRA do not take place in high latitude conditions. Kamp and Stief have exclusively shown diatom DNRA to be associated with benthic conditions: complete anoxia, and induced nitrate storage. It is unlikely that these conditions are met in epipelagic, high latitude water.
Thank you for pointing this out. As for the previous comment, references to diatoms have been removed as the DNRA genomic potential using the new marker genes you and reviewer #1 suggested does not present a hotspot in high latitudes. The related statements have been removed from the revised manuscript.
21. L480 -- this argument, that genes present don’t necessarily well represent expression, seems to undercut the base assumptions of this paper.
Thank you for pointing this out. We agree that, as originally written, this statement could be interpreted as undercutting the rationale of the study. As clarified in our response to your comment #4, the genomic potential is not a direct proxy for absolute gene expression or nitrogen transformation rates under all environmental conditions. Rather, it is used as an indicator of where a given nitrogen-cycling function is most strongly represented in the prokaryotic metagenome, and therefore most likely to be expressed at large spatial scales.
L120-124: Although genomic potential should not be interpreted as a direct proxy for absolute metabolic rates, previous studies have reported similar environmental responses between nitrogen-cycle gene abundance, transcription, and process rates across marine (Pogoreutz et al., 2010; Salazar et al., 2019) and terrestrial (Morales et al., 2017) prokaryotic systems, supporting the use of genomic potential as an indicator of large-scale functional biogeography for the nitrogen cycle.
Our point in this section was not to weaken that assumption, but to acknowledge that the relationship between gene abundance, expression, and realized rates is not expected to be one-to-one everywhere, because local expression and activity remain modulated by local factors such as environmental conditions or physiological regulation, as well as the constraint of having relative abundance. We have therefore revised the text to make clear that this limitation concerns the quantitative interpretation of genomic potential as a predictor of local absolute rates, not its use as an indicator of the large-scale biogeography of nitrogen-cycling functions.
L502-513: We also estimated a genomic potential for epipelagic denitrification in the eastern tropical Pacific upwelling system, where gene expression (norB, nosZ) has previously been reported, although the mechanisms sustaining this remain unresolved (Ganesh et al., 2015). We hypothesize that this signal may reflect episodic vertical mixing that brings low-oxygen waters and associated prokaryotic assemblages from the mesopelagic into the epipelagic zone. In contrast, in more strongly stratified regions like the Arabian Sea and Bay of Bengal, where ODZs remain isolated from the epipelagic layer, we do not observe a comparable epipelagic genomic potential. This suggests that physical connectivity (e.g., vertical mixing) may contribute to local differences between the denitrification genomic potential compared to model-based denitrification rate estimates in the eastern tropical Pacific, although both approaches identify this region as a denitrification hotspot (Härri et al., 2026; Wang et al., 2019). More generally, while genomic potential captures the large-scale distribution of the functional capacity for denitrification, its realized activity likely remains sensitive to local forcing factors such as connectivity or oxygen conditions (Dalsgaard et al., 2014; Zhang et al., 2023).
L537-545: Our approach is based on genomic potential, which captures the large-scale patterns of metabolic capabilities encoded in prokaryotic communities. However, local scales, other factors such as physical connectivity or environmental regulation (i.e., not climatologies) may modulate the relationship between genomic potential and the magnitude of realized metabolic activity or process rates, although their spatial patterns concur (Dalsgaard et al., 2014; Salazar et al., 2019). Metatranscriptomic data provide complementary information on gene expression and may help refine the interpretation of local variations in genomic potential by (Dalsgaard et al., 2014). However, metatranscriptomic observations exhibit a high spatial and temporal heterogeneity and are currently limited in sampling coverage (Paoli et al., 2022). Therefore, increasing both in-situ observations density and model resolution are needed to better capture the relationships between genomic potential, expression, and metabolic rates across scales.
22. Summary: I like the idea of putting the ample amount of metagenomic data we have to use in models like CEPHALOPOD, and I’m glad to see efforts being put toward this. Novel results will only be credible if they are based on fine-tuned inputs developed from real world biogeochemistry and nitrogen cycling patterns, which I think are a bit lacking in this paper. If the paper pivoted toward using known patterns to validate the model input, I would support this. I would also like to note that some of the conclusions/introductions read as extremely text-generated, and I would encourage the authors to use the results to drive their own conclusions in correlation with existing literature.
Thank you for this constructive final comment. In response, we first refined the selection of pathway-specific marker genes following your comments and those of reviewer #1 and reviewer #2 and clarified the definition and interpretation of genomic potential (see response to comment #4, 9 and 10). All reviewer’s suggestions also helped to clarify the introduction and to better anchor our discussion in the results and in the existing nitrogen-cycle literature.
Importantly, we also added an explicit comparison between the projected genomic potential hotspots for nitrogen fixation and denitrification and independent model-based rate estimates from both mechanistic and inverse biogeochemical approaches (Fig. S4; see our response to your comment #4). We show that the spatial overlap between genomic potential and modeled rate hotspots is of the same order of magnitude as the overlap between the two modeling approaches themselves. In addition to reference to existing literature, we use this comparison to support the interpretation that genomic potential captures the large-scale biogeography of microbial functional capacity for these pathways, while recognizing that local regulation and environmental forcing can modulate realized activity and rates. Most relevant additions are stated below:
L397-405: The dominance of nitrogen fixation in epipelagic oligotrophic gyres and of denitrification in oxygen-deficient and upwelling regions is consistent with independent model-based estimates of these rates (Fig. S4). Using the 75th quantile as a threshold, the spatial overlap between genomic potential and model-based rate hotspots is of the same order of magnitude as the overlap between mechanistic and inverse model estimates themselves: 59 to 47% compared to 44% for nitrogen fixation, and 38 to 50% compared to 52% for denitrification (Fig. S4). In conclusion, metabolic requirements similarly shape the spatial distribution of the genomic potential for nitrogen transformation pathways for both bio-unavailable and bioavailable forms, that are consistent with the spatial distribution of independently modelled biogeochemical rates.
L462-467: We showed elevated abundance of ANRA and nitrogen fixation encoding-genes in the epipelagic layer of tropical and subtropical regions, strongly associated with high temperature and low nitrate concentrations, consistent with known environmental controls on nitrate assimilation and nitrogen fixation (Luo et al., 2014; Zhang et al., 2020) and the identified nitrogen fixation hotspots from biogeochemical models (Fig. S4)(Härri et al., 2026; Wang et al., 2019). This agreement supports the interpretation of genomic potential as an indicator of the large-scale biogeography of nitrogen fixation capacity in prokaryotes.
L508-513: This suggests that physical connectivity (e.g., vertical mixing) may contribute to local differences between the denitrification genomic potential compared to model-based denitrification rate estimates in the eastern tropical Pacific, although both approaches identify this region as a denitrification hotspot (Fig. S4)(Härri et al., 2026; Wang et al., 2019). More generally, while genomic potential captures the large-scale distribution of the functional capacity for denitrification, its realized activity likely remains sensitive to local forcing factors such as connectivity or oxygen conditions (Dalsgaard et al., 2014; Zhang et al., 2023).
L560-570: In this study, we show that the marine microbial metagenome captures distinct spatial patterns in nitrogen transformations associated with either cellular biosynthesis or energy metabolism. Genes encoding for metabolic pathways supporting cellular biosynthesis (ANRA, Nitrogen Fixation) dominate in warm, oligotrophic, epipelagic waters, and are mainly associated with Cyanobacteria. In contrast, those supporting energy metabolism (nitrification, DNRA, denitrification) dominate in deeper waters, associated with Nitrosophaeria in oxic conditions, or facultative anaerobes such as Gammaproteobacteria in eastern boundary upwelling systems or oxygen-deficient environments. This is consistent with known prokaryotic metabolic strategies and modelled biogeochemical rates. Therefore, we suggest that metagenomics can complement conventional biogeochemical measurements by explicitly linking microbial communities to their biogeochemical functions. This positions metagenomic data as a promising basis to develop observation-based Essential Ocean Variables (EOVs; Muller-Karger et al., 2018) describing microbial biogeochemical and ecosystem functions, and to extend this approach to other elemental cycles, such as sulfur and methane metabolism (Salazar et al., 2019).
-
AC2: 'Reply on RC2', Alexandre Schickele, 24 Jul 2026
Interactive computing environment
Data and code to map the genomic potential for nitrogen metabolism in prokaryotes A. Schickele et al. https://doi.org/10.5281/zenodo.17407277
Viewed
Since the preprint corresponding to this journal article was posted outside of Copernicus Publications, the preprint-related metrics are limited to HTML views.
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 504 | 0 | 2 | 506 | 0 | 0 |
- HTML: 504
- PDF: 0
- XML: 2
- Total: 506
- BibTeX: 0
- EndNote: 0
Viewed (geographical distribution)
Since the preprint corresponding to this journal article was posted outside of Copernicus Publications, the preprint-related metrics are limited to HTML views.
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
General comments:
The authors examine global patterns in key nitrogen cycling pathways in the ocean using genetic datasets. The tool used -- the statistical model CEPHALOPOD -- has been previously published and provides a thorough and statistically rigorous analysis of how the key genes vary with environmental factors. This is a step up from simple DIY correlations between datasets. Most of the results themselves are not novel -- we know that assimilation of N by cyanobacteria should dominate in oligotrophic gyres, and nitrification at depth, for example, but it is still useful to see how these known patterns emerge using this tool. It is very interesting to see the widespread but low abundance coverage of the denitrification genes. However, the conclusion of high relative abundance of DNRA across the ocean, especially at high latitudes in the surface, needs more work to be convincing. Not only is it very surprising (why would an anaerobic process show up in the most oxygenated waters at high latitudes?), but the genetic marker may not be appropriate here. I am not a geneticist or omics expert, but my understanding is that nir B and nir D can be involved in ANRA as well as DNRA (see for example this review: https://doi.org/10.1016/j.synbio.2025.12.017), and ANRA is actually what we expect in high latitude surface waters. I suggest the authors look for nrfA instead and reevaluate their conclusions about DNRA.
Specific scientific comments:
Use of "aerobic" (ex: L. 17): I would not use the term aerobic for autotrophic processes/photosynthesis. I realize the ANRA may be associated with aerobic heterotrophy as well in the surface, which may explain why a small section of the pie in Fig. 6c for ANRA is not cyanobacteria. I and many others consider aerobic to mean O2-consuming, not O2-producing, so may be confusing.
Genes for nitrification: Again, I am not an expert, but why not look for amoA and nxr genes for ammonia and nitrite oxidation? This is what is used for this database on nitrifiers: Tang et al "Database of nitrification and nitrifiers in the global ocean," https://doi.org/10.5194/essd-15-5039-2023 -- which should be cited and with which results here could be compared.
Nitrification taxonomy: It is reassuring that the taxonomy associated with nitrification is Nitros*, but it is perhaps a bit circular to pitch it as a conclusion (considering the literature) since for nitrification, taxonomy and function are so related. It would make the paper stronger overall to pitch this more as a "validation" that the model can find known patterns such as this one.
Fig. 1: There are many wide black arrows, but then you say you did not consider them in this study. Perhaps there is another way to say this that is not so disserving? You are using methods that just identify correlations among datapoints, so probably OK to just not even mention that you cannot consider transport in the model. Same for anammox: you DID consider it, but it was just not sufficiently abundant to include in the statistical analysis.
L. 120: "Habitat modeling" seems like it should represent a much broader category of modeling, rather than just simply this statistical approach. I would assume that any ecosystem model that aims to identify the ecological niche of a population or functional type is a habitat model.
Cutoff (50 observations; line 169): You say that a gene must be present in at least 50 observations "to ensure sufficient spatial coverage" -- this does seem incompatible with many of the anaerobic N cycling transformations in anoxic zones unless you have 50 samples IN anoxic zones. This meant that you couldn't show anammox patterns. Perhaps for these types that didn't make that cut (anammox, and perhaps nrf as a replacement for nirBD for DNRA), it's useful to explain that they do "show up where they should," i.e., the small volumes of anoxic zones, rather than saying "we didn't consider them." You did consider them, but they are just not candidates for global-scale representation which is reasonable!
L. 225: "outliers in the biological inputs were removed." Is this reasonable, given that we do think there is much heterogeneity in microbial communities and BGC fluxes in the environment? What are the implications of this? Is it more like you are finding "average" patterns and ignoring the known deviations? (This seems reasonable to do.)
DNRA results: As I state above, I'm deeply skeptical that DNRA is actually this abundant in oxygenated surface waters. It would be really amazing if these results WERE true, and so, you need to convince me and the other readers that I can believe this. See above major comment about nrf vs. nirBD.
L. 497: Levine et al 2025 is a review, not a model. Cite the models you are referring to directly.
Smaller presentation comments:
L. 85: "ammonium that is being released" -- note that this may also be considered a key N transformation process! Perhaps also mention urea somewhere, as urea cycling is becoming more recognized currently.
L. 107: do you mean "such genes" instead of "such enzymes"?
L. 113: "the one[s] associated" (not "the one")
L. 145: "Doing so, [we] explore"
L. 272: "The three features" -- could you name them here? I am curious at this stage!
L. 302: Could you clarify/specify that "bio-unavailable forms" refers to N fixation and denitrification? This makes sense to look at these two separately from the others, but just took me a second to realize what was going on.
Figs. 4 and 5: Could you just label each of the maps with the process for much easier reference?
L. 408: Again, I would not call cyanobacteria "aerobic."
L. 410/11: Nitros* are named the way they are because they are associated with nitrification, so seems naive to state a result this way (see above comment).
L. 415: "In the lat[t]er"?
L. 416: Is it known that these are diazotrophs? Perhaps more fair to just say non-cyanobacterial groups?