Shifts in regional production as a driver of future global ocean production stoichiometry

Using a global ocean biogeochemistry model, we examined three drivers of global ocean production C:N:P ratio: flexible phytoplankton stoichiometry, phytoplankton community composition, and regional production shifts. For a middle-of-the-road warming scenario (SSP2), the model predicts a substantial increase in the global export C:P ratio from 113:1 to 119:1 by the year 2100. The most important physiological driver of this stoichiometric change is the effect of the worldwide warming on cyanobacteria, followed by the effect of phosphate depletion on eukaryotes in the Southern Ocean. Also, there is a modest global shift in the phytoplankton community in favor of cyanobacteria at the expense of eukaryotes with a minimal effect on the global production stoichiometry. We find that shifts in the regional production, even in the absence of any change in phytoplankton stoichiometry or taxonomy, can change the global production C:N:P ratio. For example, enhancing the production in the polar waters, which typically have low C:N:P ratios, will have the effect of lowering the global ratio. In our model, the retreat of Antarctic sea ice has this very effect but is offset by production changes downstream and elsewhere. This study thus provides an understanding of how regional production changes can affect the global production C:N:P ratio. However, the current literature indicates substantial uncertainty in the future projections of regional production changes, so it is unclear at this time what their net effect is on the global production C:N:P ratio. Finally, our model predicts that the overall increase in the carbon content of organic matter due to flexible C:N:P ratio helps to stabilize carbon export in the face of reduced nutrient export (i.e. the decrease in C export is ~30% smaller than expected from the decrease in P export by 2100) but has a minimal effect on atmospheric CO2 uptake (~1%).


Introduction
Traditionally, the export production of carbon (C) in global ocean biogeochemical models has been directly linked to the export of macronutrients phosphorus (P) and nitrogen (N), assuming a fixed elemental stoichiometry in organic matter (e.g. Redfield 1934). For example, a ratio adopted by the Ocean Carbon Cycle Model Intercomparison Project (Najjar et al 2007) is C:N:P = 117:16:1. Today the ocean modeling community is in the midst of adopting alternate formulations of phytoplankton stoichiometry that allow for flexibility in light of observations that the C:N:P ratio actually varies with discernable patterns spatially and by species (e.g. Martiny et al 2013). Those observations indicate that C:P and N:P ratios are much higher in smaller phytoplankton like cyanobacteria and in oligotrophic waters, well above the canonical Redfield ratios of C:P = 106 and N:P = 16, and much lower in larger cells like diatoms in eutrophic waters. In state of the art Climate Model Intercomparison Project 5 (CMIP5) (Bopp et al 2013) and CMIP6 models (Arora et al 2020), some models have fully flexible C:N:P (e.g. CMCC/PELAGOS), while others have partially flexible C:N:P (CESM/-MARBL, HadGCM/Diat-HadOCC, GFDL/TOPAZ), or have fixed C:N:P (CESM/BEC, IPSL/PISCES, MPI/HAMOCC, CanESM5/CMOC).
In the context of climate change, the flexibility in the phytoplankton stoichiometry is potentially an important driver of the global ocean carbon cycle, as C export can be decoupled from N and P export.
Under warming, the expected strengthening of water column stratification in much of the global ocean will likely diminish the vertical supply of P and N to the euphotic layer and thus reduce primary production (e.g. Boyd et al 2015, Hutchins andFu 2017). A fixed C:N:P ratio would predict a reduction in C export proportional to the P and N export reductions, leading to higher atmospheric pCO 2 . However, flexible stoichiometry might reasonably predict that organic matter become more C-rich under reduced nutrient supply (Galbraith andMartiny 2015, Garcia et al 2018). This may in turn help maintain C export even as P and N export declines. Thus, flexible stoichiometry can potentially buffer or reduce the simulated sensitivity of C export to changes in N or P supply compared to estimates obtained by fixed stoichiometry.
Two previous studies investigated such buffering effect under future warming conditions using global 3D ocean models Matsumoto 2017, Kwiatkowski et al 2018). In these studies, the formulation of flexible stoichiometry is different. In our previous work, we introduced an empirically based, power law formulation of C:P as a function of only P in a model with two phytoplankton types (Tanioka and Matsumoto 2017). It recognized that C:P is much more variable than C:N, and that P is a major driver of the C:P ratio (Galbraith and Martiny 2015). In Kwiatkowski et al (2018), the formulation of flexible stoichiometry is theoretical and based on the chain model of Pahlow and Oschlies (2009). It calculates phytoplankton growth and nutrient uptake by optimizing the allocation of N and P to specific cellular functional compartments so as to maximize growth. Despite these substantial differences, the two studies both conclude that the buffering effect is small for both atmospheric CO 2 uptake (Kwiatkowski et al 2018) and C export (Tanioka and Matsumoto 2017).
In this study, we are motivated by two reasons to revisit the buffering effect and changes in the global ocean production stoichiometry under global warming. The first motivation is that previous studies focused on the phytoplankton physiology and taxonomy as the main drivers of global ocean production stoichiometry. By physiology, we refer to the flexibility in the phytoplankton C:N:P ratio that responds to environmental changes. By taxonomy, we refer to the fact that the phytoplankton community C:N:P ratio is a weighted average of the constituent phytoplankton stoichiometry. Thus the C:N:P ratio of the most abundant phytoplankton type will have the strongest imprint on the community ratio. Here we are motivated to quantitatively examine a third driver: the possibility that shifts in regional ocean productivity can drive the global mean production C:N:P ratio even if phytoplankton physiology and community composition remain unchanged. For example, the global mean ratio would rise if the oligotrophic subtropical gyres, which typically have high C:N:P ratios, were to become more productive. All else being the same, the high ratios of the gyres would make a larger contribution to the global mean in that case.
The second motivation for revisiting the future global ocean C:N:P ratio is that we now have a more complete power law model of flexible stoichiometry so that the C:N:P ratio is a function of a fuller set of environment variables. In our previous study (Tanioka and Matsumoto 2017), the lack of C:N flexibility and temperature as an environmental driver were noticeable deficiencies. Also, this expanded model of stoichiometry is applied to three phytoplankton types instead of two, as in the previous study. Thus, we are able to more credibly quantify the contributions that physiology and community composition make to the community stoichiometry. Finally, as noted by Kwiatkowski et al (2018), theirs was the only study so far that has investigated the effects of the full C:N:P ratio of multiple phytoplankton under global warming, so a revisit of the same issue, especially with a different formulation of flexible stoichiometry, is needed and timely.

