Microbial ecosystem dynamics drive fluctuating nitrogen loss in marine anoxic zones.

The dynamics of nitrogen (N) loss in the ocean's oxygen-deficient zones (ODZs) are thought to be driven by climate impacts on ocean circulation and biological productivity. Here we analyze a data-constrained model of the microbial ecosystem in an ODZ and find that species interactions drive fluctuations in local- and regional-scale rates of N loss, even in the absence of climate variability. By consuming O2 to nanomolar levels, aerobic nitrifying microbes cede their competitive advantage for scarce forms of N to anaerobic denitrifying bacteria. Because anaerobes cannot sustain their own low-O2 niche, the physical O2 supply restores competitive advantage to aerobic populations, resetting the cycle. The resulting ecosystem oscillations induce a unique geochemical signature within the ODZ-short-lived spikes of ammonium that are found in measured profiles. The microbial ecosystem dynamics also give rise to variable ratios of anammox to heterotrophic denitrification, providing a mechanism for the unexplained variability of these pathways observed in the ocean.

Bioavailable nitrogen (n class="Chemical">N) is a key macronutrient that limits the rates of biological activity. In the ocean, the concentration of nitrate (NO3−), the major form of bioavailable N, is reduced by anaerobic reduction to biologically inert N2 gas within small subsurface O2-deficient zones (ODZs) (1). The volumetric rate of N removal within these zones is limited by the downward flux of organic matter from sinking particles (2). In turn, ODZ volumes are strongly dependent on the regional O2 content of the thermocline in which they reside (3). Variations in climate have major impacts on the supply of O2 and organic matter to the ODZ, driving changes in the magnitude of N removal across a wide spectrum of timescales, from months to millennia (4–6).

Microbial community structure also plays a major role in regulating N2 gas production. Anaerobic processes, such as anammox and heterotrophic denitrification, can tolerate up to micromolar amounts of n class="Chemical">O2, allowing them to coexist with aerobic nitrifying microbes, which become limited by O2 only at nanomolar concentrations (2, 7–14). Because both anaerobic and aerobic metabolisms utilize the key N-cycle intermediates ammonium (NH4+) and nitrite (NO2−) as substrates, their coexistence results in resource competition whose outcome is determined by nanomolar variations in O2 (15). When nitrification is dominant, the reoxidation of partially denitrified NO2− to nitrate (NO3−) reduces the magnitude of N2 production and increases O2 consumption; when aerobic nitrifiers are excluded by O2 scarcity, NO3− is efficiently reduced all of the way to N2 (15, 16). Here we demonstrate that resource competition between aerobic nitrifiers, anaerobic denitrifiers, and anammox bacteria can also lead to regional-scale temporal variability in the rates of N and O2 cycling, even with constant physical fluxes of O2 and organic matter into the ODZs.

To examine the role of microbial interactions in the dynamics of fixed N loss, we analyzed a microbial ecosystem model (15) embedded within an ocean general circulation model (17, 18). The steady three-dimensional ocean circulation is optimized to fit tracer observations (temperature, salinity, n class="Chemical">radiocarbon, and CFC-11), implying realistic ventilation rates and pathways of the ODZs (19). We focus on the world’s largest ODZ, in the eastern tropical North Pacific (ETNP) (20), by restricting the boundaries of the model from the equator to 35° N, the coast to 180° W, and the surface ocean to 2,000-m depth. Observed annual mean concentrations of O2 and NO3− (21) are transported into the domain at its open boundaries to ensure their realistic supply to the ODZ region. The circulation does not vary over time, leaving microbial ecosystem dynamics as the sole source of temporal variability.

The microbial ecosystem model simulates the biomass of four microbial functional groups and the biogeochemical cycles of N and n class="Chemical">O2 (15). In the surface ocean, phytoplankton produce dissolved organic nitrogen (DON) and sinking organic particles from inorganic N (NH4+, NO2−, NO3−). DON is remineralized by heterotrophic bacteria using O2, or multistep denitrification (reduction of NO3− to NO2−, then to N2) below a critical O2 threshold (O2crit). The NH4+ and NO2− released by heterotrophs is used by autotrophs: slow growing and O2-inhibited anammox bacteria, or aerobic archaea and bacteria that either perform NH4+ or NO2− oxidation with nanomolar O2 sensitivities. Autotrophs assimilate NH4+ from seawater for growth. Because the C:N ratio of bacterial biomass (6.8 ± 1.2) (22) matches that of organic matter within the ODZ (6.8) (2) heterotrophs satisfy their nutrient demand via NH4+ remineralized from DON (23). DON is released by phytoplankton, sinking particles, and all microbial populations during mortality.

