Global nutrient cycling by commercially-targeted marine ﬁsh

. Throughout the course of their lives ﬁsh ingest food containing essential elements, including nitrogen (N), phosphorus (P) and iron (Fe). Some of these elements are retained in the ﬁsh body to build new biomass, which acts as a stored reservoir of nutrients, while the rest is excreted or egested, providing a recycling ﬂux to water. Fishing activity has modiﬁed the ﬁsh biomass distribution worldwide and consequently may have altered ﬁsh-mediated nutrient cycling, but this possibility remains largely unassessed, mainly due to the difﬁculty of estimating global ﬁsh biomass and metabolic rates. Here we quantify the 5 role of commercially-targeted marine ﬁsh between 10g and 100kg ( CTF 100 kg 10 g ) in the cycling of N, P and Fe in the global ocean, and its change due to ﬁshing activity, by using a global size-spectrum model of marine ﬁsh populations calibrated to observations of ﬁsh catches. Our results show that the amount of nutrients potentially stored in the global pristine CTF 100 kg 10 g biomass is generally small compared to the ambient surface nutrient concentrations but might be signiﬁcant in the nutrient-poor regions of the world: the North Atlantic for P, the oligotrophic gyres for N and the High Nutrient Low Chlorophyll (HNLC) 10 regions for Fe. Similarly, the rate of nutrient removal from the ocean through ﬁshing is globally small compared to the inputs, but can be important locally especially for Fe in the equatorial Paciﬁc and along the western margin of South America and Africa. We also estimate that the cycling rate of elements through CTF 100 kg 10 g biomass was on the order of 3% of the primary productivity demand for N, P and Fe globally, prior to industrial ﬁshing. The corresponding export of nutrients by egestion of fecal matter by CTF 100 kg 10 g was 2.3% (N), 3.0% (P) and 1-22% (Fe) of the total particulate export ﬂux, and was generally more 15 signiﬁcant in the low-export oligotrophic tropical gyres. Our study supports a signiﬁcant, direct role of the CTF 100 kg 10 g fraction of the ichthyosphere in global nutrient cycling, most notably for Fe, which has been substantially modiﬁed by industrial ﬁshing. Although we were not able to estimate the roles of smaller species such as mesopelagic ﬁsh because of the sparsity of observational data, ﬁshing is also likely to have altered their biomass signiﬁcantly through trophic cascades, with impacts on biogeochemical cycling that could be of comparable magnitude to the changes we assess here.


Introduction
Nutrient elements are used by organisms to construct the molecules they need to grow and metabolise, but are often scarce, limiting growth rates (Moore et al., 2013;Sterner and Elser, 2002). Plankton and bacteria dominate the cycling of nutrients from literature to quantify the amount of Fe in the global fish biomass, and the amount of Fe cycled by this biomass each year.
However, the total inventories and cycling rates have remained unquantified due to the lack of reliable global fish biomass and metabolism estimates.
Recently Bianchi et al. (2021) used a dynamical global model of marine commercially-targeted fish (CTF), constrained by observed fish catches, to estimate the total CTF biomass and cycling rates, and their distribution in the world's oceans. Our 70 study builds on Bianchi et al. (2021) by using an updated version of the model to estimate the role of CTF in the global cycling of the three most important growth-limiting nutrients: N, P and Fe (Fig. 1). Section 2 lays out the method we use to investigate the nutrient dynamics for the total CTF biomass for body sizes between 10 g and 100 kg, hereafter CT F 100kg 10g , using model simulations for both a reconstructed state prior to industrial fishing ("pristine") and a simulated global peak catch.
We then present and discuss the results as a series of relatively self-contained sections by topic. Section 3 details the total 75 nutrient content in CT F 100kg 10g , Section 4 focuses on nutrient cycling rates by CT F 100kg 10g , Section 5 discusses an extension of the CT F 100kg 10g estimates to all fish from 1g to 1000kg, and Section 6 presents the impact of fishing on CT F 100kg 10g (Fig. 1).
Section 7 concludes the paper. 80 We used an ensemble of simulations from the BiOeconomic mArine Trophic Size-spectrum (BOATS) model (Carozza et al., 2017), also used by Bianchi et al. (2021).

