Global nutrient cycling by commercially targeted marine fish

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

1 Introduction

Nutrient elements are used by organisms to construct the molecules they need to grow and metabolize but are often scarce, limiting growth rates (Moore et al., 2013; Sterner and Elser, 2002). Plankton and bacteria dominate the cycling of nutrients in the ocean (Sarmiento and Gruber, 2006), but an increasing number of studies recognize the contribution of animals to biogeochemical cycles (e.g., Wilson et al., 2009; Bianchi et al., 2013; Saba et al., 2021). Locally, it has been shown that marine animals can significantly impact the supply and storage of nutrients with consequences for primary production (Cavan et al., 2019; Roman et al., 2014) and interact with the cycling of elements through direct and indirect pathways (Vanni, 2002; Atkinson et al., 2017; Allgeier et al., 2017). For instance, Leroux and Schmitz (2015) described a theoretical framework in which animals control the flux of nutrients up the trophic chain through predation and release of waste products while also affecting the cycling of nutrients through non-consumptive effects (e.g., prey selection and stress induced in prey). In addition, since animals can swim and move in the water column, they are also able to transport nutrients from one place to another, over distances that increase with animal size (Hall et al., 2007; Vanni, 2002; Roman et al., 2014). But thus far there has been little effort to estimate how the global fish population, which we term “ichthyosphere”, influences large-scale nutrient cycling in the ocean.

During their life cycle, fish assimilate, store, and recycle essential elements that they need to build their body tissues. This storage of nutrients within fish biomass is important for human nutrition as wild-caught fish globally provide essential proteins and other micronutrients (Hicks et al., 2019). Apart from a dietary interest for humans, the reservoir of nutrients comprised by fish biomass can play a role in ecosystem function as a nutrient stockpile, the importance of which depends on how much nutrients are available otherwise (Allgeier et al., 2016). For example, the nutrients embedded in fish biomass are not directly available for primary producers, which can appear as a type of competition for resources (Hjerne and Hansson, 2002). But at the same time, the cycling of elements by fish acts as a source of nutrients to primary producers, as fish recycle elements through the excretion of dissolved bioavailable components.

The cycling of N and P by fish has often been studied in freshwater systems, but little is known for the global ocean (Schindler and Eby, 1997; Vanni, 2002; Vanni et al., 2006; Griffiths, 2006). For instance, McIntyre et al. (2008) showed that fish are able to create hotspots of recycled nutrients in streams that could meet more than 75 % of the algal and microbial requirement for N, and Allgeier et al. (2014) showed that fish are a fundamental nutrient source to coral reefs. Additionally, fish egest particulate products that can be mineralized and enhance productivity or that can sink to depth and export elements (Davison et al., 2013; Saba and Steinberg, 2012), so that they are no longer available for primary producers. As such, fish biomass can act as a bank account for nutrients, depositing nutrients when fish are feeding rapidly and withdrawing nutrients to the water column when metabolism exceeds predation. Finally, the egestion of particulate products by fish has been shown to modify the stoichiometry of biogenic particles, including dramatic changes of Fe : C, implying that egested material may also modify the relative availability of nutrients through the water column (Le Mézo and Galbraith, 2020).

The amount of nutrients stored and cycled by fish can vary with different environmental and physiological factors, in space and time (e.g., Halvorson and Small, 2016; Prabhu et al., 2016; Francis and Côté, 2018; Czamanski et al., 2011; Allgeier et al., 2014). In addition to natural variations, anthropogenic activities, mostly fishing, modify the storage and cycling of nutrients by the ichthyosphere. For instance, Layman et al. (2011) and Allgeier et al. (2016) analyzed the cycling of N and P by fish in fished and un-fished coastal sites of the Bahamas and the Caribbean, respectively. Layman et al. (2011) showed lower recycling rates of nutrients by fish in fished sites, due to biomass reduction and habitat fragmentation. Beyond fish biomass reduction, Allgeier et al. (2016) stressed the role of community size structure that, influenced by fishing, also led to reduction in nutrient storage and cycling.

