Ecological control of nitrite in the upper ocean.

Microorganisms oxidize organic nitrogen to nitrate in a series of steps. Nitrite, an intermediate product, accumulates at the base of the sunlit layer in the subtropical ocean, forming a primary nitrite maximum, but can accumulate throughout the sunlit layer at higher latitudes. We model nitrifying chemoautotrophs in a marine ecosystem and demonstrate that microbial community interactions can explain the nitrite distributions. Our theoretical framework proposes that nitrite can accumulate to a higher concentration than ammonium because of differences in underlying redox chemistry and cell size between ammonia- and nitrite-oxidizing chemoautotrophs. Using ocean circulation models, we demonstrate that nitrifying microorganisms are excluded in the sunlit layer when phytoplankton are nitrogen-limited, but thrive at depth when phytoplankton become light-limited, resulting in nitrite accumulation there. However, nitrifying microorganisms may coexist in the sunlit layer when phytoplankton are iron- or light-limited (often in higher latitudes). These results improve understanding of the controls on nitrification, and provide a framework for representing chemoautotrophs and their biogeochemical effects in ocean models.

Introduction

Nitrogen proximally limits primary production in much of the surface ocean, and the n class="Chemical">nitrogen cycle exerts a strong influence on the coupled cycles of carbon and other elements[1]. Most fixed nitrogen in the ocean is in the oxidized form of nitrate (). Primary producers and other microorganisms in the surface ocean consume and reduce to build organic molecules. Microorganisms then oxidize detrital organic nitrogen back to in a series of steps, with intermediates ammonium (, here considered interchangeable with ammonia, NH3) and nitrite ()[2].

Despite its typically low concentration, plays a central role in global nitrogen and n class="Chemical">carbon cycles, providing a key resource for significant microbial metabolisms. is an intermediate of nitrification, the two-step, microbially mediated oxidation of to () that occurs in association with distinct clades of metabolically diverse chemoautotrophic archaea and bacteria in many environments[2-5].

The primary nitrite maximum (Pn class="Chemical">NM), an accumulation of at the base of the euphotic (sunlit) zone at concentrations of 10–1000 nmol L−1, is ubiquitous in oxygenated subtropical oceans[2] (Fig. 1). In a typical subtropical vertical profile, the PNM is located at the onset of the nitricline, just below the deep chlorophyll maximum (DCM) and an associated maximum[6-11]. The relative magnitude of the maximum varies, but peak [] is consistently lower than [] at the PNM in oligotrophic environments[10,12,13]. These subsurface maxima form in strongly stratified water columns, i.e., when the euphotic zone is deeper than the mixed layer[14,15]. Nitrification rates and the biomass of nitrifying microbes are often observed to peak at, or just below, the PNM[10,13,16-22]. In contrast, in subpolar regions, concentrations are elevated to similar magnitudes throughout all of the upper ocean, including the surface (Fig. 1). In low-oxygen environments, a secondary maximum forms below the PNM at much higher concentrations (order 10 μmol L−1) due to anaerobic activity[2]; however, we focus here only on the dynamics of aerobic environments.

Fig. 1

Observed and modeled nitrite () concentrations along three transects in the ocean. a Map of the GLODAPv2 transect locations. b Transects from the GLODAPv2 database[60, 61]. c Transects from the global simulation. Map generated with Python version 3.5.1, Matplotlib version 1.5.1, and Basemap version 1.0.7[86]

Despite the widespread occurrence of the PNM, the mechanisms that determine the locations and magnitude of accumulated are still not fully resolved. What is the source of the accumulated in n class="Chemical">oxygenated waters? Two hypotheses, not mutually exclusive, have been advanced[23]: (i) excretion of due to incomplete assimilatory reduction of by phytoplankton[14,24,25], and (ii) chemoautotrophic oxidation, the first step of nitrification[23,26]. Here, we focus on outstanding questions of the latter hypothesis, since isotopic evidence suggests that oxidation is a major source of at the PNM[10,21].

First, why does peak nitrification occur at the base of the euphotic zone? Photoinhibition of nitrifying microorganisms has been documented[26-30] and incorporated into ecosystem models[31,32], but observations of nitrification in the euphotic zone[16,18,33] and close to the sea surface[34-36] suggest that this is not universally the case. An alternate (but not exclusive) hypothesis is that phytoplankton outcompete slow-growing, chemoautotrophic nitrifiers for (and ) in the euphotic zone[37,38], but not deeper where light limits photoautotrophy.

Second, while both and can accumulate in the oxygenated thermocline, why does consistently accumulate to a higher concentration in oligotrophic envn class="Chemical">ironments? Hypotheses include differential effects of photoinhibition of and oxidation[26,39,40] and differential temperature sensitivity of their rates[41], though some observations conflict with the latter[42]. Here, we will suggest that this reflects differences in the metabolisms of the distinct clades of nitrifiying microorganisms.

Third, why does nitrification sometimes occur in the euphotic zone, and why is surface [] also elevated in some areas, such as in subpolar regimes? Possible explanations include accumulation due to rapid entrainment of from below the euphotic zone, chemoautotrophic nitrification within the euphotic zone, or the incomplete reduction of by phytoplankton.

Here we synthesize and address these questions and hypotheses using a hierarchy of mathematical models and simulations that are grounded in theoretical and laboratory evaluations of the kinetics and efficiencies of marine nitrifying microorganisms. The models are used to interpret both existing and new water column observations in the subtropical n class="Chemical">North Pacific as well as the global distribution of in the oxygenated upper ocean.

We first present a general population dynamics model for and oxidizers. The model couples estimates of cellular substrate uptake rates with an energetically informed stoichiometry of whole-organism metabolism. We then examine point balance solutions of the model and find that the relative concentrations of and near the PNM reflect the respective subsistence concentrations for -oxidizing and -oxidizing organisms, with exact concentrations also influenced by vertical mixing. We next discuss how competitive interactions control the vertical structure of nitrification and the position of the subtropical Pn class="Chemical">NM using simulations of water column profiles of an upper ocean microbial ecosystem. Finally, we implement the ecosystem model in a global ocean simulation, and find that deep mixed layers and the light or iron limitation of primary producers can explain [] and nitrification at the surface in some regions.

Results

Microbial population dynamics

We consider the clades of microorganisms carrying out the two steps of nitrification as two metabolic functional types: ammonia-oxidizing and n class="Chemical">nitrite-oxidizing organisms (hereafter, AOO and NOO, respectively). We describe the population dynamics of each type as:where B (mol L−1) is the biomass of type i with population growth rate μ (d−1), loss rate L (d−1), and physical transport. Loss rate is a function of the population density of grazers Z and other factors (Methods). We describe population growth rate μ using Monod kinetics with limiting resource R (determined by Liebig’s Law of the minimum). Growth rate depends on the yield of biomass with respect to the limiting resource y (mol biomass mol R−1), the maximum specific resource uptake rate (mol R per mol biomass per day), and a half-saturation concentration K (mol L−1) as:

The uptake parameters together give an expression for the specific uptake affinity (), a measure of competitive strength at low resource concentrations[43]. For the two nitrifying metabolisms, we assume that DIN and n class="Chemical">oxygen are the potentially limiting resources (R). Over most of the subtropical thermocline, oxygen exists at sufficiently high concentrations (greater than nanomolar) to serve as the terminal electron acceptor for AOO and NOO[44,45].

