nature geoscience

Article https://doi.org/10.1038/s41561-022-01066-2

Global patterns in marine organic matter

stoichiometry driven by phytoplankton
ecophysiology

Received: 15 June 2021 Keisuke Inomura 1,2,3 , Curtis Deutsch 2,4, Oliver Jahn 3
,
Stephanie Dutkiewicz3 & Michael J. Follows3
Accepted: 29 September 2022

Published online: xx xx xxxx

The proportion of major elements in marine organic matter links cellular
Check for updates processes to global nutrient, oxygen and carbon cycles. Differences in the
C:N:P ratios of organic matter have been observed between ocean biomes,
but these patterns have yet to be quantified from the underlying small-scale
physiological and ecological processes. Here we use an ecosystem model
that includes adaptive resource allocation within and between ecologically
distinct plankton size classes to attribute the causes of global patterns in
the C:N:P ratios. We find that patterns of N:C variation are largely driven by
common physiological adjustment strategies across all phytoplankton,
while patterns of N:P are driven by ecological selection for taxonomic
groups with different phosphorus storage capacities. Although N:C varies
widely due to cellular adjustment to light and nutrients, its latitudinal
gradient is modest because of depth-dependent trade-offs between nutrient
and light availability. Strong latitudinal variation in N:P reflects an ecological
balance favouring small plankton with lower P storage capacity in the
subtropics, and larger eukaryotes with a higher cellular P storage capacity in
nutrient-rich high latitudes. A weaker N:P difference between southern and
northern hemispheres, and between the Atlantic and Pacific oceans, reflects
differences in phosphate available for cellular storage. Despite simulating
only two phytoplankton size classes, the emergent global variability of
elemental ratios resembles that of all measured species, suggesting that the
range of growth conditions and ecological selection sustain the observed
diversity of stoichiometry among phytoplankton.

Nearly a century ago, Alfred C. Redfield described the relationship of nutrient demand by phytoplankton6 affect microbial resource com-
between the average elemental ratios of phytoplankton and biogeo- petition, regulating the balance between denitrification and nitrogen
chemical cycles1. Since then, the ‘Redfield ratios’ (C:N:P = 106:16:1) have fixation and the global N inventory7–9. While the plasticity of elemental
become a cornerstone of oceanography. They affect carbon export and ratios in phytoplankton populations has long been recognized, fixed
storage in the deep ocean and thus atmospheric CO2 concentration2–4, Redfield ratios continue to be widely used in ecological and biogeo-
and the intensity of the oxygen minimum zones5. The elemental ratios chemical models10–12.

Graduate School of Oceanography, University of Rhode Island, Narragansett, RI, USA. 2School of Oceanography, University of Washington, Seattle,
1