In previous work (15), we assessed the model fit to observed long-term mean (climatological) fields of O2 and n class="Chemical">NO3−, and profile compilations of NH4+, NO2−, and biologically produced N2 gas (N2xs) from within the ETNP. These data reflect the characteristic vertical profiles of key chemical indicators of the metabolic status of the ODZs: subsurface maxima in NO2− and N2xs, reduced accumulations of NO3−, and nanomolar levels of NH4+ and O2 (). To constrain uncertainty in model parameters, we varied microbial growth, mortality, and nutrient affinities over two orders of magnitude, spanning values observed in laboratory cultures and process studies (), and compared the resulting simulated profiles to the observations (). Of the 90 parameter combinations tested, half reproduce all observed chemical profiles simultaneously, implying a realistic balance of physical and biological fluxes of N and O2. This ensemble of model simulations that reproduce the data are used for further analysis and to quantify the sensitivity of our main results.

The simulated rate of regional N loss aligns with geochemical estimates based on measurements of the accumulation of n class="Chemical">N2xs, the deficits of NO3−, and its isotopes (19, 24), but fluctuates strongly over time (Fig. 1) despite the steady rates of ocean circulation. These fluctuations are not caused by changes in the flux of organic matter or the physical supply of O2 to the ODZ region, which are stable (). The fluctuations persist across a wide range of physiological and ecological assumptions: regardless of the precise O2 sensitivities of the microbial populations (yellow, red, and green lines in Fig. 1); with and without inclusion of dissimilatory NO3− reduction to NH4+ (DNRA) (25) (blue line in Fig. 1); whether heterotrophic denitrification is represented as a facultative or obligate process (1, 26) or if its steps are mediated by a single or multiple populations (27) (). Fluctuations in N loss are found under all ecosystem model parameter combinations that satisfy the available tracer data constraints (). Their amplitude is large relative to time-mean rates, averaging 43% ± 35% (SD) on regional scales and 233% ± 123% (SD) at the locations where fluctuations occur. While the regionally integrated N-loss rate lacks a characteristic frequency, local rates of N loss vary through semiregular oscillations (Fig. 1). The complex fluctuations in the regional-scale N loss (Fig. 1) thus arise from the integration of the many localized oscillators with distinct periods, phasing, and amplitudes.

Fig. 1.

Time series of unforced variability in the regional and local rates of N loss from the ODZ of the eastern tropical North Pacific. (A) Rates (1 Tg N y−1 = 1012 g N y−1) are spatially integrated across the ETNP in the standard model simulation (black) and sensitivity cases (gray and colors). Fluctuations occur regardless of physiological or ecological uncertainties (): whether O2 tolerances of anaerobes are 1 µM (yellow) or ≥10 µM (red), if the two steps of nitrification (NH4+ and NO2− oxidation) have different nanomolar O2 sensitivities (green), or if an additional metabolism [dissimilatory nitrate reduction to ammonium (DNRA)] is incorporated into the model (blue). They also hold across wide ranges in other microbial ecosystem parameters (gray, ). (B) Time series of local rates of N loss in locations with representative ecosystem oscillations (12°N, 90°W at 100 m and 25°N, 113°W at 400 m).

Time series of unforced variability in the regional and local rates of N loss from the ODZ of the eastern tropical n class="Chemical">North Pacific. (A) Rates (1 Tg N y−1 = 1012 g N y−1) are spatially integrated across the ETNP in the standard model simulation (black) and sensitivity cases (gray and colors). Fluctuations occur regardless of physiological or ecological uncertainties (): whether O2 tolerances of anaerobes are 1 µM (yellow) or ≥10 µM (red), if the two steps of nitrification (NH4+ and NO2− oxidation) have different nanomolar O2 sensitivities (green), or if an additional metabolism [dissimilatory nitrate reduction to ammonium (DNRA)] is incorporated into the model (blue). They also hold across wide ranges in other microbial ecosystem parameters (gray, ). (B) Time series of local rates of N loss in locations with representative ecosystem oscillations (12°N, 90°W at 100 m and 25°N, 113°W at 400 m).

The oscillations are also evident in aerobic metabolic rates, which together with the changes in N loss, drive large-scale fluctuations in the concentrations of n class="Chemical">O2, NH4+, and NO2−, (Fig. 2). Fluctuations are strongest at the edge of the ODZ’s anoxic core, in a “suboxic” zone, where the full diversity of simulated microbial populations coexist—aerobic NH4+ and NO2− oxidizers as well as autotrophic anammox bacteria and heterotrophic denitrifiers (15). In the anoxic core of the ODZ, where aerobic metabolisms are excluded, the chemical environment, the resident microbial populations, and their metabolic rates are relatively stable over time. The coincidence of variability in zones of long-term nitrifier–denitrifier coexistence implies that the oscillations are driven by interactions between these microbial groups. Indeed, if the nitrifiers are separated from autotrophic anammox and heterotrophic denitrifiers by imposing nonoverlapping O2 thresholds, oscillations do not arise in the model simulations ().