The model
BOATS provides a global, size-based numerical simulation of commercially-targeted marine organisms (including molluscs and crustaceans, here collectively termed "fish") by coupling an ecological module with a fishery economics module (Carozza 85 et al., 2016(Carozza 85 et al., , 2017. The ecological model is based on processes derived from macro-ecological theory (Carozza et al., 2016), parameterized through a Bayesian Monte Carlo approach that optimizes model performance by minimizing the discrepancy between observed and simulated catch and biomass for globally-distributed LMEs (Carozza et al., 2017). Modeled fish are divided into 3 "super-species" groups defined by the asymptotic mass of fish: small (0.3 kg), medium (8.5 kg) and large (100

Nutrient cycling
Catch N, P, Fe N, P, Fe N, P, Fe N, P, Fe N, P, Fe N, P, Fe N, P, Fe N, P, Fe N, P, Fe N, P, Fe N, P, Fe N, P, F e N, P, F e N, P, F e N, P, F e N , P ,  through the global fish biomass (solid arrow), the nutrients contained within the biomass, and the removal of nutrients from the ocean through fishing (dashed arrow). fin and fout are the fluxes in and out of the global fish biomass, respectively. E indicates excretion (release of dissolved compounds) and F indicates egestion (release of solid feces). kg). Individual sizes are binned logarithmically into 30 mass classes ranging from 10 g to 100 kg, with all three size groups 90 starting at 10 g (Carozza et al., 2016). The three groups are not intended to represent the entire marine ecosystem, but rather the sum of all species that have been commercially harvested (and are therefore accounted for in harvest records and stock assessments, which are used to constrain the model). The underlying philosophy of the model is that, although these very diverse species differ widely in their biological strategies, all are competing for food energy ultimately provided by the fixation of organic carbon through photosynthesis (which has been shown to limit fish harvests (Chassot et al., 2010;Stock et al., 95 2017)), while inhabiting the same environment, which therefore makes them subject to the same basic ecological constraints.
The constraints we apply in the model are the impacts of water temperature on growth, mortality, and phytoplankton size, and net primary production. Although this biologically 'coarse-grained' approach precludes resolution of species-level dynamics, it is solidly-founded in bioenergetic principles, and is well-suited to providing a global view of the entire ecosystem on long timescales, given that it is likely to be relatively robust under any changes in the distribution, abundance or evolution of 100 commercial stocks under historical fishing pressure (Guiet et al., 2020). Figure 2 provides a schematic overview of the model structure.
The BOATS model was deliberately developed to represent marine organisms over the size range most heavily targeted by fisheries, since this is the range for which fisheries data can offer useful constraints on the ecosystem function (Carozza et al., 2016(Carozza et al., , 2017. The starting point of 10 g coincides roughly with the weight of mature anchovy (approximately 11 cm in length 105 Net primary production Water T°r ecruitment harvest growth fishing effort Physiological limit: maximum possible somatic growth rate  according to Pauly and Tsukayama (1984)), while above 100 kg the growing significance of mammals in the ocean makes a strictly ectothermic model less capable of capturing the full marine animal community (Hatton et al., 2021).
Fishing effort and catch are computed assuming open-access dynamics and based on the Gordon-Schaefer fishery economics model (Gordon, 1954;Schaefer, 1954) with increasing catchability over time due to technological progress (Galbraith et al., 2017).

110
The model represents fish on a two-dimensional grid, i.e. longitude and latitude, which is divided in regular 1°x 1°grid cells. Thus, the model does not resolve the vertical dimension, but sums all ecosystem productivity and biomass within the water column at each horizontal point. Given that the model does not resolve interactions in space between individuals, this reduction in dimensionality does not -on its own -introduce any bias. The model is forced by monthly climatologies of observed net primary production and surface ocean temperature (Dunne et al., 2007).
historical catch maxima across LMEs as reconstructed by the Sea Around Us Project (SAUP) (Pauly and Zeller, 2016), while simultaneously falling within the acceptable bounds for the catch:biomass ratio as constrained by stock estimates (Ricard et al.,125 2012). We then perform simulations with our Fe-limited verision of the model using these 31 model parameter combinations.
Each simulation includes 200 years without catch to estimate pristine biomass at equilibrium, followed by 220 years with fishing driven by the only increase of the catchability of biomass at 7%.y −1 to reproduce the historical progression of the global fishery (Galbraith et al., 2017). The simulations slighlty underestimate the observed LME catch peak of ∼ 110M t.y −1 and the variation between LME peaks corresponds to observation with a squared Pearson product-moment correlation coefficient 130 r 2 = 0.42, p − value < 10 −9 , when averaged across all ensemble members (Supp. Figs. A1 and A2). More details on the observational constraints and on the Monte Carlo approach used here can be found in Bianchi et al. (2021). From the 31 simulations of the ensemble, we analyze the global biomass and cycling rates at pristine state and at the time of the global peak catch.

135
To estimate the quantity of a nutrient X stored in fish biomass, we convert our modeled biomass estimates from wet weight to carbon (C) weight and multiply by an average mass ratio of nutrient to carbon, X:C. Prior work has suggested that body nutrient concentration of fish may be affected by several factors such as body size, ontogeny, species, sex, diet, temperature or water nutrient concentration (e.g., Halvorson and Small, 2016;Prabhu et al., 2016;Allgeier et al., 2017), with species appearing to be the most important factor (Allgeier et al., 2020). Among these factors, our model could potentially account for change during 140 ontogeny as organisms grow in size, but analysis of the available data (see supplement) shows little to no systematic variation of specific nutrient content with size. We thus assumed constant nutrient proportions throughout food webs.
Although the model implicitly includes molluscs and crustaceans, they represent only a small proportion of the commercial catch between 10g and 100kg (from SAU data invertebrates comprised about <14% of total catch). Additionally, the measured nutrient content of molluscs and crustaceans falls within the uncertainty range around the mean value for fish (Supp. 145 Tables A2-A3). As a result, we did not attempt to account for invertebrates separately, but applied the fish nutrient ratios to all CT F 100kg 10g biomass. The Fe content of fish is poorly constrained (few whole body measurements), so we rather use the 95% confidence-interval from available data, which ranges between 10 and 200 µmolFe/molC (Galbraith et al., 2019).

Nutrient cycling by fish 150
First, we define the terms we use by considering the fate of a nutrient element when ingested by an individual fish (Fig. 3a). A fish i ingests a mass flux of a nutrient, I i (e.g. gN.d −1 ), of which a fraction α (between 0 and 1) is absorbed across the gut to produce an absorption flux, A i , while the remainder is egested in the form of feces, F i = (1−α)I i (Fig.3). Note that absorption efficiency here is defined as the difference between ingestion and egestion. Published estimates of absorption efficiencies of fish are listed in Table 1. A fraction of the element absorbed across the fish gut is then used for growth and reproduction,  of the production of new biomass to the ingestion rate, and is equal to γα. The remainder is used for maintenance and excreted back to the water; E i = (1 − γ)αI i (Fig. 3a,b). We define the cycling rate of a nutrient through the fish as which can be thought of as the biologically-processed outputs of the fish.
Based on these terms, we follow Bianchi et al. (2021) to estimate the total flux of carbon through the CT F 100kg 10g based on the 160 community average τ and growth rates G i (C) of all simulated fish within the size spectrum. We calculate the ingestion for each fish within the size bins by dividing G i (C) by τ (C) (which is 0.17 when averaged across the ensemble). We then subtract the growth G i (C), to arrive at the carbon cycling rate, and sum over all fish (Fig. 3b). We note that, at steady state for the community level there is no net growth since new production is balanced by predation, but by subtracting G i (C) we ensure that there is no double-counting of internal cycling within the simulated size spectrum through predation on simulated fish.

165
Because the model mean predator-prey mass ratio for the ensemble is 0.6 10 4 , there is very little predation within the resolved size spectrum (there are only four orders of magnitude across the smallest to the largest size bin), so that the subtraction of G i (C) could lead to a small systematic underestimate in our cycling estimates. However, we also consider that predation by large predators not resolved by the model, including very large fish, marine mammals or birds, will cause a portion of the resolved fish growth to be cycled by other organisms outside of the CT F 100kg 10g . Any consequent underestimate is therefore 170 likely to be significantly less than the value of τ (C) (i.e. significantly less than 17%).
Resolved size-spectrum 10g 100kg Predation from outside the spectrum Recruitment We can then use the input rate of carbon to the fish population, f in (C), to calculate the cycling terms for any nutrient element based on the stoichiometry of the ingested organic matter and element-specific absorption. We do so by taking the nutrient X 175 to carbon mass ratio within the prey, (X : C) prey , and the absorption efficiency of nutrient X, α X , to compute: where f in (X) is the ingestion flux of nutrient X into all fish, and E(X) and F(X) are the excreted and egested fluxes of nutrient X, respectively, for all fish.
Note that we separate the ingested fraction, I, and egested fraction, F, using a constant mean absorption coefficient for each 180 nutrient, α X (Table 1).
As indicated by eq.2, the cycling rate we calculate for an element X will be proportional to the average value of (X : C) prey .
To maintain tractability, we assumed that zooplankton provide a representative indication of the stoichiometry of the fish food source and we used mean zooplankton N and P nutrient content to compute the cycling of these elements through the fish biomass (Table 1). For zooplankton Fe content, spatial variability appears to be more important than for N and P, given 185 significant differences between Fe-rich and Fe-poor regions (Table 1, Galbraith et al. (2019)). In order to be thorough, we computed Fe cycling in three different ways based on the various computation of the Fe:C distribution of zooplankton: 1) we used a mean Fe:C in zooplankton of 30.6 µmolFe/molC, 2) we used a spatial variation between Fe-rich and Fe-poor regions and used the Fe:C mean values in Fe-poor and Fe-rich regions, respectively, from Table 1, 3) we used the same spatial variation, but with the low and high Fe:C estimates of zooplankton from Table 1. For the spatial variations of Fe:C, we assumed that  conditions are encountered in HNLC regions, which are determined by a concentration of surface nitrate [N O − 3 ] larger than 5mmolN.m −3 . In order to take into account the gradient between these regions, we locally weighted zooplankton Fe content using a Michaelis-Menten function: 2.4 Primary producers demand, nutrient concentrations, export and atmospheric deposition 195 We compare the nutrient cycling by fish to the nutrient demand by primary producers, in order to provide a rough characterization of its magnitude. We use an averaged satellite-based primary productivity (PP) (Dunne et al., 2007) to compute the PP demand for N, P and Fe. We predict the C:P and C:N ratios in phytoplankton using empirical relationships with P O 3− 4 and N O − 3 surface concentrations as described in Galbraith and Martiny (2015): where nutrient concentrations are in µmol/L.
Similarly to zooplankton, the Fe:C of phytoplankton is computed by allowing variation in stoichiometric ratios between the mean value found in Fe-poor conditions and the mean value found in Fe-rich conditions, or between the high and low phytoplankton Fe:C estimates (Table 1), using a Michaelis-Menten equation analogous to equation 3.
We then use the phytoplankton nutrient ratios to compute the export of nutrients from a satellite-based estimate of total C 205 export (Dunne et al., 2007).
We use the World Ocean Atlas observed P O 3− 4 and N O − 3 water concentrations (Locarnini et al., 2010), and dissolved Fe concentrations simulated by the biogeochemical model TOPAZ2 (Dunne et al., 2013) to compare the fish biomass nutrient content to the surface ocean ambient concentrations of nutrients, and to compute the stoichiometric ratios (e.g. in eq.4). The Finally, we compare the rate of nutrients removal by fishing at the time of global peak catch to current atmospheric deposition fields of soluble N (Brahney et al., 2015) and Fe (Mahowald et al., 2009) (Supp. Fig. B1). The catch rate and its spatial distribution are simulation by the coupled economic-ecological model, and systematically differs from the actual historical 215 peak in that the model ensemble simulates higher catch rates in the open ocean than observed.

Fish biomass: a living pool of nutrients
Our results show that the nutrients contained within CT F 100kg 10g biomass represent a non-negligible proportion of the ambient dissolved concentrations in areas where these concentrations in seawater are low and/or where CT F 100kg 10g biomass is high.
The highest amounts of N, P and Fe in the pristine fish biomass are located in the most productive regions along the coasts, 220 where most simulated fish biomass occurs (Fig. 4a). Globally, the estimated pristine biomass of CT F 100kg 10g , which represents 2.5 ± 0.8 Gt of wet biomass, contains 69 ± 31 Tg of N, 15 ± 14 Tg of P and 0.012 − 0.23 Tg of Fe, of which about half is found in the Large Marine Ecosystems (LMEs) ( Table 2). Note that, in our computations, fish biomass is the only term that varies spatially since the nutrient contents of fish, (X:C), are held globally constant (see Methods section).

225
The N content of CT F 100kg 10g biomass is locally on the same order of magnitude as ambient surface in the oligotrophic gyres where nutrient concentrations are low, and in coastal shelves where large fish biomass accumulates prior to industrial fishing ( Fig. 4b, Table 2). The amount of P in CT F 100kg 10g biomass represents a high proportion of available P O 3− 4 in the North Atlantic ocean, which is relatively P-poor, exceeding 30% in some areas (Fig. 4b). The ratio also exceeds 20% in the western North Pacific and in a few locations such as in the Arabian sea. CT F 100kg 10g biomass 230 stores much higher Fe compared to dissolved surface concentrations in the highly productive subtropical front waters of the North Pacific, North Atlantic, and Southern Oceans, with relative values exceeding 50% for the high-end estimate (Fig. 4d).
Despite low surface Fe concentrations in the Southern Ocean, the proportion in CT F 100kg 10g is particularly low due to the low modeled fish biomass, while in the tropical Atlantic the input of Fe from dust greatly overwhelms the iron in fish (Mahowald et al., 2009;Myriokefalitakis et al., 2016).

235
In summary, the storage of nutrients by CT F 100kg 10g biomass can be quite significant, compared to the dissolved nutrient inventory of the euphotic zone, but only where ambient dissolved nutrient concentrations are low and/or fish biomass is high. In these areas, CT F 100kg 10g biomass could be expected to have a greater potential as a source (if stored nutrients are made available) or a sink (if the nutrients cannot be used by primary producers) of nutrients.

Comparison to previous estimates 240
Our estimates for the N and Fe contents of the ichthyosphere differ from previous studies to some degree, which can be explained by differences in the estimates of global fish biomass and/or uncertainty regarding fish nutrient contents (we were not able to find prior estimates for P).
The amount of N stored within the global fish biomass has been previously estimated to be about 23 Tg (Allgeier et al., 2017)   which is about 66% less than our global pristine estimate of 68.7 ± 30.5 TgN (Table 2). Their computation is based on an 245 estimation of fish biomass of 0.9 Gt (Jennings et al., 2008) while our ensemble pristine CT F 100kg 10g biomass is 2.5 ± 0.8 Gt (a smaller biomass estimate than that of Bianchi et al. (2021) due to Fe-limitation of fish), a difference of biomass of 64%.
Additionally, we used a mean N content of 2.8%, slightly higher than their value of 2.6%, because we used measurements made only on wild fish and did not try to account for all the catch diversity in organisms. To reflect this uncertainty, we indicate the species-related uncertainty around the mean value we use, ±0.4%, which covers the value used in Maranger et al. (2008), 250 for our calculations in Tables 2 and 4. For Fe, Moreno and Haffa (2014) estimated that the global fish biomass stored between 0.07 and 0.7 Tg, which is a roughly threefold higher than our range of 0.012-0.23 Tg of Fe (Table 1). To compute these values, Moreno and Haffa (2014) used an estimated fish biomass of 0.9-2 Gt, which is lower than our modeled estimate of 2.5 Gt for commercially-targeted fish only due to a conservative maximum estimation (Wilson et al., 2009). However, they used a range of Fe:C values, 0.073-0.324 gFe per kg of wet weight (ww) for ray-finned fish, that is about 3-12 times larger than our 10-200 µmolFe/molC range, equivalent to 0.006-0.12 gFe/kgww (assuming 12.5% C in ww). We are more confident in our compilation of Fe:C values, which is updated, more complete, and only based on peer-reviewed studies (Galbraith et al., 2019). Nonetheless, the differences highlighted here, and the large range of estimates (low-high values) highlight the large uncertainty on the Fe content of whole fish and the great need for more measurements in this domain.

260
Note that our modeled estimates are likely to be more reliable in LMEs since fish biomass in these regions is better constrained by fish catch data.

Nutrient content variations in fish
Many factors contribute to variations of fish body nutrient concentrations, which we cannot model explicitly but instead capture within our uncertainty estimates. Among the different influencing factors, fish species is important, for instance bony fish 265 contain larger quantities of P compared to other fishes (El-Sabaawi et al., 2016). Species can also vary in terms of the size of storage components (e.g., Czamanski et al., 2011) and the number and size of bones if vertebrates, which generally increases with adult size (Sterner and Elser, 2002). Variations in body nutrient contents can also be explained by sex and life stage, e.g.
decreased Fe body content occurs in female rainbow trout during sexual maturation (Shearer, 1984), though these would tend to average out at the population level. Studies have also shown that ontogeny may affect nutrient content, for example juveniles 270 have less P than adults (Pilati and Vanni, 2007) while Fe content varies throughout the life cycle of salmon (Shearer et al., 1994).
Because our model uses size to differentiate between fish, we analyzed aquatic animal body nutrient and body size data from Vanni et al. (2017) for any existing relationship between size and nutrient content. We found no relationship between body N content and size, and only a weak but significant relationship between body P and size (Supp. Fig. C1). In this data set, 275 the changes of P content with size seem to be more related to the difference between vertebrates and invertebrates, and to be significant for benthic organisms more than pelagic organisms. A recent study as indeed showed that the taxonomic identity is prominent in driving nutrient content variations compared to size (Allgeier et al., 2020). In addition, Hjerne and Hansson (2002) found no significant changes in the N and P content of fish with species (sprat and herring), fish size, seasons or different areas of the Baltic Sea, as did Griffiths (2006) for the P content of different fish species in lakes. The scant available data for 280 Fe, does not allow us to draw conclusions on the variations in Fe content with size.
Although nutrient ratios show many fascinating variations between species, these are small relative to variations of fish biomass density in the ocean. Thus, most of the spatial variations of the nutrient content of fish are likely to be due to variations in fish abundance rather species assemblages, and the uncertainty related to nutrient ratios is likely to be small compared to the uncertainty on fish biomass.  (Table   3). Like nutrient storage, modeled cycling by fish is larger where the biomass is higher (Fig. 5 and Supp. Fig. D1). The three different ways of computing Fe cycling by the commercial fish biomass show similar spatial patterns, but Fe cycling is larger 295 when using the weighted spatial variation between the low and high Fe:C estimates in zooplankton (Fig. 5). The spatially weighted computations suggest the possibility that Fe cycling by CT F 100kg 10g might be reduced in HNLC regions (principally in the Southern Ocean and the subarctic Pacific Ocean).

Nutrient cycling by commercial fish relative to primary production
The modeled N cycling by pristine CT F 100kg 10g biomass contributes on average to 2.2 ± 1.2% of the N demand of primary pro-300 ducers (PP) in LMEs (Table 3), and generally accounts for less than 5% of the demand, except in some coastal areas where it can be as high as 14% (Fig. 6a). Similarly, the modeled pristine P cycling by CT F 100kg 10g represents 1.2% of the total P demand in the LMEs (Table 3), less than 4% of the global P demand with a larger contribution in the North and equatorial Atlantic coastal regions, and larger than 6% contributions in some coastal areas (Fig. 6b). The high-end estimate of the Fe cycling by relative to PP demand for Fe is slightly more significant than the ratios for N and P, as it represents up to 13% of which is about 0.6% of our global pristine estimate (0.8% compared to global cycling at peak catch), or 1.2% of the LMEs pristine estimate (2.1% of our N cycling by CT F 100kg 10g biomass at global peak catch) (Table 3)). This is consistent with the cycling could provide for primary producers. Given that the biomass distribution of all marine organisms has a slope of 0 (Hatton et al., 2021),(Supp. Fig. E1), and assuming that 315 zooplankton mass ranges over 12 orders of magnitude, zooplankton biomass would be expected to be roughly 3 times the biomass of the total fish biomass in the 10 g to 100 kg size range. In addition, due to their higher metabolic and ingestion rates (Maldonado et al., 2016;Griffiths, 2006), zooplankton cycling rates are likely to be much higher than fish cycling rates, also exemplified by our modeled cycling size spectrum which has a negative slope (Supp. Fig. E2). McIntyre et al. (2008) showed that fish excretion can be important in supplying N and P to primary producers when conditions of high fish biomass and high 320 PP demand or low ambient nutrient concentrations are combined, and when nutrient inputs from anthropogenic sources are low. Our results indeed suggest an increased contribution of fish cycling when these conditions are combined. However, the strength of this contribution also depends on the timing between the release of nutrients by fish and the primary producers demand for these nutrients. For Fe, CT F 100kg 10g cycling represents a more important fraction of the PP demand compared to N and P, but large uncertainties remain in its computation. At the global scale, Moreno and Haffa (2014) estimated that the amount of Fe excreted by the commercial marine fish biomass ranged between 0.4-1.5 TgFe per year. Our modeled estimate range of 0.12-0.77 TgFe/yr, based on the high and low estimates ( Fig. 6e and Supp. Fig.D2), is lower but overlapping. The difference can in part be attributed to the lower Fe:C values we used for zooplankton (Table 1), and highlights the uncertainty on the Fe cycling 330 computation (Fig. 5). In the Southern Ocean, whales have been shown to contribute to a maximum of 0.2-0.3% of the phytoplankton demand for Fe in a pre-whaling ecosystem and no more than 0.03-0.04% in a post-whaling ecosystem, making their contribution negligible compared to that of zooplankton (>70% for microzooplankton only) (Maldonado et al., 2016).
With a modeled contribution of about 0.05-0.5% of the phytoplankton demand for Fe in the Southern Ocean, modeled pristine CT F 100kg 10g cycling coherently might be able to sustain a larger part of primary productivity than the current whale population, 335 but still far less than zooplankton as discussed before for N cycling.
Finally, fish movements also allow the transport of nutrients laterally in the ocean, constituting a sink of nutrients where the fish forage and a source of nutrients where the fish excrete, egest or die (Vanni et al., 2013;Francis and Côté, 2018), an effect we do not explicitly include here.