WA, USA. 3Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA. 4Department of
Geosciences and High Meadows Environmental Institute, Princeton University, Princeton, NJ, USA. e-mail: inomura@uri.edu

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Regional variation of elemental ratios in organic matter produc- termed MITgcm-CFM; see Methods, Extended Data Figs. 1 and 2, and
tion and export has been observed in the ratios of surface nutrient Extended Data Table 1). The cellular model relates macromolecular and
drawdown13 and detected in large-scale geochemical tracer distribu- elemental allocation to light intensity and growth rate, and has been
tions5,14,15. The latitudinal gradients inferred from geochemical data calibrated and validated with laboratory data from several species of
have been confirmed by global compilations of direct measurements phytoplankton21. The ocean model simulates nutrient and carbon cycles
of bulk organic matter, supporting a strong latitudinal variation of N:P driven by multiple processes, including the growth and mortality of
and P:C in organic matter, and weaker trends in N:C16. Together, these phytoplankton, the formation and decay of particulate and dissolved
observations reveal substantial correlations between elemental stoichi- organic matter10,38, and transport by the ocean’s general circulation35–37.
ometry and phytoplankton community structure, with high N:P and P:C The model’s biological component, CFM-Phyto, predicts the ele-
in subtropical regions dominated by small phytoplankton14, and lower mental stoichiometry of phytoplankton based on resource allocation
ratios elsewhere. While the stoichiometries of bulk organic matter and under different light intensities and nutrient concentrations. Essential
residual surface nutrients probably originate from phytoplankton, the elements are apportioned to major groups of macromolecules (for
potential underlying physiological and ecological mechanisms have example, proteins, carbohydrates, lipids, DNAs and RNAs) with distinct
not been elucidated. stoichiometric ratios. The model predicts relationships between the
Laboratory studies have revealed substantial variations in elemen- abundance of these molecules, growth rate and light intensities (see
tal proportions of phytoplankton within and across taxa17–20, both of Methods). It reproduces laboratory-measured relationships between
which could underlie large-scale patterns of stoichiometric variation. these factors and cellular stoichiometry that are shared across multiple
Intra-taxonomic variations reflect the balance of major macromolecu- species of both eukaryotes and prokaryotes23,39,40.
lar pools, having distinct elemental ratios19,21,22. While carbon occurs The implementation of CFM-Phyto in the MITgcm is adapted
in most macromolecules, nitrogen is predominantly associated with from its original form in two ways. First, inter-specific variation is
proteins, and phosphorus is largely contained in RNA and storage represented by two size classes of phytoplankton: ‘small’ to represent
molecules19,21. The allocation to protein and RNA increases with growth prokaryotes and ‘large’ to represent eukaryotes. Informed by empirical
rate20–23, enabling faster biosynthesis and leading to higher N and P data, large (eukaryotic) phytoplankton are assumed to have a higher
relative to C. Thus, phytoplankton in a fast-growing environment (for P storage capacity26 and higher maximum growth rates41,42 than small
example, high nutrient) are expected to have high N:C and P:C24,25. (prokaryotic) phytoplankton (see ‘Differentiating small and large phy-
Intracellular storage can also have a large impact on elemental ratios, toplankton’ in Methods). Second, the elemental allocation is expanded
especially for P, for which the structural pools are much smaller than to include two intracellular pools of Fe: Fe in photosystems and Fe in
those of N and C21. Differences in resource allocation at the cellular scale storage (see ‘Relating Fe quota and growth rate’ in Methods).
represent a physiological mechanism by which global stoichiometric
patterns can arise from depth and/or spatial variations in the growth Cellular-scale variations
conditions of phytoplankton. Global patterns of nutrient uptake yield distributions of surface macro-
There is also a distinct pattern in inter-taxonomic variation of nutrient concentrations and of the growth-limiting nutrient that are
the elemental ratios. Measurements across multiple taxa show that consistent with observations (Extended Data Figs. 3 and 4)43. The model
N:P ratio of eukaryotes is on average lower than that of prokaryotes26. predicts substantial variation in element ratios between small and
The difference can be partly explained by the capacity to hold excess large plankton, which are consistent with the observed differences
phosphorus24. Ocean regions with a higher fraction of larger eukary- between eukaryotes and prokaryotes (Fig. 1). These differences are
otic phytoplankton may lead to lower N:P than regions dominated most pronounced in ratios involving P27. Differences in P storage explain
by smaller prokaryotes27, provided enough P is available. Large-scale most of the P:C and N:P variability in CFM-Phyto (turquoise arrows in
differences in the size structure of phytoplankton communities intro- Fig. 1a,b). To a smaller degree, the differences also reflect the struc-
duce additional ecological mechanisms that may generate global tural composition (molecules other than P storage) of the cell, which
stoichiometric patterns26, modulating the physiological factors that is dominated by proteins. However, N storage has a relatively minor
arise from cellular-scale allocation. effect on the size-based differences, as implied by the similarity of N:C
The biological causes of observed large-scale distributions of distribution between prokaryotes and eukaryotes (Fig. 1a,b). Varia-
organic matter C:N:P remain poorly understood. In particular, the tions in N:C are instead driven by allocation to protein in response to
relative contribution of plasticity within phytoplankton groups ver- environmental factors.
sus ecological selection among groups with systematic differences To quantify the role of P storage in generating large differences in
in stoichiometric ratios is not known. Theoretical studies of variable N:P among plankton size classes, we estimated the level of the excess
stoichiometry typically employ the internal-stores model27,28, which phosphorus uptake for both prokaryotic20,39 and eukaryotic44 phyto-
is empirically informed but does not resolve the macromolecular plankton in laboratory experiments. Under N-limited conditions, cells
allocation or its acclimation to changing environmental conditions. accumulate storage P due to excess P availability and increasing P:C.
Models that relate macromolecular allocation to the physiological Under P limitation, cells maintain a minimum, necessary P content in,
status of the organism provide more mechanistic detail (for example, for example, nucleic acids with a lower P:C. This difference manifests
refs. 25,29–32, with reviews by refs. 33,34 and the supplement of ref. 21) but as stored P (per C) associated with luxury uptake24. Laboratory stud-
their implications for global stoichiometric patterns have yet to be fully ies reveal substantially higher P storage capacity in eukaryotic cells
analysed. Here we address this gap by implementing a cellular resource (Fig. 1c), contributing to the lower overall N:P observed among large
allocation model within a global model of ocean circulation and bio- cells (compare Fig. 1a and b).
geochemistry. Model simulations reproduce observed variability at the
relevant scales, from cells to biomes, allowing an empirically validated Spatial patterns of elemental ratios
diagnosis of key physiological and ecological factors. The observed latitudinal variation in the N:C of organic matter is distinct
from that of N:P16 (Fig. 2). The N:C ratio has relatively small variation (the
Simulating cellular macromolecules in a global coefficient of variation (CV) = 0.11) but increases slightly towards higher
ocean model latitudes. In contrast, N:P varies strongly (CV = 0.33) and systematically
We incorporated an explicit representation of macromolecular alloca- with latitude across all ocean basins; low in the high latitudes, high in
tion by phytoplankton (CFM-Phyto21) into a global general circulation the subtropical gyres and intermediate values near the Equator. The
and biogeochemical model, MITgcm35–37 (here the combined model is values are slightly higher in the Northern Hemisphere, especially in the

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

a Small phytoplankton b Large phytoplankton c

0.030 0.030 0.040
Contour: Contour:
Model: 5 Model: 5 N:P (mol mol–1)
N:P (mol mol–1) 0.035
0.025 Structure 0.025 Structure
Storage (P) Storage (P) 0.030
10 10

PStore:C (mol mol–1)

0.020 0.020
P:C (mol mol–1)

P:C (mol mol–1)

Data: Data: 0.025
Prokaryote Eukaryote
0.015 15 0.015 15 0.020
Redfield Redfield
20 20 0.015
0.010 25 0.010 25
30 35 30 35 0.010
0.005 40 0.005 40
45 45 0.005

0 0 0
0 0.05 0.10 0.15 0.20 0.25 0 0.05 0.10 0.15 0.20 0.25 Prokaryote Eukaryote
–1
N:C (mol mol ) –1 N:C (mol mol )

Fig. 1 | Observed and modelled elemental ratios in phytoplankton. a,b, vectors (modelled P storage). Larger, eukaryotic cells in b are associated with
Variations in elemental ratios in ‘small’ (prokaryotic) phytoplankton (a) and higher storage contributions (longer blue vector) than smaller prokaryotic cells
‘large’ (eukaryotic) phytoplankton (b). Colour shading indicates N:P, computed in a. c, Differences in P storage between small and large plankton size classes in
as the ratio of N:C (x axis) and P:C (y axis). Laboratory data for small prokaryotic the model are based on empirical estimates derived from laboratory studies of
cells (white points, a) and eukaryotic cells (white points, b) at a variety of growth prokaryotic20,39 and eukaryotic44 phytoplankton (n = 43 and 28 for prokaryotes
rates and light intensities (excluding a few outliers) (see Supplementary Data and eukaryotes, respectively. A few data with excess growth-limiting nutrients at
and references there). Arrows indicate the stoichiometric ratios predicted by the steady state39 are not included). The box represents median (centre line) and
the allocation model decomposed into structural and storage components first and third quartiles (box). The whiskers represent the value range without
based on average nutrient and light conditions from the surface ocean at 50° S, outliers (those outside the box by 1.5 times the interquartile range). The P storage
where N and P nutrients are largely replete. Lilac arrows indicate the modelled is estimated based on the differences in P:C under N and P limitations for closest
contribution from acclimation in the absence of P storage. Observed points fall growth rates.
above those lines due to P storage, the sense of which is indicated by the light blue

a b
0.30 40

Pacific
0.25 Atlantic
Indian 30
0.20
N:C (mol mol–1)

N:P (mol mol–1)

0.15 20

0.10
10
0.05

0 0
50° S 0° 50° N 50° S 0° 50° N
Latitude Latitude

Fig. 2 | Latitudinal variation of N:C and N:P ratios of bulk organic matter in with the highest uncertainty is not included. The CV based on the data and model
each basin. a, N:C ratio. b, N:P ratio. Cyan points are averaged data52,53 (data are 0.11 and 0.11 for N:C and 0.33 and 0.29 for N:P, respectively, implying stronger
distribution in Extended Data Fig. 5). Curves are model predictions: blue, Pacific latitudinal variations for N:P.
Ocean; orange, Atlantic Ocean; green, Indian Ocean. In each panel, a data point

Atlantic Ocean. These broad-scale, meridional patterns are reproduced driven by acclimation (Extended Data Fig. 6). The N:C ratio varies with
by the simulations (Fig. 2). These model–data comparisons are based on the primary environmental factors that govern growth rate: nutrient
particulate organic matter (POM), including phytoplankton biomass, availability and light. Low light requires a high investment in photosyn-
which accounts for a substantial fraction of the modelled standing stock thetic proteins to support growth. In the surface polar oceans, relief
of POM. Thus, the model fidelity to data suggests that the elemental of nitrogen stress enables higher growth rates and allocation to bio-
ratios of total organic matter are largely controlled by phytoplankton synthetic and photosystem proteins21, leading to high N:C of primary
with only limited alteration from organic matter recycling through the producers. In surface subtropical waters, nitrogen availability limits
microbial food web. growth rates, reducing investment in biosynthetic protein, while high
The N:C of modelled organic matter pools is strongly influenced by light intensity reduces investment in light-harvesting proteins, leading
the physiological acclimation of phytoplankton, which are the ultimate to low protein allocation and low N:C of primary producers. However,
source and substantial component of the standing stock of the total deeper in the subtropical water column, lower light and higher nutrient
particulate detritus. Thus, we seek to interpret the patterns of POM concentrations favour greater investment in protein and modelled N:C
by considering the controls on phytoplankton stoichiometry. Differ- increases with depth.
ences in N:C between the size classes are small (Fig. 1 and Extended Despite the strong variation of modelled N:C with latitude and
Data Fig. 6), whereas variations within each size class are large and depth (Fig. 3), the depth-averaged trends in N:C across latitude are

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

a N:C (mol mol–1) b NO3– (mmol m–3)

0 0
30
0.20
25
0.18
50 50
20
0.16
Depth (m)

Depth (m)
100 0.14 15
100

0.12 10

150 0.10 150

5

0.08

50° S 0° 50° N 50° S 0° 50° N

Latitude Latitude

Fig. 3 | Depth–latitude variation of predicted N:C and NO3−. a,b, Longitudinally which dominate the oligotrophic regime. The values are zonally averaged.
averaged N:C of phytoplankton (a) and NO3− concentrations (b). The black curve Though few in number, such observed growth rate profiles do reveal subsurface
in a indicates the depth of phytoplankton biomass maximum. The white curve maxima45,54–56.
in b indicates the depth of the highest growth rate of the small phytoplankton,

relatively small (Fig. 2) because of vertical trends in N:C and biomass. north–south asymmetry. The variation that cannot be explained by
A substantial fraction of the modelled, depth-integrated biomass is plankton size classes is here quantified by δN:P = N:P − N:PSize (Fig. 4b,d).
associated with a subsurface maximum at low latitudes (Fig. 3), as has This component of N:P variability is inversely related to the distribution
been observed in some subtropical profiles45–48. The subsurface chlo- of PO43− (Extended Data Fig. 3), because high PO43− concentrations lead
rophyll and biomass maxima arise in part from the trade-off between to enhanced P storage (equations (32)–(34) and (25) in the Methods).
opposing vertical gradients of nutrients and light49. At the depth of the Surface PO43− is more depleted in the northern than the southern sub-
emergent phytoplankton biomass maximum, phytoplankton N:C thus tropical gyres (Extended Data Fig. 3b) and limits cellular storage of
exhibits little variation with latitude (Fig. 3), except at some subpolar phosphorus regardless of size class (Fig. 4e,f). This results in higher N:P
latitudes where Fe is limiting (Extended Data Fig. 4), and the modelled in the Northern Hemisphere subtropics (particularly North Atlantic)
latitudinal variation in the depth-averaged N:C (CV = 0.11) is weaker than in the southern (Figs. 2b and 4a).
than at the surface (CV = 0.17). Cellular P storage also explains the asymmetry between polar
The latitudinal variation in N:P predicts higher values in the oceans of the northern and southern hemispheres. Even though large
low latitudes and lower N:P in high latitudes with stronger varia- phytoplankton are dominant in high latitudes in both hemispheres
tion (CV = 0.29) than for N:C (CV = 0.11), consistent with particulate (Fig. 4c), surface PO43− concentrations are higher in the Southern Ocean
observations (Fig. 2) and inferences from nutrient distributions14. The than in the North Atlantic (Extended Data Fig. 3b). This leads to a higher
model also predicts substantial N:P differences between ocean basins accumulation of plankton P storage (Fig. 4e,f) and lower particulate
(Fig. 4a). In contrast to N:C, there are substantial differences in N:P N:P in the southern high latitudes (Figs. 2b and 4a). The effect of P
between size classes, caused primarily by the distinct P storing capac- storage also explains the observed correlation between organic P:C
ity between small and large phytoplankton (Fig. 1). In the model, phy- and PO43− concentration4. In contrast, model simulations do not pre-
toplankton traits are guided by allometric constraints so selective dict a hemispheric N:C asymmetry, nor is it evident in observations16.
pressures in the oligotrophic regimes favour the smaller size class with Thus, we hypothesize that hemispheric asymmetry in N:P is created
lower P storage capacity (Fig. 1). At the same time, even though P is in mostly by the level of P storage per cellular C rather than variations of
excess over much of the subtropical ocean, the low concentrations N per cellular C. Further investigation is needed to clarify the form of
cause accumulation in P storage to be lower than the full potential. P storage; a large part of it may be polyphosphate, which can account
Thus, the patterns of total biomass N:P are caused by both physiological for a substantial fraction of cellular P in diatoms19, although RNA or
acclimation of each size class to its local light and nutrient levels, and by P-containing lipids may also contribute.
the ecological selection of the dominant size class in each environment.
The model also predicts a high fraction of total P in cellular stor- Global diversity of elemental ratios
age (Fig. 1a,b), which thus plays an important role in setting the overall On a global scale, modelled phytoplankton span a wide range of stoi-
stoichiometry. In turn, storage capacity is linked to cell size. How much chiometric ratios (Fig. 5), despite the explicit representation of only
of the N:P variation can be explained by the distribution of plankton size two plankton size classes. While the preponderance of emergent
classes? The answer to this question can be estimated from the local plankton biomass occurs with a stoichiometry near the canonical
fraction of total biomass and respective stoichiometries associated Redfield proportions (N:P = 16), each of the elemental ratios exhib-
with each size class (Fig. 4c), according to: its an approximately fourfold range across modelled populations
(Fig. 5). This variability across populations closely resembles the
N ∶ PSize = fS MS + (1 − fS ) ML (1) observed patterns of stoichiometric ratios sampled across plankton
species in laboratory data (Fig. 1). Moreover, the stoichiometric dif-
where MS and ML are global mean N:P ratios within each size class ferences between small and large plankton classes are similar to the
(Fig. 4b) and fS is the fraction of phytoplankton biomass in the small differences in the median observed species traits between small pho-
size class. Variations in N:PSize reflect the global scale pattern of N:P tosynthetic prokaryotes and larger eukaryotes. The peak modelled
(Fig. 4a,c) and account for about half of the total difference between biomass of small plankton occurs at a P:C ratio of ~0.05, less than half
subtropical and subpolar regimes. that of large plankton. A similar difference is also observed between
Intra-taxonomic variation of N:P due to acclimation also con- the median values of prokaryotic versus eukaryotic plankton species
tributes substantially to the subtropical enhancement and drives a (Fig. 5). In contrast, the ranges of highest biomass for N:C are

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

a b c
N:P (mol mol–1) 50 fS and N:PSize (mol mol–1) fS N:PSize
50 1.0
22
60° N 40 60° N
40 2
0.8 20
40° N 10 40° N

N:P (mol mol–1)

20° N 30 20° N 18
0.6
Latitude

Latitude
30 MS
0° 0°
16
20
20° S δN:P 10
1
20° S 0.4
20 14
40° S 40° S
10 0.2
12
60° S ML N:PSize 60° S
10
0
0 10 0 10
0° 120° E 120° W 0° 0.2 0.4 0.6 0.8 0° 120° E 120° W 0°
Longitude fS Longitude

d δN:P (mol mol–1) e f PStore:N (mol mol–1)

30 0.12
0.12
60° N 25 0.10 Large 2 60° N
PStore:N (mol mol–1) 10
40° N 20 40° N 0.10
0.08
20° N 15 20° N 0.08
Latitude

Latitude
0° 10 0.06 0°
1 0.06
10
20° S 5 20° S
0.04 0.04
40° S 0 40° S
0.02 Small 0.02
60° S –5 60° S
0
–10 10 0
0° 120° E 120° W 0° 0.5 1.0 1.5 2.0 2.5 0° 120° E 120° W 0°
3– –3
Longitude PO4 (mmol m ) Longitude

Fig. 4 | Global distribution of modelled N:P in phytoplankton and its that is linearly related to fS (see equation (1)). d, The residual N:P variation (δN:P)
ecological and physiological origins. a, N:P of total phytoplankton biomass. b, represents the physiological acclimation within size classes, including variations
The relationship between total N:P and the fraction of phytoplankton biomass in in P storage. e, Relationship between PStore:N for different PO43− concentrations. f,
the small size class (fS). c, The global mean N:P for small and large phytoplankton Spatial pattern of P storage per cellular N. All mapped values are ratios of biomass
(MS and ML, labelled on y axis in b) are used to estimate the N:P variation, N:PSize, N and P integrated between 0 and 260 m depth.
due to ecological selection for plankton size classes, which has a global pattern

a b
Large plankton Small plankton
0.025 0.025
14
14.5

0.020 13 0.020 14.0

13.5
12
0.015 0.015
13.0
P:C
P:C

11
12.5
0.010 0.010
10 12.0

0.005 9 0.005 11.5

11.0
8
0.05 0.01 0.15 0.20 0.25 0.05 0.01 0.15 0.20 0.25
N:C N:C

Fig. 5 | Stoichiometric variability among phytoplankton in modelled against data for eukaryotic and prokaryotic phytoplankton, respectively. The
populations and species observations. a,b, Model output of total model structure of the stoichiometric variability is an emergent property of the model,
phytoplankton biomass (shaded, log10 scale in mmol C) is binned according to but broadly consistent with the observed range of traits and the differences
its local N:C and P:C ratios (mol mol−1) alongside the laboratory observations between large and small size classes.
of elemental ratios (points). Large (a) and small (b) phytoplankton compared

similar between small and large as in the data (Fig. 5) because nevertheless predicts a wide range of stoichiometric ratios within and
N storage is small relative to the structural N pool in the cell 21. between distinct biomes. This stoichiometric diversity and its associ-
Similar to P:C ratios, and again consistent with observations, ated spatial patterns emerge from the physiological acclimation to
the N:P ratios of model phytoplankton populations have diver- different environmental conditions, and the ecological selection for
gent median values between large and small plankton (Extended populations that are well adapted to those conditions. The similarity
Data Fig. 7). between the observed and emergent elemental ratios suggests that the
Although the model only coarsely resolves phytoplankton size combination of physiological and ecological selection represented in
classes, and does not explicitly represent variation at a species level, it this simple model may be important selective pressures that generate

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

and sustain the diversity of stoichiometric ratios observed among References

global phytoplankton species. 1. Redfield, A. C. In James Johnstone Memorial Volume (ed. Daniel, R.
There are finer-scale taxonomic variations of elemental stoichiom- J.) 177–192 (Univ. Press of Liverpool, 1934).
etry18 than represented here, for example, the high N:P of diazotrophs, 2. Matsumoto, K., Rickaby, R. & Tanioka, T. Carbon export buffering
particulate inorganic carbon associated with coccolithophores and and CO2 drawdown by flexible phytoplankton C:N:P under glacial
the intimate connection of diatoms (which contribute to the large, conditions. Paleoceanogr. Paleoclimatol. 35, e2019PA003823
high-storage class) with the silica cycle. Given the broad qualitative (2020).
success of the simple model presented here, we hypothesize that 3. Broecker, W. S. Ocean chemistry during glacial time. Geochim.
additional taxonomic resolution and the impact of top-down con- Cosmochim. Acta 47, 1689–1705 (1982).
trols have next-order impacts on C:N:P. This remains to be tested in 4. Galbraith, E. D. & Martiny, A. C. A simple nutrient-dependence
future studies. mechanism for predicting the stoichiometry of marine
ecosystems. Proc. Natl Acad. Sci. USA 112, 201423917 (2015).
Implications for plankton stoichiometric 5. DeVries, T. & Deutsch, C. Large-scale variations in the
diversity stoichiometry of marine organic matter respiration. Nat. Geosci. 7,
Our results provide a bridge between the recognized diversity 890–894 (2014).
and plasticity of plankton stoichiometry at cellular scales, and the 6. Karl, D. M. et al. Ecological nitrogen-to-phosphorus stoichiometry
coherent large-scale patterns of stoichiometry documented in at station ALOHA. Deep Sea Res. II 48, 1529–1566 (2001).
nutrient distributions and bulk marine organic matter. The consist- 7. Lenton, T. M. & Watson, A. J. Redfield revisited 1. Regulation of
ency between phytoplankton C:N:P predicted by a physiologically nitrate, phosphate, and oxygen in the ocean. Regulation 14,
based ecosystem model and particulate observations supports the 225–248 (2000).
hypothesis that the resource allocation by phytoplankton is the pri- 8. Tyrrell, T. The relative influences of nitrogen and phosphorus on
mary influence on the composition of the material that is ultimately oceanic primary production. Nature 400, 525–531 (1999).
removed from the photic zone, and exerts long-term controls on 9. Weber, T. & Deutsch, C. Oceanic nitrogen reservoir regulated by
the coupling of biogeochemical cycles. Evaluating the potential for plankton diversity and ocean circulation. Nature 489, 419–422
heterotrophic processes to decouple the stoichiometry of exported (2012).
organic matter from that of phytoplankton50,51 will ultimately 10. Dutkiewicz, S. et al. Capturing optically important constituents
require simultaneous stoichiometric data on the various living and and properties in a marine biogeochemical and ecosystem
detrital pools. model. Biogeosciences 12, 4447–4481 (2015).
Our results suggest that phytoplankton N:C and N:P are controlled 11. Somes, C. J., Schmittner, A., Muglia, J. & Oschlies, A. A
by distinct physiological and ecological factors, indicating that elemen- three-dimensional model of the marine nitrogen cycle during the
tal ratios cannot simply be modelled with the conventional approach last glacial maximum constrained by sedimentary isotopes. Front.
that considers all the elements as uniform pools in a similar way across Mar. Sci. 4, 108 (2017).
species33. In particular, the modest latitudinal gradient of N:C, despite 12. Stock, C. A. et al. Ocean biogeochemistry in GFDL’s Earth System
its high variability at cellular and population scales, reflects physiologi- Model 4.1 and its response to increasing atmospheric CO2. J. Adv.
cal acclimation to local environmental conditions through resource Model. Earth Syst. 12, e2019MS002043 (2020).
allocation. In contrast, the large latitude gradients in P:C and N:P high- 13. Arrigo, K. R. et al. Phytoplankton community structure and the
light the key role of phosphorus storage and PO43− availability27, and the drawdown of nutrients and CO2 in the Southern Ocean. Science
ecological selection for or against plankton communities whose size 283, 365–367 (1999).
structure can accommodate large P storage. 14. Weber, T. S. & Deutsch, C. Ocean nutrient ratios governed by
While our model considers only two explicit size classes of plank- plankton biogeography. Nature 467, 550–554 (2010).
ton, it nevertheless captures a wide range of stoichiometric ratios 15. Teng, Y. C., Primeau, F. W., Moore, J. K., Lomas, M. W. & Martiny,
that approximates the range observed among species in laboratory A. C. Global-scale variations of the ratios of carbon to phosphorus
measurements. This stoichiometric diversity arises naturally from in exported marine organic matter. Nat. Geosci. 7, 895–898
the range of growth conditions that occurs in the ocean model. The (2014).
similarity of the simulated and observed ranges of stoichiometries and 16. Martiny, A. C. et al. Strong latitudinal patterns in the elemental
the measured differences among distinct taxa suggests that similar ratios of marine plankton and organic matter. Nat. Geosci. 6,
selective processes may be responsible for generating and sustaining 279–283 (2013).
stoichiometric diversity exhibited by phytoplankton species in the 17. Geider, R. J. & La Roche, J. Redfield revisited: variability of C:N:P
modern ocean. in marine microalgae and its biochemical basis. Eur. J. Phycol. 37,
The model presented here can be used to predict and interpret 1–17 (2002).
macromolecular distributions in phytoplankton biomass in the ocean 18. Finkel, Z. V. et al. Phylogenetic diversity in the macromolecular
and provides a physiological framework for predicting biological composition of microalgae. PLoS ONE 11, e0155977 (2016).
and ecological responses to climate change. Given that the elemental 19. Liefer, J. D. et al. The macromolecular basis of phytoplankton
ratios of phytoplankton influence the export of nutrients and carbon C:N:P under nitrogen starvation. Front. Microbiol. 10,
to the deep ocean, mechanistic representations of plankton resource 763 (2019).
allocation may be essential in understanding long-term responses 20. Garcia, N. S., Bonachela, J. A. & Martiny, A. C. Growth-dependent
to and feedbacks between marine elemental cycles, the carbon cycle cell size controls interactions between nutrient supply and
and climate. cellular elemental stoichiometry of marine Synechococcus. ISME
J. 10, 2715–2724 (2016).
Online content 21. Inomura, K. et al. A mechanistic model of macromolecular
Any methods, additional references, Nature Research reporting sum- allocation, elemental stoichiometry, and growth rate in
maries, source data, extended data, supplementary information, phytoplankton. Front. Microbiol. 11, 86 (2020).
acknowledgements, peer review information; details of author contri- 22. Omta, A. W. et al. Quantifying nutrient throughput and DOM
butions and competing interests; and statements of data and code avail- production by algae in continuous culture. J. Theor. Biol. 494,
ability are available at https://doi.org/10.1038/s41561-022-01066-2. 110214 (2020).

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

23. Felcmanová, K. et al. Carbon use efficiencies and allocation 44. Elrifi, I. R. & Turpin, D. H. Steady-state luxury consumption and the
strategies in Prochlorococcus marinus strain PCC 9511 during concept of optimum nutrient ratios: a study with phosphate and
nitrogen-limited growth. Photosynth. Res. 134, 71–82 (2017). nitrate limited Selenastrum minutum (Chlorophyta). J. Phycol. 21,
24. Sterner, R. W. & Elser, J. J. Ecological Stoichiometry: The Biology of 592–602 (1985).
Elements from Molecules to the Biosphere (Princeton Univ. Press, 45. Furuya, K. Subsurface chlorophyll maximum in the tropical
2002). and subtropical western Pacific Ocean: vertical profiles of
25. Klausmeier, C. A., Litchman, E., Daufresne, T. & Levin, S. A. phytoplankton biomass and its relationship with chlorophyll a and
Optimal nitrogen-to-phosphorus stoichiometry of phytoplankton. particulate organic carbon. Mar. Biol. 107, 529–539 (1990).
Nature 429, 171–174 (2004). 46. Li, W. K. W., Dickie, P. M., Irwin, B. D. & Wood, A. M. Biomass of
26. Deutsch, C. & Weber, T. Nutrient ratios as a tracer and driver of bacteria, cyanobacteria, prochlorophytes and photosynthetic
ocean biogeochemistry. Annu. Rev. Mar. Sci. 4, 113–141 (2012). eukaryotes in the Sargasso Sea. Deep Sea Res. A 39, 501–519
27. Sharoni, S. & Halevy, I. Nutrient ratios in marine particulate (1992).
organic matter. Sci. Adv. 6, eaaw9371 (2020). 47. Pollehne, F., Klein, B. & Zeitzschel, B. Low light adaptation and
28. Sauterey, B. & Ward, B. A. Environmental control of marine export production in the deep chlorophyll maximum layer in the
phytoplankton stoichiometry in the North Atlantic Ocean. Proc. northern Indian Ocean. Deep Sea Res. II 40, 737–752 (1993).
Natl Acad. Sci. USA 119, e2114602118 (2022). 48. Cornec, M. et al. Deep chlorophyll maxima in the global ocean:
29. Moreno, A. R., Hagstrom, G. I., Primeau, F. W., Levin, S. A. & occurrences, drivers and characteristics. Glob. Biogeochem.
Martiny, A. C. Marine phytoplankton stoichiometry mediates Cycles 35, e2020GB006759 (2021).
nonlinear interactions between nutrient supply, temperature, and 49. Cullen, J. J. Subsurface chlorophyll maximum layers: enduring
atmospheric CO2. Biogeosciences 15, 2761–2779 (2018). enigma or mystery solved? Annu. Rev. Mar. Sci. 7, 207–239 (2015).
30. Shuter, B. A model of physiological adaptation in unicellular 50. Pourtois, J., Tarnita, C. E. & Bonachela, J. A. Impact of lytic phages
algae. J. Theor. Biol. 78, 519–552 (1979). on phosphorus- vs. nitrogen-limited marine microbes. Front.
31. Bonachela, J. A., Allison, S. D., Martiny, A. C. & Levin, S. A. A model Microbiol. 11, 221 (2020).
for variable phytoplankton stoichiometry based on cell protein 51. Talmy, D., Martiny, A. C., Hill, C., Hickman, A. E. & Follows, M. J.
regulation. Biogeosciences 10, 4341–4356 (2013). Microzooplankton regulation of surface ocean POC:PON ratios.
32. Daines, S. J., Clark, J. R. & Lenton, T. M. Multiple environmental Glob. Biogeochem. Cycles 30, 311–332 (2016).
controls on phytoplankton growth strategies determine adaptive 52. Martiny, A. C., Vrugt, J. A. & Lomas, M. W. Concentrations and
responses of the N : P ratio. Ecol. Lett. 17, 414–425 (2014). ratios of particulate organic carbon, nitrogen, and phosphorus in
33. Follows, M. J. & Dutkiewicz, S. Modeling diverse communities of the global ocean. Sci. Data 1, 140048 (2014).
marine microbes. Annu. Rev. Mar. Sci. 3, 427–451 (2011). 53. Martiny, A. C., Vrugt, J. A. & Lomas, M. W. Data from:
34. Bonachela, J. A., Klausmeier, C. A., Edwards, K. F., Litchman, E. & Concentrations and ratios of particulate organic carbon,
Levin, S. A. The role of phytoplankton diversity in the emergent nitrogen, and phosphorus in the global ocean. Dryad https://doi.
oceanic stoichiometry. J. Plankton Res. 38, 1021–1035 (2016). org/10.5061/dryad.d702p (2015).
35. Marshall, J., Hill, C., Perelman, L. & Adcroft, A. Hydrostatic, 54. Vaulot, D., Marie, D., Olson, R. J. & Chisholm, S. W. Growth of
quasi-hydrostatic, and nonhydrostatic ocean modeling. J. Prochlorococcus, a photosynthetic prokaryote, in the equatorial
Geophys. Res. 102, 5733–5752 (1997). Pacific Ocean. Science 268, 1480–1482 (1995).
36. Marshall, J., Adcroft, A., Hill, C., Perelman, L. & Heisey, C. A 55. Vaulot, D. & Marie, D. Diel variability of photosynthetic
finite-volume, incompressible Navier Stokes model for studies of picoplankton in the equatorial Pacific. J. Geophys. Res. Oceans
the ocean on parallel computers. J. Geophys. Res. 102, 5753–5766 104, 3297–3310 (1999).
(1997). 56. Liu, H., Nolla, H. A. & Campbell, L. Prochlorococcus growth rate
37. Follows, M. J., Dutkiewicz, S., Grant, S. & Chisholm, S. W. and contribution to primary production in the equatorial and
Emergent biogeography of microbial communities in a model subtropical North Pacific Ocean. Aquat. Microb. Ecol. 12, 39–47
ocean. Science 315, 1843–1846 (2007). (1997).
38. Ward, B. A., Dutkiewicz, S., Jahn, O. & Follows, M. J. A
size-structured food-web model for the global ocean. Limnol. Publisher’s note Springer Nature remains neutral with regard to
Oceanogr. 57, 1877–1891 (2012). jurisdictional claims in published maps and institutional affiliations.
39. Healey, F. P. Interacting effects of light and nutrient limitation on
the growth rate of Synechococcus linearis (Cyanophyceae). J. Open Access This article is licensed under a Creative Commons
Phycol. 21, 134–146 (1985). Attribution 4.0 International License, which permits use, sharing,
40. Sakshaug, E., Andersen, K. & Kiefer, D. A. A steady state adaptation, distribution and reproduction in any medium or format,
description of growth and light absorption in the marine as long as you give appropriate credit to the original author(s) and the
planktonic diatom Skeletonema costatum. Limnol. Oceanogr. 34, source, provide a link to the Creative Commons license, and indicate
198–205 (1989). if changes were made. The images or other third party material in this
41. Litchman, E., Klausmeier, C. A., Schofield, O. M. & Falkowski, article are included in the article’s Creative Commons license, unless
P. G. The role of functional traits and trade-offs in structuring indicated otherwise in a credit line to the material. If material is not
phytoplankton communities: scaling from cellular to ecosystem included in the article’s Creative Commons license and your intended
level. Ecol. Lett. 10, 1170–1181 (2007). use is not permitted by statutory regulation or exceeds the permitted
42. Marañón, E. et al. Unimodal size scaling of phytoplankton growth use, you will need to obtain permission directly from the copyright
and the size dependence of nutrient uptake and use. Ecol. Lett. 16, holder. To view a copy of this license, visit http://creativecommons.
371–379 (2013). org/licenses/by/4.0/.
43. Moore, C. M. et al. Processes and patterns of oceanic nutrient
limitation. Nat. Geosci. 6, 701–710 (2013). © The Author(s) 2022

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Methods tuting this equation into equation (3) gives:

We incorporated a macromolecular model of phytoplankton
(CFM-Phyto) into the global ocean model (MITgcm). This combined QN = QPro
C
YN∶C
Pro
+ (APRNA μQPro
C
+ QRNA
P,min
)
(5)
model predicts cellular growth rate based on the macromolecular allo- YN∶P + QDNA YN∶C + QChl YN∶C + QNsto
RNA C DNA C Chl N
cation, which in turn is used to determine the elemental stoichiometry
of phytoplankton for the next model time step.
The phytoplankton component of the model is executed using In CFM-Phyto, we resolve three types of protein, photosynthetic,
the following algorithm, which is illustrated graphically in Extended biosynthetic and other:
Data Fig. 2: (1) relate the growth rate and elemental stoichiometry of
Pro_Pho Pro_Bio Pro_Other
phytoplankton based on the macromolecular allocation; (2) evalu- QPro
C
= QC + QC + QC (6)
ate the possible growth rates under four different limiting nutrient
assumptions and select the lowest rate: Liebig’s Law of the Minimum; Photosynthetic proteins represent those in chloroplasts
(3) evaluate storage of non-limiting elements; (4) evaluate excess of largely responsible for light harvesting and electron transport. The
non-limiting elements relative to maximum quotas; (5) based on that model assumes a constant composition of chloroplasts; thus, the
excess, evaluate effective nutrient uptake rate; and (6) evaluate the amount of photosynthetic protein is proportional to the amount of
change in the elemental stoichiometry based on the balance between chlorophyll21:
the growth rate and effective nutrient uptake rate. We describe the
Pro_Pho
procedural details in the following text. Parameter values are listed in QC = APho QChl
C
(7)
Extended Data Table 1. See ref. 21 for further details and justification
of each equation in CFM-Phyto; here we repeat equations essential to Biosynthetic proteins represent proteins related to producing
explain the model used in the current study. new material such as proteins, carbohydrates, lipids, RNAs, DNAs and
other molecules. The models use the following empirically derived
Connecting the elemental stoichiometry and the growth rate relationship21:
The first step of the algorithm is to obtain the relationship between the
Pro_Bio
current elemental stoichiometry and the growth rate (μ). To do that, QC = ABio μ (8)
we use CFM-Phyto21 (Extended Data Fig. 1). The model is based on the
assumption of pseudo-steady state with respect to macromolecular Substituting equations (6)–(8) (in this order) into equation (5)
allocation; in other words, the cellular-scale acclimation occurs rap- leads to the following equation:
idly relative to environmental changes. Laboratory studies show that
Pro_Other
macromolecular re-allocation occurs on the timescale of hours to QN = (APho QChl
C
+ ABio μ + QC ) YN∶C
Pro
days19. This is fast relative to the rates of environmental change in our Pro_Other
+ (APRNA μ (APho QChl + ABio μ + QC ) + QRNA ) YN∶P (9)
coarse-resolution ocean simulations and so steady state solutions21 are C P,min RNA

used to relate growth rate, macromolecular allocation and elemental +QDNA YN∶C + QChl YN∶C + QSto
C DNA C Chl N
stoichiometry, as described in detail below. We first describe the case
of N quota (here defined as QN; moles cellular N per mole cellular C)
in detail, and then we briefly explain the case of P and C quotas as the Empirically, chlorophyll depends on the growth rate and equation
overall procedures are similar. After that, we describe the case with Fe (10) accurately describes the relationship between the growth-rate
quota, which extends the previously published model21 for this study. dependences of chlorophyll under different light intensities21:

Relating N quota and growth rate QChl

C
= AChl μ + BChl (10)
CFM-Phyto describes the allocation of N quota as follows, focusing on
the quantitatively major molecules: with AChl = (1 + E) /vI and BChl = m/vI with E (dimensionless) as a cons­
tant representing growth-rate-dependent respiration, and m (d−1)
QN = QPro
N
+ QRNA
N
+ QDNA
N
+ QChl
N
+ QSto
N
(2) describing maintenance respiration. vI (mol C (mol C in Chl)−1 d−1)
represents chlorophyll-specific photosynthesis rate based on an
where QN is total N quota (per cellular C: mol N (mol C)−1), the terms on established function of light intensity I (μmol m−2 s−1)21,57:
the right-hand side are the contributions from protein, RNA, DNA, chlo-
rophyll and N storage. We use empirically determined fixed elemental vI = vmax
I
(1 − eAI I ) (11)
stoichiometry of macromolecules21 (Extended Data Table 1) to connect
the macromolecular contributions of different elements (here C and P): where vmax
I
is the maximum chlorophyll-specific photosynthesis rate,
e is the natural base and AI is a combined coefficient for absorption
QN = QPro
C
YN∶C
Pro
+ QRNA
P
YN∶P
RNA
+ QDNA
C
YN∶C
DNA
+ QChl
C
YN∶C
Chl
+ QNsto
N
(3) cross-section and turnover time. Substitution of equation (10) into
equation (9) leads to the following quadratic relationship between
Here Yj∶k
l
represents the imposed elemental ratio (elements j and QN and μ:
k) for each macromolecular pool (l). QxC and QxP describe the contribu-
tions of macromolecule x to the total C quota (mol C (mol C)−1) and P QN = aN μ2 + bN μ + cN + QSto
N
(12)
quota (mol P (mol C)−1), respectively.
CFM-Phyto uses the following empirically supported relationship where
to describe QRNA
P
(ref. 21):
aN = APRNA (APho AChl + ABio ) YN∶P
RNA
QRNA
P
= APRNA μQPro
C
+ QRNA
P,min
(4) Pro_Other
bN = (APho AChl + ABio ) YN∶C
Pro
+ AChl YN∶C
Chl
+ APRNA (APho BChl + QC ) YN∶P
RNA

Pro_Other
where APRNA is constant (below, A values represent constant except AChl; cN = BChl YN∶C
Chl
+ (APho BChl + QC ) YN∶C
Pro
see below), μ is growth rate (d−1) and QRNA
P,min
represents the minimum
+QRNA YN∶P + QDNA YN∶C
amount of RNA in phosphorus per cellular C (mol P (mol C)−1). Substi- P,min RNA C DNA

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Relating P quota and growth rate we do not resolve nitrogen-fixing organisms here, Fe allocation (mol
Similarly, CFM-Phyto describes the relationship between the current P Fe (mol C)−1) represents only the first two:
quota QP and μ. P is allocated to its major molecular reservoirs:
QFe = QPho
Fe
+ QSto
Fe
(19)
Thy
QP = QRNA
P
+ QDNA
C
YP∶C
DNA
+ QP + QOther
P
+ QSto
P
(13)
As for equation (7) and equation (14), we relate the allocation of
Similar to equation (7), with the assumption of the constant com- Fe to photosystems to the investment in chlorophyll, QChl
C
:
position of photosynthetic apparatus, the model connects the amount
of the chlorophyll to phosphorus in thylakoid membranes: QPho
Fe
= AFe QChl
Pho C
(20)

Thy
QP = AP∶Chl
Pho
QChl
C
(14) This is a strong simplification because the pigment to photosys-
tem ratio is observed to vary59, but enables an explicit, mechanistically
As for N allocation, substitution of equations (14), (4), (6), (7), (8) motivated representation of Fe limitation, which, a posteriori, results in
and (10) (in this order) into equation (13) leads to a quadratic relation- global scale regimes of iron limitation that resemble those observed43
ship between QP and μ: (Extended Data Fig. 4). With equations (10), (19) and (20), we obtain the
following relationship between QFe and μ:
QP = aP μ2 + bP μ + cP + QSto
P
(15)
QFe = AFe A μ + AFe
Pho Chl
B + QSto
Pho Chl Fe
(21)
where

aP = APRNA (APho AChl + ABio ) Evaluating the growth rate

bP = APRNA (APho BChl +
Pro_Other
QC ) YN∶P + AP∶Chl AChl We assume that the cellular growth rate is constrained by the most
RNA Pho
limiting element within the cell (and its associated functional mac-
cP = QRNA
P,min
+ QDNA
C
YP∶C
DNA
+ AP∶Chl
Pho
BChl + QOther
P romolecules). Thus, at each time step and location, and for each cell
type, the evaluation of growth rate is based on the following two steps:
(1) computation of the growth rate for each element without storage;
Relating C quota and growth rate that is, the case when all of the elemental quotas are allocated
Similarly, CFM-Phyto describes C allocation as follows: to functional macromolecules; and (2) selection of the lowest growth
rate among these; Liebig’s Law of the Minimum. For the first step, we
Plip−Thy
QC = 1 = QPro
C
+ QRNA
C
+ QDNA
C
+ QOther
C
+ QC define μi (i = C, N, P, Fe) as the growth rate, assuming that nutrient i is
(16)
+QCsto + QNsto limiting. Under this condition, QSto i
should be small as element i is
C C
allocated to other essential molecules. We assume that QSto N
is also
small under C limitation because N storage molecules are rich in
where Plip−Thy indicates P lipid in thylakoid membranes. The equation carbon. With these assumptions, the solution for μi is obtained by
represents the allocation per total cellular C in mol C (mol C)−1, so the solving the standard quadratic relationships of equations (12),
sum of the macromolecules in C (QC) becomes 1. Using the imposed (15) and (18) for N, P and C, respectively, neglecting any
elemental ratios of macromolecular pools (Yj∶kl
) we relate the elemental QSto
i
terms:
contributions:
−bi + √b2i − 4ai (ci − Qi )
Thy
QC = QPro + QRNA YC∶P + QDNA + QOther + QP YC∶P + QSto + QSto YC∶N (17) μi = (22)
C P RNA C C Plip C N Nsto 2ai

where QC = 1. For μFe, equation (21) without QSto

Fe
leads to
Following the steps similar to those for the N and P allocations,
substituting equations (14), (4), (6), (7), (8) and (10) (in this order) into QFe − AFe B
Pho Chl
μFe = (23)
equation (17) leads to the following quadratic relationship between AFe A
Pho Chl
cellular C quota QC (=1 mol C (mol C)−1) and μ:
Once the μi values are obtained, we determine the effective growth
QC = aC μ2 + bC μ + cC + QSto
C
+ QSto
N
YC∶N
Nsto
(18) rate, μ, based on the lowest value, which identifies the limiting element
based on current intracellular quotas:
where
μ = min (μN , μP , μC , μFe ) (24)
aC = APRNA (APho AChl + ABio ) YC∶P
RNA
Pro_Other
bC = AChl (1 + APho + AP∶Chl
Pho
YC∶P
Plip
) + ABio + APRNA (APho BChl + QC ) YC∶P
RNA
Evaluating nutrient storage
Pro_Other
cC = (1 + APho + AP∶Chl
Pho
YC∶P
Plip
) BChl + QC In CFM-Phyto, non-limiting nutrients can be stored in an intracellular
reserve21, reflecting commonly observed luxury uptake. Storage is
+QRNA YC∶P + QDNA + QOther
P,min RNA C C evaluated as the difference between the total elemental quota (updated
later) and the functionally allocated portion of that element:
Non_Sto
Relating Fe quota and growth rate QSto
i
= Qi − Qi (25)
In order to capture global scale biogeochemical dynamics including
the iron-limited high-nitrogen, low chlorophyll regimes, CFM-Phyto21 Here QNon_Sto
i
represents the contribution to element i by functional,
is extended to resolve Fe quota and allocation. The model is guided by non-storage molecules. For N, P and C, QNon_Sto
i
is represented by the
a laboratory proteomic study58 in which the major Fe allocations are to non- QSto
i
terms on the right-hand side in equations (12), (15) and (18),
photosystems, storage and nitrogen-fixing enzymes (nitrogenase). As respectively:

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Non_Sto phytoplankton63. Transporter proteins could be represented in a model

Qi = ai μ2 + bi μ + ci (26)
with a finer-scale resolution of the proteome64.

Similarly, for Fe, from equation (21): Differentiating small and large phytoplankton
In this model, ‘small’ phytoplankton broadly represent picocyanobac-
Non_Sto
QFe = AFe A μ + AFe
Pho Chl
B
Pho Chl
(27) teria, which have high nutrient affinities and low maximum growth
rates (for example, Prochlorococcus), whereas ‘large’ phytoplankton
When there is N storage, QSto
C
must be recomputed to consider the represent eukaryotes with higher maximum growth rates (for example,
allocation of C to it: diatoms). The former are associated with a gleaner strategy adapted
to oligotrophic regimes, while the latter are opportunistic, adapted to
Non_Sto
QSto
C
= QC − QC − QSto
N
YC∶N
Nsto
(28) variable and nutrient-enriched regimes. To encapsulate this, the large
phytoplankton have overall higher imposed Vmax i
(~µmaxQi), Ki and vmax
I
than for the small phytoplankton (Extended Data Table 1), consistent
Evaluating the excess nutrient 10
with the previous models (for example, ref. ). In addition, the larger
Storage capacity for any element is finite and we define excess nutrient cells are assigned a higher Qmax
P
following the observed trends (Fig. 1
as a nutrient (N, P, Fe) that is in beyond an empirically informed, and Extended Data Table 1).
imposed maximum phytoplankton storage capacity. Excess nutrient
is assumed to be excreted (see below). Excess of element i ( QExc i
) is Computing the change in the elemental stoichiometry
computed: The computation of the change in the elemental quotas is done based
on the balance between the effective nutrient uptake rate and the loss
QExc
i
= max (Qi − Qmax
i
, 0) (29) of nutrient to the new cells:

dQi
where Qmax
i
is maximum capacity for nutrient i. For N, CFM-Phyto com- = VEff − μQi (34)
dt i
putes Qmax
i
as a sum of non-storage molecules and prescribed maximum
nutrient storing capacity according to model–data comparison21:
This change in the elemental quotas based on the cellular pro-
Non_Sto Sto_max
Qmax
i
= Qi + Qi (30) cesses and the passive transport of elements in phytoplankton by the
flow field created by MITgcm governs the elemental stoichiometry of
Laboratory studies suggest that when P is not limiting, the phos- the next time step at a specific grid box, as in other versions of ecologi-
phorus quota maximizes to a value that is almost independent of cal models with MITgcm10.
growth rate21,39,44. Storage of each element is finite and the upper limit
to storage is imposed by specifying the maximum cellular quotas ( Qmax P
Calculation of CV values
(ref. 21) and also Qmax
Fe
) with size and taxonomic dependencies (for exam- We computed the CV values based on the following equation:
ple, refs. 27,41). Thus, the maximum storage is represented by the differ-
σ
ence between the prescribed maximum quota and the actual quota21: CV = (35)
x̄
Sto_max
Qi = Qmax
i
− Qi (31) where σ is the standard deviation and x̄ is the mean. The purpose of this
computation is to quantify the latitudinal variation of the averaged
elemental stoichiometry. Thus, we used the averaged values for each
In the case where QSto
i
computed in the previous section exceeds latitude (as plotted in Fig. 2) for the calculation of σ and x̄ .
Sto_max
Qi the value of QSto
, i
is replaced by QSto_max
i
and the difference is
placed in the excess pool, QExc
i
. MITgcm-CFM
The biogeochemical and ecological component of the model resolves
Computing effective nutrient uptake rate the cycling of C, P, N and Fe through inorganic, living, dissolved and
One factor that influences the cellular elemental quota is the effective particulate organic phases. The biogeochemical and biological trac-
nutrient uptake rate (mol i (mol C)−1 d−1) of N, P and Fe, which we define ers are transported and mixed by the MIT general circulation model
as follows: (MITgcm)35,36, constrained to be consistent with altimetric and hydro-
graphic observations (the ECCO-GODAE state estimates)65. This
QExc
i three-dimensional configuration has a coarse resolution (1° × 1° hori-
VEff
i
= Vi − (32)
τExu
i zontally) and 23 depth levels ranging from 5 m at the surface to 5450 m
at depth. The model was run for three years, and the results of the third
where Vi (mol i (mol C)−1 d−1) is nutrient uptake rate and the second term year were analysed, by which time the modelled plankton distribution
represents the exudation of the excess nutrient based on the timescale becomes quasi-stable. Equations for the biogeochemical processes
τExu
i
(d−1). For Vi, we use Monod kinetics60,61: are as described by equations and parameters in previous studies10,38.
Here, however, we include only nitrate for inorganic nitrogen, and do
[i]
Vi = Vmax
i
(33) not resolve the silica cycle. We simulated eukaryotic and prokaryotic
[i] + Ki
analogues of phytoplankton (as ‘large’ and ‘small’ phytoplankton). The
eukaryotic analogue has a higher maximum C fixation rate for the same
where Vmax
i
is maximum nutrient uptake, [i] (mmol m−3) is the environ- macromolecular composition and higher maximum nutrient uptake
mental concentration of nutrient i and K i (mmol m −3) is the rates, but also has overall higher half-saturation constants for nutrient
half-saturation constant of i. Previous models have resolved the rela- uptake. We used light absorption spectra of picoeukaryotes, which sits
tionship between nutrient uptake and allocation to transporters31,62. in-between small prokaryotes and large eukaryotes10. In MITgcm, the
Here we do not explicitly resolve allocation to transporters, as prot- mortality of phytoplankton is represented by the sum of a linear term
eomic studies indicate that it is a relatively minor component of the (ml), representing sinking and maintenance losses, and quadratic terms
proteome compared with photosystems and biosynthesis in representing the action of unresolved next-trophic levels66,67, implicitly

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

assuming that the higher-trophic-level biomass scales with that of its 66. Steele, J. H. & Henderson, E. W. The role of predation in plankton
prey. We assumed that the latter term is small to avoid introducing models. J. Plankton Res. 14, 157–172 (1992).
additional uncertainties. Similarly, we do not resolve the stoichiometric 67. Aumont, O., Ethé, C., Tagliabue, A., Bopp, L. & Gehlen, M.
effects of prey selection due to the nutritional status of prey, or viral PISCES-v2: an ocean biogeochemical model for carbon and
partitioning of nutrients in the environment50. Atmospheric iron depo- ecosystem studies. Geosci. Model Dev. 8, 2465–2513 (2015).
sition varies by orders of magnitude around the globe and has a large 68. Mahowald, N. M. et al. Atmospheric iron deposition: global
margin of uncertainty, including the bio-availability of the deposited distribution, variability, and human perturbations. Annu. Rev. Mar.
iron, which in turn depends on the source and chemical history of the Sci. 1, 245–278 (2009).
deposited material68. To realize a realistic global net primary produc- 69. Garcia, H. E. et al. World Ocean Atlas 2018, Volume 4: Dissolved
tion, we doubled the atmospheric iron input from ref. 10; this factor is Inorganic Nutrients (Phosphate, Nitrate and Nitrate+Nitrite, Silicate)
well within the uncertainty of the iron supply estimates. Each of the (NOAA Atlas NESDIS 84, 2018); https://archimer.ifremer.fr/
two phytoplankton groups has variable C:N:P:Fe as determined by doc/00651/76336/
the component macromolecules at each time step. The pools of C, N,
P and Fe are tracked within the modelled three-dimensional flow fields. Acknowledgements
We thank H. Frenzel and K. Bryan for computational support, E. J.
Data availability Zakem for sharing a Bash script for efficient MITgcm runs, and E.
The model output from this study is available at https://doi. M. Howard, A. J. Margolskee, H. Frenzel and J. L. Penn for helpful
org/10.6084/m9.figshare.21197578. feedback. This study was supported by the US National Science
Foundation (OCE-2048373, C.D. and K.I.; OCE-1558702, M.J.F.), the
Code availability Gordon and Betty Moore Foundation (GBMF no. 3775, C.D.; GBMF
The physical model and ecosystem code are available at http://mitgcm. no. 3778, M.J.F.) and the Simons Foundation (Simons Postdoctoral
org and https://gitlab.com/darwinproject/gud, respectively, and the Fellowship in Marine Microbial Ecology, Award 544338, K.I.; Simons
specific modification for this study can be downloaded from https:// Collaboration on Ocean Processes and Ecology, Award 329108,
zenodo.org/record/4730684. M.J.F.; CBIOMES, Award 549931, M.J.F.). We thank these foundations
for their support.
References
57. Cullen, J. J. On models of growth and photosynthesis in Author contributions
phytoplankton. Deep Sea Res. A 37, 667–683 (1990). K.I. and O.J. developed a model with advice from S.D. and M.J.F. K.I.
58. Saito, M. A. et al. Iron conservation by reduction of performed numerical simulations with advice from C.D., O.J., S.D.
metalloenzyme inventories in the marine diazotroph and M.J.F. K.I. and C.D. analysed the results and wrote the original
Crocosphaera watsonii. Proc. Natl Acad. Sci. USA 108, 2184–2189 manuscript, which was revised by K.I., C.D., S.D. and M.J.F.
(2011).
59. Strzepek, R. F., Boyd, P. W. & Sunda, W. G. Photosynthetic Competing interests
adaptation to low iron, light, and temperature in Southern Ocean The authors declare no competing interests.
phytoplankton. Proc. Natl Acad. Sci. USA 116, 4388–4393 (2019).
60. Monod, J. The growth of bacterial cultures. Annu. Rev. Mar. Sci. 3, Additional information
371–394 (1949). Extended data is available for this paper at https://doi.org/10.1038/
61. Healey, F. P. Slope of the Monod equation as an indicator of s41561-022-01066-2.
advantage in nutrient competition. Microb. Ecol. 5, 281–286
(1980). Supplementary information The online version contains
62. Smith, S. L., Yamanaka, Y., Pahlow, M. & Oschlies, A. Optimal supplementary material available at https://doi.org/10.1038/s41561-
uptake kinetics: physiological acclimation explains the pattern of 022-01066-2.
nitrate uptake by phytoplankton in the ocean. Mar. Ecol. Prog. Ser.
384, 1–12 (2009). Correspondence and requests for materials should be addressed to
63. McKew, B. A., Metodieva, G., Raines, C. A., Metodiev, M. V. & Keisuke Inomura.
Geider, R. J. Acclimation of Emiliania huxleyi (1516) to nutrient
limitation involves precise modification of the proteome to Peer review information Nature Geoscience thanks C. Mark Moore and
scavenge alternative sources of N and P. Environ. Microbiol. 17, the other, anonymous, reviewer(s) for their contribution to the peer
4050–4062 (2015). review of this work. Primary Handling Editors: Tamara Goldin and Xujia
64. Zavřel, T. et al. Quantitative insights into the cyanobacterial cell Jiang, in collaboration with the Nature Geoscience team.
economy. Elife 8, e42508 (2019).
65. Wunsch, C. & Heimbach, P. Practical global oceanic state Reprints and permissions information is available at
estimation. Phys. D 230, 197–208 (2007). www.nature.com/reprints.

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Extended Data Fig. 1 | Schematic of CFM-Phyto component of this study. and P thyla represent precursor carbohydrates and P in thylakoid membranes,
Orange, red, blue, and black arrows represent C, N, P, Fe fluxes, respectively. respectively. Green space and orange outline represent the intracellular space
Yellow, red, blue, gray squares are molecules that mainly influence C, N, P and Fe and cell membrane layers, respectively.
budgets, respectively. Dashed squre represents photosynthetic molecules. CH

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Extended Data Fig. 2 | Solution flow of CFM-Phyto part of the modeling. In CFM-Phyto, cellular elemental quotas Qi and growth rate μ are linked by macromolecular
allocation.

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Extended Data Fig. 3 | Annual mean NO3− and PO43− concentrations at the surface. (a)(b) are model output. (c)(d) climatological average from World Ocean Atlas69.

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Extended Data Fig. 4 | Model-emergent and observed nutrient limitations observed primary nutrient limitation from field data compilation43, which did
for phytoplankton growth. (a) Small phytoplankton. (b) Large phytoplankton. not include an assessment for C limitation.
Color field shows the limiting nutrient diagnosed from MITgcm. Circles are

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Extended Data Fig. 5 | Distribution of the elemental ratios measured in particulate organic matter from global field data. (A) N:C. (B) N:P. A small fraction of data
(<1%) lies outside of the plotted range. Data are based on the global compilation52,53.

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Extended Data Fig. 6 | Differences in N:C between small phytoplankton and large phytoplankton. (A) The difference in percent. (B) N:C in small phytoplankton.
(C) N:C in large phytoplankton.

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Extended Data Fig. 7 | N:P variability among phytoplankton in modeled Distinct curves of biomass versus elemental ratios are separately plotted for
population and species observations. Model output of total model small (blue) and large (red) size classes, and compared against the median values
phytoplankton biomass (Gmol C) is binned according to its local N:P (mol mol−1). of traits observed for small (prokaryote, blue) and large (eukaryote, red).

Nature Geoscience
Article https://doi.org/10.1038/s41561-022-01066-2

Extended Data Table 1 | Parameter values

-: Dimensionless (for example, mol mol−1), *Values from the original CFM-Phyto study21, #Empirical stoichiometric parameters21. The rest of the parameters were chosen to reproduce realistic
nutrient limitation, reasonable net primary production, and observed magnitude of elemental stoichiometry. We have given overall higher Vmax and K values for large phytoplankton as typically
modeled (for example, ref. 10), and higher Qmax
P
for large phytoplankton based on observations26.

Nature Geoscience