Microbial feedbacks optimize ocean iron availability.

Iron is the limiting factor for biological production over a large fraction of the surface ocean because free iron is rapidly scavenged or precipitated under aerobic conditions. Standing stocks of dissolved iron are maintained by association with organic molecules (ligands) produced by biological processes. We hypothesize a positive feedback between iron cycling, microbial activity, and ligand abundance: External iron input fuels microbial production, creating organic ligands that support more iron in seawater, leading to further macronutrient consumption until other microbial requirements such as macronutrients or light become limiting, and additional iron no longer increases productivity. This feedback emerges in numerical simulations of the coupled marine cycles of macronutrients and iron that resolve the dynamic microbial production and loss of iron-chelating ligands. The model solutions resemble modern nutrient distributions only over a finite range of prescribed ligand source/sink ratios where the model ocean is driven to global-scale colimitation by micronutrients and macronutrients and global production is maximized. We hypothesize that a global-scale selection for microbial ligand cycling may have occurred to maintain "just enough" iron in the ocean.

Iron, nitrogen, carbon, and other elements in the surface ocean are consumed by microbes, leading to a flux of sinking and subducted organic matter that is respired in deeper waters, ultimately releasing those nutrients back into dissolved (inorganic) form. Ocean circulation returns resource-rich deep waters back to the surface layer, maintaining global ocean productivity and sustaining a biologically mediated reservoir of inorganic carbon at depth, drawing down atmospheric .

Oceanic primary production may be constrained by one of several key resources. Nutrient addition experiments have shown that most of the surface ocean is locally limited by the availability of either fixed nitrogen or iron (1–6). On a global scale, the strength of biological ocean carbon storage can be determined by the efficiency of macronutrient use, the fraction (0 to 1) of dissolved inorganic nutrients in the abyss resulting from remineralization of sinking and subducted organic matter (7, 8). Low efficiency indicates that macronutrients in the deep ocean arrive there by physical transport, while high efficiency indicates that biological export dominates. Evaluated from observations, modern macronutrient use efficiency is 0.36 (9), confirming that the biosphere is prevented from fully consuming macronutrients at the surface. Indeed, iron availability is the primary limiting factor for production over about half of the surface ocean (10–12), and therefore the global ocean might be broadly considered colimited by iron and macronutrients.

Although plentiful on Earth, ferric iron, Fe(III), has trace concentrations in the ocean because of its low solubility in oxygen-rich waters [0.01 nM (13)]. Due to rapid precipitation or scavenging, iron has a short residence time compared to macronutrients, around 100 y (14, 15), so older upwelling waters are typically depleted in iron but macronutrient-rich, leading to iron limitation in the southern, equatorial, and sub-Arctic Pacific “high-nutrient, low-chlorophyll” oceans.

The conclusion that oceanic production is colimited by macronutrients and iron is also supported by modeling studies exploring responses to changing iron supply (e.g., refs. 10 and 16–21). Stimulated by interest in iron fertilization to mitigate anthropogenic emissions or as a driver of past climate, models have generally found net global productivity responds only moderately or weakly to changing iron inputs. High latitudes appear to be close to local limitation by iron and light, while, globally, relieving iron limitation in one region leads to compensating macronutrient limitation in adjacent areas (22). Thus, while observations often suggest local iron limitation (1–3), the complex dynamics of global simulations highlight an interplay between different resources that shape the patterns of productivity in the oceans.

Why is the global ocean poised at a state of colimitation between iron and macronutrients? We argue that it is because marine organisms can modify the residence time of iron in the oceans. Most oceanic dissolved iron (<0.2-m filter) is associated with organic chelating ligands (23–31), organic molecules that bind to iron and act as a refuge from loss by scavenging and precipitation. The specific nature of ligands is difficult to characterize (32), with bulk abundance and conditional stability coefficients typically inferred electrochemically (e.g., refs. 23–26 and 28–38). These studies reveal multiple classes of ligand with different binding affinities, concentrations, and distributions. This spectrum of organic chelators includes, but is not limited to, siderophores, which are strong-binding molecules secreted by microbes to competitively capture iron (39–41), and weaker-binding molecules that have an affinity for iron, like exuded polysaccharides, porphyrins, degraded protein remnants, and humic substances (23–27, 38, 42–46). Ligand production is likely heterogenous in space and time: Siderophore synthesis occurs under iron stress (40), but elevated ligand levels have been documented following experimental iron addition (24, 44, 47, 48), suggesting other ligand production pathways.