Equations (1) and (2) provide a general description of microbial population dynamics. We next consider idealized approximations and numerical solutions to interpret the controls on and in the subsurface ocean.

Point balance solution for []:[] at the PNM

Why does often accumulate to a higher degree than at the PNM[10,12,13]? Consider the local dynamics of the subsurface oligotrophic envn class="Chemical">ironment, assuming the simplest approximation of the model, in which physical transport is negligible and the system is close to steady state (). Using Eqs. (1) and (2), we evaluate the concentration of resource R that limits the growth of nitrifier type i as:

is the subsistence concentration of type i[46,47]. In theory, at a steady state, the population with the lowest R* excludes all others limited by the same resource (if its maximum growth rate, , is larger than loss rate L), and the environmental concentration is set to this . We hypothesize that in the vicinity of the Pn class="Chemical">NM, and are the respective limiting resources for AOO and NOO, respectively, and that the ratio of the environmental concentrations will therefore reflect the ratio of their respective subsistence concentrations.

What are the differences in for the nitrifier guilds? We first estimate the yields of the AOO and n class="Chemical">NOO functional types. Following established methodology[48], we derive stoichiometrically balanced equations to describe catabolic and anabolic processes as a function of the redox chemistry that underlies each step of nitrification. Constraining thermodynamic efficiency with published laboratory observations of marine AOO and NOO growth, we approximate the population-level growth stoichiometries of ammonia-oxidizing biomass BAOO and nitrite-oxidizing biomass BNOO as:where biomasses are assumed to be synthesized from one mole of for both types[4] (Methods). The predicted stoichiometry of and O2 demand for the AOO is 1:1.45, close to that observed in culture[37] (1:1.52) indicating that the model captures key aspects of real metabolisms. In Eqs. (4) and (5), the yields for AOO and NOO are and moles biomass N synthesized per mole DIN used, respectively. In other words, the NOO type consumes three times as much to produce a unit of biomass, relative to the consumption of by the AOO type.

Second, we speculate that NOO specific uptake affinity may be lower than that of n class="Chemical">AOO if (as is observed) -oxidizing bacteria are larger than -oxidizing archaea[37,49,50]. Established empirical and theoretical allometric relationships[51,52] suggest that specific affinity decreases with increasing cell radius r as r−2 (Methods). If we assume a 10-fold larger cell volume for NOO, this estimates a 4.6-fold lower specific affinity.

Using these differences in yield (y) and affinity (), we evaluate the ratio of the subsistence concentrations of and for the AOO and n class="Chemical">NOO, respectively, in Eq. (3). If AOO and NOO are of similar size/affinity, and if loss rates (L) are low relative to maximum growth rates (and thus negligible in the denominator), then we predict that []:[] ~  ~ 1:3 in the vicinity of the PNM. This is consistent with the observed ratio at the PNM in oligotrophic environments[10,12]. If yields are similar, but affinity is 4.6-fold different, we predict a similar order of enhanced accumulation. Together, differences in both yield and affinity would result in an even lower ratio, underestimating []: []. Below, we examine the independent and additive effects of both distinctions between the nitrifier metabolisms in the water column model.

Understanding vertical structure with a water column model

What sets the vertical structure of the PNM? To address this, we simulate an idealized oligotrophic n class="Chemical">water column in one dimension, in which the photon flux attenuates with depth and vertical mixing of biomass and resources is described by turbulent diffusion. Enhanced vertical mixing close to the surface simulates the mixed layer.

Our ecosystem model resolves the interactions and mixing of inorganic and n class="Chemical">organic nitrogen and nitrogen-based biomass of five microbial metabolic functional type populations (Fig. 2). AOO and NOO are modeled as above (parameter values in Supplementary Table 1), with growth rate limited by the supply of or , the supply of oxygen[45], or internal constraints (i.e., the maximum rate). For simplicity, we neglect other limitations such as iron availability (Supplementary Note 3). A phytoplankton functional type assimilates , , and into biomass with light- and nutrient-limited growth[53,54]. The traits of the phytoplankton type are based on the cyanobacterium Prochlorococcus, since it is often the most abundant photoautotroph in oligotrophic marine environments, and perhaps the strongest photoautotrophic competitor for DIN due to its high-nutrient affinity, a consequence of its small size[52,55]. A heterotrophic bacterial type remineralizes organic detritus to . A microzooplankton grazer type consumes the other microbes and excretes . Temperature modifies the rates of all organisms except for the nitrifying types, following experimental results[56] (Supplementary Note 1). A pool of sinking organic detritus is produced from the mortality of all populations. We do not impose photoinhibition or any other direct inhibitor of nitrification, and instead examine the vertical structure that emerges as a function of the ecological interactions.

Fig. 2

Schematic of the marine ecosystem model with explicit nitrification. The model resolves three species of inorganic fixed nitrogen (ammonium (), nitrite (), and nitrate ()), particulate organic nitrogen (PON), phytoplankton (P), zooplankton (Z), heterotrophic bacteria (Bhet), ammonia-oxidizing organisms (AOO), and nitrite-oxidizing organisms (NOO)

Schematic of the marine ecosystem model with explicit nitrification. The model resolves three species of inorganic fixed n class="Chemical">nitrogen (ammonium (), nitrite (), and nitrate ()), particulate organic nitrogen (PON), phytoplankton (P), zooplankton (Z), heterotrophic bacteria (Bhet), ammonia-oxidizing organisms (AOO), and nitrite-oxidizing organisms (NOO)

For comparison with the model, we show relevant observations collected at several stations in the North Pacific subtropical gyre in 2014 (Fig. 3; station locations in Supplementary Fig. 1). Peak [], detectable at two out of four stations, is positioned below the n class="Disease">DCM, reaching 110 and 130 nmol L−1. Detectable oxidation rates peak at this apparent PNM at these locations, reaching 3 and 6 nmol L−1 d−1. Gene abundances indicating -oxidizing archaea (amoA) also correlate with []. In general, the resulting vertical structures are consistent with other observations of the PNM in oligotrophic systems[10,12].

Fig. 3

Stratified water column model solutions and observations from the subtropical North Pacific. Model solutions are indicated with lines, with shaded areas denoting 1 s.d. of the model ensemble. Observations from the four stations are indicated with marked points in a–d, with error bars denoting 1 s.d.: a [] and [] (the shaded region of [] reaches about 1000 nM), b ammonia-oxidizing and nitrite-oxidizing organism (AOO and NOO) abundances (observed amoA gene abundances), c nitrification rates, and d [Chl a], PAR (scaled to fit plot), and []. Observations below the detection limit are indicated with open (vs. filled) markers. Also shown are solutions to additional model experiments (e–g), with only the difference in yield (reflecting the threefold energetic difference), with only the difference in affinity (here, as the 4.6-fold decrease related to a 10-fold larger NOO cell volume), and with no quantitative difference in the parameters describing AOO and NOO metabolisms