Although numerous works have identified significant local effects, little is known about the contribution of the ichthyosphere to nutrient budgets at the global scale. Maranger et al. (2008) used global fish catch data to estimate the total removal of N by commercial fisheries. They integrated a spatial component in their analysis by computing N removal in 58 large marine ecosystems (LMEs) and compared it to fertilizer runoff in each off these LMEs. Moreno and Haffa (2014) took a similar approach to estimate the amount of Fe removed each year by fishing from 1950 to 2010. They also used biomass estimates 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 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$, 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 nutrient content in ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$, Sect. 4 focuses on nutrient cycling rates by ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$, Sect. 5 discusses an extension of the ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ estimates to all fish from 1 g to 1000 kg, and Sect. 6 presents the impact of fishing on ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ (Fig. 1). Section 7 concludes the paper.

Figure 1 The ichthyosphere role in global nutrient cycling. The diagram schematically indicates the cycling of nutrients (N, P, and Fe) 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). f_in and f_out 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).

2 Methods

2.1 Model description and simulations

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).

2.1.1 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 et al., 2016, 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 three “super-species” groups defined by the asymptotic mass of fish: small (0.3 kg), medium (8.5 kg), and large (100 kg). Individual sizes are binned logarithmically into 30 mass classes ranging from 10 g to 100 kg, with all three size groups 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., 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 commercial stocks under historical fishing pressure (Guiet et al., 2020). Figure 2 provides a schematic overview of the model structure.

Figure 2 Schematic overview of the BOATS model. The red, green, and black arrows indicate dependencies of model components on external forcings. The top panel indicates the energetic limits of growth as a function of fish size, while the bottom panel illustrates the three size spectra of fish groups (for simplicity only the large group is represented), their internal dynamics, and the link to economics via harvest and the interactive effort.

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 fishery data can offer useful constraints on the ecosystem function (Carozza et al., 2016, 2017). The starting point of 10 g coincides roughly with the weight of mature anchovy (approximately 11 cm in length 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).

The model represents fish on a two-dimensional grid, i.e., longitude and latitude, which is divided in regular 1^∘ × 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).

Galbraith et al. (2019) hypothesize that fish growth is reduced under iron scarcity in the wild and demonstrate that the implementation of a simple form of iron limitation of fish in BOATS improves the simulated fish catch in oceanic regions known to have low iron concentrations. We thus used a different version of the model from Bianchi et al. (2021) that includes an Fe limitation of fish growth, using surface nitrate concentrations as a proxy for iron limitation as described in Galbraith et al. (2019).

2.1.2 Simulations

We use the parameter sets selected in Bianchi et al. (2021) to run our updated Fe-limited model version. Briefly, from the Monte Carlo ensemble of 10 000 simulations, Bianchi et al. (2021) selected the parameter combinations (31 in total) that best match the 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-to-biomass ratio as constrained by stock estimates (Ricard et al., 2012). We then perform simulations with our Fe-limited version 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 in the catchability of biomass at 7 % yr⁻¹ to reproduce the historical progression of the global fishery (Galbraith et al., 2017). The simulations slightly underestimate the observed LME catch peak of ∼ 110 Mt yr⁻¹, and the variation between LME peaks corresponds to observation with a squared Pearson product-moment correlation coefficient r²=0.42, p value $<{\mathrm{10}}^{-\mathrm{9}}$, when averaged across all ensemble members (Figs. S1 and S2 in the Supplement). 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.

2.2 Nutrient content of fish

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 ontogeny as organisms grow in size, but analysis of the available data (see Supplement) shows little to no systematic variation in 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 10 g and 100 kg (from SAUP 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 (Tables S2–S3 in the Supplement). As a result, we did not attempt to account for invertebrates separately, but we applied the fish nutrient ratios to all ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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 µmol Fe (mol C)⁻¹ (Galbraith et al., 2019).

2.3 Nutrient cycling by fish

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., g N d⁻¹), 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=(\mathrm{1}-\mathit{\alpha })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, $G_i=\mathit{\gamma }{A}_i=\mathit{\gamma }\mathit{\alpha }{I}_i$, where γ is the somatic assimilation efficiency. The trophic efficiency, τ, is commonly defined as the ratio 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=(\mathrm{1}-\mathit{\gamma })\mathit{\alpha }{I}_i$ (Fig. 3a, b). We define the cycling rate of a nutrient through the fish as $I_i-G_i=E_i+F_i$, which can be thought of as the biologically processed outputs of the fish.