Nutrient export by feces
Our results show that fish fecal material has the potential to affect the distribution of nutrients within the water column, especially in regions of low export intensity. Egested nutrients are integrated into fecal pellets that sink out of the surface layer and are recycled at greater depths than if bound to smaller particles (Wotton and Malmqvist, 2001;Turner, 2015), especially fish fecal pellets that can sink faster and deeper than marine snow and phytodetritus (Saba and Steinberg, 2012). Figure 6c-d,f 345 quantifies how much CT F 100kg 10g egestion may contribute to the export of N, P and Fe to depth using the mean absorption efficiencies in Table 1 and assuming exported nutrient ratios, without fish, are on average equal to those of phytoplankton. For all the nutrients, CT F 100kg 10g -mediated export accounts for a larger part of the export in the warm, low export regions of the world oceans, i.e. the tropical gyres, where it can contribute up to 50% of the exported Fe for the high-end estimate (Fig. 6f), 6% of the exported N and 10% of the exported P (Fig. 6c,d). Globally, modeled pristine CT F 100kg 10g biomass egests 29.4 ± 15.9

350
Tg of N, 4.5 ± 2.3 Tg of P and 0.009-0.59 Tg of Fe each year, which on average roughly accounts for 2.3 ± 1.2%, 3.0 ± 1.5% and 1.1-22% of the export of N, P and Fe, respectively, out of the euphotic zone (Table 3).
These results are in agreement with Davison et al. (2013) who showed that the contribution of mesopelagic fish to the carbon export (via respiration, excretion, egestion and death) is higher in regions where the total export is small. However, Davison et al. (2013) also showed that locally, in the California Current, the active transport of C by mesopelagic fish alone, which we do 355 not model, represents about 15-17% of the total carbon export at depth. In their modeling study, Aumont et al. (2018) estimated that, globally, diurnal vertical migration of epipelagic organisms (all migrating fish and zooplankton) contributes to the flux of carbon to depth of about 18% of the passive flux. Fish egesting and respiring at depth transport significant amounts of carbon, and thus also transfers nutrients from the surface to deeper layers, a process that would have increased the contribution of fish to the export of nutrient if represented in the model.