Stratified water column model solutions and observations from the subtropical n class="Chemical">North Pacific. Model solutions are indicated with lines, with shaded areas denoting 1 s.d. of the model ensemble. Observations from the four stations are indicated with marked points in a–d, with error bars denoting 1 s.d.: a [] and [] (the shaded region of [] reaches about 1000 nM), b ammonia-oxidizing and nitrite-oxidizing organism (AOO and NOO) abundances (observed amoA gene abundances), c nitrification rates, and d [Chl a], PAR (scaled to fit plot), and []. Observations below the detection limit are indicated with open (vs. filled) markers. Also shown are solutions to additional model experiments (e–g), with only the difference in yield (reflecting the threefold energetic difference), with only the difference in affinity (here, as the 4.6-fold decrease related to a 10-fold larger NOO cell volume), and with no quantitative difference in the parameters describing AOO and NOO metabolisms

The emergent PNM

With the subtropical configuration, a PNM emerges consistently in the n class="Chemical">water column as a consequence of the ecological interactions. We illustrate the model equilibrium solution that corresponds to our observed profile from the North Pacific (Fig. 3; full solutions in Supplementary Fig. 2). The simulation qualitatively reproduces the vertical structures of [], rates of oxidation, and AOO abundances. We explored the sensitivity of the model solution with an ensemble of ~1000 model realizations in which ecosystem model parameters were sampled randomly over reasonable ranges of uncertainty (see Methods and Supplementary Table 1). The s.d. of the ensemble of solutions (Fig. 3, shaded region) shows that the key results outlined below are indeed qualitatively robust despite the uncertainties.

With the steady-state solutions, we calculate the subsistence concentrations (R*s) to diagnose ecological control (Fig. 4). We use Eq. (3) for AOO and n class="Chemical">NOO, with maximum growth rate , and Eq. (28) for phytoplankton, incorporating their maximum light-limited growth rate (Methods). We then use the resulting loss rate L, a function of the dynamic grazer population Z (Eq. (29)), to calculate the depth-varying R*s for and .

Fig. 4

Water column model solutions and subsistence concentrations. a, b [] and [] with associated R*s of ammonia-oxidizing and nitrite-oxidizing organisms (AOO and NOO) and picophytoplankton (P). The R*s of P here are calculated with the maximum light-limited growth rate (Eq. (28)), and thus become negative at depth once light renders photoautotrophy unsustainable (see Supplementary Fig. 3 for further details)

Water column model solutions and subsistence concentrations. a, b [] and [] with associated R*s of n class="Chemical">ammonia-oxidizing and nitrite-oxidizing organisms (AOO and NOO) and picophytoplankton (P). The R*s of P here are calculated with the maximum light-limited growth rate (Eq. (28)), and thus become negative at depth once light renders photoautotrophy unsustainable (see Supplementary Fig. 3 for further details)

We find that the modeled phytoplankton exclude the chemoautotrophic nitrifiers in the euphotic zone, where both are limited by nitrogen. This is because the photoautotrophs have lower subsistence concentrations for and at the surface (Fig. 4). Theory and observations show that -oxidizing archaea and picophytoplankton have similar uptake affinities for [37,51,52]. However, for picophytoplankton, the effective half-saturation concentrations for n class="Chemical">nitrogen uptake with respect to growth rate are significantly lower[57], and maximum growth rate is higher. This reflects that nitrifiers use DIN for energy production, which is relatively inefficient, while the phytoplankton use it as a nitrogen source for synthesis with a much higher biomass yield.

Deeper in the water column, once light limits photoautotrophy, the maximum light-limited growth rate decreases significantly, leading to sharply increasing R*s (the green dashed lines in Fig. 4). With still increasing depth, phytoplankton losses become larger than their growth rate (Supplementary Fig. 3c), and the phytoplankton can no longer survive (where the green dashed lines end in Fig. 4). Once phytoplankton are excluded, the nitrifiers become competitive for and , and the ambient nutrient concentrations are set by the nitrifier R*s instead. This accounts for the large increase in [] and [] at about 100 m (in this example).

R*s decrease with depth as grazing pressure is reduced (see blue dashed line in Supplementary Fig. 3c), and so, after the increase, and then decline with depth. Thus, subsurface maxima in both (the PNM) and emerge in the simulations, controlled by a combination of top–down and bottom–up processes. As in the point balance, accumulates to a higher maximum concentration than because of the yield and affinity distinctions between n class="Chemical">AOO and NOO.

We further demonstrate the control of []:[], and its sensitivity to the nitrifier parameters, using model experiments in which we isolate the differences in yields and affinities of the NOO and n class="Chemical">AOO (Fig. 3e–g). With both the yield and affinity differences between the nitrifiers included, []:[] is about 1:10, and [] is overestimated. When yields and affinities are assumed identical for AOO and NOO, using AOO parameters for both (giving them identical R*s), peak [] and [] are identical, and [] is underestimated. With only a difference in yield or the uptake affinity, []:[] is about 1:3 or slightly lower, as quantified above with the point balance, and consistent with observed [] and observed ratios[10,12,13]. This suggests that either the yield or the affinity difference, or a smaller combination of both, may best represent natural assemblages of nitrifiers.

Even deeper in the water column, [] decreases, resembling the observed deep “tail” of the Pn class="Chemical">NM[6]. Observed pelagic [] has been shown to be much lower (nanomolar)[13] than the deep concentration in the model here (tens of nanomolar). Assuming no other sinks for [], we hypothesize that deep [] concentrations still reflect the R*s of the NOO, but of a diverse community. Affinities and efficiencies may vary with clades, such as those with maximum growth rates too low for survival at the PNM, or may be plastic in nature. A slower-growing, efficient, deep NOO clade could deplete to lower levels.

Nitrifier abundances and nitrification rates

Modeled nitrifier abundances also reflect the distinctions between the AOO and n class="Chemical">NOO (Fig. 3b, f). With both yield and affinity differences, NOO abundance is 30-fold lower than AOO. In the model experiment with only the yield difference, NOO abundance is threefold lower than AOO. This threefold difference in abundance is consistent with observations of AOO and NOO cell abundances: a fourfold difference in the abundances of -oxidizing marine group 1 (MG1) (Thaumarchaea) and -oxidizing Nitrospina in the California Current[19], a 1–5-fold lower abundance of Nitrospina-like bacteria compared to amoA gene and MG1 abundances in Monterey Bay[58], and a 1–4-fold difference at Station ALOHA[58]. Thus, observed abundances are quantitatively consistent with observed []:[] in oligotrophic environments, providing compelling evidence for the differences between AOO and NOO metabolisms in natural assemblages.

Rates of and oxidation are identical below 150 m depth for all model experiments (Fig. 3c, g). This matches the lack of consistent differences in observed rates[2], though differences in coastal waters have been documented[41,42]. Thus, differences in the modeled AOO and n class="Chemical">NOO emerge in nutrient distributions and nitrifier abundances, but not in subsurface nitrification rates. This reinforces our understanding that exported organic matter determines the rate of all steps in the sequence of remineralizing metabolisms below the euphotic zone[2,12,13,22].

Vertical mixing affects the magnitude of the PNM

Modeled [] is higher than R* from the peak of the PNM to about 175 m depth. This indicates that vertical mixing is non-negligible there. In other model simulations, particularly in less stratified n class="Chemical">water columns, this effect is stronger, and DIN accumulates to concentrations much higher than R* (Supplementary Fig. 4). In Fig. 3, vertical mixing sweeps cells away from their location of favorable growth so that the NOO cannot sustain a population large enough to draw down to their R*, resulting in accumulation to a higher concentration (see quantitative analysis in Supplementary Note 2).