Figure 3 Schematic of the flow of elements (a) within an individual and (b) through the size spectrum. I: ingestion; G: growth and reproduction; A: absorption; E: excretion; F: egestion (feces); α: absorption efficiency, i.e., assimilation efficiency net of egestion; γ: somatic assimilation efficiency, i.e assimilation allocated to growth; f_in: flux entering the spectrum; f_out: flux going out of the spectrum.

Galbraith and Martiny (2015)

Table 1 Mean nutrient content in fish, zooplankton, and phytoplankton and mean absorption efficiencies (A) of N, P, and Fe for fish used in this study.

Data from ^a Czamanski et al. (2011) (percent converted back to wet weight using a 25 % dry weight in wet weight (Galbraith et al., 2019), geometric mean, and P absorption efficiency computed from the linear regression between predator and prey C $/$ P ratio). ^b Schindler and Eby (1997). ^c Galbraith et al. (2019). ^d Thodesen et al. (2001). ^e Prabhu et al. (2016). ^f Moore et al. (2013). “ww” is for wet weight, Fe-poor and Fe-rich refer to the conditions in which the organisms lived, and low/high are the low and high estimates from gathered data in Galbraith et al. (2019). Standard deviations are the arithmetic standard deviations, associated with the arithmetic mean. Ratio mean values are geometric means, and the ranges are the 95 % confidence interval, except for fish Fe : C for which the range is estimated from Galbraith et al. (2019).