Fig. 2.

Spatial distribution of ecologically driven oscillations within the ODZ. The spatial distribution of oscillation amplitudes (colors) is shown along a zonal cross-section through the model ODZ (20−28°N). Oscillation amplitudes are computed as the difference between maximum and minimum values over a 10-y simulation for (A) the N loss rate, (B) the O2 consumption rate, and the concentrations of (C) NH4+ and (D) NO2−. Variability is overlain by time-mean concentrations of O2 (in µM; black contours). Gray shading denotes the western coastline of North America.

Spatial distribution of ecologically driven oscillations within the ODZ. The spatial distribution of oscillation amplitudes (colors) is shown along a zonal cross-section through the model ODZ (20−28°N). Oscillation amplitudes are computed as the difference between maximum and minimum values over a 10-y simulation for (A) the n class="Chemical">N loss rate, (B) the O2 consumption rate, and the concentrations of (C) NH4+ and (D) NO2−. Variability is overlain by time-mean concentrations of O2 (in µM; black contours). Gray shading denotes the western coastline of North America.

The mechanism of these oscillations derives from a fundamental ecosystem dynamic: consumption of O2 by aerobic microbes provides an advantage for anaerobes, but their niche cannot be sustained against the physical n class="Chemical">O2 supply without intermittent dominance of the aerobes. In the model ODZ, the consumption of O2 by NO2− oxidation (∼41 Tg O2 y−1) vastly outweighs NH4+oxidation (∼4.9 Tg O2 y−1); NH4+oxidation thus plays little role in the oscillatory dynamic. The complete ecological sequence of the oscillation is illustrated by the phase diagram of NH4+ and O2 at a single point in space (Fig. 3). When O2, NH4+, and NO2− are plentiful, NO2− oxidizing bacteria experience net population growth (location Fig. 3 , i). Their metabolic rate exceeds the physical O2 supply and depletes the available O2 and NO2− (Fig. 3 , ii). The loss of O2 promotes anaerobic metabolisms, but the loss of NO2− also depletes the energy available for heterotrophic growth fueled by denitrification. The NO2− oxidizers can short circuit complete heterotrophic denitrification to N2 because of their higher efficiency of NO2− utilization (15), which is required by the model to reproduce the observed distribution of NO2− within oxic and anoxic waters (). In contrast, the depletion of NO2− has little effect on anammox bacteria because they are generally limited by NH4+ in the model, consistent with rate measurements from the ODZs (1, 28). Thus, as NO2− is drawn down by oxidation, the decline of heterotrophic denitrification relative to anammox (Fig. 3, iii) depletes NH4+ to levels that, in turn, limit the NO2− oxidizers, slowing their rate of O2 and NO2− utilization (Fig. 3, iii). The cessation of O2 consumption allows its concentration to be gradually replenished by the physical supply, while NO2− simultaneously accumulates due to an excess of NO3− reduction over NO2− oxidation (Fig. 3, iv). O2 accumulation selects against anammox (Fig. 3, iv), while NO2− accumulation fuels a rapid burst of N loss through heterotrophic denitrification (Fig. 3, i). The NH4+ liberated from DON during this denitrification pulse restores it to levels that sustain NO2− oxidizer growth, a condition that again favors net O2 consumption, and the oscillation starts anew.

Fig. 3.

Dynamics of the ecosystem oscillation. The oscillation of key ecosystem variables is shown in the phase space of NH4+ and O2, from a representative location at the suboxic boundary between the anoxic zone and the oxic ocean (i.e., same as in Fig. 1 at 400 m). Time proceeds in the counterclockwise direction, indicated by spiraling arrows. NH4+ and O2 levels are colored by (A) the concentration of NO2− (µM), (B) the rate of O2 consumption by NO2− oxidation (µM O2 y−1), (C) the contribution of anammox to total N2 production (ƒamx), and (D) the rate of total N2 production (µM N y−1). Light colors are always either low concentrations or low rates of activity. Straight arrows in A identify the dominant process driving changes in NH4+ and O2 during each phase of the cycle. Locations i–iv marked on the phase diagrams are described in the text.