360
More than the fish contribution to total export, their effect on particles may be most relevant for stoichiometric changes, especially for Fe. Indeed, since the absorption efficiency of Fe is smaller than that of C, the Fe:C in feces will be greater than in the ingested particles, and exported fecal material will have a greater Fe:C than biogenic sinking particles made of phytoplankton aggregates or dead organisms (Le Mézo and Galbraith, 2020). This is potentially important for mesopelagic organisms feeding on sinking material in light of the possible Fe limitation of marine animals (Le Mézo and Galbraith, 2020; Until now we have considered CT F 100kg 10g as represented explicitly by the BOATS model, which ignores non-commercially targeted fish as well as the small and large ends of the fish size range. We provide a rough estimate of how the total fish biomass within a more inclusive marine size-spectrum spanning 1g larvae to 10 6 g sharks (F 1000kg 1g ) would compare to our estimates in 370 two steps. First, we take the Bianchi et al. (2021) estimate that the CT F 100kg 10g is supported by roughly half of the total primary production that could be available to this size range (best estimate 58%), from which we calculate that if the non-commercial fraction were included the total biomass estimate would be about 5.2 Gt. We then consider that the expanded size range includes 6 orders of magnitude compared to the 4 orders of magnitude of the standard BOATS size range. In BOATS, the size spectrum of abundance has a slope of about -1, giving the biomass size spectrum a slope of 0 (Supp. Fig. E1) consistent with 375 the observed Sheldon spectrum (Hatton et al., 2021). Thus, by extension, we estimate that the biomass contains about 1.5 times (6 orders of magnitude versus 4 orders of magnitude) more biomass, giving a total of 7.8 Gt of wet biomass. Overall, we estimate that F 1000kg 1g exceeds CT F 100kg 10g by a factor of 3.1.
The extended size spectrum contains more biomass, but it also contains more small fish than the standard size range of 380 BOATS, a combination that would change the cycling rates of the total fish biomass. Since smaller fish are shown to have higher cycling rates than large ones (Clarke and Johnston, 1999), we expect higher cycling of nutrients by the extended size spectrum than would be implied by the larger biomass alone. In the pristine ocean modeled here by BOATS, the C cycling rates decrease with size following a slope of −0.37molC.m −2 .g −1 (Supp. Fig. E2). Based on this slope, and the biomass slope of 0, we would expect that the CT F 1000kg 1g would have cycled about 2.4 times more C in the pristine ocean compared 385 to the CT F 100kg 10g , with the 1g-10g size range cycling 40% more C than the rest of the spectrum combined. Taken all together, the total cycling by F 1000kg 1g would have been a factor of 4.8 greater than the CT F 100kg 10g rates discussed throughout the paper above. Assuming the additional organisms eat prey with similar body nutrient contents as the organisms already modeled, all of the nutrient cycling rates presented above would be increased proportionally.
6 Fish catch: anthropogenic extraction of nutrients from the ocean 390 As fishing activity represents a direct removal of nutrients from the ocean, we estimated how much nutrients were extracted at the global peak catch and how these extraction rates would compare to nutrient inputs to the ocean. Globally, we estimate that modeled fishing activity removes 5.4 ± 0.7 T gN.yr −1 , 1.2 ± 0.3 T gP.yr −1 and 0.09 − 1.8 10 10 gF e.yr −1 from the ocean at the time of global peak catch, of which a little less than 50% in LMEs (Table 4). Although our model is calibrated to agree reasonably well with observed catch maxima in LMEs over the years 1950 to 2010 (Pauly and Zeller, 2016;Bianchi et al., 395 2021), the total catch varies between ensemble members and tends to be overestimated in the open ocean, even with the Felimitation we added in this study, for reasons that remain unclear. In addition, the catch estimate we present is at the time of the global peak catch in idealized simulations, rather than being historically accurate. Consequently, the simulated average catch at global peak (196.2 ± 57 Tg wet weight) is higher than the peak catch estimated from fishery observations (130 ± 65 Tg), so our estimates may exceed actual wild capture extractions by 50%, with the most significant overestimates in the open ocean 400 (see Bianchi et al. (2021) for more details).
low high

Nitrogen
Our estimate of N extraction from the CT F 100kg 10g catch is consistent with previous computations, and the differences mostly reflects the fish biomass estimates. For example, using catch data, Maranger et al. (2008)  Our study framework allows to spatially compare extracted nutrients to nutrient inputs, which was not the case in previous 410 work. For N, even though N extraction by fishing can be significant locally compared to N deposition at the surface, N extraction is negligible compared to the other sources of N to the surface layers. Figure 8a compares the amount of N extracted by fishing to the modeled soluble atmospheric N deposition from Brahney et al. (2015). Globally, fishing removal of N is smaller than current modeled atmospheric deposition of soluble N, with the higher values, up to more than 60% of the N deposition, in the southern equatorial Pacific, along the western margins of Africa and South America and in the Arabian Sea, where catch is 415 high and deposition is low (Supp. Fig. B1). However, most of N supply to the surface ocean occurs through vertical diffusion and mixing of the upper layers (Sarmiento and Gruber, 2006), which likely accentuates the fact that N extraction by fishing is insignificant at the global scale.

Phosphorus
Similarly to N, our estimate of P extraction by fishing is coherent with previous work and shows that it is very small compared 420 to inputs of P to the ocean and resupply from vertical mixing and diffusion in the water column. We estimate that the amount of P removed by fishing at the global scale amounts to 1.2 ± 0.3 TgP/yr, of which 0.6 ± 0.2 TgP/yr occur in the LMEs (Table   4). Huang et al. (2020) estimated that wild and aquaculture fisheries, including finfish, crustaceans and molluscs, represented 1.1 Tg of P in 2016, which is superficially similar to our estimate. However, our estimated catch of 196.2 ± 57 Tg of wet weight, is larger than the global amount of catch they used of 169 Tg, even though they considered aquaculture in addition to 425 wild captures. Our catch estimate at global peak catch overestimates the high seas catch as discussed at the beginning of this section and inBianchi et al. (2021). Additionally, our estimate is solely based on fish P content (Table 1), which may slightly overestimate the amount of P extracted by fishing activity since crustaceans and molluscs have lower P content than finfish (Huang et al., 2020).
Contrary to N, the modeled removal of P from harvest would largely exceeds the atmospheric deposition of soluble P as P

430
inputs to the ocean mostly occur through riverine inputs (Table 4). Consequently, catch transfers P from the ocean to land where as P supply to the ocean is mostly occurring in the coastal areas, with possible impacts on the P budget of the open ocean (Huang et al., 2020). But similarly to N, vertical diffusion and mixing of the upper layers supplies P to the surface ocean in quantities that most likely render P extraction by fishing relatively insignificant.

435
Fe extraction by fishing activity is within the range of previous estimates but large uncertainties remain due to uncertainty regarding the Fe:C of fish. Moreno and Haffa (2014) (Fig. 8b), where modeled Fe deposition is small and harvest is high (Supp. Fig. B1). Contrary to N and P, Fe has a much shorter residence time and thus is subjects to local perturbations, among which Fe extraction by fishing could be important.

Local and time-dependent nutrient budgets
Nutrient budgets are subject to perturbations in space and in time that can modify the relative strength of the nutrient extraction by fishing activity. Some local nutrient budgets have been investigated to compare the amount of nutrient extracted by fishing to the nutrient loads (e.g. Hjerne and Hansson (2002)). If we were to do similar budgets, at the global scale, assuming all P inputs come from rivers and atmospheric deposition, which represents 48.5 TgP/yr (Table 4 and  Note that fish extracted from a given area may have foraged elsewhere, especially large fish able to undertake long-distance migrations like tuna, salmon or sharks (e.g., Afonso et al., 2017;Gresh et al., 2000). Consequently, the ratios between extracted nutrients and nutrient deposition may be over-or under-estimated, thus over-or under-estimating the role of fishing as a local 455 sink of nutrients (Vanni et al., 2013). In addition, the relative timing of fishing effort along with phytoplankton growth, nutrient input seasonality and residence times may also modify the importance of fishing activity as a sink of nutrients (e.g., Francis and Côté, 2018;Vanni et al., 2006).

Reductions in nutrient cycling caused by fishing
Fishing has had a dramatic influence on nutrient cycling by CT F 100kg 10g , as it has permanently removed a large amount of 460 biomass, especially in the large size classes. In our ensemble of simulations, the cycling rates decrease by about 30% for the three elements considered at the time of the global peak catch (Table 3), due to the global reduction in fish biomass of about 60%. Since large fish are heavily targeted by fisheries, the size-spectrum of CTF appears truncated at larger size class at the time of the global peak catch compared to a pristine state (Supp. Fig. E1). This reduction of the mean community size enhances the cycling of elements as smaller animals tend to have higher metabolic rates (e.g. Schmitz et al., 2010;Vanni, 2002), so that 465 the reduction in overall cycling rate is roughly half the reduction of biomass.
Finally, as the size classes below 10g are not resolved (Fig. 2) we are not able to account for biomass changes that fishing might induces through trophic cascades, or how it reverberates up to fish through food supply (Dupont et al., submitted), which has the potential to further modify fish-mediated nutrient cycling.

470
In this study, we estimate the amount of N, P and Fe contained in and cycled by the global CT F 100kg 10g biomass, both in its pristine state and at the time of the global peak catch. The overall contribution of this commercial ichthyosphere to oceanic nutrient cycling is relatively small, but is more significant in regions of low ambient nutrient concentrations, high fish biomass and low export production. The industrial fish catch generally represents a small extraction of nutrients globally compared to external inputs, though it removes significant amount of P from the open ocean compared to external inputs (mainly riverine).

475
In general, local cycling of N and P by fish is less significant than Fe cycling by fish because N and P are resupplied glob-ally through large-scale circulation processes, while Fe cycling is more local and susceptible to perturbations through rapid scavenging for example. In addition, poor absorption of Fe by fish leads to an enrichment of the Fe content of fecal matter.
Globally, nutrient cycling by the modeled CT F 100kg 10g biomass is small compared to primary producers demand for these nutrients. The highest contributions are found close to the coasts where fish biomass and productivity demand are high. Fish 480 egestion of nutrients via faecal pellets is likely to comprise the largest fraction of the sinking flux in regions of low export production, i.e. the subtropical gyres. Fecal pellets may also significantly impact the stoichiometry of sinking particles, especially for Fe, with consequences for mesopelagic organisms.
Our study provides a first global glimpse of nutrient cycling by the ichthyosphere. Although the contribution of CT F 100kg 10g simulated by our model tends to be on the order of only a few percent of total surface nutrient budgets and fluxes, we estimate 485 that these would be a factor of 3 or more larger for the 1 g to 1000 kg range, including all fish. Furthermore, the role of fish in shaping the ecosystem processes through top-down pressure on their prey may be of similar or greater magnitude than the quantities estimates here (Frank et al., 2005;Baum and Worm, 2009;Hessen and Kaartvedt, 2014;Kavanagh and Galbraith, 2018). Fish can also be highly relevant at the local scale, for example by changing the vertical distribution of Fe in the water column as mentioned above, and as highlighted by many studies on fish in coral reefs. Unresolved factors, such 490 as fish migrations, would alter our results and the sensitivity of fish-mediated nutrient cycling to warming and deoxygenation due to climate change (e.g., Lefort et al., 2015;Lotze et al., 2019). There remains much to be learned about the role of fish in global nutrient cycling.
Appendix A: Model-data comparison The Bayesian Monte Carlo approach was not performed for the updated model with Fe limitation that we use in this study.

495
Instead, we took the 31 combinations of selected parameters from the Monte Carlo approach on the standard version of the model performed in Bianchi et al. (2021) and run 31 simulations with the updated Fe-limited version of the model. The Fe limitation decrease global fish biomass slightly (Galbraith et al., 2019) so that we underestimate the global LME peak catch of ∼ 110M t.yr −1 (Fig. A1). The comparison of modeled versus observed peak catch across LMEs gives an r 2 of 0.42 (p − value < 10 − 9) for the ensemble average, while the ensemble meber with the best fit to observations has an r 2 of 0.57 500 (p − value < 10 − 9) (Fig.A2) (Bianchi et al., 2021).     Figure C1. Body N (% of dry weight) and body P (% of dry weight), as a function of body mass (log(g)) for pelagic (purple) and benthic    (Table 1). Cycling at global peak catch (10g-100kg) small medium large All y=13.34-0.57x Figure E2. Size spectrum of modeled CT F 100kg 10g C cycling a) in the pristine state and b) at global peak catch.