Power law model of C:N:P
In this study, the power law formulations of flexible P:C and N:C respectively take the form: These are the uptake P:C and N:C ratios by phytoplankton during organic matter production. From these two ratios, the N:P uptake ratio can be derived. The environmental drivers are [PO 4 ] and [NO 3 ] concentrations, temperature (T), and light (I). The exponents s ratio driver indicate the sensitivity and sign of the drivers in controlling the P:C or N:C ratio. The reference values for the environmental drivers and the base ratios are indicated with subscript 0. (See supplemental information table S1 (available online at stacks.iop.org/ERL/15/124027/mmedia) caption for further details) A full description of the stoichiometry model represented by equations (1) and (2) is given elsewhere . Briefly, we determined s ratio driver through an exhaustive meta-analysis of laboratory studies of marine phytoplankton stoichiometry (Tanioka and Matsumoto 2020), which we translated to the three phytoplankton functional types (PFTs) represented in our global ocean model (table S1). The three PFTs are eukaryotes, cyanobacteria, and diazotrophs. Note that non-zero s ratio driver exists only for those combinations of ratios and drivers where the relationship was statistically significant in the meta-analysis. For example, s P:C [PO4] =0.58 for eukaryotes means that a 1% increase in [PO 4 ] leads to a 0.58% increase in P:C. The eukaryotes have the largest s P:C [PO4] of all the PFTs, which could be understood in terms of the larger storage capacity of eukaryotes that allows luxury storage of P. Also, s P:C T = −8.0 for cyanobacteria and diazotrophs indicates that a 1% increase in temperature (in terms of • K) leads to an 8% drop in P:C. As discussed in the meta-analysis by Tanioka and Matsumoto (2020), one of multiple reasons why phytoplankton might become more C rich under warmer conditions is increased metabolic stimulation of C fixation. The interested reader is referred to the metaanalysis for a discussion of the underlying biochemical mechanisms.
Note that s ratio driver =0 when there is no environmental dependence on a ratio (e.g. eukaryote N:C ratio is not dependent on [PO 4 ], so s N:C [PO4] =0). When a ratio has no dependence on any of the environmental drivers, all terms with the exponents in equations (1) and (2) would become 1, and the ratio would be given by the base ratio (e.g. [N:C] = [N:C] 0 ).