This influence of mixing at the model PNM is consistent with residence time analyses that find that older exists at the peak Pn class="Chemical">NM[10,21]. Biological turnover will be slower if nitrifiers are unable to metabolize all of the at a particular location. We may thus consider two regimes within the PNM: a transport-influenced upper regime, and a nitrifier-controlled lower regime. This provides an interpretation of the double-peaked character of the PNM in some locations, as an alternative to the standing hypothesis that the upper PNM forms because of reduction by photoautotrophs[9,18].

Predicting global distributions of nitrite and nitrification

Why does accumulate at the subsurface in the subtropics, but also throughout the surface in subpolar regions in Fig. 1? To understand this distribution, we resolve the population dynamics of nitrifiers in a three-dimensional ocean circulation and biogeochemistry model. Building upon previous work[54,59], the model couples the cycles of n class="Chemical">nitrogen, carbon, phosphorus, iron, silica, and oxygen, includes both particulate and dissolved organic pools, and resolves multiple populations of phytoplankton and zooplankton functional types. We introduce the microbial nitrifying and heterotrophic types as described above. In the illustrated global simulation, the affinity of the NOO was increased (though still lower than AOO affinity; see Supplementary Table 1) to match the distribution of the maximum water column [] in the GLODAPv2 compilation (Fig. 5c).

Fig. 5

Observed and simulated maximum nitrite concentration in the oxygenated water column. a From the GLODAPv2 database[60, 61], b from the global model, and c the distributions of both. Only locations that co-occur with O2 concentrations greater than 10 μM are plotted. (Note: a threshold of 50 μM O2 results in visually indistinguishable versions of a and b.) Map generated with Python version 3.5.1, Matplotlib version 1.5.1, and Basemap version 1.0.7[86]

Observed and simulated maximum nitrite concentration in the n class="Chemical">oxygenated water column. a From the GLODAPv2 database[60, 61], b from the global model, and c the distributions of both. Only locations that co-occur with O2 concentrations greater than 10 μM are plotted. (Note: a threshold of 50 μM O2 results in visually indistinguishable versions of a and b.) Map generated with Python version 3.5.1, Matplotlib version 1.5.1, and Basemap version 1.0.7[86]

The global model simulates the distribution of in the transects from the GLODAPv2 database (Fig. 1)[60,61], predicting a PNM throughout the subtropics and the accumulation throughout the upper part of the n class="Chemical">water column poleward of about ±45°. What gives rise to this latter feature? At high latitudes, deep mixing transports fixed nitrogen (including ) to the surface. If photoautotrophs are limited by light or iron, they do not deplete this DIN. Since the R* of the NOO is higher than that of the phytoplankton, surface reflects this higher concentration.

In this case, surface and may be accessible to the - and -limited chemoautotrophs, and they can coexist with light- or iron-limited phytoplankton. The global simulation predicts this coexisting nitrification in the near-surface at high latitudes (Fig. 6b). In the model, this is mostly due to vertical mixing, as inferred from observations in the Southern Ocean[62,63]: DIn class="Chemical">N and nitrifier cells are swept up to the euphotic zone because of deep mixed layers, and the nitrifiers coexist transiently with the phytoplankton. We diagnose that the nitrifiers can also coexist stably (locally) with the phytoplankton at some locations (Supplementary Fig. 5). Seasonal diagnostics of the global model reveal enhanced surface nitrification rates in the winter in both hemispheres as mixed layers deepen and phytoplankton growth becomes more light-limited (Supplementary Fig. 6). This is consistent with observations of increased amoA gene abundances and potential nitrification rates in the winter in surface waters in the coastal Arctic Ocean[36].

Fig. 6

Simulated nitrite oxidation rate. a Water column maximum rate and b surface rate (0–55 m mean). Map generated with Python version 3.5.1, Matplotlib version 1.5.1, and Basemap version 1.0.7[86]

The simulated maximum water column nitrification rates correlate broadly with primary production (Fig. 6, Supplementary Fig. 7), since subsurface nitrification depends on surface production for substrate[2,10,11,13,22]. The range of maximum rates is consistent with the 10–100 nmol n class="Chemical">N L−1 d−1 range of compiled marine nitrification rate measurements[2]. Globally integrated and oxidation in the model are 3190 and 2460 TgN per year, respectively, and integrated primary production is 7150 TgN per year (40.4 PgC per year). As a global average, 17% of oxidation is above 100 m depth in the model, providing a shallow source of for phytoplankton. This is consistent with the understanding that nitrification within the euphotic zone fuels a non-negligible fraction (here, 17%) of primary production[64].

In Fig. 5, we illustrate the maximum [] in the water column where n class="Chemical">oxygen concentration is high (i.e., where O2 > 10 μM, which excludes the domain of anaerobic activity and the secondary maximum). The model captures the lowest values in the subtropical gyres, and the higher values in high latitude and equatorial upwelling regions. The model does overestimate the values in the North Atlantic and Indian Equatorial regions. However, this overestimate disappears when nitrifiers are assumed to have the same temperature sensitivity as the other microbes (Supplementary Note 1). The model also underestimates the highest values in the Pacific Equatorial region. The coarse resolution of the physical model does not resolve the sharpness of equatorial circulation, although the climatology might also be biased by aggregating the effects of eddies, showing high [O2] and [] at the same location that were not measured at the same time. Spatial patterns of mismatches between modeled and observed magnitudes could also indicate unaccounted-for diversity among the nitrifying community, which we discuss below.

Discussion

We have developed a theoretical framework that predicts the locations of nitrification and accumulation broadly. Competition with photoautotrophy explains why nitrification rates often peak at depth. When photoautotrophic phytoplankton and chemoautotrophic nitrifiers compete for or , nitrifiers lose the competition because their metabolisms are much less efficient at using DIN for growth. This set of dynamics should characterize much of the n class="Chemical">N-limited surface ocean.

When phytoplankton growth is limited by something other than nitrogen (e.g. light or n class="Chemical">iron), nitrifiers may locally sustain growth if sufficient DIN is supplied. This characterizes the base of the euphotic zone, and subpolar surface environments during winter, where phytoplankton are light-limited. Surface nitrification rates may also be significant in high-nutrient, low-chlorophyll regions where iron limits phytoplankton growth, if nitrifiers are less iron-limited. Observations suggest that this may be the case for AOO, which require copper rather than iron for redox machinery[3]. Additional model experiments that include an iron limitation to AOO and NOO growth show similar solutions (Supplementary Note 3), though a more thorough analysis of the iron requirements of nitrifying microorganisms is needed.

Future observations could test these hypotheses. While timescales of marine nitrifier growth are too long for nitrification to become competitive at the surface at night (Supplementary Note 4, Supplementary Fig. 10), this could be tested in areas where phytoplankton have been light- or n class="Chemical">iron-limited for weeks or months at a time. Resource ratio theory[47] suggests that the degree of coexisting nitrification and primary production in the euphotic zone reflects the degree of limitation of phytoplankton growth. Less iron or less light availability should allow for a proportionally higher rate of coexisting nitrification.

Additionally, continual supply of , , and nitrifying biomass prevents competitive exclusion by phytoplankton and allows nitrifiers to sustain a transiently coexisting population. This may characterize a steady state near the surface in areas with sufficient vertical transport, such as when mixed layers are deeper than the euphotic zone or when upwelling is significant. This also may characterize coastal areas with high-nutrient injection from river runoff.