The Ligand–Iron–Microbe Feedback

Given the organic origin of iron-binding chelators, we hypothesize that a positive feedback between microbial activity, ligand abundance, and iron availability could emerge (Fig. 1). Iron is supplied by dust deposition, sediment mobilization, and hydrothermal activity. In a hypothetical watermass where ligand abundance is initially very low, Fe(III) is largely insoluble, but a small population of microbes subsist. Their turnover produces ligands such as siderophores, excreted organic carbon, or chelating detritus. Greater ligand abundance retains more iron in seawater, incrementally relieving iron limitation, promoting further biological production, and so on. Eventually, other requirements, such as macronutrients, become limiting (Fig. 1), and additional iron no longer increases productivity. At appropriate time and space scales, the global ligand pool is regulated, supporting “just enough” iron to match availability of other resources redistributed by ocean circulation, maximizing overall nutrient consumption and global productivity.

Fig. 1.

Schematic of the “ligand–iron–microbe” feedback. Free iron () is added to the ocean from external sources like dust and sedimentary mobilization, and lost due to scavenging, precipitation, and burial. Iron retention is increased by complexation with organic ligands (; and ). Microbial production of biomass () is dependent on iron, and is a source of ligands, for example, by siderophore production, excretion of organic carbon, and release of cell detritus during remineralization. The production of ligands retains a greater concentration of iron, fueling more microbial production, and so on, until the other resources that microbial production requires, such as macronutrients, become limiting. Thus the “ligand–iron–microbe” feedback maintains “just enough” iron in the ocean to match the availability of other resources, resulting in colimitation.

In the following, we show that such a feedback does emerge in an idealized three-box ocean biogeochemistry model that resolves a dynamic pool of iron-binding ligands.

Model Framework

The model (Fig. 2) represents an ocean with an overturning circulation linking a low iron input, upwelling regime, analogous to the “high-nutrient, low-chlorophyll” regions of the ocean, to a high iron input, deep water formation regime. In this idealized system, we consider two basins, which we will refer to as the “Southern Ocean” and “Atlantic Ocean,” respectively, but the overall concepts could equally apply to other regions. Further details are given in .

Fig. 2.

Schematic of the idealized three-box ocean biogeochemistry model. See text and for further details.

Microbial net production (, moles N per cubic meter per second) in the upper ocean consumes macronutrients (, dissolved inorganic nitrogen, moles per cubic meter), and total dissolved iron (, moles per cubic meter) using Michaelis–Menten kinetics (Eq. ). The maximum rate, , is limited by , , or light (, watts per square meter) levels (with coefficients , , and ).This fixed organic material has Redfield stoichiometry (). Coupled macronutrient and iron cycles (Eqs. and ) are solved following ref. 49. Defining the substantial derivative including ocean transport () and mixing () of tracer, , as in moles A per cubic meter per second, the governing equations in each model box (see text below for definitions) are and are taken up by biological activity at the surface, exported to the deep box, and completely remineralized ( is a step function: +1 at the surface and −1 at depth). The global inventory is assumed constant in this idealized model.

Iron is supplied to the surface (, moles Fe per cubic meter per second), representing aeolian (50) or shelf sediment (51) sources, with 100 times less entering the “Southern Ocean” than “Atlantic Ocean” boxes (50, 52). A hydrothermal iron source is also added to the deep box (53).

For illustrative purposes, we consider a single, representative pool of total ligand (, moles per cubic meter; Eq. ) with a conditional stability coefficient that determines the partitioning between free and complexed iron (49). Free iron (, moles per cubic meter) is scavenged and precipitated when in excess of its low solubility with rate constant (per second).

