Model scenarios<\/p>\n
To explore the impact of global change (fishing and climate change) on carbon export and sequestration by marine macrofauna of commercial interest, we force BOATS with two contrasting scenarios of climate change: the respect of the Paris Agreement (SSP 1\u20132.6) and the business-as-usual scenarios (SSP 5\u20138.5). Two variables were used to force the model of biomass prediction (BOATS, see the next section), the NPP and the SST. In BOATS, following macroecological theory, NPP drives the accumulation of biomass in marine macrofauna based on factors such as trophic efficiency, predator-prey mass ratio, and phytoplankton size, while temperature affects physiological rates, including growth, mortality, and reproduction24<\/a>,26<\/a>. These two variables were obtained from two ESMs from the Couple Model Intercomparison Project 6 (CMIP6), namely the Geophysical Fluid Dynamics Laboratory (GFDL) and the Institut Pierre-Simon Laplace (IPSL). They cover the historical period 1950\u20132014 and future projections from 2015 to 2100. To explore the cumulative impact of climate change and fisheries, we also ran climate change scenarios together with contrasting scenarios of fishing effort: one without exploitation and another one where exploitation evolves dynamically to match the historical harvest dynamic until 2014, and where fishing effort is kept constant in the future until 2100. This protocol for simulations is the protocol of the Fisheries and MEM Intercomparison Project (FishMIP) for the ISIMIP3b29<\/a>. These simulations were contributed, along with other MEMs, to the intercomparison of global fish biomass projections22<\/a>. To compare BOATS biomass outputs in the fishing scenario, we also used global reconstruction data based on stock assessments from Worm and Branch31<\/a>.<\/p>\n Alongside BOATS, to acknowledge structural and parameter differences between MEMs in simulating the impact of climate change on marine biomass22<\/a>, we also used simulation outputs of a set of five other MEMs (Apecosm, DBPM, EcoOcean, FEISTY, and ZooMSS), referred to as \u201cFishMIP simulations\u201d in the main text. The simulations and model characteristics are described in the global fish biomass projection intercomparison by Tittensor et al.22<\/a> (see Supplementary Table\u00a04<\/a> for model references and acknowledgments of the contributing modeling groups).<\/p>\n Description of the model BOATS<\/p>\n We used the BOATS model (Carozza et al.26<\/a>,27<\/a> to estimate the change in biomass of marine macrofauna from 10\u2009g to 100\u2009kg targeted by fisheries (hereafter called macrofauna for simplicity). The main ecological model parameters and variables are listed in the Supplementary Tables\u00a01<\/a>\u20133<\/a>. Macrofauna biomass is divided into three species groups k, the small, the medium, and the large species groups, and each group is itself divided into size classes. For each group, the life history of an individual starts in the first size class (10\u2009g) and can grow until it reaches its asymptotic mass, so 0.3\u2009kg, 8.5\u2009kg, and 100\u2009kg for the small, medium, and large species, respectively. The main model equation that estimates the variation through time (t) of macrofauna biomass per mass class m in each group k, fk(m,t) corresponds to:<\/p>\n $$\\frac{\\partial }{\\partial t} \\, {f}_{k}\\left(m,t\\right)=\t \\frac{-\\partial }{\\partial m}{\\gamma }_{S,k}\\left(m,t\\right)\\bullet {f}_{k}\\left(m,t\\right)+\\frac{{\\gamma }_{S,k}\\left(m,t\\right)\\bullet {f}_{k}\\left(m,t\\right)}{m} \\\\ \t -[{\\Lambda }_{k}\\left(m\\right)+{\\Lambda }_{F}\\left(m,t\\right)]\\bullet {f}_{k}\\left(m,t\\right)$$<\/p>\n \n (1)\n <\/p>\n The first two terms on the right-hand side represent, respectively, the somatic growth of an individual across size classes, at a rate \u03b3S, and the accumulation of macrofauna biomass within the same size class. The third and last term of the equation represents losses from natural mortality that occur at a rate \u039bk and from fishing mortality that occurs at a rate \u039bF. Both rates vary with time, with the average temperature in the upper ocean (top 75\u2009m) and\/or the vertically integrated NPP. The somatic growth rate (\u03b3S,k) corresponds to the net energy input needed for growth (somatic\u2009+\u2009reproduction) (\u03beI,k) minus the fraction of energy allocated to reproduction (\u03b3R,k):<\/p>\n $${\\gamma }_{S,k}={\\xi }_{I,k}-{\\gamma }_{R,k}$$<\/p>\n \n (2)\n <\/p>\n Here, \u03beI,k is the minimum value between the theoretical maximum temperature-dependent growth rate in ad libitum food conditions derived from the assumption of a Von Bertalanffy growth (\u03beVB,k), and a growth rate limited by NPP (\u03beP,k), so based on the amount of energy available from primary production. This means that the net growth rate cannot exceed the biologically determined maximum rate, even under high NPP:<\/p>\n $${\\xi }_{I,k}=\\min ({{\\xi }_{{VB},k},\\xi }_{P,k})$$<\/p>\n \n (3)\n <\/p>\n with<\/p>\n $${\\xi }_{{VB},k}=A\\bullet {m}^{b}-{k}_{a}\\bullet m$$<\/p>\n \n (4)\n <\/p>\n where \\(A={A}_{0}\\bullet {a}_{A}\\left(T\\right)\\), A0 being the allometric growth constant69<\/a>, modulated by aA(T) the van\u2019t Hoff-Arrhenius exponential temperature dependence of growth, \\({a}_{A}\\left(T\\right)=\\exp \\left(\\frac{{\\omega }_{a,A}}{{k}_{B}}\\bullet \\left(\\frac{1}{{T}_{{ref}}}-\\frac{1}{T}\\right)\\right)\\), with \u03c9a,A the growth activation energy of metabolism, kB the Boltzmann constant while T and Tref are the local temperature and a reference temperature, respectively. The term \\({k}_{a}\\bullet m\\) accounts for the energy dissipation through activity, it is also proportional to aA(T) (see Carozza et al.26<\/a>). When primary production is limiting:<\/p>\n $${\\xi }_{P,k}=\\frac{{\\phi }_{C}\\bullet \\pi \\bullet m}{{f}_{k}}$$<\/p>\n \n (5)\n <\/p>\n $$\\pi=\\frac{{\\Pi }_{\\psi }}{{m}_{\\psi }}\\bullet {\\left(\\frac{m}{{m}_{\\psi }}\\right)}^{\\tau -1}$$<\/p>\n \n (6)\n <\/p>\n where \u03a0\u03c8 is the vertically integrated NPP, \u03c4 is the trophic scaling exponent (i.e., how efficiently energy is transferred thought the trophic web), and m\u03c8 is the representative size of phytoplankton at which NPP takes place, estimated as the geometric mean of phytoplanktonic groups mass, \\({m}_{\\psi }={m}_{L}^{{{\\Phi }}_{L}(t)}+{m}_{S}^{{1-{\\Phi }}_{L}(t)}\\), \u03a6L being the proportion of large cells and m the mass of large and small cells (as in Dunne, et al.70<\/a>).<\/p>\n Finally, the energy allocated to reproduction (\u03b3R,k) corresponds to the fraction (\u03a6k) of net energy input (\u03beI,k) allocated to reproduction:<\/p>\n $${\\gamma }_{R,k}={\\phi }_{k}\\bullet {\\xi }_{I,k}$$<\/p>\n \n (7)\n <\/p>\n where \u03a6k is the mass-dependent fraction of input energy that is allocated to reproduction. Rewriting Eq. (2<\/a>) provides:<\/p>\n $${\\gamma }_{S,k}={(1-{\\phi }_{k})\\bullet \\xi }_{I,k}$$<\/p>\n \n (8)\n <\/p>\n and assuming that \\({\\gamma }_{R,k}\\) is proportional to mass and scales with a size-dependent rate sk(m), so that \\({\\gamma }_{R,k}={k}_{r}^{\\max }\\bullet {s}_{k}\\bullet m\\), we can estimate \u03d5k(m) as:<\/p>\n $${\\phi }_{k}=\\frac{{k}_{r}^{\\max }\\bullet {s}_{k}\\bullet m}{{\\xi }_{I,k}}{;} \\, {\\mbox{with}} \\, {\\xi }_{I,k}={\\xi }_{{VB},k}$$<\/p>\n \n (9)\n <\/p>\n Where krmax, the normalizing constant, is estimated by solving the equation \\({\\phi }_{k}({m}_{\\infty })=1\\), so all energy input goes to reproduction when m\u221e is reached:<\/p>\n $${k}_{r}^{\\max }=\\frac{{\\xi }_{{VB},k}({m}_{{\\infty},k})}{{s}_{k}({m}_{{\\infty},k})\\bullet {m}_{{\\infty},k}}$$<\/p>\n \n (10)\n <\/p>\n Thus, Eq. (8<\/a>) can be written as:<\/p>\n $${\\gamma }_{S,k}=\\left(1-\\frac{{\\xi }_{{VB},k}\\left({m}_{{\\infty},k}\\right)\\bullet {s}_{k}\\bullet m}{{s}_{k}\\left({m}_{{\\infty},k}\\right)\\bullet {m}_{{\\infty},k}\\bullet {\\xi }_{{VB},k}}\\right)\\bullet \\min ({{\\xi }_{{VB},k},\\xi }_{P,k})$$<\/p>\n \n (11)\n <\/p>\n For the mortality term in Eq. (1), \u039bF corresponds to the fishing mortality rate. It equals:<\/p>\n $${\\Lambda }_{F,k}\\left(m,t\\right)={q}_{k}\\left(t\\right)\\bullet {\\sigma }_{k}\\left(m\\right)\\bullet {E}_{k}\\left(t\\right)$$<\/p>\n \n (12)\n <\/p>\n It represents the catch of macrofauna group k at the size m, and is a function of the catchability q, a measure of the technological efficiency of fishing gear to capture macrofauna biomass, the selectivity function \u03c3(m), which varies between 0 and 1 and expresses the intrinsic ability of a given gear to target a given size class m, and the fishing effort E.<\/p>\n In the simulations without fishing, there is no fishing effort E and no fishing mortality \u039bF,k. With fishing, over the historical period (1950\u20132014), fishing effort evolves dynamically under an open-access dynamic proportional to the local net profit for the difference between fishing revenues and costs27<\/a>. Effort is held constant at the 2015 level for projections to the year 2100. By design, BOATS is tuned to reproduce key aspects of global fisheries catch, including total catch, size structure of the catch, catch-to-biomass ratio, and maximum catch potential of large marine ecosystems (LME), when forced with observed satellite-based NPP and temperature29<\/a>. Here, we use the same model configuration and parameters, but force the model with historical ESM simulations and future projections. To account for parameter uncertainty, an ensemble of 10,000 simulations at the LME level was conducted during model tuning, from which five representative parameter sets were selected (see Supplementary Table\u00a03<\/a>)27<\/a>,28<\/a>. Here, each simulation is replicated five times using the five parameter sets, resulting in a total of 20 simulations: 5 parameter sets x 2 ESMs (GFLD and IPSL) x with or without fishing. The model equations are solved numerically with a monthly time-step on a global 1\u00b0\u2009\u00d7\u20091\u00b0 grid, for 3 species groups (small, medium, and large asymptotic size), with a discretisation of the biomass size spectrum into 50 size classes. We report the mean and standard deviation of the five replicate parameter sets for each combination of ESM and fishing scenario.<\/p>\n Modeling carcass production<\/p>\n In BOATS, the natural mortality rate (\u039bk) represents all forms of natural mortality, including predation and non-predation mortality:<\/p>\n $${\\Lambda }_{{{{\\rm{k}}}}}={\\Lambda }_{{pred}}+{\\Lambda }_{{non}-{pred}}$$<\/p>\n \n (13)\n <\/p>\n In BOATS, \u039bk depends on temperature, mass, and asymptotic mass and is derived from Gislason et al.71<\/a>, where all mortality causes are nested in a general formula, used for each species group k:<\/p>\n $${\\Lambda }_{{{{\\rm{k}}}}}=\\lambda \\left(T\\right)\\bullet {m}^{-h}\\bullet {m}_{{\\infty},k}^{h+b-1}$$<\/p>\n \n (14)\n <\/p>\n where m, m\u221e, h, and b are respectively the mass, the asymptotic mass, the allometric mortality scaling exponent, and the allometric growth scaling exponent. \u03bb(T) is the temperature-dependent term estimated as:<\/p>\n $$\\lambda \\left(T\\right)=\\frac{{e}^{{\\zeta }_{1}}\\bullet {A}_{0}\\bullet {a}_{\\lambda }(T)}{3}$$<\/p>\n \n (15)\n <\/p>\n with \u03b61 a mortality constant estimated from Gislason et al.71<\/a>, A0 is the allometric growth constant from Andersen and Beyer69<\/a> and \\({a}_{{\\lambda }}\\left(T\\right)=\\exp \\left(\\frac{{\\omega }_{a,{\\lambda }}}{{k}_{B}}\\bullet \\left(\\frac{1}{{T}_{{ref}}}-\\frac{1}{T}\\right)\\right)\\) the van\u2019t Hoff-Arrhenius exponential function for mortality.<\/p>\n The dead biomass of a group k (BDead,k) therefore corresponds to the portion of the biomass in that group (fk) that dies naturally:<\/p>\n $${B}_{{Dead},k}={\\Lambda }_{{{{\\rm{k}}}}}\\bullet {f}_{k}$$<\/p>\n \n (16)\n <\/p>\n To estimate the proportion of that biomass that sinks to the deep-sea floor (i.e., the biomass not dead from predation), we estimated the contribution of \u039bnon-pred,k to the total mortality \u039bk, noted \u0394contrib,k, as:<\/p>\n $${\\Delta }_{{contrib},k}=\\frac{{\\Lambda }_{{non}-{pred},k}}{{\\Lambda }_{k}}$$<\/p>\n \n (17)\n <\/p>\n We used the formula in Rogers, et al.72<\/a> to estimate \u039bnon-pred,k, while adding a temperature dependency term a\u03bb(T):<\/p>\n $${\\Lambda }_{{non}-{pred},k}=\\left({\\mu }_{R}*{m}^{{-\\beta }_{R}}+{\\mu }_{S} * {\\left(\\frac{m}{{m}_{S}}\\right)}^{{\\beta }_{S}}\\right)\\bullet {a}_{\\lambda }\\left(T\\right)$$<\/p>\n \n (18)\n <\/p>\n Here, the first term (\\({\\mu }_{R}*{m}^{-{\\beta }_{R}}\\)) corresponds to the extrinsic mortality, so mortality mainly from starvation, which decreases exponentially with an increase in body mass (m). The second term \\(\\left({\\mu }_{S}*{\\left(\\frac{m}{{m}_{S}}\\right)}^{{\\beta }_{S}}\\right)\\) mimics the mortality from senescence, so that larger individuals have a higher probability of dying once the mass at senescence (mS) is reached. The addition of both terms leads to the well-known U-shaped mortality function.<\/p>\n Although this function is derived from a model for coral-reef fish, U-shaped age-dependent mortality curves are well-documented for several species across the tree of life, including fish73<\/a>. U-shaped age-dependent mortality functions are also used to model offshore species population dynamics such as tuna and tuna-like species (see Lehodey et al.74<\/a> with an application to carcass production in Mouillot et al.23<\/a>). Although refining the parametrization of this function with species-specific data would further improve the accuracy of our estimate, the current knowledge, as well as the structure of the BOATS model, does not allow us to use a species- or group-specific (inshore vs offshore species) formula. We believe our approach, using a general U-shaped function, provides an ecologically grounded first-order estimate of carcass production across diverse fish groups.<\/p>\n The total dead biomass that sinks to the deep-sea floor (Bsink,k) therefore equals to:<\/p>\n $${B}_{{sink},k}={\\Delta }_{{contrib},k} * {B}_{{Dead},k}$$<\/p>\n \n (19)\n <\/p>\n For the largest macrofauna only (i.e., size classes 50 of group 3) (Supplementary Figs.\u00a04<\/a>\u20136<\/a>), it happens that the contribution of the non-predation mortality to the total mortality exceeds 1. In that case, we estimated that these organisms cannot die from predation but only from senescence, by limiting the value to 1.<\/p>\n Modeling fecal pellet production<\/p>\n We estimated carbon loss by macrofauna (respiration + fecal pellets) at the individual level (Ci,k) by first computing total energy input or gross ingestion (\u03beTOT), i.e., the energy needed for reproduction, growth, plus loss19<\/a>. \u03beTOT is estimated using the net biomass production \u03beI (see Eq. (3<\/a>)), divided by the assimilation efficiency (\u03b1):<\/p>\n $${\\xi }_{{TOT}}=\\frac{{\\xi }_{{I}}}{\\alpha }$$<\/p>\n \n (20)\n <\/p>\n As carbon loss corresponds to the difference between total ingestion (\u03beTOT) and the part used for growth and reproduction (\u03beI,k), Ci,k equals to:<\/p>\n $${C}_{i,k}={\\xi }_{{TOT}}-{\\xi }_{{I,k}}=\\frac{{\\xi }_{{I,k}}}{\\alpha }-{\\xi }_{{I,k}}$$<\/p>\n \n (21)\n <\/p>\n Factorizing Eq. (21<\/a>) leads to:<\/p>\n $${C}_{i,k}=\\frac{1-\\alpha }{\\alpha }\\bullet {\\xi }_{{I,k}}$$<\/p>\n \n (22)\n <\/p>\n To scale carbon loss (respiration\u2009+\u2009fecal pellets) to the entire population (Ctotal), we integrated the individual carbon loss Ci,k over the entire size spectrum for each species group k, and summed it:<\/p>\n $${C}_{{total}}={\\sum }_{k}{C}_{k}$$<\/p>\n \n (23)\n <\/p>\n Where Ck, so the entire carbon loss from group k equals:<\/p>\n $${C}_{k}=\\int _{{m}_{0}}^{{m}_{{\\infty},k}}{C}_{i,k}\\left(m\\right)\\bullet {n}_{k}\\left(m\\right){dm}$$<\/p>\n \n (24)\n <\/p>\n Where nk(m) corresponds to the abundance of macrofauna in a given group k and size class, \\({n}_{k}\\left(m\\right)=\\frac{{f}_{k}(m,t)}{m}\\).<\/p>\n As in Bianchi et al.19<\/a>, we assume that about \u03d5E\u2009=\u200920% of the biomass cycling rate is returned to the environment as fecal pellets, the remaining 80% being respired.<\/p>\n Modeling carbon export from fecal pellets and carcasses<\/p>\n We model the export of fecal pellets and carcass at depth z (FFPz) in each cell using an exponential function75<\/a> with a temperature-dependent microbial remineralization rate (kremin) and a sinking speed of macrofauna fecal pellets or carcass (Wsink). Once reaching the bottom of the ocean, we considered that the carbon is entirely remineralized in the last layer and is not buried in the sediment.<\/p>\n $${FF}{P}_{z+1}={FF}{P}_{{{{\\rm{z}}}}}\\times {e}^{\\left(-\\frac{{k}_{{remin}}}{{W}_{{sink}}}\\times \\left((z+1)-{{{\\rm{z}}}}\\right)\\right)}$$<\/p>\n \n (25)\n <\/p>\n To take into account the impact of climate change on the microbial activity (remineralization), we used a temperature coefficient (Q10) of 2, so the microbial activity increases by a factor of 2 for each 10\u2009\u00b0C increase, so that kremin equals:<\/p>\n $${k}_{{remin}}={k}_{1}\\times {Q}_{10}^{\\frac{{T}_{2}-{T}_{1}}{10}}$$<\/p>\n \n (26)\n <\/p>\n where k1 is the initial remineralization coefficient at a reference temperature T1, and T2 is the temperature in each cell at a given depth. The parameter Wsink is set to 787\u2009m\u2009d\u22121 and 1150\u2009m\u2009d\u22121 for fecal pellets and carcass, respectively12<\/a>,76<\/a>. We assumed that all fecal pellets were produced above the euphotic zone, at a depth z of zeu\u2009=\u200975\u2009m. As the temperature of the ocean varies with depth, we divided the ocean into coarse layers, 0\u2013200\u2009m, 200\u2013500\u2009m, 500\u20131000\u2009m, 1000\u20132000\u2009m, and 2000\u2009m to the bottom, and used the temperature fields from the CMIP6 GFDL and IPSL model forcings. Since there is a strong drift in potential temperature in the two last layers, we identified and removed the trend (linear detrending) using preindustrial control simulations. Finally, we ran the export model in a three dimensional ocean of dimensions 180\u2009\u00d7\u2009360\u2009\u00d7\u200924, the last dimension corresponding to the 24 depth layers used later in the OCIM (next section).<\/p>\n Modeling carbon sequestration<\/p>\n We estimate the amount of carbon stored in the ocean interior by marine macrofauna using the OCIM30<\/a>. The spatial resolution of the OCIM is 2\u00b0 horizontally with 24 vertical depth layers. We use it here to track DIC injected (Cinj) by carbon remineralization in the interior ocean. Cinj was estimated using the result of Eq. (25<\/a>) for each depth layer as:<\/p>\n $${C}_{{in}{j}_{z}}={FF}{P}_{z}-{FF}{P}_{z+1}$$<\/p>\n \n (27)\n <\/p>\n The tracer equation used to estimate the concentration and distribution of the injected DIC (C) corresponds to:<\/p>\n $$\\frac{{dC}}{{dt}}=T * C+{C}_{{inj}}$$<\/p>\n \n (28)\n <\/p>\n Here, T corresponds to the transport matrix operator that accounts for ocean currents (i.e., advection-diffusion matrix), C is the concentration of DIC, and Cinj is the source of DIC from the remineralization of macrofauna-derived organic carbon converted in \u00b5mol\/kg\/yr. For simplicity, we assume a boundary condition of C\u2009=\u20090 at the sea-surface, which assumes instantaneous air\u2013sea equilibration62<\/a>. Given that a new steady-state resulting from changes in carbon export will be reached in several centuries (~1000 years), we ran both transient and steady-state simulations. Transient simulations illustrate the real-time impact of changing exports on carbon sequestration, while steady-state simulations assume that transient changes will be sustained indefinitely over time, illustrating the potential maximum long-term effects of climate change and fisheries on carbon sequestration by macrofauna. An initial pre-industrial baseline was found by solving Eq. (27<\/a>) at steady-state and using C_inj for 1950. Then, a transient simulation with yearly changing export was performed using the backward Euler method, starting at 1950 and ending in 2100. Finally, we found another steady-state assuming that the C_inj value for 2100 extends indefinitely into the future. We obtained the total carbon sequestered by macrofauna in the ocean interior by globally integrating C to obtain a value of carbon sequestration (Cseq) in GtC. The time of sequestration (in years) was obtained by dividing Cseq by the globally integrated carbon export (in GtC\u2009yr\u22121). Finally, we calculated the carbon sequestration origin, or regions where carbon sequestered by macrofauna in the ocean interior comes from, as the product of carbon injected at each depth and the mean sequestration time at that depth, also known as the mean first-passage time (e.g., Siegel et al.77<\/a>, Primeau78<\/a>). The OCIM was run on MATLAB.<\/p>\n Estimation of KMCAs<\/p>\n KMCAs are defined using an upper quartile method as areas where both the commercial macrofauna biomass and their contribution to carbon sequestration are high (i.e., value of biomass and contribution to carbon sequestration in the upper quartile). To quantify the impact of fishing and climate change on the spatial distribution of the KBCA, we estimated the spatial distribution of the KBCA for two decades: in the 1950s in a scenario without fishing, and the 2010\u2019 in a scenario with historical fishing. All analyses were done in the R statistical computing environment (version 4.1.0)79<\/a>. All outputs are available at Mariani et al.80<\/a>.<\/p>\n Reporting summary<\/p>\n