Potential increasing dominance of heterotrophy in the global ocean

Autotrophy is largely resource-limited in the modern ocean. Paleo evidence indicates this was not necessarily the case in warmer climates, and modern observations as well as standard metabolic theory suggest continued ocean warming could shift global ecology towards heterotrophy, thereby reducing autotrophic nutrient limitation. Such a shift would entail strong nutrient recycling in the upper ocean and high rates of net primary production (NPP), yet low carbon export to the deep ocean and sediments. We demonstrate transition towards such a state in the early 22nd century as a response to business-as-usual representative concentration pathway forcing (RCP8.5) in an intermediate complexity Earth system model in three configurations; with and without an explicit calcifier phytoplankton class and calcite ballast model. In all models nutrient regeneration in the near-surface becomes an increasingly important driver of primary production. The near-linear relationship between changes in NPP and global sea surface temperature (SST) found over the 21st century becomes exponential above a 2–4 ◦ C ?> global mean SST change. This transition to a more heterotrophic ocean agrees roughly with metabolic theory.


Introduction
The metabolic theory of ecology (Brown et al 2004) predicts that the exponential relationship existing between temperature and metabolism (Gillooly et al 2001) of the individual is scalable to ecosystems. Marine planktonic communities are no exception. It has been demonstrated that the differential temperature dependencies of phytoplankton community respiration (a heterotrophic process) and primary production (an autotrophic process), of which the dependency in respiration is stronger (Duarte and Agusti 1998, López-Urrutia et al 2006, Regaudie-de Gioux and Duarte 2012, determines the spatial and temporal metabolic balance of the global surface ocean (Duarte and Agusti 1998, Hoppe et al 2002, López-Urrutia et al 2006. Growth of phytoplankton, however, is also determined by nutrient limitation, and most phytoplankton communities in the modern ocean are resource and not temperature limited (Marañón et al 2014). In the modern surface ocean, regions of net heterotrophy (where community respiration exceeds gross primary production) are compensated for by regions of high primary production (Duarte andAgusti 1998, Hoppe et al 2002). However, whether zones such as the oligotrophic subtropical gyres are truly net heterotrophic remains an ongoing debate (e.g., Ducklow et al 2013, Williams et al 2013. The attenuation of the vertical flux of particulate organic carbon is likewise strongly influenced by temperature, which regulates remineralization (bacterial respiration) rates below the euphotic zone (Marsay et al 2015).
Regions with high rates of community respiration tend to be classified as heterotrophic, and also act as net sources of carbon to the atmosphere (Duarte and Agusti 1998). There are obvious implications for the global ocean carbon sink in a warmer world, if higher near-surface remineralization rates lead to a reduction in oceanic carbon uptake, producing a positive feedback on greenhouse warming and 'short-circuiting' of the carbon cycle (López-Urrutia et al 2006, Regaudiede Gioux and Duarte 2012, Marsay et al 2015. Indeed, isotopic reconstructions of remineralization profiles during the warmer Eocene epoch (55.5-33.7 Ma) reveal both shallower remineralization and higher metabolic rates than found in the modern ocean (Olivarez Lyle andLyle 2006, John et al 2013), and are suspected of contributing to the maintenance of higher global temperatures over the period (Olivarez Lyle and Lyle 2006). Earth system modelling of the Eocene using temperature-dependent remineralization rates by John et al (2014) recently achieved good agreement with the isotopic reconstructions, supporting the hypothesis that enhanced metabolic rates played a critical role in reducing observed carbon fluxes.
Major changes are currently underway in open ocean biogeochemistry owing to anthropogenic additions of greenhouse gases to the atmosphere. These changes are affected through three independent but synergistic phenomena; ocean warming, acidification and deoxygenation (summarized by Gruber 2011). Recent interest in the potential net biogeochemical impacts of these phenomena has resulted in several studies (e.g., Schmittner et al 2008, Steinacher et al 2010, Bopp et al 2013 examining the combined biogeochemical response in coupled Earth system models. The common biogeochemical metrics in these studies are change in net primary production (NPP), carbon export to the deep ocean, and seawater oxygen content. There is general agreement across models and modelling studies that increasing atmospheric CO 2 concentrations over the 21st century will increase sea surface temperatures (SSTs) and stratification, and decrease carbon export production (Schmittner et al 2008, Steinacher et al 2010, Bopp et al 2013, Yool et al 2013. Simulated 21st century NPP also generally decreases, although Schmittner et al (2008) found an increase using the University of Victoria Earth System Climate Model (UVic ESCM). Although model structures vary and regional differences exist, the main reasons for the decline in NPP and export production are similar and related to increasing restrictions on autotrophic primary production. Increased water column stratification in the low and middle latitudes diminishes resupply of nutrients to the euphotic zone from the deep ocean (Bopp et al 2001). Global NPP declines as nutrient limitation increases in these regions (with declining carbon export as a consequence), even though high latitudes might show increases in NPP as light and temperature limitation recedes (e.g., Schmittner et al 2008, Steinacher et al 2010, Yool et al 2013. Longer-term ocean biogeochemical trends under a high carbon emissions scenario have been speculated to follow those of the 21st century (e.g., Moore et al 2013). However, Schmittner et al (2008) demonstrated the combination of reduced stratification and warming-enhanced biological processes after the 21st century result in large increases in global NPP beyond 2100, though carbon export still declines in their study. Taucher and Oschlies (2011) also demonstrated in the UVic ESCM an increase in temperature-dependent microbial fast recycling by 2100 that increases temperature-dependent NPP, which does not occur when these processes are considered temperature-independent. A recent update of the UVic ESCM biogeochemistry model by Keller et al (2012) refined the phytoplankton growth rate and added new formulations for grazing and iron limitation. These modifications have resolved the discrepancy between the sign of 21st century NPP trends in the Schmittner et al (2008) study and other model studies, though the longer-term increases in NPP remain (described in section 3). Kvale et al (2015) recently implemented a calcifying phytoplankton functional type and prognostic particulate calcium carbonate (CaCO 3 ) in the model, which increased the ecological dynamism and improved the mechanistic realism of the carbon export processes. These modifications impact the model's behavior under anthropogenic forcing but do not eliminate the long-term trend in NPP (also described in section 3). Here we argue that the apparently robust simulated trend reversal in NPP in the UVic ESCM must be given serious consideration in light of evidence presented by paleoproxy and modern observations that phytoplankton metabolic processes have the capacity to short-circuit the marine carbon cycle. ). The UVic ESCM is a coarse-resolution (1.8°× 3.6°× 19 ocean depth layers) ocean-atmosphere-biosphere-cryosphere-geosphere model. The standard UVic ESCM version 2.9 includes a NPZD submodel with 'mixed phytoplankton' and 'diazotroph' phytoplankton functional types and one 'zooplankton' functional type (described in Schmittner et al 2008). The MIXED model is the standard UVic ESCM NPZD submodel, updated with the changes made by Keller et al (2012) and referred to as 'NOCAL' in Kvale et al (2015). The CAL model builds upon MIXED and additionally includes a 'small and calcifying phytoplankton' functional type (referred to throughout this paper simply as 'calcifiers') and a prognostic CaCO 3 tracer (Kvale et al 2015). The NOCACO3TR model only includes the additional 'small and calcifying phytoplankton' functional type but does not include the prognostic CaCO 3 tracer. Table 1 lists relevant biogeochemical model tracers and functional types for the three model configurations used here. These configurations were chosen for this study because of hypothesized nutrient-driven shifts in future community composition between 'small' (commonly assumed as a proxy for calcifiers) and 'large' (commonly assumed as a proxy for diatoms and other non-calcifiers) phytoplankton under anthropogenic climate change forcing (e.g., Cermeno et al 2008, Marinov et al 2010, 2013. Shifts in relative abundance of phytoplankton functional types have the potential to alter carbon export to the deep ocean (Bopp et al 2005), for which carbonate from calcifiers is a significant vector (Jin et al 2006). Please refer to Kvale et al (2015) for a complete description of the CAL and MIXED ecosystem and biogeochemical models.

Methods
The UVic ESCM phytoplankton biomass equations utilize a modified Eppley curve (Eppley 1972) to calculate maximum possible growth rate as a function of seawater temperature and iron availability. Growth limitation by nitrate and phosphate are calculated using half saturation constants and nutrient concentration, where calcifiers in the CAL and NOCACO3TR models are assigned nutrient affinities higher than mixed phytoplankton (Le Quéré et al 2005). Actual growth rates are determined by the most limiting factor, given the maximum possible growth rate and availability of nutrients and light. Biomass loss terms include a temperature-independent linear mortality term, grazing by zooplankton, and a loss term referred to as 'microbial fast recycling'. There are three temperature-dependent heterotrophic loss terms in the UVic ESCM: microbial fast recycling, zooplankton excretion, and detrital remineralization. The microbial fast recycling term accounts for nearsurface losses to the microbial loop and dissolved organic matter cycling. Similar processes govern detrital remineralization and both are parameterized using the Eppley curve. Zooplankton excretion is a proportion of zooplankton grazing, which is subject to prey availability and thus phytoplankton growth rate, as well as the zooplankton growth rate which is temperature dependent until capped at 20 C° (Keller et al 2012). Respiration due to metabolic maintenance is not explicitly accounted for, but is implicitly included in primary production. Redfield stoichiometry is used throughout.
Organic carbon export in all model configurations occurs through a prognostic detrital tracer made up of dead plankton that sinks at a rate that increases linearly with depth and remineralizes at a rate dependent on temperature. Any organic detritus that reaches the sediments is assumed to dissolve instantly back into the water column. In the CAL model, a portion of this organic detritus is 'protected' from remineralization by a fixed fraction of the associated particulate inorganic CaCO 3 as a rough parameterization of ballasting (Klaas and Archer 2002). This protected detritus loses its protection at the rate of CaCO 3 dissolution, which is a function of the degree of carbonate saturation. CaCO 3 is produced at a fixed ratio based on dead plankton (mixed phytoplankton and zooplankton in MIXED, calcifiers and zooplankton in CAL and NOCACO3TR). It is exported and dissolved instantaneously in MIXED and NOCACO3TR (Schmittner et al 2008, Keller et al 2012, but it is traced and assigned prognostic sinking and dissolution rates in CAL (Kvale et al 2015). In all model configurations, any CaCO 3 reaching the sediments is buried.
The three model versions were first integrated for ten thousand years in fully coupled mode using a constant pre-industrial CO 2 atmospheric concentration of 283.8 ppm to establish equilibrium. Historical CO 2 forcing, as well as historical agricultural, volcanic, sulphate aerosol and CFC emissions, and changes to land ice and solar forcing were then applied from year 1800 to 2005 using the PMIP3 (Paleoclimate Modelling Intercomparison Project Phase 3) data compilation (Machida et al 1995, Battle et al 1996, Etheridge et al 1996, Flückiger et al 1999, Ferretti et al 2005, Meure et al 2006. From year 2005 to 2400 the models were forced using increasing CO 2 concentrations and radiative forcing from all non-CO 2 greenhouse gases, fractions of the land surface devoted to agricultural uses, and the direct effect of sulphate aerosols as an alteration of the surface albedo following 'business-as-usual' RCP scenario 8.5 (RCP8.

Results and discussion
The physical response (i.e., changes in temperature, salinity, ocean circulation, sea ice) to radiative forcing is the same across all three model configurations. Atlantic meridional overturning rapidly weakens from about 22 Sv (10 6 m 3 s −1 ) before year 2000 to less than 8 Sv by year 2300 (not shown). Antarctic Bottom Water formation increases slightly after year 2000 by 1 Sv (not shown). Zonally averaged SSTs warm up to 11 C • relative to the year 1800 value by year 2400, and strong near-surface stratification occurs on a global scale (figure 1). Differences in model response to climate forcing demonstrate the effect of biogeochemical model structure only, such as presence or absence of calcifiers (compare MIXED to NOCACO3TR) and  Figure 1 shows changes in zonally integrated plankton concentration with respect to the year 1800 values. Changes in concentration are relatively small over the 21st century in all models, but show large changes after year 2100. Mixed phytoplankton decrease in the low latitudes due to increasing nutrient limitation in all models, and increase in the high latitudes due to increasing temperature and light availability. Diazotroph concentrations increase in the low latitudes, though the scale of their response is small in comparison to the other phytoplankton types. Calcifiers in the NOCACO3TR and CAL models have a lower nutrient half saturation constant (a higher nutrient affinity) than mixed phytoplankton, so decreasing surface nutrient concentrations in the low and middle latitudes increase their relative advantage, and the population increases in the southern Pacific and Atlantic basins (shown in zonal integration). Changes in zooplankton concentrations broadly follow the spatial pattern of total phytoplankton biomass.

Changes in biogeochemical pathways
Carbon and nutrient recycling fundamentally shifts after the 20th century in all model configurations. The ratio of NPP to near-surface (to 130 m depth) heterotrophic respiration declines rapidly (figure 2). While all three respiration terms increase over time (figure 2), only fast recycling increases for all three models during the 21st century. While the overall contribution of fast recycling to respiration is small (not shown), it is the most sensitive to 21st century increases in SSTs and dominates the shift towards increasing respiration rates until zooplankton excretion and total water column remineralization rapidly increase in the 22nd century. Fast recycling rates of mixed phytoplankton never decline in spite of an overall decline in NPP (figure 3) and detritus production (figure 4) in the first half of the MIXED simulation. The accelerating effect of microbial fast recycling on increasing global NPP and (decreasing the NPP): (heterotrophic respiration) ratio is most apparent in the NOCACO3TR simulation. Calcifiers replace mixed phytoplankton in the lower latitudes, which leads to more fast recycling and higher nutrient availability. More nutrients allow for higher NPP (including more calcifiers), which maintains the zooplankton population in these regions throughout the simulation. This allows zooplankton excretion rates to increase faster in NOCACO3TR than in the other models. A mitigating effect by CaCO 3 ballast is apparent in the CAL simulation when compared to NOCACO3TR, where eventual increases in all heterotrophic rates are not as large because of the protection from remineralization provided to a portion of sinking detritus. Figure 3 shows changes in globally integrated NPP from pre-industrial values, as well as the % change as a function of global mean SST for all model configurations. Our simulations show that the relationship between global mean NPP and SST is complex and not simply linear as previously discussed by Bopp et al (2013) and Moore et al (2013). In MIXED, the relationship between NPP and SST anomalies is negative and quasi-linear until 3 to 4 C • of warming, whereupon the relationship becomes positive and exponential. At global mean SST changes below 3 C • , reduced nutrient supply dominates the global NPP trend (a physically-driven feedback on autotrophy), but a transition occurs at about 4 C • whereupon Figure 1. Hovmöller diagrams of variable anomalies relative to year 1800 for the three UVic ESCM model configurations used in this study, and atmospheric CO 2 concentration timeseries.   warming-enhanced heterotrophic respiration and a small global reduction of stratification (figure 1) increases the supply of near-surface nutrients and increases NPP (Schmittner et al 2008).

Changes in NPP
The inclusion of calcifiers in the NOCACO3TR model alters the initial negative linear NPP:SST relationship in MIXED to an initial positive linear relationship. Global NPP never declines in NOCACO3TR and CAL in spite of globally declining surface nutrients because of the replacement of mixed phytoplankton populations by calcifiers.
The addition of a prognostic calcite tracer reduces model NPP sensitivity to SST warming relative to the NOCACO3TR configuration. In the CAL configuration trends in all plotted variables share the same temporal and spatial patterns, but not the magnitude, with the NOCACO3TR configuration. This is due to the added effect of calcite-ballasted detritus, which increases over the simulation alongside the calcifier population and effectively removes nutrients from the surface (thereby mitigating the heterotrophy feedback). Eventually, however, even the mitigating effect of ballasting cannot override the effect of warming-enhanced remineralization, and the NPP:SST relationship becomes exponential.
Choice of parameterization of particle export and remineralization can control surface nutrient resupply and resultant NPP trends (e.g., Taucher andOschlies 2011, Segschneider andBendtsen 2013). For instance, low latitude increases in implicit calcifiers in Moore et al (2013) were not able to offset declines in global NPP driven by the other phytoplankton because of differences in the export efficiency across the phytoplankton types. In the NOCACO3TR model there are no differences between calcifier and mixed phytoplankton carbon export efficiency, and thus calcifier replacement of mixed phytoplankton serves as an equivalent substitution. In CAL, however, calcifiers (including zooplankton) are more efficient exporters because of their ability to ballast detritus (particulate organic carbon) in protective CaCO 3 . Increased global carbon and nutrient export efficiency as phytoplankton calcifier populations expand could potentially reduce global NPP. However, our simulations show the relevance of calcifier biogeography, as regions where calcifiers expand are also the regions where warming-enhanced biological processes increase the earliest (figure 1). The effect of increasing carbon and nutrient export efficiency in CAL does therefore not offset warming-enhanced remineralization.

Changes in carbon export
Eventual increases in global NPP relative to preindustrial values in MIXED do not result in an increase in carbon export away from the surface (figure 4). Global CaCO 3 and particulate organic carbon export from the near-surface reproduce the NPP trend over the historical period and 21st century in the MIXED configuration, but while NPP and CaCO 3 flux eventually increase beyond the pre-industrial level (figures 3 and 4), particulate organic carbon export production does not, and deep water particulate organic carbon flux declines continuously over the simulation. The decline in particulate organic carbon flux in MIXED is mostly caused by the warmingenhanced detrital remineralization rate in the nearsurface, which removes particulate organic carbon (none of which is protected by CaCO 3 ) from export even as NPP and CaCO 3 export increase.
In the NOCACO3TR simulation, replacement of mixed phytoplankton by low and middle latitude calcifiers results in a more modest decline in globally integrated particulate organic carbon export production at 130 m depth, but a larger sensitivity to increasing temperature-dependent remineralization rates results in a similar decline in deep ocean particulate organic carbon flux by 2400 ( figure 4). CaCO 3 flux, on the other hand, increases more than in other model configurations because of the large NPP:SST sensitivity ( figure 3).
The moderating effect of the ballast model is apparent in the CAL results, which show the smallest total decline in global deep ocean particulate organic carbon flux and the smallest overall increase in CaCO 3 export production (figure 4). Increasing globally integrated CaCO 3 export production as a response to increasing NPP agrees with Schmittner et al (2008), and is not reproduced in models that project declining NPP (e.g., Moore et al 2013, Yool et al 2013. The addition of small phytoplankton and calcifiers to the model in the NOCACO3TR and CAL configurations decreases the early particulate organic carbon export production response (nearsurface particulate organic carbon fluxes are reduced less in these models over the first part of the simulations) by maintaining NPP in the warmest (and lower nutrient) regions. The application of a CaCO 3 tracer mitigates the longer-term response at 2 km depth because explicit CaCO 3 ballasting protects some particulate organic carbon from remineralization, so deep ocean particulate organic carbon fluxes do not decline as dramatically in the CAL model as in the other two models. Kriest et al (2010) have previously demonstrated strong equilibrium nutrient distribution sensitivity to parameterization of particle export flux and remineralization length scale, both of which are modified between the NOCA-CO3TR and CAL simulations. Addition of calcifiers (but no ballasting) to the PISCES-T model further reduced global export production and ocean carbon uptake relative to the control experiment in Manizza et al (2010) because expansion of calcifiers reduced the zooplankton population (who prefer silicifiers in their model), which reduced carbon export production and export efficiency more than if no calcifiers were included in the model. This zooplankton response is not seen in CAL because zooplankton grazing is not preferential between calcifiers and mixed phytoplankton, and reinforces the conclusion that ecosystem model parameterization choice critically determines model export production and export efficiency behavior in climate transitions (Taucher andOschlies 2011, Segschneider andBendtsen 2013).
Despite the large differences in carbon export responses between models described above, only small differences exist in the accumulated anthropogenic emissions diagnosed in the simulations. From years 1800-2400, total diagnosed CO 2 emissions are 5125 (MIXED), 5085 (NOCACO3TR), and 5082 Pg C (CAL). This result underscores that while large shifts in marine ecology may occur in the short term, these shifts must be maintained over long timescales if global carbon air-to-sea fluxes are to be substantially impacted.

Potential caveats
There are several potential shortcomings in the simulations presented here that might affect our results. Circulation change due to melting of continental ice is not simulated by the model and would likely affect the biogeochemistry. Replacement of mixed phytoplankton by small phytoplankton and calcifiers is dependent on shifts in relative competitive advantage, at least some of which is artificial to the model and would change by including more phytoplankton types. Specifically, the additions of diatoms and silicate ballasting could mitigate biogeochemical changes associated with temperature-dependent processes. The experimental set-up does not take transient adjustment in iron limitation into account, which might affect phytoplankton community composition. Effects of ocean acidification (e.g., Gao et al 2012) on photosynthesis are not included. None of the models simulate globally decreasing net calcification with increasing CO 2 concentration (Riebesell et al 2000), which would increase the long-term model sensitivity to warming-enhanced remineralization, but is also debatable (Lohbeck et al 2012, Engel et al 2014. Furthermore, using a single dissolution parameterization for zooplankton and phytoplankton CaCO 3 ignores the contribution of aragonite dissolution (Gangstø et al 2008).
It is also important to note that the UVic ESCM does not explicitly account for autotrophic respiration and so a common metric for calculating net heterotrophy, gross primary productivity (e.g., Regaudie-de Gioux and Duarte 2012), cannot be calculated. New production in the model implicitly includes the effect of autotrophic respiration, which can be expected to increase with increasing temperature.

Conclusions
The experiments described here highlight the importance of biogeochemical model parameterizations for the determination of the behavior of the ocean carbon cycle during climate transitions. Though each model analysed in this study (MIXED, NOCACO3TR, and CAL) has very similar equilibrium states, their respective biogeochemical behaviors under transient forcing are unique both in global and regional trends. Overall, the models share a common response of a transition between two dominant regimes. The first regime is physically-dominated and regulates autotrophy for lower changes in SSTs. Reductions in nutrient resupply from the deep ocean cause mixed phytoplankton populations to decline in the low and middle latitudes, while increased warming and reduced sea ice cover increase the mixed phytoplankton population in high latitudes. In the MIXED model, global NPP, particulate organic carbon and CaCO 3 export production decline over the first half of the simulation as a result of the low latitude trend. Models including an additional plankton functional type representing small phytoplankton and calcifiers (NOCACO3TR and CAL) show a partial replacement of mixed phytoplankton by calcifiers in the stratifying low latitudes due to their higher nutrient affinities. Consequently, global NPP, particulate organic carbon and CaCO 3 export production change less. All models start to transition into the second regime around the start of the 21st century, when temperature-dependent heterotrophic respiration rates increase relative to NPP.
During the second regime a biologically-driven feedback on NPP dominates the physically-driven one that defined the 21st century. It is driven by warmingenhanced heterotrophic respiration rates, which compensate for reductions in nutrient resupply by rapidly recycling carbon and nutrients in the near-surface. This feedback increases global NPP and CaCO 3 export production in all model configurations and represents a fundamental ecological shift towards increasing heterotrophy. Microbial fast recycling in the near-surface is the most sensitive to warming and accelerates the transition to this second regime. The simulated transition to this second regime in all model configurations raises the possibility trends observed in NPP over recent decades cannot be extrapolated far into the future. The net carbon export into the deep ocean during this second regime is dependent on whether the model accounts for detrital ballasting by calcite. The model that contains a ballasting parameterization (CAL) shows a mitigated increase in near-surface particulate organic carbon production, and a mitigated decline in deep ocean particulate organic carbon export because the calcite protects particulate organic carbon from near-surface remineralization and continues to effectively remove it to the deep ocean throughout the simulations. No CMIP5 model currently accounts for the ballasting effect, which suggests