Thus, we suggest that light is an indirect control on nitrification in the water column. This can be reconciled with laboratory data showing direct light inhibition[26-28,30,39,40], if long-term exclusion from the surface has reduced the value of photoprotection and facilitated the evolution of photoinhibition in many clades of nitrifiers. We find that the base of the euphotic zone is an optimal location for nitrification in stratified n class="Chemical">water columns: nitrifying chemoautotrophy is outcompeted by photoautotrophy above and limited by the availability of reduced DIN, sourced from the remineralization of organic matter by heterotrophs, below.

We present two reasons for a distinction between the two steps of nitrification that allow for higher [] than [] in oligotrophic environments: (1) energetic constraints suggest that n class="Chemical">NOO require significantly more to sustain their population relative to the amount of required by AOO, and (2) allometry suggests that the specific nutrient affinity of NOO may be lower if they are on average larger than the AOO in a particular community, as recent work has suggested[50]. Both distinctions reduce the competitive strength of the modeled NOO, which cannot deplete to as low of a concentration.

How well do these theoretical distinctions represent real nitrifying populations? The observed abundances of AOO:n class="Chemical">NOO and []:[] in oligotrophic environments are predicted by the distinction in yield alone. However, we know that organisms can develop enzymatic machinery to overcome physiological or energetic constraints. We expect a diversity of yields and affinities to characterize the marine nitrifying community. Here, growth efficiencies were estimated from marine batch cultures grown in initially nutrient-rich conditions, but efficiencies of cells may be significantly different across natural environments. Though the estimates included some mixotrophic cultures (Supplementary Note 5), other metabolic diversity may contribute to energy production[3,4], increasing yields. A contribution from “comammox” bacteria that completely oxidize to is now fathomable, given recent evidence that it may be oligotrophically adapted[65]: measureable -oxidation rates may be lower if comammox metabolizes a portion of the pool. Even more speculatively, a yet-to-be-discovered small, high-affinity NOO population may also exist in the open ocean. Comparisons of marine AOO and NOO affinities and efficiencies could test these hypotheses, and further connect the dots between the proposed mechanisms and natural assemblages.

Though the model here resolves only one bulk type of AOO and n class="Chemical">NOO each with fixed parameters, the framework developed may be useful for linking patterns of diversity among multiple types to biogeochemical patterns. Our model shows that DIN concentrations are sensitive to the assumptions of the hypothesized yield and affinity differences (Fig. 3e–g), and this sensitivity may be exploited to gain insight into the controls on nitrifier diversity. For example, observations show a shift in dominance from ammonia-oxidizing bacteria (AOB) to ammonia-oxidizing archaea (AOA), as well as to different clades of AOA, from coastal to more oligotrophic environments and with an increase in depth[10,13,19]. Since larger AOB have lower affinity for than AOA[37], AOB should have a larger R*, which may contribute to an explanation for the observed decrease in []:[] across the productivity gradient of the California Current[10]. We note, however, that time-varying fluxes of from faster growing organisms closer to the dynamic mixed layer make it less likely that the steady-state approximation should hold for , and so [] may be less predictable than []. Eddy circulation, not resolved in this global model, will also contribute spatial and temporal heterogeneity to the basin-scale patterns simulated here[66].

Even at steady state, factors other than the nitrifier metabolic parameters cause variation in the magnitudes of [] and []. In Eq. (3), the subsistence resource concentration depends also on the loss rate, which varies in space and time according to the population density of predators (or viral lysis). This top–down control, one of the most uncertain parameterizations in ecosystem models, could cause additional differences between subsistence concentrations. Different grazing parameterizations do not qualitatively affect model solutions, though they change the specific depth at which nitrification becomes sustainable. In addition, we found that the modeled PNM was higher than predicted because of vertical mixing. High rates of mixing can drive the actual resource concentration away from the subsistence resource concentration predicted by Eq. (3)[67].

Thus, AOO and n class="Chemical">NOO that have a potential to draw down DIN to very low concentrations (e.g., in batch cultures) can be associated with much higher DIN concentrations where losses and mixing are significant. Our models suggest that this is often the case. This highlights the utility of a dynamic ecosystem model in synthesizing complex interactions. Resolving the fluxes of all of the components results, at times, in unintuitive relationships between the standing stocks of nutrients and biomass.

In the current model, nitrification is sufficient to simulate global distributions of without phytoplankton excretion of . Previous modeling suggests that excretion by phytoplankton at depth can also form the Pn class="Chemical">NM[14]. Laboratory observations[23-25] and observed associations of Prochlorococcus with the PNM[68] support this mechanism. Such excessive reduction by phytoplankton may also reflect the energetic constraints examined here, since the reverse of the sequence of nitrification redox reactions happens within phytoplankton cells. Here we do not present an argument against phytoplankton excretion also contributing to accumulation except to point out that isotopic analysis suggests otherwise in some environments[10,21]. Rather, we provide explanations for accumulation where nitrification contributes significantly to the pool. Observations of nitrifier abundances and ambient nitrification rates are consistent with our hypotheses.

To conclude, we have hypothesized the main controls on the distributions of in oxygenated waters and on the biogeography of nitrifying microorganisms. Additionally, we have presented a dynamic parameterization of chemoautotrophic metabolisms suitable for global ocean biogeochemical models. The model articulates the rapid n class="Chemical">nitrogen cycling by microbial respiration at the base of the euphotic zone, which should aid in sharpening descriptions of export production. The model can also be extended to examine the interactions of anaerobic as well as aerobic nitrogen-cycling metabolic functional types within an aerobic–anaerobic microbial ecosystem, which may provide insight into the fate of fixed nitrogen in anoxic zones.

Methods

Theoretical nitrifier functional type development

Following established methodology[48], three half-reactions combine to form the catabolic and anabolic full reactions for each nitrifier metabolism: (1) the oxidation of or , (2) the reduction of oxygen, and (3) biomass synthesis. The parameter f represents the fraction of electrons from the electron n class="Species">donor that are channeled into biomass synthesis vs. respiration. Combined, these reactions inform a “whole organism” stoichiometry. The resulting correlation of yield y and growth rate yV (Eq. (2)) thus characterizes the overall efficiency of particular metabolism (though tradeoffs between y and growth rate may be important in modifying further characteristics).

We assume a fixed C:N biomass of 5 ± 1, as observed for heterotrophic marine bacteria[69], and that both n class="Chemical">AOO and NOO use reduced nitrogen (i.e., the same oxidation state as and organic nitrogen) as the source of the elemental nitrogen for synthesis. This latter assumption represents the NOO population most realistically[4] and increases our burden of proof to distinguish the two metabolisms energetically, since the electron donor reaction is the sole energetic difference between them. For simplification, and following previous methodology[48], we here neglect the formation of N2O as a byproduct of oxidation, which should have a negligible impact on AOO stoichiometry, given N2O yields per mol N nitrified of < 1%[70].

For the oxidizer (here considering and NH3 interchangeably), the three half-reactions, for generic biomass CHOn class="Chemical">N, and their electron-partitioning coefficients, are:where d normalizes the biomass synthesis reaction to one electron (see definition below). The sum gives the full metabolism for -oxidizing biomass BAOO (ignoring water and lumping bicarbonate into the CO2 pool for simplification), as a function of f:

For the oxidizer, the three half-reactions are:which when summed gives the full metabolism -oxidizing biomass BNOO as:where the requirement of one mole of per mole n class="Chemical">NOO biomass is effectively negligible in all model simulations.

The resulting nitrifier yields y, defined as moles biomass N synthesized per mole DIn class="Chemical">N used, are:

Since we assume that the oxidizer also uses reduced nitrogen, we use the same estimate for d for both functional types. Following previous methodology, d represents the number of electron equivalents that corresn class="Chemical">pond to the oxidation states of the inorganic constituents of that synthesis[48]. Assuming generic microbial biomass composition of C5H7O2N and d = 4(5) + 1(7) − 2(2) − 3(1) gives d = 20.

Estimate of nitrifier efficiency from data

The whole-organism stoichiometries require a constraint on the fraction of electrons donated to biomass synthesis vs. respiratory energy production (f). Smaller f equates to a lower growth efficiency and a lower yield, y. While thermodynamics can provide some theoretical constraints on f, the efficiency of energy production in marine microbes is not well known. Thermodynamics and wastewater treatment studies suggest that f is very similar for the two nitrifying types[48]—the difference in the free energy of the oxidation of to compared to to (+41.65 and +32.93 kJ e− eq−1, respectively, at standard state with pH = 7.0[48]) does not significantly impact theoretical estimates—so we employ published laboratory observations of marine n class="Chemical">AOO and NOO growth to constrain a common f for both metabolisms[37,49,71-73] (Supplementary Table 2).

Yields and were estimated from observations of cell growth on or , and the value of f was inferred for each using Eqs. (6) and (7) (Supplementary Table 2). Some of the observed growth was mixotrophic, which exhibited yields about 10–20% higher than obligate chemoautotrophic growth (Supplementary Note 5). When required, the yield calculations assumed n class="Chemical">nitrogen cell quotas of 0.12 and 1.2 fmol N per cell for the AOO and NOO groups, respectively, with a range of 0.07–0.16 and 0.7–1.6 fmol N per cell contributing to uncertainty in the yields. These quotas are computed from the 10.2 ± 1.1 fg protein per cell content of AOA, as measured[37], and an assumption of a 10-fold larger quota for the NOO, based on the measured minimum 10-fold difference in protein content between AOA and AOB[37] and an assumed similar protein content for AOB and NOB. The nitrogen content of protein was assumed to be 16% by weight, and additional uncertainty was incorporated by considering a range of 10–20% (Supplementary Table 2), giving the above nitrogen quotas. The NOO nitrogen quota can also be independently estimated from the sphereoidical volume of a marine strain of Nitrospina of size 0.3–0.4 μm × 1–3 μm[49]: converting from an average bacterial carbon quota[74] of 0.22 g C cm−3 with a C:N of 5 gives a quota of order 1 fmol N per cell, consistent with our estimate.

Our analysis supports the assumption that f is similar for marine NOO and n class="Chemical">AOO populations in aggregate, though it ranges widely for NOO, and results in an average value on the order of 0.03. For the AOO and the NOO groups, the average yields are (112 ± 32)−1 mol biomass N per mol oxidized, and (310 ± 320)−1 mol biomass N synthesized per mol oxidized, respectively. The average f values corresponding to these yields are 0.030 and 0.032, respectively. A value of 0.03 is about one-fifth of the value of f inferred for wastewater bioreactors[48], perhaps explained by the need for marine organisms in oligotrophic environments to have higher-affinity, energetically-expensive nutrient transport systems.

The stoichiometries presented assume f = 0.03, d = 20 ± 4, with uncertainty in d from the 1 mol/mol s.d. of the C:N[69], where changes in H and O stoichiometry were neglected to give the largest impact on the resulting stoichiometric uncertainty. We examine the sensitivity of the yields to variation in f and d in Supplementary n class="Chemical">Note 6 (Supplementary Fig. 11).

Nitrifier uptake kinetics

Kinetics experiments with cultured -oxidizing archaea Nitrosopumilis provide values for the parameters for the uptake of by n class="Chemical">AOO[37] (Supplementary Table 1), including the information needed to convert to maximum specific uptake rate Vmax (in mol per mol biomass N per day). The model incorporates a conversion from the reported maximum rate[37] of uptake of 24.2 ± 2.23 μmol per mg protein per hour at 30 °C (Supplementary Table 1) with an associated half-saturation concentration of 133 ± 38 nM . Although a maximum rate of 51.9 μmol  per mg protein per hour was reported in batch culture, the growth of the cells was impaired by agitation, and so we do not expect this maximum rate to represent the maximum rate of cells in the ocean. The s.d. of the maximum rate was estimated from the s.d. of the corresponding oxygen uptake (36.29 ± 3.35 μmol O2 per mg protein per hour: (3.35 × 24.2)36.29−1 = 2.23). Assuming a 16% N content of protein gives a specific maximum uptake rate of 50.8 ± 4.68 mol per mol biomass N per day at 30 °C, which is the value of Vmax we use in the model. (See Supplementary Note 1 for discussion of temperature sensitivity.)

This maximum uptake rate corresponds to a maximum per cell nitrification rate of 5.92 ± 0.84 fmol per cell per day (at 30 °C), assuming the measured cell quota for AOA[37] of 10.2 ± 1.1 fg protein per cell per day. This maximum rate is consistent with the measured cellular nitrification rates of enriched cultures of AOA[72,75] of about 2 and 2–4 fmol N per cell per day.

Measurement of the kinetics of natural assemblies of marine oxidizers[29] show a lower half-saturation constant than Nitrosopumilis (27.2 ± 4.4 vs. 133 ± 38 nM ) with respect to a bulk oxidation rate of 24.9 ± 1.3 nM n class="Chemical">N d−1. We would not expect all measured half-saturation constants to be identical, because they vary with their associated maximum rate, which is why the specific affinity () is the relevant trait[43]. However, we can infer from this comparison that natural assemblages may have a lower maximum rate, and thus that the model here may overestimate nitrification rates in some locations.

Allometric theory and affinity

Known -oxidizing bacteria are significantly larger in size than marine -oxidizing archaea[37,49,50]. We assume the 10-fold difference in volume discussed above. (Note: Recently, observations have suggested almost a fourfold larger cell diameter between dominant n class="Chemical">AOO and NOO types, equating to a 50-fold larger cell volume[50]). We then use established empirical and theoretical allometric relationships[51,52] to predict the kinetic parameters for the NOO relative to those of the AOO. Allometric theory predicts that though the cellular uptake rate should increase with cell size, the specific uptake rate should decrease due to a decrease in the surface to volume ratio: the cellular rate scales with surface area as cell radius r2, volume increases as r3, and so the specific rate scales as r−1. Theory also predicts that the half-saturation concentration increases with cell size. The diffusion-limited cellular uptake rate (which explains the steep slope of the Michaelis–Menton form) increases as r. The cellular affinity equates to the quotient of the cellular uptake rate and K. This suggests that K scales as r2r−1 = r. Together, allometry thus suggests that the specific affinity decreases with cell size as r−2, which is supported by more detailed analysis of nutrient uptake models[43].