Ligands are produced as a fraction (, moles per moles C) of organic matter turnover in both surface and deep ocean associated with net production in the surface layer, portraying exudate release (27) and siderophore production (54), and deep ocean remineralization of sinking particulate organic matter (46), such as release of porphyrins and degraded protein chelators. Ligand loss is a first-order process with rate (per second), associated with microbial degradation. The prescribed rate of loss is modified to be 100 times longer in the deep box than in the surface boxes (55) to account for the gradient introduced by enhanced depletion through photochemical degradation (24, 56) and microbial uptake (57) in the upper ocean, slower losses at depth due to an order of magnitude decrease in bacterial abundance (58), and slower metabolism resulting from cooler temperatures.

In reality, ligands are diverse and only partially characterized (e.g., refs. 32 and 34), so we have little direct constraint on their source and sink rates (i.e., and values). At steady state, when transport terms are constant (), Eq. implies a relationship between the representative ligand source/sink rates ( and ) and measurable bulk characteristics of the whole ocean,We can approximate a likely range of from an observational perspective by substituting in Eq. with an order of magnitude estimate of global oceanic net production [10 GtC y−1 to 100 GtC y−1 (59)], and observed ligand concentrations () between 0.6 mol and 6 mol (e.g., refs. 29–31, 52, and 60), in the open ocean. First, converting net production into concentration units by distributing the integrated estimate (59) over either mixed layer volume or full ocean volume (61) gives s. An equivalent range is found using the same mixed layer volume over which organic carbon production occurs and accounting for the gradient in , assuming deep ocean values are 100 times smaller than prescribed surface values.

The ligand to organic carbon ratio () has been estimated between mol (mol C)−1 and mol (mol C)−1 (46, 62–64). Combining these values with our approximate range of leads to ligand decay rates () of to , or ligand turnover times () between 20 d and 960 y. These bulk values fall within a range of degradation timescale estimates for partially labile and semirefractory dissolved organic carbon [1 y to 1,000 y (64)], and have an upper bound (e.g., ref. 31) consistent with estimates of the residence time of ligands in deep waters of the Atlantic [779 y to 1,039 y (30)].

Illustration of the Ligand–Iron–Microbe Feedback

The emergent “ligand–iron–microbe” feedback can be cleanly illustrated with a single model equilibration time series (Fig. 3). The model is initialized with globally uniform macronutrient concentrations (33 mmol N) and no dissolved iron or ligands, then integrated for 10,000 y. Ligand characteristics are set such that s. Initially, the model is iron-limited globally, with elevated macronutrients in both “Atlantic Ocean” and “Southern Ocean” surface boxes (Fig. 3). Relatively high iron delivery to the “Atlantic Ocean” box leads to an initial rise in productivity (Fig. 3) and depletion of surface macronutrients (Fig. 3). This drives ligand production (Fig. 3), allowing accumulation of a standing stock of deep ocean iron (Fig. 3). In the following centuries, macronutrients stay depleted with elevated productivity, and “North Atlantic” ligand levels converge toward steady state due to transport and loss processes.

Fig. 3.

Illustration of the “ligand–iron–microbe” feedback: 10,000-y time series (note the log-scale axis) of (A) macronutrient, (B) dissolved iron, and (C) ligand concentrations in the “Southern Ocean” box (yellow), “Atlantic Ocean” box (green), and Deep Ocean box (blue), and (D) globally integrated export production. Steady state is reached after 1,000 y. The model is initialized with uniform 33 mmol N of macronutrients and no iron and ligands, with 4,500 s ( mol (mol C)−1 and y from the observationally constrained range).

In contrast, lower iron input to the “Southern Ocean” box cannot support rapid macronutrient drawdown. On longer timescales, as ligand levels increase throughout the ocean, upwelled chelated iron drives a gradual incomplete reduction of surface macronutrients.

Steady state is reached after 1,000 y. The “Southern Ocean” box is iron-limited (Fig. 3), with incomplete macronutrient use (Fig. 3) fueled by iron supply from the deep ocean, while the “Atlantic Ocean” box is macronutrient-limited, with sufficient iron to fully consume macronutrients primarily delivered by the overturning circulation.

Therefore, using ligand source and sink rates guided by present-day observations (i.e., constraining ), we demonstrated an emergent positive feedback that promotes global-scale iron and macronutrient colimitation with uniform ligand distribution (2 mol) and slightly elevated ligand abundance in the upwelling box, reminiscent of the modern ocean (31, 60, 64).

How robust is this result? The next section presents a suite of simulations over a wide range of and initial conditions.