Based on these terms, we follow Bianchi et al. (2021) to estimate the total flux of carbon through the ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ based on the 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. Because the model mean predator–prey mass ratio for the ensemble is 0.6 × 10⁴, there is very little predation within the resolved size spectrum (there are only 4 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$. Any consequent underestimate is therefore likely to be significantly less than the value of τ(C) (i.e., significantly less than 17 %).

Thus, our carbon cycling rate equation for ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ is given by the output flux equation:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{f}_{\mathrm{out}}\left(\mathrm{C}\right)& =\sum _i\left(I_i\right(\mathrm{C})-G_i(\mathrm{C}\left)\right)=\sum _i\left(E_i\right(\mathrm{C})+F_i(\mathrm{C}\left)\right)\\ & =f_{\mathrm{in}}\left(\mathrm{C}\right)-\sum _i\left(G_i\right(\mathrm{C}\left)\right).\end{array}\end{array}$$

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-to-carbon mass ratio within the prey, (X : C)_prey, and the absorption efficiency of nutrient X, α_X, to compute

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& f_{\mathrm{in}}\left(X\right)=f_{\mathrm{in}}\left(\mathrm{C}\right)(X:\mathrm{C}{)}_{\mathrm{prey}},\\ & E\left(X\right)=\sum _i\left({\mathit{\alpha }}_X{I}_i\right(X)-G_i(X\left)\right),\\ & F\left(X\right)=f_{\mathrm{in}}\left(X\right)(\mathrm{1}-{\mathit{\alpha }}_X),\end{array}\end{array}$$

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 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 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 µmol Fe (mol C)⁻¹, (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, and (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 in Fe : C, we assumed that Fe-poor conditions are encountered in HNLC regions, which are determined by a concentration of surface nitrate $\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$ larger than 5 mmol N m⁻³. In order to take into account the gradient between these regions, we locally weighted zooplankton Fe content using a Michaelis–Menten function.

$$\begin{array}{}\text{(3)}& \begin{array}{rl}(\mathrm{Fe}:\mathrm{C}{)}_{\mathrm{zoo}}& =(\mathrm{Fe}:\mathrm{C}{)}_{\mathrm{zoo}}^{\text{Fe-rich or high}}+\left((\mathrm{Fe}:\mathrm{C}{)}_{\mathrm{zoo}}^{\text{Fe-poor or low}}\right.\\ & \left.-(\mathrm{Fe}:\mathrm{C}{)}_{\mathrm{zoo}}^{\text{Fe-rich or high}}\right){\displaystyle \frac{[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{surf}}}{\mathrm{5}+[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{surf}}}}\end{array}\end{array}$$

2.4 Primary producer demand, nutrient concentrations, export, and atmospheric deposition

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 ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ surface concentrations as described in Galbraith and Martiny (2015):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& (\mathrm{N}:\mathrm{C}{)}_{\mathrm{phyto}}=\mathrm{12.5}\mathit{\%}+\mathrm{3}\mathit{\%}{\displaystyle \frac{[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{surf}}}{\mathrm{0.32}+[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{surf}}}},\\ & (\mathrm{P}:\mathrm{C}{)}_{\mathrm{phyto}}=\mathrm{0.6}\mathit{\%}+\mathrm{0.69}\mathit{\%}[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}{]}_{\mathrm{surf}},\end{array}\end{array}$$

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 Eq. (3).

We then use the phytoplankton nutrient ratios to compute the export of nutrients from a satellite-based estimate of total C export (Dunne et al., 2007).

We use the World Ocean Atlas-observed ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ water concentrations (Garcia et al., 2013) 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 TOPAZ2 model represents the cycles of different elements from carbon to calcite and Fe with 30 different tracers and the dynamics of three groups of phytoplankton. Surface concentrations are computed using the 2002–2019 annual mean euphotic depth from Aqua MODIS (2018).

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) (Fig. S3). The catch rate and its spatial distribution are simulated by the coupled economic–ecological model and systematically differs from the actual historical peak in that the model ensemble simulates higher catch rates in the open ocean than observed.

3 Fish biomass: a living pool of nutrients

Our results show that the nutrients contained within ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass represent a non-negligible proportion of the ambient dissolved concentrations in areas where these concentrations in seawater are low and/or where ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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, where most simulated fish biomass occurs (Fig. 4a). Globally, the estimated pristine biomass of ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$, 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).

Figure 4 Modeled commercial pristine fish biomass mean nutrient content and relative to surface nutrient concentrations. (a) P content (mmol P m⁻²), N content (mmol N m⁻²), and Fe content (low and high estimates, µmol Fe m⁻²) of the global pristine ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass, and (b) N content relative to surface ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations (%), (c) P content relative to surface ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations (%), and (d) Fe content (low and high estimates) relative to modeled surface dissolved Fe (%) from the TOPAZ model (Dunne et al., 2013). All surface concentrations are integrated over the 2002–2019 annual mean euphotic depth (Aqua MODIS, 2018). The N, P, and Fe content of fish is represented on a single map given that we used a globally constant nutrient ratio for each element, so that all spatial variation is caused by fish biomass.

Table 2 Table of values in LMEs from the model ensemble simulations in the pristine state and at the global peak catch. This table contains integrated values of nutrient content in ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass (Tg) and the ratio of nutrient content in ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass with surface nutrient concentrations (%).

The N content of ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass is locally on the same order of magnitude as ambient surface ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations, exceeding 50 % of $\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$ 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass represents a high proportion of available ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ 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. ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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).

In summary, the storage of nutrients by ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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, ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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.