The 10-fold larger volume equates to a cell radius of NOO larger than that of AOO by 101/3 = 2.2, and so for the NOO, we estimate a specific uptake rate of about half, and a K of about double that of the AOO. This gives a 4.6-fold lower specific affinity () of the NOO relative to the AOO.

Empirical results linking cell size to affinity[52] show that the affinity of AOA for is of the same order of magnitude for that of picophytoplankton such as Prochlorococcus. For the AOA, using the above specific uptake rate of 50.8 mol per mol biomass N per d and the 133 nM half-saturation concentration, the specific affinity is 382 L per μmol biomass n class="Chemical">N per d. In comparison, the specific affinity for inorganic nitrogen for a cell of diameter 0.6 μm, the average diameter of the Prochlorococcus[55], is of order 100 (calculated from allometric relationships[52] as Vmax(QminK)−1 = 88.2 L per μmol biomass N per d). This literature[52] also suggests that picophytoplankton affinity for in particular may be up to an order of magnitude higher than this: of order 1000. In sum, the theoretical specific affinity for a Prochlorococcus-sized cell for inorganic nitrogen is in the range of 100–1000 L per μmol biomass N per d, bracketing that of the nitrifiers. Thus, couched only in terms of affinity for with respect to uptake (not growth), Prochlorococcus and AOA may be close competitors.

Measurements

Measurements were made in the subtropical North Pacific on n class="Chemical">NEMO (“Nutrient Effects on Marine Microorganisms”) Cruise NH1417 in August and September 2014. Using satellite data, stations thought to be at productive locations were chosen for nitrification rate measurements to increase chances of detectable concentrations and rates. Measurements were taken at stations 44, 58, 70, and 99 (Supplementary Fig. 1; Supplementary Fig. 8). and concentrations were measured using standard colorimetric techniques[76], using a Varian Cary 100 Bio ultraviolet-visible spectrophotometer. Ammonia oxidation was measured using 15 as a tracer. Triplicate samples were spiked with enriched tracer (15NH4Cl, 99%, Cambridge Isotope, 100 nM) incubated for 24 h in 500 mL polycarbonate bottles in the dark and at close to in situ temperatures. A 30 mL sample at the beginning and at the end of the incubation (T0 and T24) was extracted from each bottle, filtered through 0.2 μm pore-sized nylon filters, and frozen. Thawed aliquots were treated with sodium azide to convert all of the to N2O gas[77]. The 14N:15N ratio of the N2O gas was analyzed on an isotope ratio mass spectrometer by the University of California Davis Stable Isotope Facility. Unlabeled carrier (to 1 μM) from laboratory stock solution was added to reach instrument detection limits. A second treatment using to measure the rate of oxidation was unsuccessful because incomplete reduction of the remaining in situ spiked inhibited accurate measurement of the enriched pool (which does not mean that oxidation rates were zero). The isotopic composition of the original samples was computed from the measured isotope ratio using a mass balance with the known carrier concentration and isotopic composition. Ammonia oxidation rate (AOR; nM d−1) was calculated as a function of the atom % 15N (a) of initial and final samples as:where is the atom percent of the pool after spiking (, where is the measured average of background samples), V is the sample volume (L), and Δt is the incubation time (d).

For the Chlorophyll a (Chl a) measurements, sean class="Chemical">water was filtered onto 25 mm Whatman GF/F filters (0.7 μm nominal pore size). Filters were extracted in the dark in 5 mL of 90% acetone for 24 h at +3 °C prior to measurement[78] on a Turner Designs TD-700 fluorometer calibrated with pure Chl a (Sigma-Aldrich). Abundances of thaumarchaeal ammonia monooxygenase aubunit A (amoA) genes were measured with quantitative PCR.

Ecosystem model

Nine state variables are resolved as concentrations of n class="Chemical">nitrogen: the biomass of five functional type populations (ammonia-oxidizing organisms BAOO, nitrite-oxidizing organisms BNOO, phytoplankton P, heterotrophic bacteria Bhet, and microzooplankton grazer Z), three inorganic nutrients (, , and ), and organic detritus D. Total nitrogen is conserved in sum over the domain. Supplementary Table 1 lists all parameters, their dimensions, and the default values used in the model. The equations for the nine state variables are as follows, with the substantial derivative notation D/Dt including the diffusive flux as function of the diffusive coefficient K as ∇⋅ for tracer C, and advective fluxes as functions of velocity as ∇ ⋅ :

Phytoplankton growth

Phytoplankton grow as a function of a maximum growth rate μmax (d−1), with limitation by nutrients (γ), modification by temperature (γ), and limitation by light. Light limitation was parameterized using an exponential form as a function of the instantaneous photosynthetic rate Γ (d−1) and the Chl a to n class="Chemical">Carbon ratio θ (g/g)[53,79] as:

Photosynthetic rate Γ was computed as a function of photosynthetically active radiation I, the maximum quantum yield of carbon fixation ϕ (mol C mol−1 photons), and the absorption of light by phytoplankton (m2 (mgChl)−1) representing a mean value over all wavelengths, as:

The Chl:C θ varies with photoacclimation, and is computed using a steady-state solution[79] as:where θmax is a maximum ratio. In the water column model, θ is allowed to reach an arbitrary minimum value of 0.1θmax, which does not affect solutions of any state variables, but does set the minimum modeled Chl a concentration in Fig. 3.

Nutrient limitation is a function of the total concentration of all species of DIN:

The inhibition of and assimilation in the presence of had a negligible effect in the water column model solutions, and so was not included (though it is included in the global model). The specific rates of uptake V (d−1) of each DIn class="Chemical">N species by the phytoplankton type are resolved as:

Values for the maximum growth rate and the half-saturation constants were computed as functions of cell size following data-based allometric relationships[52]. Cell volume v was converted from diameter assuming a spherical cell, and traits calculated with generic form:

For the picophytoplankton functional type (P in the water column model), a and b are chosen for small cells[80], with v computed assuming the average Prochlorococcus cell diameter of 0.6 μm[55]. For the larger additional phytoplankton types in the global model, a and b are representative of larger phytoplankton cells using published values[52,54].

The effective half-saturation constants for DIN uptake with respect to μmax were calculated by conversion from the allometric relationships[52] for half-saturation constants K with respect to maximum uptake rate Vmax and minimum cell quota Qmin as[57]:with published values of a and b for uptake and Qmin[52]. The half-saturation constant for was assumed to be half that of the more oxidized species: . Final values for all model parameters are listed in Supplementary Table 1.

Heterotrophic bacterial growth

The bacterial heterotrophic functional type grows as a function of organic matter (detritus D) as:where y partitions consumption of D into biomass synthesis (as y) and its remineralization into (as y − 1). The growth efficiency y is assumed as the average bacterial growth efficiency of 0.14 (s.d. 0.14) (mol biomass synthesized mol D−1) for the open ocean[81]. The maximum uptake rate Vmax and half-saturation constant K are best estimates that constrain the model heterotrophic bacterial growth rate to about 0.1 d−1, matching the average bulk bacterial growth rate[82].

Grazing