Dynamics of the ecosystem oscillation. The oscillation of key ecosystem variables is shown in the phase space of NH4+ and n class="Chemical">O2, from a representative location at the suboxic boundary between the anoxic zone and the oxic ocean (i.e., same as in Fig. 1 at 400 m). Time proceeds in the counterclockwise direction, indicated by spiraling arrows. NH4+ and O2 levels are colored by (A) the concentration of NO2− (µM), (B) the rate of O2 consumption by NO2− oxidation (µM O2 y−1), (C) the contribution of anammox to total N2 production (ƒamx), and (D) the rate of total N2 production (µM N y−1). Light colors are always either low concentrations or low rates of activity. Straight arrows in A identify the dominant process driving changes in NH4+ and O2 during each phase of the cycle. Locations i–iv marked on the phase diagrams are described in the text.

The ecosystem oscillations predicted here arise in a completely steady physical environment, but the supply of organic matter and O2 to the ODZ exhibit strong temporal variations in the real ocean. We tested the impact of physical variability on the intrinsic ecosystem oscillations by first imposing empirically derived seasonal fluxes of organic particles and then stochastic changes in the rates of ocean circulation (). Ecosystem-driven variability persists, and is even amplified, in the presence of these external forcings, suggesting the oscillations would act as a strong source of variability in the natural environment.

Top-down ecological controls on microbial populations also have the potential to limit fluctuations caused by resource competition. We represented grazing losses in the model by applying a quadratic mortality term to all microbial populations, assuming predation is unselective (). We varied the intrinsic grazing rate by an order of magnitude and find that while the variance in regional N loss is unchanged under weak grazing, under strong grazing the variance is decreased by an order of magnitude (). However, adding this strong grazing term causes an unrealistic build up of n class="Chemical">NH4+ concentration in the anoxic core of the ODZ (). By reducing the biomass of the slow-growing anammox bacteria, grazing lessens the main sink of NH4+ in these zones, allowing it to accumulate to persistently high concentrations. The observed distribution of NH4+ therefore does not support the grazing rates needed to stabilize the ecosystem oscillations.

The distribution of NH4+ within the ODZ provides a unique and detectable geochemical signature of the microbial oscillations (Fig. 4). Over the course of the oscillation, shifts in the balance of n class="Chemical">NH4+ sources and sinks lead to its temporary accumulation within the ODZ, at levels up to ∼10 times the measurable detection limit of the most sensitive technique (∼10–15 nM)—the orthophthalaldehyde (OPA) method. These NH4+ spikes are short lived, however, occurring only ∼5% of the time throughout the ODZ (O2 < 5 µM), such that the average model concentration of NH4+ remains below detection. We looked for this potential signature of the oscillation by analyzing 18 depth profiles of NH4+ from the ETNP measured using the OPA method (). Consistent with the predicted time-mean NH4+ concentrations, the average observed concentration of NH4+ in waters with O2 < 5 µM falls below detection. However, in ∼8% of these measurements, we find NH4+ concentrations exceeding this detection limit, consistent with the frequency predicted by ecological oscillations. NH4+ measurements made with less sensitive conventional techniques suggest spikes are also present in the eastern tropical South Pacific and Arabian Sea ODZs outside the model domain (e.g., refs. 11, 29, and 30), but a quantitative analysis of these features will require more high precision data.

Fig. 4.

NH4+ depth profiles from the ODZ of the eastern tropical North Pacific in model simulations and observations. Depth profiles were sampled monthly over the course of a year in the standard model simulation (pink circles) and measured on a cruise to the ETNP in 2012 (black diamonds, ). NH4+ exceeds the detection limit (∼15 nM) ∼5% of the time in the model simulation and in ∼8% of observations at O2 < 5 µM, but on average is below detection in both. Model and observed NH4+ values below 15 nM are set to this detection limit. Diel vertical migration depth for the ETNP is plotted (mean indicated by line, SD by shading) (32). The time-dependent NH4+ profiles are also shown from a model simulation with a data derived seasonal cycle of net primary production (NPP), but weak internal oscillations (violet). Seasonal fluctuations in the supply of organic matter to the ODZ cannot produce the magnitude of NH4+ spikes implied by the observations.

NH4+ depth profiles from the ODZ of the eastern tropical n class="Chemical">North Pacific in model simulations and observations. Depth profiles were sampled monthly over the course of a year in the standard model simulation (pink circles) and measured on a cruise to the ETNP in 2012 (black diamonds, ). NH4+ exceeds the detection limit (∼15 nM) ∼5% of the time in the model simulation and in ∼8% of observations at O2 < 5 µM, but on average is below detection in both. Model and observed NH4+ values below 15 nM are set to this detection limit. Diel vertical migration depth for the ETNP is plotted (mean indicated by line, SD by shading) (32). The time-dependent NH4+ profiles are also shown from a model simulation with a data derived seasonal cycle of net primary production (NPP), but weak internal oscillations (violet). Seasonal fluctuations in the supply of organic matter to the ODZ cannot produce the magnitude of NH4+ spikes implied by the observations.

