An iron cycle cascade governs the response of equatorial Pacific ecosystems to climate change

Earth System Models project that global climate change will reduce ocean net primary production (NPP), upper trophic level biota biomass and potential fisheries catches in the future, especially in the eastern equatorial Pacific. However, projections from Earth System Models are undermined by poorly constrained assumptions regarding the biological cycling of iron, which is the main limiting resource for NPP over large parts of the ocean. In this study, we show that the climate change trends in NPP and the biomass of upper trophic levels are strongly affected by modifying assumptions associated with phytoplankton iron uptake. Using a suite of model experiments, we find 21st century climate change impacts on regional NPP range from −12.3% to +2.4% under a high emissions climate change scenario. This wide range arises from variations in the efficiency of iron retention in the upper ocean in the eastern equatorial Pacific across different scenarios of biological iron uptake, which affect the strength of regional iron limitation. Those scenarios where nitrogen limitation replaced iron limitation showed the largest projected NPP declines, while those where iron limitation was more resilient displayed little future change. All model scenarios have similar skill in reproducing past inter‐annual variations in regional ocean NPP, largely due to limited change in the historical period. Ultimately, projections of end of century upper trophic level biomass change are altered by 50%–80% across all plausible scenarios. Overall, we find that uncertainties in the biological iron cycle cascade through open ocean pelagic ecosystems, from plankton to fish, affecting their evolution under climate change. This highlights additional challenges to developing effective conservation and fisheries management policies under climate change.


| INTRODUC TI ON
Ocean net primary production (NPP) by marine phytoplankton is a primary catalyst for ocean ecosystem services. By removing inorganic carbon from surface waters, NPP primes the biological carbon pump, which contributes to the ability of the ocean to regulate atmospheric carbon dioxide levels (Falkowski, Barber, & Smetacek, 1998).
NPP also introduces energy into the base of the food chain and therefore plays a key role in supporting pelagic marine ecosystems and provisioning services, such as fisheries (Lotze et al., 2019). The magnitude of NPP and its variability is driven by a range of environmental factors, all sensitive to future climate change. In the presence of sufficient light, nutrient resources regulate the specific rates of NPP, while the size of the phytoplankton standing stock is also governed by losses due to grazing and mortality. The eastern equatorial Pacific region is a particular hotspot for NPP (Behrenfeld & Falkowski, 1997;Silsbe, Behrenfeld, Halsey, Milligan, & Westberry, 2016;Westberry, Behrenfeld, Siegel, & Boss, 2008), the biomass of upper trophic levels and fish catch (Watson, 2017). The high productivity of this region supports livelihoods and directly benefits coastal communities (Boyce, Lotze, Tittensor, Carozza, & Worm, 2020;FAO, 2018). The key resource regulating phytoplankton activity in this region is iron (Moore et al., 2013), a trace micronutrient that is essential to major cellular processes such as photosynthesis, respiration and assimilation of nitrate (Morel & Price, 2003). The links between the ocean iron cycle and NPP in iron-limited regions are complex , with many potential interactions with climate change (Hutchins & Boyd, 2016).
Climate change is projected to have a negative impact on NPP and pelagic ecosystems in the tropical oceans in general and in the eastern equatorial Pacific in particular (Bindoff et al., 2019).
The equatorial upwelling strength, which supplies nutrients from deeper waters, is projected to decline in magnitude and extent in experiments performed with Earth System Models using the high emissions RCP8.5 scenario (Terada, Minobe, & Deutsch, 2020). A reduction in the equatorial upwelling, coupled with greater density stratification of the upper ocean, tends to result in an end of 21st century reduction in depth integrated NPP by 4%-11% globally under climate change (Bindoff et al., 2019) or by around 20% (over 10 mol C m −2 year −1 ) in the eastern tropical Pacific, with strong inter-model disagreement . The combination of reduced NPP and ocean warming is also projected to reduce marine animal biomass and maximum catch potential in the eastern tropical Pacific significantly (Lotze et al., 2019).
The ocean iron cycle components of Earth System Models are known to display substantial disagreement between models and between models and observations (Tagliabue et al., 2016), which, as iron is the main limiting nutrient in this region, raises new questions regarding the robustness of NPP and animal biomass future projections. A key area of uncertainty that is fundamental to the response of NPP to climate change is the strength of biological iron uptake , which is known to display substantial variability (Boyd, Ellwood, Tagliabue, & Twining, 2017;Marchetti & Maldonado, 2016). State-of-the-art biogeochemical models used in climate change experiments account for variability in phytoplankton iron uptake in a simple manner (Tagliabue et al., 2016). However, recent field observations and laboratory experiments find a wide variety in the strength of biological iron uptake that results from the ability of algae to both remove and store iron Lampe et al., 2018), suggesting this important process is poorly constrained in ocean models. Moreover, unlike macronutrients like nitrogen and phosphorus, a significant amount of dissolved iron is also removed by particle aggregation and scavenging, which operates alongside biological removal to deplete upper ocean iron levels and contribute to surface ocean iron limitation in the eastern equatorial Pacific, north Pacific, north Atlantic and Southern Oceans Tagliabue et al., 2017). Earth System Models also represent resource limitation of NPP in a simple manner, diagnosing the most limiting resource in a given space and time, which then controls the growth rate.
In this work, we take one single complex ocean general circulation and biogeochemistry model-PISCES (Aumont, Ethé, Tagliabue, Bopp, & Gehlen, 2015), routinely used for IPCC-type climate change simulations, and conduct a range of parallel climate change experiments using output from the IPSL Earth System Model to address how varying assumptions regarding the strength of biological iron removal affect future projections of NPP and animal biomass under the RCP8.5 scenario. We then used the results from our ocean model to conduct additional scenarios with two upper trophic level models, APECOSM (Maury, 2010) and EcoTroph (Gascuel, Guénette, & Pauly, 2011), to quantify the impact on upper trophic levels. Typically, inter-model comparisons of NPP changes in response to climate change compare different physical and biogeochemical models, which makes it difficult to disentangle the role of differences in the physical and biogeochemical parameterizations. Here, we focus specifically on the biological iron cycle within one model with a common ocean physical configuration and identify an iron cycle cascade, whereby uncertainties in biological iron removal by phytoplankton cascade up the food chain to impact upper trophic levels and their response to climate change in the eastern equatorial Pacific.