3.1 Comparison to previous estimates

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 Tg N (Table 2). Their computation is based on an estimation of fish biomass of 0.9 Gt (Jennings et al., 2008) while our ensemble pristine ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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), 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 roughly 3-fold 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 g Fe kg⁻¹ of wet weight (ww) for ray-finned fish, that is about 3–12 times larger than our 10–200 µmol Fe (mol C)⁻¹ range, equivalent to 0.006–0.12 g Fe kg⁻¹ ww (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.

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.

3.2 Nutrient content variations in fish

Many factors contribute to variations in 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 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 in 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 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 (Fig. S4). In this data set, 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 has indeed shown 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 Fe do 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 in fish biomass density in the ocean. Thus, most of the spatial variations in the nutrient content of fish are likely to be due to variations in fish abundance rather than species assemblages, and the uncertainty related to nutrient ratios is likely to be small compared to the uncertainty on fish biomass.

4 Nutrient cycling by the commercial fish biomass

Because they are capable of large-scale movement and alter the stoichiometry of particles, fish can have impacts on nutrient cycling that differ from single-celled heterotrophic plankton. In this section, we gauge the potential scale of these impacts by estimating the rates at which fish cycle nutrients, what fraction of primary productivity (PP) this cycling represents, and how much it can contribute to the export of nutrients from the euphotic zone as sinking egested materials. Note that this is intended only to illustrate the potential magnitude, which could be built upon with coupled fish biogeochemistry modeling.

Table 3 Table of values from the model ensemble simulations in the pristine state and at the global peak catch at both the global scale and in LMEs. This table contains integrated values of the amount of nutrients cycled by the ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass (Tg yr⁻¹), the ratio of this cycling with the global primary producers' demand for these nutrients (%), the amount of nutrients egested by the ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass (Tg yr⁻¹), and its ratio with the exported nutrient quantities (%).

The global cycling, i.e., excretion plus egestion, of nutrients by the pristine ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass represents about 210±113 Tg N yr⁻¹, 15.6±8.0 Tg P yr⁻¹, and 0.12–0.77 Tg Fe yr⁻¹, of which about half is in the LMEs (Table 3). Like nutrient storage, modeled cycling by fish is larger where the biomass is higher (Figs. 5 and S5). The three different ways of computing Fe cycling by the commercial fish biomass show similar spatial patterns, but Fe cycling is larger 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ might be reduced in HNLC regions (principally in the Southern Ocean and the subarctic Pacific Ocean).

Figure 5 Iron cycling by ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass computed using (a) a mean Fe : C in zooplankton, (b) a spatial variation between zooplankton Fe : C mean values in Fe-poor and Fe-rich areas and based on [NO₃] concentrations, and (c) the same spatial variation but using the low and high estimates of zooplankton Fe : C from data. (For more details see Methods.)

4.1 Nutrient cycling by commercial fish relative to primary production

The modeled N cycling by pristine ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass contributes on average to 2.2±1.2 % of the N demand of primary producers (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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ relative to PP demand for Fe is slightly more significant than the ratios for N and P, as it represents up to 13 % of the PP demand for Fe in some coastal areas and less than 10 % everywhere else (Fig. 6e).

Figure 6 Total nutrient cycling and egestion relative to export. Ratio (%) between the amount of nutrients cycled by the modeled pristine global ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass and estimated primary producers' demand for (a) N and (b) P, ratios of the modeled amount of nutrients egested by the pristine ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass and export at the base of the euphotic zone of (c) N and (d) P, and high-end estimate of the ratio (%) between the modeled amount of (e) Fe cycled by the modeled pristine global ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass and estimated primary producers' Fe demand, and (f) Fe egested by the pristine ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass and Fe export at the base of the euphotic zone. The high-end estimates are obtained using Fe cycling computed from the weighted spatial variation between the low and high Fe : C values of zooplankton and the weighted spatial variation between the averaged Fe : C ratios of phytoplankton in Fe-poor and Fe-rich conditions.

Our global estimates are broadly consistent with the order of magnitude influence of fish estimated by prior local studies. For all the coral reefs in the ocean, Allgeier et al. (2014) estimated that the total fish community supplies about 1.2 Tg N yr⁻¹, which is about 0.6 % of our global pristine estimate (0.8 % compared to global cycling at peak catch) or 1.2 % of the LME pristine estimate (2.1 % of our N cycling by ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass at global peak catch) (Table 3)). This is consistent with the fact that coral reefs cover about 255 000–600 000 km² (Spalding and Grenfell, 1997), which is about 0.4 %–0.8 % of the LME area. Hernández-León et al. (2008) estimated that zooplankton supply about 1780 Tg N yr⁻¹ worldwide, representing 12 %–23 % of the requirements of phytoplankton and bacteria. This estimated zooplankton cycling is about 11 times our modeled recycling rate by the pristine ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass and 6–12 times what ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ cycling could provide for primary producers. Given that the biomass distribution of all marine organisms has a slope of 0 (Hatton et al., 2021) (Fig. S7), and assuming that 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 (Fig. S8). 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 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, ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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 Tg Fe yr⁻¹. Our modeled estimate range of 0.12–0.77 Tg Fe yr⁻¹, based on the high and low estimates (Figs. 6e and S6), 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 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ cycling coherently might be able to sustain a larger part of primary productivity than the current whale population 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.

4.2 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 quantify how much ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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, ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$-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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ biomass egests 29.4±15.9 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 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 transfer 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.

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; Galbraith et al., 2019).