Robustness of the Feedback: A Large Model Ensemble

A suite of 10,000 simulations of the three-box ocean biogeochemistry model was integrated for 10,000 y each, with independently assigned values (Fig. 4) of between mol (mol C)−1 and mol (mol C)−1 and between 0.1 y and 1,000 y. This range of extends far beyond that inferred for the modern ocean. Integrations start with random initial concentrations of iron, ligand, and macronutrients (but the same total macronutrient inventory).

Fig. 4.

Ensemble of 10,000 box model simulations. (A) Range of prescribed and values, (B) envelope of resulting surface nitrate (blue), iron (red), and ligand (green) concentrations, and (C) envelope of export production rates for the “Southern Ocean” box (blue), “Atlantic Ocean” box (green), and global total (red). Envelopes are calculated from median values of simulations with equal ratios, that is, averaged within colored contours in Fig. 4 the median absolute deviation (note the log scales). The range of informed by oceanic observations is indicated by arrows on the x axis. Guided by these data, the experiments can be partitioned into three regimes: 1) iron-replete (macronutrient limited) simulations, 2) iron-limited (macronutrient replete) simulations, and 3) iron and macronutrient colimited simulations (region shaded gray)—see text for details.

Results are summarized in Fig. 4, which shows the envelopes of globally averaged, steady-state surface nutrient concentrations as a function of , and Fig. 4, which shows envelopes of export production rate in each box, as well as the global total. It is notable that the concentrations and rates of microbial production converge to a tight curve as a function of while representing numerous combinations of individual and values, as well as a range of random arbitrary initial conditions. In other words, the outcomes are robust and predictable for any given , and independent of initial conditions: The feedback restores the degree of biological limitation between macronutrients and iron for the specific value of , in the same way as in the time series shown in Fig. 3.

The 10,000 simulations fall into three groups that express fundamentally different nutrient regimes (Fig. 4 ):

The first is an iron-limited regime with relatively low that cannot accumulate a sufficient ligand pool due to low ligand production or high rates of ligand decay. The standing stock of iron is very low, and any iron input (particularly in the “Atlantic Ocean” box) is rapidly consumed or lost to scavenging, leaving surface macronutrients incompletely used. Resulting export production, largely in the “Atlantic Ocean” box, is low.

The second is an iron-replete regime, with relatively high and a large ligand pool due to rapid ligand production or slow decay rates. This supports maximum iron abundance (set by total external iron input accumulated over 10,000 y) relieving limitation even in the “Southern Ocean” box, which receives little external iron input. Surface macronutrients are low and export production is high, predominantly in the “Southern Ocean” box where nutrients are upwelled by the overturning circulation.

The third is an iron and macronutrient colimited regime in which the “ligand–iron–microbe” feedback restores the system to global-scale colimitation, with both macronutrient- and iron-limited regions. Ligand abundance supports sufficient iron levels to obtain near-maximum global export production rates mediated by macronutrient, iron, and ligand redistribution via ocean circulation, as in Fig. 3.

In short, simulations fall into three robust regimes, governed by the characteristic ratio of ligand source and sink rates. Global colimitation, with regional limiting resource variation, is subject to a reinforcing feedback if ligand characteristics are in the appropriate regime, but it is not inevitable.

Where does the modern ocean fit into this continuum? We scored ensemble members by quantitative comparison to oceanic observations (: possible range 0 to 1; Fig. 5 and ) taking the median of scores for simulations with the same (Fig. 5). Scores were compared to a benchmark () evaluated from a model simulation where ligand concentration was held fixed at a uniform, global value of 1 mol (as employed in global biogeochemistry models; e.g., ref. 65). Dynamic ligand simulations with model–data comparison scores greater than the benchmark value are outperforming the standard parameterization. The best scores are obtained in a band of intermediate with moderate levels of macronutrients, iron, and ligands. The optimum log10(γ/λ) value is 3.643 ( = 4,398 s) with a score of 0.81 (Fig. 5).

Fig. 5.