A transient accumulation of NH4+ within the ODZ might also be expected from excretion at depth by vertically migrating zooplankton and micronekton (31, 32). Measured n class="Chemical">NH4+ spikes occur up to 100–300 m below the mean depth of diel vertical migration recorded for this region (line and shading in Fig. 4). In contrast, elevated NH4+ within the ODZ is found over a similar depth range to where ecological oscillations occur in the model. Temporary spikes of NH4+ could also arise from transitory pulses of sinking organic matter that release NH4+ into the ODZ faster than it can be consumed. We tested whether changes in the particle flux can produce NH4+ spikes, by adding the observed seasonal cycle in net primary production to a model simulation with weak internal oscillations, and thus inherently small pulses of NH4+. In this case, even with forced fluctuations in the supply of organic matter into the ODZ, the predicted time-varying concentrations of NH4+ barely exceed the measured detection limit at any depth. The measured spikes in NH4+ therefore support strong nonequilibrium ecosystem behavior.

Ecological oscillations within the ODZ have direct consequences for the fraction of total N loss that derives from anammox (ƒamx) as opposed to heterotrophic denitrification (Fig. 5). The contribution of these metabolic pathways to n class="Chemical">N loss has been observed to vary across and within the ODZs from direct rate measurements in the field, but the causes of these variations remain hotly debated (e.g., refs. 2, 4, and 32). During the course of the oscillation, when NO2− oxidizing bacteria are ascendant, the NO2− that would otherwise be reduced by heterotrophs is reoxidized to NO3−. The suppression of heterotrophic denitrification temporarily allows NH4+-limited anammox to contribute 100% of local N2 production. However, after NO2− accumulates, the rapid bursts of heterotrophic denitrification vastly exceed previous rates of anammox (Fig. 3 ) and thereby dominate total N loss over a complete oscillatory cycle (horizontal lines Fig. 5). These local variations in the balance of N loss processes can temporarily obscure the time-mean gradients in ƒamx across the ODZ (15). Because they occur over an extremely narrow range in the concentrations of O2, NH4+, and NO2−, evaluating this ecological contribution to observed variations in ƒamx will require frequent and high-precision measurements of these chemical abundances and associated metabolic rates.

Fig. 5.

The contribution of anammox to total N2 production (ƒamx) over time. Time series of ƒamx in representative locations across the ODZ, from 115 m to 450 m. At the oxic–anoxic interface (oxycline), ƒamx can vary over wide ranges that temporarily obscure its time-mean gradient (blue, cyan, and yellow lines). Within the secondary NO2− maximum, ƒamx approaches the value of 0.28 and oscillations are weak (orange, red, and gray lines). Solid lines are from the heart of the ODZ, whereas dashed lines are from its margins. Time-mean contributions of anammox to N loss are shown as colored horizontal lines on the Right axis.

The contribution of anammox to total N2 production (ƒamx) over time. Time series of ƒamx in representative locations across the ODZ, from 115 m to 450 m. At the oxic–anoxic interface (n class="Chemical">oxycline), ƒamx can vary over wide ranges that temporarily obscure its time-mean gradient (blue, cyan, and yellow lines). Within the secondary NO2− maximum, ƒamx approaches the value of 0.28 and oscillations are weak (orange, red, and gray lines). Solid lines are from the heart of the ODZ, whereas dashed lines are from its margins. Time-mean contributions of anammox to N loss are shown as colored horizontal lines on the Right axis.

Oscillatory behavior is a common feature of idealized ecosystem models with multiple interacting populations (33, 34), but is rarely shown to persist in realistic representations of the environment such as a three-dimensional ocean circulation. Intrinsic ecosystem oscillations provide a mechanism to generate variations in marine microbial community structure and N and n class="Chemical">O2 cycling, which are often ascribed to externally forced changes in physical and chemical conditions. Because these oscillations lack spatial coherence and power at decadal and longer timescales (Fig. 1), they are unlikely to explain large-scale decadal variations in N loss (5). However, dynamics such as these may be pervasive beyond the ODZs, occurring wherever the physical supply of resources selects for a microbial community that over time undermines its own ecological niche by shifting the chemical environment to temporarily favor the growth of its competitors or degrade the growth of its facilitators.