5 Extended size spectrum and total fish biomass

Until now we have considered ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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⁶ g sharks ($F_{\mathrm{1}\mathrm{g}}^{\mathrm{1000}\mathrm{kg}}$) would compare to our estimates in two steps. First, we take the Bianchi et al. (2021) estimate that the ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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 (Fig. S7) consistent with 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_{\mathrm{1}\mathrm{g}}^{\mathrm{1000}\mathrm{kg}}$ exceeds ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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 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.37 mol C m⁻² g⁻¹ (Fig. S8). Based on this slope, and the biomass slope of 0, we would expect that the ${\mathrm{CTF}}_{\mathrm{1}\mathrm{g}}^{\mathrm{1000}\mathrm{kg}}$ would have cycled about 2.4 times more C in the pristine ocean compared to the ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$, with the 1–10 g size range cycling 40 % more C than the rest of the spectrum combined. Taken all together, the total cycling by $F_{\mathrm{1}\mathrm{g}}^{\mathrm{1000}\mathrm{kg}}$ would have been a factor of 4.8 greater than the ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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

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 Tg N yr⁻¹, 1.2±0.3 Tg P yr⁻¹, and 0.09–1.8 × 10¹⁰ g Fe yr⁻¹ from the ocean at the time of global peak catch, of which a little less than 50 % is in LMEs (Fig. 7, 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., 2021), the total catch varies between ensemble members and tends to be overestimated in the open ocean, even with the Fe limitation 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 (see Bianchi et al., 2021, for more details).

Figure 7 Distribution of the modeled amount of N (g N m⁻² yr⁻¹), P (g P m⁻² yr⁻¹), and Fe ( µg Fe m⁻² yr⁻¹) extracted from the ocean at the time of the global peak catch. The two color bars for Fe represent the low and high estimates based on the 95 % confidence interval for Fe : C values in fish (10–200 µmol Fe (mol C)⁻¹).

Table 4 Table of values in LMEs from the model ensemble simulations in the pristine state and at the global peak catch. This table contains integrated values of the amount of nutrients removed by fishing in LMEs (Tg yr⁻¹) and the ranges of values of the global nutrient inputs to the ocean from the literature (Tg yr⁻¹).

^a Brahney et al. (2015). ^b Mahowald et al. (2009). ^c Okin et al. (2011). ^d Myriokefalitakis et al. (2016). ^e Fowler et al. (2013). ^f Ito (2015). ^g Wang et al. (2015). ^h Gruber and Galloway (2008). ⁱ Benitez-Nelson (2000). ^j Kanakidou et al. (2012). ^k Moreno and Haffa (2014).

6.1 Nitrogen

Our estimate of N extraction from the ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ catch is consistent with previous computations, and the differences mostly reflect the fish biomass estimates. For example, using catch data, Maranger et al. (2008) estimated the amount of N returned to land via fishing in LMEs to be about 0.9 Tg N yr⁻¹ in the 1960s and 2.3 Tg N yr⁻¹ in 2000. Our estimate of the N content of fish catch at global peak catch in LMEs is 2.8±0.4 Tg N yr⁻¹, which is slightly larger than Maranger et al. (2008) estimate in 2000. Allgeier et al. (2017) also estimated the amount of N globally harvested to be about 2.072 Tg N yr⁻¹ using the Maranger et al. (2008) 2.6 % N content and the FAO catch data. The FAO data only contain reported catches and consequently is lower than the SAUP data used in Maranger et al. (2008) and used to calibrate BOATS, which explains part of the difference.