Ensemble of 10,000 box model simulations. (A) Model–data comparison scores (; Eq. ), with gray contours of macronutrient usage efficiency (note the uneven contour interval). (B) Envelope of model–data comparison scores and (C) envelope of macronutrient use efficiency, both calculated from the median of simulations with equal the median absolute deviation. The range of informed by oceanic observations is indicated by arrows on the x axis. Best-fit simulations have an optimal (4,398 s) marked by the vertical black dashed line. The optimum ratio lies within the data-constrained range of , and within the iron and macronutient colimited regime identified in Fig. 4. The horizontal dashed line in B is a benchmark model–data comparison score (0.69) calculated from an instance where ligand concentration is a fixed, uniform, 1 mol. The horizontal dashed line in C is at the emergent value of macronutrient usage efficiency (0.29) at the optimal ratio.

Macronutrient use efficiency was diagnosed as the fraction of macronutrients transported from the surface to the deep by biological activity (7–9). Efficiency of macronutrient use is positively correlated with (Fig. 5 ), and reflects the pattern of export production in the low-iron, upwelling “Southern Ocean” box (Fig. 4; e.g., refs. 66 and 67). Model solutions exceeding the benchmark score have efficiencies of 0.22 to 0.40, with an optimum of 0.29 at the model–data best fit (Fig. 5).

Discussion

Iron availability in the modern, oxygen-rich, ocean is intimately tied to the cycling of organic chelating ligands, and so could form a reinforcing cycle: Production of ligands enhances the availability of iron, in turn increasing biological activity, leading to further production of ligands. Logically, this feedback should be halted when another resource becomes limiting (Fig. 1). In other words, the feedback drives the system toward global macronutrient and micronutrient colimitation and maximum net productivity, mediated by macronutrient, iron, and ligand redistribution via ocean circulation (Fig. 3).

We explored this emergent “ligand–iron–microbe” feedback in an idealized ocean biogeochemistry model with a dynamic ligand pool. The ratio of aggregate ligand production and loss rates () is the critical parameter that controls the biogeochemical state of the global ocean. Our model suggests three major regimes: the “modern” regime with intermediate ligand lifetime and global-scale colimitation, a long-lived ligand regime with an iron-replete global ocean, and an iron-limited short-lived ligand regime where scarce chelators cannot maintain a standing stock of iron. Guiding with observations, the hypothesized feedback drives the system into global-scale colimitation between iron and macronutrients (Fig. 3) reminiscent of the observed ocean (60), with relatively uniform ligand levels and slightly higher abundance in the upwelling regime (31, 64). Macronutrient use efficiency (0.29, range 0.22 to 0.40) is also consistent with the modern ocean value of 0.36 (9). Globally, the feedback maintains a pool of ligands that can support “just enough” iron to match the abundance of other resources, such as macronutrients, leading to global colimitation and maximized rates of production (Fig. 4).

We note that this highly idealized model represents the aggregate properties of what is a rich spectrum of chelating molecules with differing origins and cycles. Here we consider broad shifts in the overall properties of this suite of molecules, which is a practical choice that ought to be revisited as understanding of ligands and their cycles is revealed in the future. There are other properties of ligands, such as their conditional stability coefficient(s), which also vary but are not explicitly explored here. More-comprehensive approaches to modeling ligand cycles are being explored (e.g., ref. 64) and depend on determining constrained parameter values, for which the global GEOTRACES dataset (60) and target process studies identifying ligand sources and sinks and iron retention mechanisms (e.g., refs. 37, 38, and 68) should provide fundamental constraints.

There are analogies between the feedback suggested here, linking ligand and iron availability to macronutrient concentration at certain scales, and discussions of interplay between the nitrogen and phosphorus cycles (4–6). Our discussion parallels the perspective that phosphorus may be the ultimate limiting nutrient in the ocean, with nitrogen held close to Redfieldian proportions on the global scale through the balance of nitrogen fixation and denitrification (4–6). Our focus on iron means we assumed a fixed global macronutrient inventory and have not resolved this complex and subtle interplay. In practice, some interesting connections are likely, for example, through iron limitation of nitrogen fixation (5), and interactions with the oxygen cycle and redox chemistry. These could be addressed within a more complex modeling framework (69–71).

Keeping in mind simplifications of the model, a key inference is that the modern ocean falls into the intermediate nutrient regime, where ligand characteristics lead to global-scale macronutrient and iron colimitation. What sets the properties of the ligand pool? Could they reflect an evolutionary self-organization of the global ocean ecosystem?