| Biogeochemical and ecosystem models
The PISCES biogeochemical model is a core component of the IPSL coupled climate model and has been used in the past CMIP5 and ongoing CMIP6 exercises. The version of PISCES used in this study includes five nutrients (nitrate, phosphate, ammonium, silicate and iron), two phytoplankton (nanophytoplankton and diatoms), two zooplankton (micro-and meso-zooplankton), two particle size classes, dissolved inorganic carbon, dissolved organic carbon, oxygen, alkalinity, calcium carbonate, biogenic silica, calcite and represents the processes of nitrogen fixation, denitrification, calcification and ammonia oxidation, as well as a dynamic representation of iron-binding ligands . The overall biological removal in the model is a combination of a maximum uptake rate (V max ), modulated by iron availability and the degree of iron accumulation Buitenhuis & Geider, 2010), with the assumption that iron uptake is slowed as cells get closer to their prescribed maximum iron quota (Q FeCMax ). The magnitude of V max and Q FeCMax is controlled by the Bio FeMax parameter, which is assigned a value of 80 μmol Fe/mol C in the control model. V max is set by the maximum rate of carbon production multiplied by Bio FeMax to produce a maximum rate of iron uptake, while Q FeCMax is set at the value of the Bio FeMax parameter . PISCES has a relatively complicated iron cycle and also accounts for iron loss via scavenging, due to both adsorption and colloidal pumping processes, with variable iron-binding ligand concentrations. The model also includes variable recycling of dissolved iron by zooplankton that is driven by food quality (Richon, Aumont, & Tagliabue, 2020) and implicit bacterial iron cycling. The model includes sources from dust, sediments, rivers and hydrothermal vents and iron is lost from the ocean as particulate iron sinks out of the bottommost model grid cell. At present, the PISCES model is one of the better performing global ocean models (Tagliabue et al., 2016) and thus an ideal platform to assess how uncertainties around total biological iron removal, connected via the Bio FeMax parameter, affects climate change projections.
We conducted simulations with PISCES embedded into the NEMO framework in the ORCA2 configuration, which has mean horizontal resolution is approximately 2° by 2° cos(latitude) and the meridional resolution is enhanced to 0.5° at the equator. The model has 30 vertical levels that have a resolution that increases from 10 m at the surface, with 12 levels located in the first 125 m, to 500 m at depth. Each experiment was spun up using typical physical forcing  for 250 years before we began our climate change experiments. Our climate runs were forced with offline output from the IPSL-CM5A model at monthly resolution, using a preindustrial control simulation from 1801 to 2100 and an addi- APECOSM (Maury, 2010) is a global ecosystem model that mechanistically represents the 3D dynamics of three size-structured generic pelagic communities (epipelagic, mesopelagic and migratory). It integrates individual, population and community levels and includes the effects of life-history diversity using a trait-based approach (Maury & Poggiale, 2013). APECOSM scales dynamic energy budget (DEB)-based individual bioenergetics (Kooijman, 2000) up to the species and community level via a trait-based approach (Maury & Poggiale, 2013). The uptake and use of energy for individual growth, maintenance and reproduction are modelled according to the DEB theory (Kooijman, 2000). The model considers important ecological processes such as opportunistic size-structured trophic interactions and competition for food, key physiological aspects such as vision and respiration, as well as essential behaviours such as 3D habitat-based movements, schooling and swarming. In APECOCM, physical drivers from NEMO (3D fields of temperature and horizontal currents, vertical mixing) as well as biogeochemical drivers from PISCES (3D fields of primary and secondary production-small and large phytoplankton, small and large zooplankton-detritus, light and oxygen) control the biological and ecological dynamics at various levels. One of the main characteristics of APECOSM is that the particular structure of regional ecosystems is not prescribed a priori. On the contrary, the ecosystem structure emerges dynamically from the interactions between the global generic structure and set of rules of the model and the regional environmental drivers.
EcoTroph (Gascuel et al., 2011) is a simplified ecosystem model, where the structure of the ecosystem is summarized by biomass or production trophic spectra, which is the distribution of the biomass or production according to trophic levels, with all species combined.
The trophic functioning of the marine food web is represented by the flow of biomass, with biomass entering the system at trophic level 1 due to NPP. Then, the biomass flow reaching each trophic level is defined by the trophic transfer efficiency which connects the flow from a given trophic level to the next. Here, we used a parameterization (du Pontavice, Gascuel, Reygondeau, Maureaud, & Cheung, 2020), where the trophic transfer efficiency at lower trophic levels is derived from the plankton food web model COBALT (Stock, Dunne, & John, 2014), while for higher trophic levels, it is based on a semi-empirical model depending on the sea surface temperature (SST). The biomass residence time at every level of the food web depends on the species life expectancy and determines the productivity of the ecosystem. This is a key parameter in EcoTroph, linking the production of a given trophic level to its biomass, and which is also forced by SST (Gascuel, Morissette, Palomares, & Christensen, 2008).