A single zooplankton grazer consumes picophytoplankton (solely the Prochlorococcus-like type in the global model), heterotrophic bacteria, and the chemoautotrophic nitrifiers. This parameterization was chosen based on the assumption that grazing preferences are predominantly governed by size[83], and that Prochlorococcus and the microbial functional types are all roughly the same size. (Thus, there is a possibility that differences in size between AOO and n class="Chemical">NOO may influence this top–down control.) The total amount of grazing is calculated as a saturating function of total prey biomass, giving the rate of grazing g (L mol−1 d−1), which is then multiplied by the biomass of each prey in Eqs. (13–16) and summed for total consumption by Z in Eq. (17), as:with maximum grazing rate gmax and half-saturation K. Though these values are uncertain, a recent compilation[83] suggests that the two values are of the same order of magnitude, constraining their influence on g. Here, we assume gmax = K = 1, which is of the order of magnitude for these values used in previous marine ecosystem models[54].

The zooplankton grazer also excretes via respiration as governed by its growth efficiency ζ. Studies of zooplankton nitrogen growth efficiency is variable, and that 0.5 is a reasonable mid-range value[83].

Temperature

Other than the nitrifier growth[56], all microbial growth, grazing, and mortality rates are represented as a function of temperature (non-dimensional γ) using a formulation that follows the Arrhenius equation[54] as:where T is the ambient temperature, T0 is a reference temperature, A regulates the temperature modification, and τ normalizes the maximum value. (See Supplementary Note 1 for discussion of temperature sensitivity).

Calculation of R*s

R*s are calculated with Eq. (3) for AOO and n class="Chemical">NOO, where is the maximum growth rate. For phytoplankton R*s, the maximum growth rate is considered to be its maximum light-limited growth rate μlight, from Eq. (18) without the limitation by nutrients, as:which is acquired from the steady-state solutions. Thus, R* for P is calculated as:

For both phytoplankton and nitrifiers, the loss rate L is calculated as the sum of respiratory or maintenance losses m (d−1) and the resulting steady-state gZ (d−1), aswhere g (L mol−1 d−1) is the above grazing rate (equation 25) and Z (mol L−1) is the population density of predators. In the model here, the loss rate L is the same for all of the microbial prey populations, since m is assumed constant and g depends on the sum of all prey biomass. (This specific rate L is multiplied by the biomass of each prey B in Eqs. (13–16) to give biomass-dependent losses.)

1D water column physical environment

In the water column model, the mixed layer was imposed by varying the vertical diffusion coefficient K with depth, from a maximum Kmax at the surface to a minimum Kmin with a length scale of zmld. The fixed (no flux) boundary conditions result in some accumulation of D at the bottom of the 2000 m domain, conceptually representing a sediment layer. To smooth over numerical error, vertical mixing was allowed to increase there with a 100 m length scale, simulating a bottom boundary mixed layer. K (m2 s−1) is thus calculated at cell faces as:where H is the height of the domain (2000 m).

Light energy I decreases with depth according to the attenuation coefficients for water kw and for n class="Chemical">chlorophyll kChl, following a previous approach[84] as:where TChl is the sum of the concentrations of chlorophyll, and kChl is an upper estimate of the absorption by chlorophyll to account for additional absorption by colored dissolved organic matter[54]. Iin is the incoming irradiance, which is 0.5Imax, or for resolution of the daily cycle, Iin(t) = 0.5Imax(cos(ωt) + 1) where ω = 2π d−1.

A temperature curve was fit to the mean observations from the four stations sampled on Cruise NH1417 (Supplementary Fig. 9) with an exn class="Chemical">ponential form:with temperature T in °C.

The illustrated domain was 2000 m in height, with 5 m vertical resolution. Equations were integrated forward in time using the 4th order Runge-Kutta method until equilibrium (i.e., solutions independent of initial conditions). Advection for the sinking organic matter pool was carried out using the QUICK advection scheme, consisting of a linear interpolation between points weighted by an upstream second order curvature, resulting in third order accuracy. Fluxes were calculated at the faces of each grid cell, and concentrations at the centers.

1D model ensemble with parameter uncertainty

For the ~1000 model solutions in Fig. 3, we vary the values of the biological parameters in the illustrated model to communicate their uncertainty (see Supplementary Table 1 for the ranges of values). Values for the kinetic parameters of the nitrifiers were assigned randomly from the normal distributions from results of the oxygen kinetics experiments[37]. The relative larger size of the n class="Chemical">NOO as compared to AOO was varied from a factor of 1 to a factor of 20, thus, impacting the relative values of Vmax and K by a factor of 11/3 to 201/3. Because of their large uncertainty, the following parameters were sampled randomly over the linear range of ±50% of the default value: the maximum growth rate μmax and nutrient half-saturation constants for phytoplankton (with computed accordingly); the parameters governing the consumption and respiration organic matter by heterotrophic bacteria: Vmax, K, and y; the grazing parameters gmax and K, mortality rates m and m, and zooplankton efficiency ζ.

Global model

The present configuration of the 3D MITgcm biogeochemical model resolves a total of six phytoplankton populations with parameters that represent the traits of the following six functional types: diatoms, picoplankton, n class="Chemical">diazotrophs, coccolithophores, and other large and other small phytoplankton. Four zooplankton types graze on the phytoplankton: one each specifically on the picoplankton (which also consumes the six introduced microbial types below), the other small phytoplankton, and coccolithophore types, and the fourth grazes on the diatom, the other large phytoplankton, and diazotroph type. The three-dimensional ocean circulation state estimate (the ECCO-GODAE state estimate) is from the configuration of the MITgcm as constrained by observations[85], and has a horizontal resolution of 1° × 1° and 23 levels of vertical resolution, from 10 m at the surface to 500 m at depth. The model was numerically integrated until rates of microbial activity equilibrated throughout the thermocline (200 years in the illustrated model).

In addition to the six phytoplankton types, six microbial metabolic functional types are included in the global ecosystem model, and are responsible for all n class="Chemical">organic matter remineralization, nitrification, and denitrification. Particulate and dissolved organic matter (POM and DOM) are consumed and subsequently remineralized by the aerobic heterotrophic bacterial functional type (as in the water column model), an anaerobic nitrate-reducing () heterotrophic type, and an anaerobic denitrifying () heterotrophic type. The growth of each heterotrophic type is limited by the sum of POM and DOM, and POM and DOM are both taken up, weighted by the fraction of the limitation imposed by each as a function of the local concentration (analogous to the uptake of the three species of DIN by phytoplankton). Redfieldian C:N:P:Fe stoichiometries of bacterial types and demands are constant. The two aerobic nitrifier types (AOO and NOO) are included as in the water column model. A chemoautotrophic anammox functional type () is also included. The stoichiometries for the three anaerobic metabolic functional types are described in Supplementary Note 7.

The depletion of oxygen and switch from aerobic to anaerobic respiration occurs dynamically in the model, following previous parameterization[45]. The lower n class="Chemical">organic matter yield for the anaerobic types results in the competitive exclusion of the anaerobic heterotrophs in oxygenated environments. The anammox type is likewise excluded from oxygenated environments by the aerobic nitrifiers. The assumed stoichiometries of anaerobic metabolisms results in [] greater than 1 μM, simulating a secondary maximum, for which analysis is beyond the scope of the investigation here.

Code availability

Fortran code for the ecosystem model and the one-dimensional water column configuration is available at https://github.com/emilyzakem/eco-nitrify. The global model (the MITgcm) is available at http://mitgcm.org and the ecosystem comn class="Chemical">ponent including the bacteria is available from git://gud.mit.edu/gud1.

Data availability

Supplementary Data 1 contains all data presented from Cruise NH1417.