Proterozoic oceans were likely rich in dissolved ferrous iron (Fe(II)) due to anoxia (72–74), but production was limited by low concentrations of phosphate and nitrogen (75–79). It has been argued that nutrients were limiting due to several intertwined factors including slow remineralization in the absence of oxygen (80, 81), phosphate scavenging by abundant Fe(II) ions (78, 82), and suppression of nitrogen supply due to scavenging of molybdenum (83, 84), an essential nitrogen cycle cofactor, by sulfide ions in coastal environments (85, 86).

Around the Neoproterozoic (87, 88), ancient marine cyanobacteria expanded into the nutrient-poor oceans as evolutionary selection drove internal stoichiometry toward lower subsistence nutrient ratios and also increased exudation of dissolved organic carbon such as polysaccharides and carboxylic acids (89). The iron-chelating properties of organic carbon exudates may have initially promoted the release of phosphorus from dust by dissolution of iron minerals, which is a strategy some microbes employ today under phosphorus stress (90). As the ocean became oxygenated, these organic molecules may have facilitated the aerobic availability of Fe(III) by acting as iron-stabilizing ligands (43), helping photosynthesizers overcome a key negative feedback between oxygen-producing photoautotrophy and their dependency on iron (79). Enhancing supply of iron and organic carbon in nitrogen-poor surface oceans may have created new opportunities for nitrogen fixers (87), increasing oceanic nitrogen supply (75–77, 79).

Transitioning from an anoxic ocean rich in free Fe(II) to an oxygenated ocean with organic carbon-bound Fe(III) (79) would have initiated the “ligand–iron–microbe” feedback, with implications for nutrient consumption and ocean productivity.

Phanerozoic innovations continued to alter ocean ecosystem and marine iron chemistry interactions. Siderophore production and uptake genes are widely distributed in heterotrophic bacteria (91), while microbes that cannot produce competitively binding ligands, like some eukaryotic phytoplankton, developed methods to steal siderophore-bound iron by employing ligand nonspecific iron reductases (e.g., refs. 40, 44, 47, and 92). This suggests temporal variations in .

The “ligand–iron–microbe” feedback may merit reassessment of the climatic effects of changing iron sources on atmospheric . We would anticipate little change to the global outcome of iron fertilization experiments when resolving dynamic ligands, because the ocean should already be at or near global colimitation with near-maximum export production, although local changes may occur to allow small atmospheric anomalies (e.g., refs. 16–21). On longer timescales, the “ligand–iron–microbe” feedback could potentially act by self-amplifying or self-damping climatic changes, depending on the associated nutrient dynamics in the oceans. During ice ages (93), there may have been an increase in iron supply as well as greater macronutrient availability through possible erosion of exposed shelf sediments or greater marine nitrogen fixation (e.g., refs. 5, 8, 94, and 95). We predict that, under these conditions, the “ligand–iron–microbe” feedback would act in a self-amplifying manner, by increasing production of ligands and iron abundance so as to drive the system toward colimitation with macronutrients, maximizing productivity. In a warmer, more humid, and therefore less dusty climate (e.g., ref. 20), we, in turn, predict the “ligand–iron–microbe” feedback would first act to buffer the ocean from the effects of decreased iron supply by maintaining a high rate of production and thus ligand supply that extends the lifetime of iron. Once the global ocean transitions from the colimited regime into the iron-limited regime, the feedback might then exacerbate iron scarcity through reduction in productivity and declining ligand production.

Finally, we note that our proposed feedback could apply to other trace metals that are subject to chelation by organic ligands in the ocean. There is close correspondence between the stoichiometry of organic matter and available dissolved concentrations of trace metals (figure 1 in ref. 12), suggesting that the system may have coevolved to optimize micronutrient supply and demand for several other trace elements (26, 96, 97). Perhaps there is also “just enough” of these other micronutrients in the ocean to match levels of the more abundant macronutrients, promote global resource colimitation, and maximize production.

Materials and Methods