| Evaluation of NPP and the emergent constraint from past ENSO events
Across Earth System Model ensembles, the sensitivity of NPP to El Niño-Southern Oscillation (ENSO) climate variability has been shown to be related to the projected climate impact on future NPP (Kwiatkowski et al., 2017). The tropical NPP-ENSO sensitivity of all PISCES parameterizations was calculated in the respective preindustrial control simulations, following the approach of Kwiatkowski et al. (2017). In the 249 years of preindustrial simulation for each model experiment, we diagnosed the linear regression coefficient between annual Niño 3.4 region (5°N-5°S, 120°-170°W) SST anomalies and tropical (30°N-30°S) NPP anomalies. Observationally derived tropical NPP-ENSO sensitivities were calculated using the Hadley Centre Sea Ice and Sea Surface Temperature data set (HadISST1; Rayner, 2003) and two NPP products, the Vertically Generalised Production Model (VGPM; Behrenfeld & Falkowski, 1997) and the Marra model (Marra, Ho, & Trees, 2003) over the SeaWIFS satellite record (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007). The observationally derived sensitivities range from −1.50%/°C (Marra model) to −3.16%/°C (VGPM). We chose to use the SeaWIFS satellite period for this analysis rather than the MODIS era as the SeaWIFS era encompassed a period with strong ENSO dynamics, making this the time period most suitable for this analysis. For that same reason, MODIS data are used to evaluate the modelled rates of NPP as they are not subject to strong interannual variations not reproduced in the model.