Our study framework allows the spatial comparison of extracted nutrients to nutrient inputs, which was not the case in previous 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 high and deposition is low (Fig. S3). 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.

Figure 8 Ratio (%) between simulated extracted nutrients and current aeolian soluble nutrient inputs at the surface of the ocean for (a) N and (b) Fe. “low” and “high” refer to the use of 10 or 200 µmol Fe (mol C)⁻¹ in fish, respectively.

6.2 Phosphorus

Similarly to N, our estimate of P extraction by fishing is coherent with previous work and shows that it is very small compared 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 Tg P yr⁻¹, of which 0.6±0.2 Tg P yr⁻¹ occurs 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 wild captures. Our catch estimate at global peak catch overestimates the high seas catch as discussed at the beginning of this section and in Bianchi 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 exceed the atmospheric deposition of soluble P as P inputs to the ocean mostly occur through riverine inputs (Table 4). Consequently, catch transfers P from the ocean to land whereas 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 supply P to the surface ocean in quantities that most likely render P extraction by fishing relatively insignificant.

6.3 Iron

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) investigated the extent to which commercial catch has globally translocated Fe from the ocean to land. They estimated the global rate of translocation of Fe to be between 0.007 and 0.03 Tg in 2010. Our modeled global range of Fe removal is about 0.0009–0.018 Tg Fe yr⁻¹ (Fig. 7, Table 4). Our lower estimated values can once again be explained by the difference in the Fe : C ratios used for fish. Even though our estimate is lower, it shows that locally Fe extraction can be significant compared to Fe inputs from dust deposition. Although the high-end estimate of Fe extracted is globally small relative to modeled soluble Fe deposition from Mahowald et al. (2009), it reaches values larger than 100 % in the coastal eastern equatorial Pacific and in some other coastal areas such as western South Africa, northern Europe, and Canada (Fig. 8b), where modeled Fe deposition is small and harvest is high (Fig. S3). Contrary to N and P, Fe has a much shorter residence time and thus is subject to local perturbations, among which Fe extraction by fishing could be important.

6.4 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 nutrients 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 Tg P yr⁻¹ (Tables 4 and S1), then the extracted P flux represents 2.5 % of the global input flux (1.2 % over the LMEs). For N, global catch represents about 2 % of the combined N inputs from atmospheric deposition (49.6 Tg N yr⁻¹), rivers (80 Tg N yr⁻¹), and N₂ fixation (140 Tg N yr⁻¹) (Tables 4 and S1).

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 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).

6.5 Reductions in nutrient cycling caused by fishing

Fishing has had a dramatic influence on nutrient cycling by ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$, as it has permanently removed a large amount of 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 a larger size class at the time of the global peak catch compared to a pristine state (Fig. S7). 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 the reduction in overall cycling rate is roughly half the reduction of biomass.

Finally, as the size classes below 10 g are not resolved (Fig. 2), we are not able to account for biomass changes that fishing might induce through trophic cascades, or how it reverberates up to fish through food supply (Dupont et al., 2022), which has the potential to further modify fish-mediated nutrient cycling.

7 Conclusions

In this study, we estimate the amount of N, P, and Fe contained in and cycled by the global ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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 a significant amount of P from the open ocean compared to external inputs (mainly riverine). In general, local cycling of N and P by fish is less significant than Fe cycling by fish because N and P are resupplied globally 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ 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 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 ${\mathrm{CTF}}_{\mathrm{10}\mathrm{g}}^{\mathrm{100}\mathrm{kg}}$ simulated by our model tends to be on the order of only a few percent of total surface nutrient budgets and fluxes, we estimate 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 estimated 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 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.