An idealized model (Fig. 2) is used to efficiently explore the “ligand–iron–microbe” feedback. There are two asymmetrical surface boxes as well as a global deep ocean box linked by an overturning circulation (). Deep water formation occurs in the “Atlantic Ocean” with upwelling in the “Southern Ocean”. These boxes are named for ease of reference, but the overall concepts could equally apply to other regions. A flux due to mixing () also occurs. See Table 1 for parameter values.

Table 1.

Fixed box model parameter values

Name	Symbol	Units	Value
No. of boxes	n_b		3
Latitudinal extent	y	m	3.0×10^6, 13.0×10^6,
			16.0×10^6
Box thickness	z	m	50.0, 50.0, 5050.0
Box volume	v	m^3	3.4×10^15, 1.0×10^16,
			1.4×10^18
Overturning rate	ψ	m^3⋅s^−1	20.0×10^6
Mixing rate	κ	m^3⋅s^−1	1.0×10^6
Maximum productivity	α	mol m^−3 y^−1	5.1×10^−4
N half-saturation	k_N	mmol m^−3	0.1
Fe half-saturation	k_Fe	μmol m^−3	0.1
Light half-saturation	k_I	W m^−2	45.0
Fe source rate	s_Fe	mol m^−3⋅s^−1	2.1 ×10^−15,
			2.1×10^−13,
			2.3×10^−17
Fe′ scavenging rate	k_s	s^−1	1.0×10^−7
L_T stability coefficient	β	m^3mol^−1	1.0×10^6
N reference mean	N_0¯	mmol m^−3	24.0, 3.0, 31.6
Fe reference mean	Fe_0¯	μmol m^−3	0.2, 0.3, 0.7
L_T reference mean	L_T0¯	μmol m^−3	1.9, 2.0, 2.4
N reference SD	σ_N_0	mmol m^−3	3.7, 5.1, 11.5
Fe reference SD	σ_Fe_0	μmol m^−3	0.4, 0.3, 2.3
L_T reference SD	σ_L_T0	μmol m^−3	1.0, 0.9, 3.4

Triplets indicate values in the “Southern Ocean,” “Atlantic Ocean,” and “Deep Ocean” boxes, respectively.

Net production () is modeled using Michaelis–Menten kinetics (Eq. ), with maximum rate, , and potentially limited by macronutrients (Eq. ), iron (Eq. ), or light, with half-saturation coefficients, , , and .

Iron is asymetrically supplied to the surface (), representing aeolian dust (50) or sediment (51) sources. A hydrothermal iron source [1 Gmol Fey−1 (53)] is also added to the deep ocean.

For illustrative purposes, we consider a single dynamic, representative pool of ligand (Eq. ) with a conditional stability coefficient () that determines partitioning between free ( and ) and complexed () iron and ligands (49): . Free iron is scavenged and precipitated when in excess of its low solubility with rate constant . Ligand cycling is described by Eq. .

Our model suite comprises 10,000 simulations of 10,000 y, each with different rates of ligand production () between 1 and 1 and ligand loss () between 0.1 y and 1,000 y. The rate of ligand loss is modified to be 100 times smaller in the deep box than in the surface box (55) to account for the gradient resulting from enhanced depletion by photochemical degradation (24, 56) and microbial ligand uptake (57) in the upper ocean, slower losses due to an order of magnitude decrease in bacterial abundance (58), and slower metabolism resulting from cooler temperatures, at depth. The global macronutrient inventory is conserved, with initial concentrations randomly assigned a fraction of the total inventory. Iron and ligand levels were randomly assigned between 0 mol and 100 mol.

Simulations are scored by their resemblance to oceanic observations (Eq. ) for equally weighted, steady-state macronutrient, iron, and ligand concentrations. Reference values (, , and ) are averages and standard deviations from representative geographic areas for each box using World Ocean Atlas 2013 (98) for nitrate and using the GEOTRACES intermediate data product 2017 (60) for iron and ligands. Owing to paucity of ligand data in the Southern Ocean, and comprise data from the entire Southern Hemisphere.where , , and .

Our overall characterization of the “ligand–iron–microbe” feedback is insensitive to geometry, circulation, or biogeochemistry parameters, although the exact details vary for choices of and .

Model code, output, and processing routines can be accessed via GitHub (http://bit.ly/lauderdale-ligand-iron-microbe-feedback) (99).