| Response of NPP to climate change
The evolution of global NPP during our experiments is not significantly altered across our model experiments and NPP reduces by around 8% from a value of around 48 Pg C/year during the period 1852-2005 ( Figure S1), similar to prior work with CMIP5 models Frölicher, Rodgers, Stock, & Cheung, 2016;Laufkötter et al., 2015). The spread in global NPP among our model experiments is much smaller than that across the CMIP5 suite ( Figure S1). Spatially, the control model produces a change in NPP averaged over 2091-2100, relative to 1986-2005, that shows general decreases in the low latitudes and increases in the polar seas ( Figure 1a; Figure S2) typical of Earth System Model projections under the RCP8.5 scenario (Bindoff et al., 2019;Bopp et al., 2013).
In the equatorial Pacific in particular, NPP is projected to show declines that exceed 5 mol C m −2 year −1 (Figure 1a), similar to previous modelling exercises .
The broad coherence between our experiments at the global scale masks substantial inter-experiment differences in the eastern equatorial Pacific. Differences in the climate change impact on NPP between our experiments can be difficult to visualize ( Figure S2), but emerge more clearly when we calculate the anomaly in the climate impact on NPP, relative to the control model signal (Figure 1; Figure S3). This shows that when overall phytoplankton iron uptake rate declines, either by lowering both uptake and storage or uptake alone, the trend in NPP due to climate change in the control model  Figure 1d). Interestingly, the general trend of the results is retained when only the maximum uptake rate (V max ) is modified and the maximum storage quota is left at the control value of 80 (dashed blue and red lines in Figure 1d and Figure S3 respectively) with NPP declines of 5.0 ± 5.2% and 0.3 ± 3.0% for V max = 40 and V max = 20 respectively. This suggests that the rate of iron uptake dominates the impact on NPP from the overall change in biological iron uptake (Bio FeMax ) with changes to the maximum iron storage potential playing little role in this low iron region.
Altering the strength of biological iron uptake leads to opposite responses in the magnitude of iron removal by scavenging over the 1986-2005 period. When the magnitude of Bio FeMax is reduced to 40 and 20 μmol/mol, relative to the control value of 80 μmol/mol, we find that total biological iron uptake in the PEQD declines as expected by 25% and 45% respectively. In parallel, the total rate of iron loss