Global ocean model
The new power law model of flexible C:N:P is implemented in MESMO (Matsumoto et al 2008, an earth system model of intermediate complexity (EMIC). The model components used here represent global ocean physics and biogeochemistry, an energy-moisture balanced atmosphere, and a thermodynamic and dynamic sea ice. The ocean model domain has an equal-area grid with 10 • increments in longitude and is uniform in the sine of latitude. It has 16 vertical levels with the top two levels representing the surface 100 m where biological production takes place. Export production is calculated at 100 m. The biogeochemistry model computes the net primary production (NPP), which depends on nutrients, light, temperature, and the mixed layer depth. Nutrients are P, N, iron (Fe), and silicic acid (H 4 SiO 4 or simply Si). Both P and Si have a closed system. In contrast, Fe and N have external sources, aeolian Fe flux and N-fixation, which are balanced by Fe burial and denitrification, respectively. The limiting nutrient during uptake is determined by the Liebig's rule of the minimum relative to the uptake stoichiometry. Remineralization of the sinking particulate organic matter depends on the sinking speed, temperature, and dissolved oxygen concentration. We recently described the power law C:N:Penabled version of MESMO and its equilibrium state . Briefly, equations (1) and (2) are applied to all ocean grid points in the top 100 m to all PFTs represented. One of the PFTs is eukaryotes, which include diatoms and are fast and opportunistic growers under nutrient replete conditions (e.g. largest maximum growth rate). Cyanobacteria represent Prochlorococcus and Synechococcus, which are dominant in oligotrophic waters (i.e. smaller halfsaturation constants in the Michaelis-Menten uptake kinetics). Diazotrophs are nitrogen fixers with slow growth rates but their growth is not limited by nitrate.
Here, the phytoplankton community stoichiometry is calculated as the average of the three PFT C:N:P ratios, weighted by their NPP in terms of P (NPPp). Therefore, the community C:N:P at any one grid point within the production layer is most closely aligned with the C:N:P ratio of the most productive PFT. The 2D maps of community stoichiometry we present below are phytoplankton community stoichiometry integrated over the 100 m thick production layer: (3) where i is the PFT index (i.e. i = 1 for eukaryotes, i = 2 for cyanobacteria, i = 3 for diazotrophs), k is the vertical level index in the production layer. Within the production layer, the top level has lower nutrient concentrations, higher temperatures, and more light than the deeper level. According to the power law model dependence (table S1), the combination of these conditions leads to higher C:N:P ratios in the top level than the deeper level everywhere. NPPp can be higher or lower in the top level versus the deeper level. Thus, C:N:P com , given by equation (3), is an NPPp-weighted, vertically averaged C:N:P ratio. The global C:N:P com is obtained by extending this averaging to the east-west and north-south directions to the global ocean.

Simulations under a future warming scenario
We carried out seven model experiments: five that use the power law formulation of C:N:P, and two 'Redfield' experiments with a fixed C:N:P = 117:16:1 (supplemental information table S2). The equilibrium model states for both formulations of C:N:P were obtained by running the model for 5000 years under preindustrial boundary conditions until the model reached a steady state. All seven runs are 750 years long, starting from their respective equilibrium runs in 1750 and terminate in 2500.
Experiments Control and C_Redfield are the control runs for each C:N:P formulation, where the equilibrium runs were simply continued for 750 years without drift.
For the global warming experiments, we use the Shared Socioeconomic Pathway 2, a new 'middleof-the road' scenario of future warming (Riahi et al 2017). The radiative forcing stabilizes at~6.5 Wm −2 with pCO 2 at~600 µatm by the year 2100. Throughout the runs, there is no change in either the aeolian Fe input or the external supply of N attributable to anthropogenic sources.
Experiment GW is our standard global warming experiment using the power law model of C:N:P. There are three sensitivity experiments GW_CNP, GW_COM, and GW_ICE of the standard experiment. As explained more in detail below, masks of the equilibrium seasonal cycles of phytoplankton C:N:P ratio, phytoplankton community composition, and sea ice were applied to each grid point in experiments GW_CNP, GW_COM, and GW_ICE, respectively. Thus, they are designed to isolate the effects of phytoplankton stoichiometry, community composition, or sea ice respectively in GW. Finally, experiment GW_Redfield is the global warming run with a fixed C:N:P ratio.

Results
We use a description of the Control run to demonstrate the three main controls on the global mean C:N:P com and establish the performance of the model in terms of global averages. Figure 1 shows that the phytoplankton community C:P (C:P com ) uptake ratio at equilibrium is elevated in the subtropical gyres and reduced in more eutrophic waters due to both the physiological and taxonomic controls. In the polar regions, eukaryotes with the lowest C:P ratios (~60) are the most dominant PFT, thus C:P com is well below 106. PFT fraction indicates the contribution that the PFT makes to total NPP in terms of C. In the subtropical gyres, the cyanobacteria with very high C:P ratios (~250) are dominant, thus elevating C:P com . The community composition shown in figure 1(a) thus indicates the taxonomic control on C:P com . Figure 1(b) shows that the equilibrium C:P becomes lower towards the poles for all PFTs as P becomes more abundant. This P effect is largest for eukaryotes, which have the largest s P:C [PO4] (see table S1). For cyanobacteria and diazotrophs, there is an additional temperature effect (see s P:C T , table S1) that elevates C:P towards the warmer, lower latitudes. The zonal mean C:P for each PFT shown in figure 1(b) thus illustrates the physiological control on C:N:P com . The third control is the regional production shifts. In figure 1(c), the global distribution C:P com at equilibrium is shown with a black line that indicates the global mean C:P com = 146:1. Even if the local phytoplankton stoichiometry and their community composition remain fixed everywhere, any NPPp increase in waters where C:P com > 146 would tend to raise the global mean C:P com ratio, because there would be a greater weighting of the elevated C:P com ratios to the global mean. Conversely, an increase in NPPp where C:P com < 146 would tend to lower the global mean C:P com ratio.
At equilibrium, the global mean C:N:P com is 146:19:1 (table S2), which is nearly identical to the observed biomass-weighted ratio of 146:20:1 (Martiny et al 2013). In comparison, the global mean export C:N:P in the model is 113:16:1, similar to the Redfield ratio, which still serves as a good global mean stoichiometry. The overall equilibrium model state is thus consistent with our current understanding of the large scale C:N:P variability across the global ocean (Martiny et al 2013, Teng et al 2014, Wang et al 2019 and serves as a good starting point for the future simulations. See Matsumoto et al (2020) for more details of the Control run including PFT nutrient limitation maps.
The reason that the export C:N:P ratio is lower than the C:N:P com ratio (table S2) is that the two levels within the production layer have different C:N:P ratios and contribute differently to export production and to water column NPPp. The water column NPPp is the sum of production in the top two levels. In contrast, export production would reflect a relatively smaller contribution from the top level than the deeper level, because the top level's organic matter would undergo a more extensive remineralization than the deeper level's organic matter before sinking below 100 m. In the model, there is no preferential remineralization of nutrients over carbon in organic matter, so the uptake and export C:N:P ratios are identical in the Redfield experiments.
The physical changes simulated in experiment GW are illustrated in figure 2. The pCO 2 forcing of the SSP2 scenario shows a rapid rise during the 2000-2100 period and peaks in about year 2200 (figure 2(a)). With pCO 2 approximately doubling, the model simulates a global mean surface temperature rise of 3 • C-4 • C ( figure 2(b)). The simulated amplitude of warming is consistent with the generally accepted equilibrium climate sensitivity of 1.5 • C-4.5 • C (Flato et al 2013). As expected, this surface warming promotes both ice melting, water column stratification, and reduced vertical mixing. Global annual mean sea ice coverage thus decreases from approximately 10% to 6.5% ( figure 2(c)). There is a substantial loss of sea ice in the Southern Ocean (SO) (figure 2(f)). The 32% loss of Antarctic sea ice cover by 2100 is consistent with a 33% mean loss projected by the CMIP3 model ensemble using a middle-of-therange SRES A1B forcing (Bracegirdle et al 2008) and a 29% mean loss projected by the most recent CMIP6 models using the same new SSP2 forcing with a large range (Roach et al 2020). Stronger stratification leads to a weakening of the Atlantic meridional overturning circulation (MOC) by from 12 Sv to 8 Sv (~30%), Corresponding to these physical changes are biogeochemical changes in experiment GW (table S2, figure 3). Due to enhanced ocean stratification and reduced vertical mixing, the global mean surface [PO 4 ], [NO 3 ], and [Si] become lower ( figure 3(a)). The largest drop is seen in [Si], which reaches~30% of the preindustrial levels. Concentrations of P, N, and Si decline most prominently in SO with little change elsewhere. Total dissolved iron (FeT), the sum of free dissolved iron and ligand-bound iron, also declines in SO but increases elsewhere. Globally, there is thus a slight surface FeT accumulation, which has been predicted to occur in the future as aeolian input of Fe continues while biological production declines (Boyd et al 2015, Hutchins andBoyd 2016).
Consistent with the reduced availability of surface P, N, and Si, export production declines globally but differs between C and P. By the year 2100, P export declines by 14% to 86% of the preindustrial level, while C export declines by only 10% from 9.4-8.5 Pg-C yr −1 in experiment GW ( figure 3(b), table S2). In the equivalent experiment GW_Redfield, the fractional change in C and P export is exactly the same. Thus, the smaller reduction in C export compared to P export in GW is a direct consequence of flexible stoichiometry, namely that organic matter being exported becomes more C-rich. By the year 2100, C:N:P com increases from 146:19:1 in Control to 157:20:1 in GW, a change of +11 in C:P com (table S2).
As shown in figure 3(c), the global mean increase in C:P by 2100 is much greater for cyanobacteria and diazotrophs (~11 in molar ratio) than for eukaryotes (~7). Maps of C:P changes for each PFT and the community are provided in supplemental information figure S1. The overall larger C:P response in the smaller PFTs is largely attributed to their dependence on temperature (i.e. s P:C T ), which is absent in eukaryotes (table S1). The P effect on C:P cannot explain the higher C:P in the smaller PFTs, because s P:C [PO4] is in fact greatest for the eukaryotes.
Taxonomically, there is relatively little change in the global mean PFT fractions ( figure 3(d)).
The small decrease in the fraction of eukaryotes (~3%) is largely balanced by an increase in the cyanobacteria fraction, with diazotrophs changing little. However, local changes in the PFT fractions can be quite large, as high as 20%-40% in some polar waters with cyanobacteria replacing eukaryotes (supplemental information, figure S2).
The physiological and taxonomic controls on the global C:N:P com were further investigated by conducting the two experiments GW_CNP and GW_COM. In GW_CNP, we repeated GW, except we applied seasonal masks of C:N:P for each PFT created from the Control run throughout the production layer. Thus at every point, the annual cycle of C:N:P for each PFT remains the same throughout the postindustrial run. Similarly, we applied seasonal masks of community composition in experiment GW_COM, so that the annual cycle of the PFT fractions remain the same, while allowing NPP to change. These two runs are thus intended as controlled experiments of the GW run, where GW_CNP is different from GW in only the physiological control on local C:N:P, and where GW_COM is different from GW in only the taxonomic control on local C:N:P. Interpreting these simulations require some caution in that while the masks work largely as intended (supplemental information, figure S3), GW_CNP and GW_COM are not perfect controlled experiments of GW, especially GW_COM . Phytoplankton growth, stoichiometry, nutrients, and production are all interdependent, so fixing one of them by the masks can lead to changes elsewhere (e.g. export production, table S2).
The C:N:P com ratio is 157:20:1 in GW, 145:19:1 in GW_CNP, and 159:21:1 in GW_COM for the year 2100 (table S2). At face value, this indicates that the physiological response of phytoplankton under warming changes the global C:P com by +12 (=157-145), nearly the same as the overall C:P com change predicted in experiment GW. The sign of this change is consistent with our understanding of the environmental drivers (table S1) in that both future warming and increased nutrient limitation would elevate uptake C:P. On the other hand, allowing for local shifts in the community composition under warming changes C:P com by −2 (=157-159). This suggests that the local shifts in the community composition have a modest effect on the global phytoplankton uptake stoichiometry.
The importance of the third driver of the global mean C:N:P com ratio, the effect of regional production changes, appears to be also modest in our experiment GW (figure 4). Figure 4(a) is a translation of the C:P com map for the Control run ( figure 1(c)) to a map of the sign of the effect of production changes on the global stoichiometry. As noted above, an increase in production in waters where C:P com < global mean (sign = −1, generally polar regions) would lower the global mean C:P com . Conversely, an increase in production in waters where C:P com > global mean (sign = +1, generally lower latitudes) would raise the global mean C:P com . The effect of local production changes on the global mean C:N:P com can be visualized in figure 4(c), which is given by the product of the local sign of the effect ( figure 4(a)) and the local ∆NPPp by the year 2100 ( figure 4(b)).
As indicated by figure 4(b), NPPp becomes lower in the North Atlantic but has a dipole pattern in SO with an increase in the poleward waters and a reduction in the equatorward waters. We investigated this ∆NPP dipole pattern in SO using the sensitivity experiment GW_ICE (figure 5). As noted above, seasonal masks of sea ice created from the equilibrium run were applied to GW_ICE. Physical changes associated with the presence of sea ice, such as warming and stratification in SO, are unaffected by the mask, so the difference between GW and GW_ICE isolates the effect of sea ice capping. Figure 5(a) shows that NPPp increases by 6 Tg-P yr −1 (=54-48) in Figure 4. Diagnosis of the effect of regional production changes on the global mean C:Pcom ratio for the year 2100 in experiment GW. (a) Anomaly of the Control run C:Pcom ratio from the global mean, which is 146 and indicated by the black line, (b) ∆NPPp in experiment GW (2100-1750) in unit of mole-P m −2 yr −1 , and (c) diagnosed effect of the regional production changes. In (a), the sign of the effect of the regional production changes is positive (+1) for locations with local C:Pcom > global mean (red) and negative (−1) for locations with C:Pcom < global mean (blue). The diagnosed effect (c) is given by the product of the sign of the effect (a) and the production changes (b). experiment GW between 1750 and 2100 in the poleward waters of SO (>55 • S). For convenience, hereafter we refer to >55 • S as 'Antarctic.' In experiments Control and GW_ICE, the Antarctic NPPp remains unchanged at~48 Tg-P yr −1 . Thus, the +6 Tg-P yr −1 change in GW is attributable to the Antarctic sea ice retreat. After 2100, the Antarctic NPPp continues to increase in GW but diverges in experiments GW_ICE and Control. Therefore, the NPPp increase in GW after 2100 is also driven by factors other than sea ice (e.g. nutrients, warming, and increased stratification).
As shown in figure 4(c), the effect of local production changes on the global mean C:N:P com is negative in the Antarctic waters, indicating that the Antarctic sea ice retreat lowers the global C:P com ratio. In contrast, changes in the subantarctic production have a positive effect, because of the dipole pattern in ∆NPPp. Interestingly, the regional production effect of the North Atlantic on the global C:P com ratio also has a dipole pattern (figure 4(c)) even while ∆NPPp is uniformly negative ( figure 4(b)). The effect is positive in the subpolar North Atlantic (south of Greenland) and negative further south at~40 • N, because the global mean C:P com (i.e. black line, figure 4(a)) runs between the two regions We estimated the net effect of these regional production shifts on the global C:P com ratio by calculating the NPPp-weighted ratio following equation (3) using the local C:P com from Control and NPPp from GW. The global C:P com ratio thus calculated is 148:1 for the year 2100. It is higher by 2 than 146:1 from the Control run (table S2). Note that these two global ratios were calculated using the same local ratios shown in figure 1(c) but had different NPPp weighting due to changes in the regional production under global warming. At face value, the result means that the regional production changes simulated by GW would change the global C:P com ratio by +2.
Finally, the global export C:N:P ratio increases from 113:16:1 in 1750 to 119:17:1 in 2100 in experiment GW (table S2). Kwiatkowski et al (2018) report a similar increase in the global export C:P ratio (4.3%) but a much smaller change for C:N ratio (0.4%). In GW, as a result of organic matter becoming more Crich, C export declines by only 10% while P export declines by 14%. The stoichiometric buffer effect on C export (BE) can be quantified with the following equation: where fPOC and fPOP are the fractional changes in the export C and P production since 1750. For the experiment GW in year 2100, fPOP = 0.86, fPOC = 0.9, so BE = 29%. For the Redfield experiments, fPOP = fPOC and thus BE = 0%, indicating no buffering. At the other extreme is perfect buffering, when POC remains unchanged (i.e. fPOC = 1) for any change in POP due to compensatory change in export C:N:P.

Discussion
Phytoplankton community stoichiometry increased in our global warming experiment GW (table S2, figure S1). Of the total change in the global C:P com of +11 by the year 2100, physiology accounted for approximately +12, while changes in the community composition and regional production both accounted respectively for −2 and +2. There is nonlinearity in the system and these quantification involved imperfect sensitivity experiments, but it is evident that phytoplankton's physiological response to the environment is the most important of the three controls on the global phytoplankton stoichiometry. The main physiological drivers of C:P increase are warming and nutrient depletion. The warming effect makes cyanobacteria more C-rich everywhere especially in the lower latitudes (figure S1(c)). The nutrient depletion effect is most prominent in eukaryotes and primarily in the polar seas (figure S1(b)), where warming and reduced sea ice promote growth. Eukaryotes do not show elevated C:P in the lower latitudes, unlike the cyanobacteria and diazotrophs, because there is no temperature dependence on eukaryote ratios. Overall, the warming effect is the stronger physiological response, a response that was previously unaccounted for in our previous study (Tanioka and Matsumoto 2017). Note that the effect of temperature on metabolism (e.g. maximum growth rate or uptake rate) is separate and already accounted for in our model. It is likely that the quota model of Kwiatkowski et al (2018) accounts implicitly for the temperature effect on stoichiometry via change in growth rate, but it is not identified as an environmental driver of stoichiometry under warming per se. Experiment GW_CNP indicates that much of the rise in the global mean uptake C:P is due to physiology. Similarly, Kwiatkowski et al (2018) attributes 95% of global change in C:N:P to physiology.
In terms of global taxonomy, the model predicts that nutrient depletion shifts the phytoplankton community in favor of cyanobacteria at the expense of eukaryotes. The shift by the year 2500 is only~3% in PFT fractions (figure 3(d)), which is consistent with a similar taxonomic shift by~2% by the year 2100 in Kwiatkowski et al (2018). Qualitatively shifts have been simulated by other models and attributed to the advantages that smaller cells have in nutrient uptake (i.e. lower half saturation constants) over larger cells under nutrient limitation (Bopp et al 2005, Marinov et al 2010. However, there are other credible projections of the future phytoplankton community composition that highlight the importance of factors other than nutrients such as temperature, light, and grazing (Laufkotter et al 2015).
We quantified the effect of shifts in the regional production on global stoichiometry in this study for the first time. The net global effect is modest. However, the regional effects are quite nuanced, because they depend on whether the regional stoichiometry is higher or lower than the global mean. For example, the decrease in NPPp in the subpolar North Atlantic (just south of Greenland) tended to raise the global mean C:P com , because the regional C:P com was lower than the global mean (figure 4). At the same time, the same decrease in NPPp in the North Atlantic further south at~40 • N tended to lower the global mean C:P com , because the regional C:P com was higher than the global mean.
Another dipole pattern in the effect of regional production on the global mean C:P com is diagnosed in SO for a different reason. There, the sign of the effect is uniformly negative but ∆NPPp has a sign flip ( figure 4(b)). The retreat of Antarctic sea ice is largely responsible for the NPPp increase in the Antarctic waters and contributes modestly to the NPPp reduction in the subantarctic waters (figure 5). To the extent that the Antarctic sea ice retreat contributes to this ∆NPPp dipole pattern, it is reminiscent of the 'downstream' effects seen in iron fertilization experiments (e.g. Dutkiewicz et al 2005, Jin et al 2008.
In these experiments, nutrients become depleted and production is enhanced at the sites of iron fertilization; however, production is reduced downstream of the sites due to nutrient-depleted waters being transported from the fertilization sites. In this study, the retreat of sea ice triggers a boost in production and reduction of all nutrients in the Antarctic waters. The advection of these nutrient depleted waters northward by Ekman transport reduces production into the subantarctic.
A somewhat analogous effect is noted in a recent future projection based on a fully coupled model (Moore et al 2018). During the period 1850-2300, the coupled model simulates elevated NPP in SO (south of 60 • S), driven by warming and ice retreat, as in this study. At the same time, the model simulates a decrease in NPP north of 30 • S (24% by 2300) due to nutrient redistribution.
It is quite possible that the NPP projection by the coupled model is more credible than ours, given that coupled models generally have more realistic physics and a higher resolution than EMICs. For example, SO NPP in the coupled model of Moore et al (2018) is also affected by variation in surface winds, which drive upwelling. In contrast, surface winds are constant in MESMO. Although the coupled model was run under a more extreme future warming scenario RCP 8.5 (van Vuuren et al 2011), it is interesting to consider how regional production shifts could affect the global stoichiometry had we used the coupled model's projected ∆NPPp instead of MESMO's. According to our calculation, if ∆NPPp were −24% north of 30 • S and +60% south of 60 • S (Moore et al 2018), the global C:P com would be 140:1. When compared to 146:1 in our Control run, the regional production effect on the global C:P would thus be −6. This is substantially larger and opposite in sign to our estimate of +2 based on MESMO.
A lesson from this simple calculation is that regional production changes could still be quantitatively important to global ocean stoichiometry in the future. However, there are a couple important caveats to this lesson. First, an accurate model simulation of Antarctic sea ice is currently one of the most challenging tasks of climate modeling. In its Fifth Assessment Report and special report on the cryosphere, the Intergovernmental Panel on Climate Change gives 'low confidence' to Antarctic sea ice projections (Collins et al 2013, IPCC 2019. Second, limited data suggest that melting sea ice serves as an important source of Fe to the surface ocean, and depending on whether Antarctic waters are more strongly light-limited or Fe-limited, it is difficult to accurately predict the response of polar Fe cycling and production to future sea ice retreat (Lannuzel et al 2016). These caveats are highlighted by an intercomparison of the 21st century projections of primary production from nine global models (Laufkotter et al 2015). The intercomparison shows that all models projected increased NPP in SO, while disagreeing on the future of Antarctic sea ice. The same study also shows substantial variability in NPP projections of lower latitudes.
Finally, BE calculated for experiment GW with equation (4) is~30% for the year 2100. It is interesting to note that the stoichiometric buffering indicated by the divergence of global C and P export (figure 3(b)) is noticeable even prior to the year 2000 even while warming is still rather slow.
One could also consider the stoichiometric buffering effect in terms of anthropogenic C uptake. Kwiatkowski et al (2018) estimate this effect to be 0.5%-3.5%, indicating a small change in C uptake due to flexible stoichiometry. The range depends on their assumptions about whether the future ocean remains largely N-limited or becomes more P-limited. Likewise, the oceanic uptake of C is larger by just 1% in experiment GW with flexible stoichiometry compared to experiment GW_Redfield. These results point to the greater importance of the air-sea pCO 2 gradient created by the massive emissions of fossil fuel CO 2 into the atmosphere as a driver of C uptake in the future.

Summary
This study examined three controls on global ocean phytoplankton stoichiometry under global warming with a focus on the last: physiology, taxonomy, and regional production. Our study predicts that marine organic matter becomes more carbon-rich, so that even under increased ocean stratification due to future warming, global carbon export production is strongly buffered from decreasing. We find that the physiological response of marine phytoplankton to warming and nutrient depletion is the primary driver of elevating the global ocean C:N:P ratio. Changes in the community composition are only modestly important in comparison. We provided a quantitative analysis of the effect of the regional production shifts on the global ocean stoichiometry for the first time and find that it is quite nuanced in that the retreat of Antarctic sea ice can trigger offsetting changes in production downstream of the Antarctic surface waters. However, a different pattern of regional production changes triggered by retreating Antarctic sea ice, as projected by a coupled model from a different study, suggests that regional production changes could still be an important driver of global ocean stoichiometry.