| Drivers of NPP trends and their evaluation
The resilience of upper ocean iron limitation in the PEQD drives the response of NPP to climate across our experiments. In the con- The performance of all our experiments over the historical period remained consistent with available data sets for NPP and dissolved iron. Specifically, we examined the magnitude of NPP, its variability and available surface ocean observations of dissolved iron across the tropics to assess whether any experiments were inconsistent. Turning first to the overall magnitude of NPP in the PEQD province, we can see that while there are slight differences in total NPP across the experiments in the present-day period, they range from 3.2 to 3.6 Pg C/year ( Figure 3a) and are much smaller than the differences derived from the three available global satellite NPP algorithms (Behrenfeld & Falkowski, 1997;Silsbe et al., 2016;Westberry et al., 2008) applied to MODIS data between 2002 and 2019 ( Figure 3a). A more direct method of evaluating the model NPP dynamics is to assess whether variations in NPP have the correct sensitivity to ENSO-driven SST changes within the Niño 3.4 region, which has been used previously to assess tropical NPP trends and reduce uncertainty across climate models (Kwiatkowski et al., 2017). Analysis of remote sensing data from past strong El Niño/La Nina transitions has found an ENSO sensitivity of tropical NPP of between −1.50% and −3.16%/°C (Kwiatkowski et al., 2017). Similarly, all our model experiments exhibit a negative relationship (p < .001) between SST anomalies in the ENSO-driven Niño 3.4 region and anomalies in tropical NPP, with sensitivities ranging from −1.69%/°C (Bio FeMax = 20) to −2.80%/°C (control model). Although our experiments with reduced iron uptake are slightly less sensitive to changes in SST, differences are small and all remain well within the observational bracket and the range from other CMIP5 climate models across their preindustrial control simulations (slopes in Figure 3b; Table 1). A final means to assess the viability of our different iron uptake experiments is to compare against available iron concentration observations. However, in this iron-limited region, iron is depleted to low levels in both the observations and our model experiments and is unable to discriminate strongly between experiments ( Figure S4). Overall, the range of projected NPP changes of -12.1% to +2.4% remains plausible future projections for the PEQD province. This implies that the sensitivity of NPP to climate oscillations such as ENSO are not always good constraints on divergent future NPP projections, at least for the PEQD under different iron limitation regimes. Overall, our model experiments did not show large differences in their NPP or nutrient limitation dynamics within the PEQD until after ~2030 (Figures 1d, 2a-c and   3a), which limits our ability to constrain their likelihood further with existing observational data.

| Impact of NPP changes on upper trophic levels
To quantify the implications of the uncertainty in the impact of climate change on NPP in the PEQD province driven by the strength of biological iron uptake on upper trophic levels, we conducted a set of specific additional experiments. We chose to use two upper trophic level models that span the typical differentiation between semi-empirical and mechanistic trait-based approaches. EcoTroph is a simplified semiempirical model of the flow of biomass through the food web, from primary producers to top predators (Gascuel et al., 2011), using SST and NPP as model inputs. APECOSM (Maury, 2010) is a more complex and mechanistic model that considers size-structured trophic

| D ISCUSS I ON
The balance between iron removal by phytoplankton uptake and scavenging controls the response of NPP to climate change in the eastern equatorial Pacific. Consistent with observations, iron is fully depleted in this region across all our experiments ( Figure S4) due to the combination of biological uptake and abiotic scavenging removal. In our model, any iron removed by phytoplankton can be retained efficiently in the upper ocean by a suite of recycling processes linked to grazer and bacterial activity , consistent with the available process observations . However, any iron that is removed by scavenging is inefficiently retained in the upper ocean and sinks into the ocean interior as particulate iron. Thus, in experiments where phytoplankton uptake is lowered, scavenging is increased.
This increase in scavenging increases the removal of iron from the upper ocean and iron limitation is intensified ( Figure 3). As assumptions regarding the strength of phytoplankton uptake control the balance between biotic and abiotic iron removal, the degree of iron limitation in the historical period is altered, which then drives the ecosystem response to climate change. It is plausible that additional processes not considered here, linked to photochemical transformation of scavenged iron, altered recycling efficiencies and bacterial dynamics, may also prove to be important in governing the strength of iron limitation and its resilience to climate change. However, more process-orientated studies will be necessary to illuminate the major iron cycle fluxes in this region, which would provide critical constraints on the model scenarios.
The role played by the strength of iron limitation in the historical period drives the future impact of climate change on marine ecosystems more broadly than the eastern equatorial Pacific province. This is clearly seen when we expand our analysis to the tropical Pacific between 30S and 30N as a whole, where there is a tight relationship between the changing extent of iron limitation and the climate change impact on NPP ( Figure 5). The largest reductions in NPP are found in regions where limitation by iron lessens to be replaced by increasing nitrogen limitation. Our results suggest that if iron remains the most deficient resource, then the sensitivity to changes in nitrogen supply linked to modifications to upwelling and stratification will be lessened. This points to the need to consider how the iron supplied to the ocean by other means, including continental margins, rivers and anthropogenic deposition linked to fossil fuel burning or fire, may alter in the future alongside physical mixing pathways to regulate regional NPP changes in response to climate change. Overall, our results highlight that shifts in resource limitation and its sensitivity to a changing climate have significant impacts on projections of how climate change impacts NPP and upper trophic levels.
Greater insight into the mechanisms of iron uptake is being provided by a combination of physiological and genomic methods.
Genomic data have shown that marine phytoplankton exploit a range of cross membrane transport pathways to maximize their ability to internalize iron (Coale et al., 2019;McQuaid et al., 2018). Limited information on iron quotas (which result from overall iron uptake) from the region is available from synchrotron x-ray fluorescence methods that specifically target phytoplankton cellular quotas, which find values ranging between around 10 and 40 μmol to the cellular iron demand (Marchetti & Maldonado, 2016). As the maximum uptake rate itself is a function of the number of transporters divided by the handling time (Aksnes & Egge, 1991), this opens up the possibility to use field or laboratory data on variations in iron transporter abundance and other biochemical responses (Caputi et al., 2019;Nunn et al., 2013) to constrain maximum uptake rates in the future, similar to work for phosphorus (Caceres et al., 2019).
It is highly likely that the maximum uptake rate for iron is also variable, which this study shows would drive significant variability in the climate trend in NPP and upper trophic levels at the regional scale.
There is a need for a new generation of observationally informed approaches to representing cellular resource limitation appropriate for inclusion in future Earth System Models to inform future climate change assessments of ocean ecosystem change.
We find a strong sensitivity of climate change impacts on eastern equatorial Pacific upper trophic levels to uncertainties in the responses of regional NPP. Consistent with previous assessments, this study projects negative impacts of ocean warming and changes in NPP on upper trophic levels in the PEQD province. However, we also find that the uncertainties regarding these projections are higher than previously suggested. Under the high emissions RCP8.5 scenario, total consumer biomass is projected to continuously decline over the century in the PEQD province, for both APECOSM and EcoTroph ecosystem models across all biological iron removal experiments (Figure 4).
Even in the scenario where NPP slightly increases (Bio FeMax = 20), any positive effect is offset by the adverse effect of warmer temperatures on the food web functioning to some degree (depending on the complexity of the ecosystem model used, Figure 4). This is because warmer temperatures are known to lead to amplified effects on upper trophic level biomass. Overall, we note that any change in consumer biomass does not necessarily cause a parallel change in catch. This is because the realized commercial catch also results from fisheries strategies and fisheries management, as well as modifications to the productivity of each remaining biomass unit, which may increase with warming (e.g. Guiet, Aumont, Poggiale, & Maury, 2016).
At present, this study is limited by the lack of in situ data on iron cycle processes in the region, nutrient limitation regimes and a need to be able to discriminate between the different tempera- tance is the role of shifts in resource limitation from one resource F I G U R E 5 Broader link between changing iron limitation and changing net primary production (NPP) in the tropical Pacific Ocean. The relationship between the percentage change in NPP (relative to 1986:2005) and the iron-limited area (diatoms) across the total area of the Pacific Ocean between 30S:30N for the control (black), Bio FeMax = 40 μmol/ mol (blue), Bio FeMax = 20 μmol/mol (red), V max = 40 μmol/mol (pink) and V max = 40 μmol/mol (green) to another that control the sensitivity of NPP to climate change.
Future work that couples field data on biogeochemical rates and pools with genomic techniques that illuminate resource limitation and the organismal responses to in situ environmental gradients has the potential to be transformative when coupled to a new generation of biogeochemical models.

| CON CLUS ION
To conclude, understanding how climate change will affect NPP