This is the accepted version of the journal article: Heneghan, Ryan F.; Everett, Jason D.; Sykes, Patrick; [et al.]. «A functional size-spectrum model of the global marine ecosystem that resolves zooplankton composition». Ecological modelling, Vol. 435 (November 2020), art. 109265. DOI 10.1016/j.ecolmodel.2020.109265 This version is available at https://ddd.uab.cat/record/274266 under the terms of the license Title A functional size-spectrum model of the global marine ecosystem that resolves zooplankton composition 5 Authors Ryan F. Heneghan1,2,3*, Jason D. Everett1,4,5, Patrick Sykes1, Sonia D. Batten6, Martin Edwards7, Kunio Takahashi8, Iain M. Suthers4, Julia L. Blanchard5, Anthony J. Richardson1,9 1 Centre for Applications in Natural Resource Mathematics, School of Mathematics and 10 Physics, University of Queensland, Brisbane, QLD, Australia 2 Institut de Ciència i Tecnologia Ambientals, Universitat Autònoma de Barcelona, Barcelona, Spain 3 School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia 15 4 School of Biological, Earth and Environmental Sciences, University of New South Wales, Sydney, NSW, Australia 5 Institute of Marine and Antarctic Studies, University of Tasmania, Hobart, TAS, Australia 6 The CPR Survey, Marine Biological Association, Nanaimo, BC, Canada 7 Marine Biological Association of the UK, Citadel Hill, The Hoe, Plymouth, United Kingdom 20 8 National Institute of Polar Research, Tokyo, Japan 9 CSIRO Oceans and Atmosphere, Queensland Biosciences Precinct, St Lucia, QLD, Australia *Corresponding Author: ryan.heneghan@gmail.com 25 Keywords Zooplankton, functional traits, size-spectrum model, marine ecosystem, modelling Highlights • A trait-based size-spectrum marine ecosystem model to resolve the global 30 zooplankton community • Zooplankton groups represented by body size range, size-based feeding characteristics and carbon content 1 • Emergent patterns of zooplankton abundance, biomass and growth rates agree well with empirical estimates 35 • Resolving zooplankton functional complexity has implications for the total biomass of fish, highlighting the critical importance of including zooplankton functional complexity in ecosystem models Abstract Despite their critical role as the main energy pathway between phytoplankton and fish, the 40 functional complexity of zooplankton is typically poorly resolved in marine ecosystem models. Trait-based approaches—where zooplankton are represented with functional traits such as body size—could help improve the resolution of zooplankton in marine ecosystem models and their role in trophic transfer and carbon sequestration. Here, we present the Zooplankton Model of Size Spectra version 2 (ZooMSSv2), a functional size-spectrum model that resolves 45 nine major zooplankton functional groups (heterotrophic flagellates, heterotrophic ciliates, larvaceans, omnivorous copepods, carnivorous copepods, chaetognaths, euphausiids, salps and jellyfish). Each group is represented by the functional traits of body size, size-based feeding characteristics and carbon content. The model is run globally at 5 resolution to steady-state using long-term average temperature and chlorophyll a for each grid-cell. 50 Zooplankton community composition emerges based on the relative fitness of the different groups. Emergent steady-state patterns of global zooplankton abundance, biomass and growth rates agree well with empirical data, and the model is robust to changes in the boundary conditions of the zooplankton. We use the model to consider the role of the zooplankton groups in supporting higher trophic levels, by exploring the sensitivity of steady- 55 state fish biomass to the removal of individual zooplankton groups across the global ocean. Our model shows zooplankton play a key role in supporting fish biomass in the global ocean. For example, the removal of euphausiids or omnivorous copepods caused fish biomass to decrease by up to 80%. By contrast, the removal of carnivorous copepods caused fish biomass to increase by up to 75%. Our results suggest that including zooplankton complexity in 60 ecosystem models could be key to better understanding the distribution of fish biomass and trophic efficiency across the global ocean. 1 Introduction 2 Marine ecosystem models are valuable tools for understanding marine ecosystem function, 65 structure and productivity under global change. However, the marine ecosystem is complex, and there remains a lack of unified understanding of key physiological and ecological processes (Tittensor et al., 2018). Modellers address this uncertainty by making assumptions about which ecosystem processes to represent in their models, depending on their objectives. A common assumption across higher trophic level modelling frameworks is to represent 70 zooplankton—the main grazers of phytoplankton and bacteria, and prey of small fish—very simply (Everett et al., 2017). There are two common approaches to representing zooplankton in marine models focussed on higher trophic levels. First, phytoplankton and zooplankton dynamics are driven externally 75 by nutrient-phytoplankton-zooplankton-detritus models (Maury, 2010; Woodworth-Jefcoats et al., 2013; Christensen et al., 2015; Carozza et al., 2016; Petrik et al., 2019) and then this is used to drive higher trophic levels in an ecosystem model. Second, phytoplankton and small zooplankton are lumped together as a resource for the smallest fish size classes in size- spectrum models (Blanchard et al., 2009, 2012; Law et al., 2009; Datta et al., 2010; Scott et 80 al., 2014; Zhang et al., 2015, 2016). These simple representations of zooplankton, together with the similar lower trophic level structure among models, precludes a deeper understanding of the role of zooplankton in ecosystem functioning. Even with a simple representation of zooplankton, many modelling studies suggest that the productivity and structure of higher trophic levels are sensitive to the representation of lower 85 trophic levels (Friedland et al., 2012; Jennings and Collingridge, 2015; Heneghan et al., 2016; Dam and Baumman, 2017). Mitra et al. (2014) demonstrated that in a modelled plankton food web, trophic dynamics were sensitive to small changes in parameterisation of zooplankton feeding rates. Similarly, Fuchs and Franks (2010) found that zooplankton with large predator- prey mass ratios (PPMR, i.e., a predator is a lot larger than its prey) gave rise to a flatter 90 plankton abundance size-spectra (relatively more large organisms), in comparison to zooplankton with small PPMRs, which led to a steeper plankton size-spectra. Moving beyond the plankton, Jennings and Collingridge (2015) demonstrated that the productivity and total biomass of the global fish community was highly sensitive to how energy moved through the lower planktonic trophic levels, from phytoplankton to zooplankton and beyond. 3 95 One of the main reasons that modellers choose to simplify zooplankton is because of their tremendous taxonomic diversity, especially at the phylum level. A typical zooplankton sample might contain organisms from 10 phyla or more, each with different body sizes, biochemical compositions, growth rates, feeding preferences and reproductive modes (Hansen et al., 1994, 1997; Kiørboe et al., 2011; Mconville et al., 2017). While it is not possible to include 100 thousands of common plankton species in ecosystem models, an alternative is to model organisms based on traits such as body size. Modelling the marine ecosystem using body size instead of species identity is rooted in over 50 years of observations that show a remarkable consistency in the distributions of organism abundance or biomass against body size (size- spectra) across a range of marine ecosystems (Blanchard et al., 2017). For zooplankton, body 105 size is strongly related to trophic position (Andersen et al., 2016a) and the size-based feeding behaviour of different zooplankton groups has been suggested to structure the community across environmental gradients (Mitra and Davis, 2010; Barton et al., 2013). Owing to their body size, zooplankton are the primary grazers of phytoplankton, which span nine orders of magnitude in body size, from picoplankton (0.2–2 m equivalent spherical diameter, ESD; 110 10−14.5–10−11.5 g wet weight) to microplankton (>20 m ESD; >10−8.4 g wet weight). For zooplankton, other traits such as predator-prey mass ratio (PPMR) and body composition are also important, since these factors also play a role in organism fitness (McGill et al., 2006; Litchman et al., 2013; Andersen et al., 2016a; McConville et al., 2017). For instance, 115 phytoplankton size structure changes across environmental gradients (Agawin et al., 2000; Brewin et al., 2010; Barnes et al., 2011). Changes in the phytoplankton have implications for the structuring of the zooplankton community, which in turn affects how energy from primary production is transported to higher trophic levels. This is because zooplankton exhibit vast diversity in their feeding behaviour, with PPMRs varying >7 orders of magnitude across 120 different functional groups, from ~10 for carnivorous copepods to ~108 for salps and larvaceans (Bone, 1997; Wirtz, 2012). Another functional trait that varies strongly across zooplankton groups is body composition, and this can have substantial effects on energy transfer to higher trophic levels. Carbon is the primary structural component of zooplankton (Kiørboe, 2013) and zooplankton vary widely in their carbon content between gelatinous 125 groups with low carbon density and non-gelatinous groups with high carbon density (Kiørboe, 2013; McConville et al., 2017). Thus, a zooplankton community dominated by more 4 gelatinous groups offers less nutritional value and growth potential for higher trophic levels. Moreover, changes in carbon content affects the relative fitness of different zooplankton groups, since critical physiological and competitive processes such as metabolism, search 130 rate, and average growth efficiency scale with carbon across zooplankton groups (Acuña et al., 2011; Kiørboe, 2011; Kiørboe and Hirst, 2014; McConville et al., 2017). The diversity of traits in the zooplankton means that global ecosystem modellers must consider more than just their body size if they are to improve the realism of zooplankton in 135 their models. In the past 15 years, trait-based modelling approaches have been applied to explain the distribution of phytoplankton groups (Follows et al., 2007; Edwards et al., 2013; Acevedo-Trejos et al., 2015), and there is a growing literature applying the approach separately to zooplankton (Fuchs and Franks, 2010; Brun et al., 2016; Heneghan et al., 2016; Schnedler-Meyer et al., 2016; Prowe et al., 2018) and fish communities (Stuart-Smith et al., 140 2013; Blanchard et al., 2014; van Denderen et al., 2018; Petrik et al., 2019). However, to our knowledge, there are no global marine ecosystem models resolving both higher and lower trophic levels that incorporate the diversity of traits in the zooplankton. Recently, Blanchard et al., (2017) proposed unifying size-spectrum ecosystem modelling (which only resolves the body-size distribution of organisms) with trait-based modelling in a new functional size- 145 spectrum framework. This framework allows modellers to take advantage of the powerful role of body size in structuring the marine ecosystem as well as other functional traits that are important, and offers a promising way forward for resolving the functional complexity of the zooplankton within a model of the entire marine ecosystem. This is the approach we adopt here. 150 Using functional traits such as body size and PPMR to resolve the diversity of the zooplankton in marine ecosystem models is not a new idea. Twenty-five years ago, Hansen et al. (1994) hypothesised that with knowledge of the size selectivity (e.g., PPMR) of different zooplankton groups, it would be possible to construct a simple size-based model of the pelagic food web. 155 More recently, zooplankton feeding traits such as PPMR were incorporated in a functional size-spectrum model (Heneghan et al., 2016), which showed that fish biomass and ecosystem stability were sensitive to zooplankton feeding characteristics. This model (Zooplankton Model of Size Spectra, ZooMSSv1) allowed a qualitative assessment of how zooplankton traits 5 such as PPMR affected energy transfer from phytoplankton to fish, and highlighted the 160 significant implications of variations in zooplankton community traits for high trophic level biomass and ecosystem stability. However, it did not represent zooplankton community composition or food quality and so could not be used to explore how those factors might affect ecosystem structure and transfer efficiency in the global ocean. Here, we follow the work of Heneghan et al., (2016) to present ZooMSSv2 (referred to 165 hereafter as ZooMSS), the first functional size-spectrum model of the marine ecosystem to resolve phytoplankton, nine zooplankton functional groups (heterotrophic flagellates and ciliates, omnivorous and carnivorous copepods, larvaceans, euphausiids, salps, chaetognaths and jellyfish) and three size-based fish groups. Zooplankton functional groups are resolved using their body-size ranges, size-based feeding characteristics and carbon content. The 170 model is run to steady-state, and the modelled zooplankton community emerges across global environmental gradients, depending on the functional traits of the different groups. The model is not fit to observational data, rather parameter values for the model’s functional trait are obtained from experimental results in the literature and the model’s ability to reproduce observed patterns in zooplankton biomass, abundance and growth is assessed. 175 The model is generally able to reproduce steady-state global patterns of zooplankton biomass and abundance, as well as empirical growth rates from flagellates to jellyfish. What is more, model output is robust to uncertainties in the boundary conditions of the smallest size class of each group. By resolving zooplankton functional diversity, the model allows us for the first time to examine the unique roles of different zooplankton groups in supporting higher trophic 180 level biomass across the global ocean. We finish by exploring the change in global fish biomass when different zooplankton groups are removed from the model, as well as the change in biomass when the zooplankton are represented simplistically as a single functional group. 2 Methods We developed the Zooplankton Model of Size Spectra version 2 (ZooMSS) based on the 185 prototype of Heneghan et al. (2016). ZooMSS uses the functional size-spectrum framework (Blanchard et al., 2017) to resolve the body size ranges, size-based feeding characteristics and carbon content of nine zooplankton groups and three fish groups. Parameter values for the functional traits are not fit to data here, but are determined using published estimates 6 available in the literature (Table 2). The model is run to steady-state on 5 grid squares across 190 the global ocean (Figure 1). For each region, the model is forced with the long-term mean satellite sea surface temperature and chlorophyll a from MODIS-Aqua. Scripts to run the model, as well as global forcings for sea surface temperature and chlorophyll a, are available for download at https://github.com/MathMarEcol/ZoopModelSizeSpectra. 2.1 Zooplankton Model of Size Spectra version 2 (ZooMSS) 195 ZooMSS represents the marine ecosystem as three communities: phytoplankton, zooplankton and fish. The zooplankton community consists of nine of the most abundant zooplankton groups, and the fish community is made up of a small, medium and large group (Figure 1; Table 1-3). 200 Figure 1: The modelled size-spectrum, over a 5 region of the global ocean. The phytoplankton spectrum is held constant, and abundances of zooplankton and fish groups are governed by size-dependent processes of growth and mortality. Total Fish is the sum of the three fish groups (small, medium and large). Dynamics of the phytoplankton are not explicitly resolved in the model, rather the size 205 structure of the phytoplankton community is estimated directly from satellite chlorophyll a 7 observations (Brewin et al., 2010; Barnes et al., 2011; Hirata et al., 2011). Abundances of the zooplankton and fish communities are driven by size-dependent processes of growth and mortality, and the abundances of each functional group governed by separate second-order McKendrick-von Foerster equations: 𝜕 𝜕 1 𝜕2 210 𝑁𝑖 (𝑤, 𝑡) = − (𝑔𝑖 (𝑤, 𝑡)𝑁𝑖 (𝑤, 𝑡)) − 𝜇𝑖 (𝑤, 𝑡)𝑁𝑖 (𝑤, 𝑡) + (𝑑 (𝑤, 𝑡)𝑁𝑖 (𝑤, 𝑡)). (E1) 𝜕𝑡 𝜕𝑤 2 𝜕𝑤 2 𝑖 The density of individuals in group 𝑖 of weight 𝑤 at time 𝑡 per m3 is given by 𝑁𝑖 (𝑤, 𝑡). Individual processes in group 𝑖 are given by growth 𝑔𝑖 (𝑤, 𝑡), mortality 𝜇𝑖 (𝑤, 𝑡) and diffusion rates 𝑑𝑖 (𝑤, 𝑡) (Figure 2). We use the McKendrick-von Foerster equation to govern zooplankton and fish communities because it is a popular choice for modelling fish-focused 215 size-spectrum models, and it is similar to governing equations in plankton-focused size-based models (Baird and Suthers, 2007; Fuchs and Franks, 2010; Zhou et al., 2010; Ward et al., 2012, Cuesta et al., 2018). Looking ahead, this means that future model developments that resolve the dynamics of the phytoplankton, or more of the functional diversity of the zooplankton, could be achieved with the same—or similar—governing equation used here. 220 Figure 2 Overview of the individual-level processes that drive the abundance of the 8 zooplankton and fish groups. From the perspective of an individual from group 𝑖 of size 𝑤, the feeding rate 𝑓𝑖 is determined by the density of suitably sized-prey, here represented as single group 𝑗. The feeding rate is calculated using a size-selection function 𝜙, centred at size 𝑤 225 𝑤 ′ = 𝛽 (𝑤) where 𝛽𝑖 (𝑤) is the individual’s predator-prey mass ratio. Captured food then 𝑖 determines the individual’s growth and diffusion rates in each time step. Finally, this individual experiences predation mortality from larger organisms 𝜇𝑝 , as well as senescence mortality 𝜇𝑆𝑖 . From the perspective of a predator from group 𝑖, the feeding rate 𝑓 (g yr-1) on prey group 𝑗 230 depends on the search rate of group 𝑉𝑖 (𝑤) (m3 yr-1; E13), and the density of suitable prey (g yr−1): 𝑤 𝑓𝑖𝑗 (𝑤, 𝑡) = 𝑉𝑖 (𝑤) ∫𝑤 𝜙𝑖 (𝑤, 𝑤 ′ )𝑁𝑗 (𝑤 ′ , 𝑡)𝑤 ′ 𝑑𝑤 ′ , (E2) 𝑝 where 𝜙𝑖 (𝑤, 𝑤 ′ ) (E11, Table 1; all subsequent equations are also found in Table 1) is the probability a predator of size 𝑤 would consume an individual of size 𝑤’ and 𝑤𝑝 is the minimum 235 size of phytoplankton – the smallest body size in the system. The growth rate of a predator from group 𝑖, of size 𝑤 at time 𝑡, is fuelled by the consumption and conversion of prey biomass to new biomass (g yr−1): 𝑔𝑖 (𝑤, 𝑡) = 𝜏 ∑𝑗 𝐸𝑖𝑗 𝑓𝑖𝑗 (𝑤, 𝑡) , (E3) where 𝐸𝑖𝑗 is the growth conversion efficiency for predators of group 𝑖 eating prey from group 240 𝑗 (E14a,b), and τ is the effect of temperature on ingestion for group 𝑖 (E4). A drawback of the Type 1 functional response used here for feeding is that it ignores the effect of satiation at high food densities. This means that maximum growth rates of zooplankton and fish groups could be unrealistically fast in areas with high food availability, i.e., high chlorophyll a concentration. However, the effect of fast growth rates would be similar across all functional 245 groups, and so would not affect their relative fitness. Further, we will compare modelled growth rates at high chlorophyll concentrations with empirical data to assess whether the model’s growth rates are reasonable. The choice of functional response not only affect feeding rates at high or low food concentrations, but can also lead to different temporal dynamics between predators and prey (Gentleman and Neuheimer, 2008). However, our 250 focus here is on long term steady-state patterns in the composition of the global zooplankton 9 community and not on the dynamics of the different groups through time. Nevertheless, ecosystem dynamics have a real impact on marine ecosystems, particularly for zooplankton, which face seasonal fluctuations in prey abundance from phytoplankton blooms and busts. The representation of zooplankton functional response will therefore be an important 255 consideration for any future work exploring the temporal dynamics of the marine ecosystem across environmental gradients. Temperature effects are represented using a 𝑄10 temperature coefficient: 𝐾−𝐾ref ( ) 10 𝜏 = 𝑄10 , (E4) where 𝐾 and 𝐾ref are the temperatures in Kelvin of each 5° grid square, and the reference 260 temperature where 𝜏 = 1, respectively. From the perspective of prey of size 𝑤, the total mortality from predation by larger size classes 𝜇𝑝 (𝑤, 𝑡) (yr−1) is given by: ̅𝑗 𝑊 𝜇𝑝 (𝑤, 𝑡) = 𝜏 ∑ ∫ 𝜙𝑗 (𝑤 ′ , 𝑤)𝑉𝑗 (𝑤 ′ )𝑁𝑗 (𝑤 ′ , 𝑡)𝑑𝑤 ′ , (E5) 𝑗 𝑤 ̅𝑗 is the maximum size of a predator from group 𝑗. Since individuals grow over time, where 𝑊 an additional source of mortality from senescence was incorporated that increased with body 265 size (yr−1): 𝜌 𝑤 𝜇𝑆𝑖 (𝑤, 𝑡) = 𝜏𝛿 (𝑊 ) , (E6) 𝑠𝑖 where 𝛿 is the coefficient of senescence mortality, 𝜌 the exponent and 𝑊𝑠𝑖 the body size after which senescence mortality rapidly increases for an individual from group 𝑖. This senescence 270 mortality term also acts as a closure term for the largest size classes, by preventing a build- up of large individuals who are not exposed to predation (Andersen et al., 2016b). For an individual of size 𝑤, at time 𝑡, from group 𝑖, total mortality 𝜇𝑖 (𝑤, 𝑡) (yr−1) is given by summing predation and senescence mortality (yr−1): 𝜇𝑖 (𝑤, 𝑡) = 𝜇𝑝 (𝑤, 𝑡) + 𝜇𝑆𝑖 (𝑤, 𝑡). (E7) 10 275 Finally, the second-order diffusion term for an individual from group 𝑖 of size 𝑤 at time 𝑡 is (g2 yr−1): 𝑤 2 𝑑𝑖 (𝑤, 𝑡) = 𝑉𝑖 (𝑤) ∑(𝜏𝐸𝑖𝑗 ) ∫ (𝑤 ′ )2 𝜙𝑖 (𝑤, 𝑤 ′ )𝑁𝑗 (𝑤 ′ , 𝑡)𝑑𝑤 ′ . (E8) 𝑗 𝑤𝑝 We incorporate the diffusion term following the work of Datta et al., (2010), who demonstrated that stochastic variation in growth at the individual level (something that the 280 more commonly used first-order McKendrick-von Foerster equation does not consider) can be approximated at the functional group level by adding a diffusion term to the first-order McKendrick-von Foerster equation. Incorporating a diffusion term means that the model incorporates both the advection of biomass from smaller to larger size classes through the processes of predation and growth (𝒈𝒊 (𝒘, 𝒕); E3) but also includes an approximation of the 285 variability in the growth rates of each community from the stochastic variation in individual- level growth rates through time (𝒅𝒊 (𝒘, 𝒕); E8). As well as increasing model realism by approximating the functional group-level effects of individual-level variability, Datta et al., (2011) demonstrated that the diffusion term improves the stability of the numerical approximation of real world predator-prey interactions and individual growth trajectories in 290 the presence of large, destabilising predator–prey mass ratios, which have been shown to cause the community size spectra to develop travelling wave attractors (Law et al., 2009; Zhang et al., 2013; Heneghan et al., 2016). This is especially important here given our model includes zooplankton predator-prey mass ratios that are many orders of magnitude larger than predator-prey size ratios used in fish-focused size-spectrum models that use the first- 295 order McKendrick-von Foerster equation (see Table 2). 2.1.1 Parameterizing the phytoplankton For each 5° region, the density of phytoplankton 𝑁𝑃 (𝑤) of size 𝑤 is given by: 𝑁𝑃 (𝑤) = 𝑎𝑤 𝑏 . (E9) 300 The slope 𝑎, intercept 𝑏, and maximum size of the static phytoplankton spectrum were derived from temporally and spatially averaged satellite chlorophyll a from MODIS-Aqua (accessed via the GIOVANNI portal: https://giovanni.gsfc.nasa.gov/giovanni/), using the synoptic model developed by Brewin et al. (2010). The Brewin model gives an estimate of the 11 percentage contribution of 3 phytoplankton size classes—pico (0.2–2 m ESD), nano (2–20 305 m ESD) and microphytoplankton (>20 m ESD)—to the total chlorophyll a concentration (mg m−3 yr−1). Picophytoplankton constitute up to 75% of the biomass in low chlorophyll a (oligotrophic) waters, declining to <10% in high chlorophyll a (eutrophic) waters as micro- phytoplankton increase from <10% in oligotrophic waters to >75% in eutrophic waters, and nano-phytoplankton increase from 20% in oligotrophic to 45% in mid-chlorophyll a waters, 310 before declining to 15% in eutrophic waters (Figure 3 a). The contribution of microphytoplankton, and the maximum size of phytoplankton, both increase with chlorophyll a concentration (Brewin et al., 2010; Hirata et al., 2011; Barnes et al., 2011). We incorporated this change in the size range of the microphytoplankton with increasing chlorophyll a by linearly increasing its maximum size from 21–60 m ESD, 315 depending on the percentage contribution of the microphytoplankton group to total chlorophyll a. We used 60 m as the maximum ESD for the phytoplankton following Barnes et al.’s (2011) finding that 90% of phytoplankton are smaller than 55–65 m across polar, tropical and upwelling environments. Total chlorophyll a concentration for each of the 3 size classes was converted to grams wet weight assuming 1 g chlorophyll a = 50 g C (Zhou et al., 320 2010), and 1 g C = 10 g wet weight (Hansen et al., 1994, Boudreau & Dickie, 1992, Woodworth- Jefcoats et al., 2013) and the three size ranges were also converted from ESD to grams wet weight (assuming 1 cm3 = 1 g wet weight; Boudreau & Dickie, 1992). Finally, slope and intercept were found analytically with E8 (see Supplementary Material A1). Phytoplankton slopes we derived ranged from −1.2 to −0.77 across the global ocean (Figure 3 b), which is 325 similar to the range reported by previous empirical studies (Huete-Ortega et al., 2011; Marañón, 2015; Moreno-Ostos et al., 2015). We use satellite chlorophyll a observations to drive the abundance and structure of the phytoplankton community in each grid cell. Observed chlorophyll a is an estimate of the surface phytoplankton biomass in a region, given the processes of predation and nutrient 330 cycling and regeneration. Thus, by using observed satellite chlorophyll a our model is implicitly incorporating these important processes in the phytoplankton, which is sufficient for our present objective of resolving steady-state patterns of zooplankton biomass and composition across the global ocean, given in situ phytoplankton biomass. In global marine 12 ecosystem models, it is common practice to use satellite or earth system model estimates of 335 phytoplankton biomass (Blanchard et al., 2009; Maury, 2010; Christensen et al., 2015; Petrik et al., 2019) or primary production (Gascuel and Pauly, 2009; Carozza et al., 2016; Jennings and Collingridge, 2015) as inputs to drive higher trophic level processes without feedback from predation or nutrient cycling. We chose to use phytoplankton biomass—calculated with chlorophyll a—over primary production here because it is the most common measure of food 340 for zooplankton (Richardson and Verheye, 1999; Hirst & Bunker 2003). 13 Figure 3 a) Proportion of phytoplankton that is picoplankton (<2 μm ESD), nanoplankton (2 μm–20 μm ESD) and microplankton (>20 μm ESD), against log10(chlorophyll a), adapted from Brewin et al., (2015); b) Phytoplankton abundance spectrum at low (10−1.5 g m−3) and high 345 (100.5 g m−3) chlorophyll a. The slope 𝒂 and intercept 𝒃 for the two spectra were calculated using satellite chlorophyll a from MODIS-Aqua, and the synoptic model from Brewin et al., (2015) (see Supplementary Material A1 for more information). 350 14 Table 1 Model equations with their units and in-text reference numbers. Description Equation Units Equation Number Master equation: 𝜕 𝜕 𝑁 (𝑤, 𝑡) = − (𝑔 (𝑤, 𝑡)𝑁𝑖 (𝑤, 𝑡)) − 𝜇𝑖 (𝑤, 𝑡)𝑁𝑖 (𝑤, 𝑡) Rate of change of the 𝜕𝑡 𝑖 𝜕𝑤 𝑖 yr −1 E1 abundance of individuals of 1 𝜕2 + (𝑑 (𝑤, 𝑡)𝑁𝑖 (𝑤, 𝑡)) group 𝑖 of size 𝑤 at time 𝑡 2 𝜕𝑤 2 𝑖 Growth and mortality: 𝑤 Feeding rate of predator g yr −1 E2 𝑓𝑖𝑗 (𝑤, 𝑡) = 𝑉𝑖 (𝑤) ∫ 𝜙𝑖 (𝑤, 𝑤 ′ )𝑁𝑗 (𝑤 ′ , 𝑡)𝑤 ′ 𝑑𝑤 ′ group i on group j 𝑤𝑝 Individual growth rate for 𝑔𝑖 (𝑤, 𝑡) = 𝜏 ∑ 𝐸𝑖𝑗 𝑓𝑖𝑗 (𝑤, 𝑡) g yr −1 E3 group 𝑖 𝑗 Temperature effect for group 𝐾−𝐾ref - E4 ( ) 10 𝑖 𝜏 = 𝑄10 Individual predation rate ̅𝑗 𝑊 𝜇𝑝 (𝑤, 𝑡) = 𝜏 ∑ ∫ 𝜙𝑗 (𝑤 ′ , 𝑤)𝑉𝑗 (𝑤′)𝑁𝑗 (𝑤 ′ , 𝑡)𝑑𝑤 ′ yr −1 E5 𝑗 𝑤 𝜌 Senescence mortality 𝑤 yr −1 E6 𝜇𝑆𝑖 (𝑤, 𝑡) = 𝜏 𝛿 ( ) 𝑊𝑠𝑖 Total mortality 𝜇𝑖 (𝑤, 𝑡) = 𝜇𝑝 (𝑤, 𝑡) + 𝜇𝑆𝑖 (𝑤, 𝑡) yr −1 E7 𝑤 Individual diffusion term for 𝑑𝑖 (𝑤, 𝑡) = 𝑉𝑖 (𝑤) ∑(𝜏𝐸𝑖𝑗 )2 ∫ (𝑤 ′ )2 𝜙𝑖 (𝑤, 𝑤′)𝑁𝑗 (𝑤 ′ , 𝑡)𝑑𝑤′ group 𝑖 𝑤𝑝 g 2 yr −1 E8 𝑗 Phytoplankton spectrum: 𝑁𝑃 (𝑤, 𝑡) = 𝑎𝑤 𝑏 m−3 E9 Functional traits: Predator-prey mass ratio for 𝛽𝑖 (𝑤) = (exp(0.02 ln(𝐷𝑤 )2 − 𝑚 + 1.832))3 - E10 group 𝑖, size 𝑤 Size selection for group 𝑖 2 - E11 𝜙𝑖 (𝑤, 𝑤 ′ ) = exp [−(ln(𝛽𝑖 (𝑤)𝑤 ′ /𝑤)) /2𝜎𝑖2 ] /(𝜎𝑖 √2𝜋) Feeding kernel width 𝜎𝑖 = 0.05log10 (𝛽𝑖̅ ) + 0.33 E12 parameter for group 𝑖 Search rate for group 𝑖 𝑉𝑖 (𝑤) = 𝛾𝑤 𝛼 −1 m3 yr −1 E13 Average growth efficiency for 𝐸𝑖𝑗 = 2.5𝐶𝑗 - E14a predator 𝑖 eating prey of group 𝑗 (for all groups except salps and larvaceans) Average growth efficiency for 𝐶𝑗 - E14b 𝐸𝑖𝑗 = 0.25 predator 𝑖 eating prey of 𝐶𝑖 group 𝑗 (for salp and larvacean predators) Lower boundary condition for 𝑁𝑖 (𝑤𝑖 , 𝑡) = 𝑃𝑖 ∑ 𝑁𝑗 (𝑤𝑖 , 𝑡) m−3 E15 group 𝑖 𝑗≠𝑖 15 355 2.1.2 Incorporating functional traits We could parameterise many of the functional traits and key processes for different zooplankton groups because of the extensive experimental work carried out on zooplankton (see Table 2). 360 2.1.2.1 Body size and predator–prey mass ratio For each zooplankton group, we established body size ranges in wet weight using measurements and conversion equations from the literature (see Table 2). Prey preference for predators of a certain size was represented using the ratio of a predator’s body size to its preferred prey body size (PPMR), rather than specifying diet preferences (Table 2; Figure 4). 365 Across zooplankton taxa, PPMR increases with predator size, due to the non-isometric scaling of feeding-related apparatus with body size (Wirtz, 2012). We used the mechanistic formulation from Wirtz (2012) to calculate the PPMR range for each zooplankton group, except for salps and larvaceans. Wirtz (2012) links PPMR to a quantitative measure of the feeding mode: raptorial, active feeding is linked to a lower PPMR because predators eat prey 370 closer to their own size. By contrast, passive, suspension feeding yields a higher PPMR. The PPMR of an individual from group 𝑖 of size 𝑤 𝛽𝑖 (𝑤) is given by 𝛽𝑖 (𝑤) = (exp(0.02 ln(𝐷𝑤 )2 − 𝑚 + 1.832))3 , (E10) where 𝐷𝑤 is an individual of size 𝑤’s equivalent spherical diameter in 𝜇m, and 𝑚 is the quantitative measure of feeding mode; a large positive 𝑚-value is linked to a more 375 carnivorous feeding strategy, with a lower PPMR. By contrast, a large negative 𝑚-value is associated with a larger PPMR and a filter feeding strategy of prey capture. For filter feeders such as salps and larvaceans, we follow Wirtz’s formulation that PPMR increases with body size, but we assume that their filter mesh does not increase with body size. This means that, for salps and larvaceans, individuals will capture the same size prey no matter their own size 380 (see Bone et al., 2003). In keeping with previous studies, the PPMR for the fish communities was held constant at 100 across their size ranges (Blanchard et al., 2009; Hartvig et al., 2011; Rochet et al., 2012; Andersen et al., 2016b). The relatively small body sizes and high PPMRs of salps, larvaceans and omnivorous copepods means that these groups feed exclusively on phytoplankton, heterotrophic flagellates and 16 385 heterotrophic ciliates (Figure 4). Euphausiids also have high a PPMR range, but their larger size means that their largest size classes can also access smaller copepods, larvaceans and euphausiids, consistent with their omnivorous diet in the oceans (Schmidt and Atkinson, 2016). The low PPMR ranges of carnivorous copepods, chaetognaths and jellyfish, coupled with their larger body size means that these groups are almost totally carnivorous; we further 390 restricted their diets so that they do not feed on phytoplankton at all, which is consistent with most current understandings of their diets (Terazaki, 2000; Purcell and Arai, 2001). 2.1.2.2 Prey size selectivity The range of available prey sizes for an individual predator of body size 𝑤 from group 𝑖 is defined by a log-normal feeding kernel. This feeding kernel is centred on the predator’s 395 preferred predator–prey mass ratio (PPMR; 𝛽𝑖 (𝑤)), with a standard deviation (𝜎𝑖 ) given by the kernel width parameter for that predator’s group (Table 2): 2 𝛽𝑖 (𝑤)𝑤 ′ ′ 𝜙𝑖 (𝑤, 𝑤 ) = exp [− (ln ( )) / 2𝜎𝑖2 ] /(𝜎𝑖 √2𝜋). (E11) 𝑤 A wider feeding kernel means a predator can feed over a larger prey size range. For zooplankton, feeding kernel width is positively correlated with PPMR (Hansen et al., 1994; Fuchs and Franks, 2010; Kiørboe, 2016); filter feeders such as larvaceans or salps with a large 400 mean PPMR feed over a wider size range than carnivorous copepods or heterotopic flagellates with a smaller PPMR. We used the empirical model developed by Fuchs and Franks (2010) to link the feeding kernel width of each zooplankton group (𝜎𝑖 ), to that group’s mean PPMR (𝛽𝑖̅ ): 𝜎𝑖 = 0.05log10 (𝛽𝑖̅ ) + 0.33. (E12) 405 The feeding kernel widths for the fish communities were held constant at 1.3, in keeping with previous studies (Andersen et al., 2016b). 2.1.2.3 Growth conversion efficiency and carbon content Previous size spectrum models have focused on fish and used a currency solely of wet mass 410 rather than considering carbon, largely because fish groups vary little in their carbon content. However, as ZooMSS resolves nine zooplankton prey groups, and the carbon content of 17 zooplankton varies by more than an order of magnitude amongst common zooplankton groups (e.g., jellyfish = 0.5%, and copepods = 12%), we have included carbon content as a measure of food quality. As Straile (1997) found that mean growth conversion efficiency—as 415 a measure of prey carbon converted to predator carbon—was ~0.25 across a large range of zooplankton taxa, we consider that prey groups with a comparatively higher carbon content (e.g., copepods, euphausiids) contribute more to the wet-weight growth of a predator compared with lower carbon (more gelatinous) groups such as jellyfish (Spitz et al., 2010; Kiørboe, 2013; Mitra et al., 2014). In terms of wet weight, the growth conversion efficiency 420 of a predator of group 𝑖 (𝐸𝑖𝑗 ) feeding on prey from group j prey is: 𝐸𝑖𝑗 = 2.5𝐶𝑗 , (E14a) where 𝐶𝑗 is the carbon–wet weight ratio of group 𝑗 (Table 2). As gelatinous filter feeders—salps and larvaceans—are the fastest growing metazoans and have growth rates an order of magnitude faster than other zooplankton of the same size (Hirst 425 et al., 2003), we included this in the model by modifying E14a for these groups: 𝐶 𝐸𝑖𝑗 = 0.25 𝐶𝑗 , (E14b) 𝑖 so that the growth conversion efficiency of predators that are salps and larvaceans is also related to their carbon content 𝐶𝑖 . This means that these groups increase their wet weight by 5–10× more than other zooplankton of the same size, when consuming a given prey. For 430 example, according to E13a, a copepod consuming heterotrophic flagellates that have a carbon-wet weight ratio of 0.15 would have a gross growth efficiency of 0.375, but according to E13b a larvacean consuming the same flagellate would have a gross growth efficiency of 1.875. However, as a corresponding trade off, we made the assumption that for salps and larvaceans the additional wet weight in excess of what other zooplankton would gain is 435 converted to detritus, not into new biomass. This is based on their unique physiology: larvaceans grow and shed mucus houses multiple times each day (Bone, 1997), and salps produce large faecal pellets that contribute disproportionately to carbon flux compared with other zooplankton (Henschke et al., 2016). In the current implementation of the model, we 18 consider the faecal pellets and mucus houses to be detritus and lost from the system. Future 440 versions of ZooMSS will track the detritus through the food web. Figure 4 Overview of the size ranges of the functional groups (solid boxes) in the model, and their prey size ranges (dashed boxes). Note that boxes for predator and prey size are further apart for groups with higher PPMRs (e.g., euphausiids, larvaceans and salps). 19 Table 2 Trait values for the nine zooplankton and three fish groups. Group Min. Size, 𝑤𝑖 ̅𝑖 Max Size, 𝑊 𝑚- log10PPMR Feeding Kernel Carbon-Wet Length ESD log10(g)* Length ESD log10(g)* value range, 𝛽𝑖 (𝑤) Width, 𝜎𝑖 Weight Ratio, 𝑪𝒊 Hetero. Flagellates - 3.3×10-4 cma –10.7a - 7×10-3 cma –6.8a 1.52 0.2–0.72 0.363^ 0.154 Hetero. Ciliates - 1×10-3 cmb –9.3b - 1×10-2 cmb –6.3a 0.042 2.5–2.92 0.473^ 0.154 Larvaceans 8×10-3 cmc 1×10-2 cmc –6.4c 3×10-1 cmc 1×10-1 cmc –3.2c - 6.8–10.87 0.73^ 0.029 Omni. Cop. - 4×10-3 cmd –7.5d 2.8×10-1 cme 9×10-2 cme –3.5e -0.52 3.6–4.62 0.573^ 0.121 Carn. Cop. - 4×10-3 cmd –7.5d 6×10-1 cme 1.8×10-1 cme –2.5e 1.52 0.8–1.92 0.43^ 0.121 Euphausiids - 6×10-2 cmf –4.2f 6 cmg 1.5 cmg 0.2g -23,15 6.6–7.83,15 0.703^ 0.121 Chaetognaths 1×10-1 cmh 1.5×10-2 cmh –5.9h 4 cmh 6×10-1 cmh –0.9h 116 1.9–3.416 0.463^ 0.041 Salps 5×10-2 cmi 5×10-2 cmi –4.7i - 3.6 cmi 1.4i - 6.8–11.719 0.73^ 0.029 Jellyfish - 1.2×10-1 cmj –3j - 6 cmj 2j 0.731 2.7–4.71 0.523^ 0.0051 Small Fish - 1.2×10-1 cmk –3k - 6 cm 2 - 222 1.322 0.123 Medium Fish - 1.2×10-1 cmk –3k - 27 cm 4 - 222 1.322 0.123 Large Fish - 1.2×10-1 cmk –3k - 125 cm 6 - 222 1.322 0.123 * g wet weight calculated from ESD, assuming 1 gram = 1 cm3, ^ Feeding kernel widths were calculated with the empirical equation derived in (3), using mean log 10(PPMR) for this group. 445 a From Table 3 in (1), b From figure 1 in (5), c Minimum and maximum larvacean trunk lengths taken from (6) and (8) respectively, and converted to ESD and wet weight using equation derived in (7), d Carbon mass obtained from supplementary material in (10), converted to wet weight and ESD using carbon: wet weight ratio from (1) e Maximum omnivorous and carnivorous copepod lengths taken from (11) and converted to ESD and then wet weight using equation derived in (12), f Euphausiid embryo ESD from figure 2 in (13), g Maximum length taken from supplementary material in (3) and converted to ESD and wet weight using equation from (14), h Minimum and maximum ESD from supplementary material in (3), lengths derived using head width: body length ratio from (16) i Minimum and maximum salp length taken from (17) and converted to ESD and wet weight using equation derived in (18), maximum salp body size taken as geometric 450 mean of Salpida and Pyrosomatida from (17), after using equation in (18) j From supplementary material in (20), k From (21). 1. Hansen et al. (1997), 2. Wirtz (2012), 3. Fuchs and Franks (2010), 4. Menden-Deuer and Lessard (2000), 5. Taylor (1978), 6. López-Urrutia (2004), 7. Deibel (1998), 8. Hopcroft et al. (1998), 9. Sato et al. (2001), 11. Kiørboe & Hirst (2014), 11. Benedetti et al. (2016), 12. Acevedo et al. (2012), 13. Kawaguchi et al. (2011), 14. Meyer and Teschke (2016), 15. Schmidt and Atkinson (2016), 16. Pearre (1980), 17. Henschke et al. (2016), 18. Heron et al. (1988), 19. Bone et al. (2003), 20. Acuña et al. (2011), 21. Heneghan et al. (2016), 22. Andersen et al., (2016b), 23. Pauly 455 and Christensen (1995). 20 Table 3 Model parameter values. Symbol Definition Value Unit Source 𝛾 Coefficient of search 𝛾 = 640 g −α m−3 yr −1 1,2 rate 𝛼 Exponent of search rate 𝛼 = 0.8 - 1,2 𝑊𝑆𝑖 Body size at which ̅𝑖 0.01𝑊 g 3,4 senescence mortality begins for group 𝑖 𝛿 Coefficient of 1 yr −1 - senescence mortality 𝜌 Exponent of 0.3 - 3,4 senescence mortality 𝑤𝑝 Smallest phytoplankton 10-14.5 g 5 body size 𝑃𝑖 Relative abundance of Flagellates = 0.1 - - smallest size class of Ciliates = 0.1 group 𝑖 Larvaceans = 0.1 Omnivorous copepods = 0.04 Carnivorous copepods = 0.06 Euphausiids = 0.1 Chaetognaths = 0.1 Salps = 0.01 Jellyfish = 0.01 𝑄10 Temperature scaling 2 - 6,7 coefficient 𝐾ref Reference temperature 303.15 K - 1. Blanchard et al., (2009), 2. Peters (1983), 3. Hall et al. (2006), 4. Heneghan et al. (2016), 5. Brewin et al., 2015, 6. Eppley (1972), 7. Goldman and Carpenter (1974). 460 2.2 Numerical implementation of ZooMSS 2.2.1 Boundary conditions for the size-spectrum model We do not explicitly resolve reproduction for any zooplankton or fish group, which affects how the boundary conditions of the smallest size class of each functional group is parameterised. Although reproduction in fish is a relatively simple process of production of 465 eggs of roughly the same size (1 mg), and being sexually fertilised (Neuheimer et al., 2016), reproduction in the zooplankton is much more diverse, from alternating generations of sexual and asexual reproduction in salps and jellyfish (Fautin, 2002; Daponte et al., 2013), to the hermaphroditism of chaetognaths (Bone, 1991), and intersexuality in copepods (Gusmão and 21 McKinnon, 2009). Therefore, for simplicity, we assume constant recruitment for each 470 zooplankton and fish group, which keeps the abundance of their smallest size classes fixed at a proportion of the total abundance in those size classes from other groups. For each zooplankton group, the abundance of the smallest size class at time 𝑡, 𝑁𝑖 (𝑤𝑖 , 𝑡), was fixed with respect to the total abundance of the other groups 𝑗 in that size class: 𝑁𝑖 (𝑤𝑖 , 𝑡) = 𝑃𝑖 ∑ 𝑁𝑗 (𝑤𝑖 , 𝑡) , (E 15) 𝑗≠𝑖 475 where 𝑃𝑖 is the relative abundance of group 𝑖 in size class 𝑤𝑖 , with respect to the total abundance of all other groups. For the smallest zooplankton group, heterotrophic flagellates, the abundance of their smallest size class was fixed to be 10% of the abundance of the phytoplankton community in the same size class. In keeping with ZooMSv1, the smallest size class abundance for the total fish community was fixed at the total zooplankton abundance 480 in that size class, divided equally among the three fish groups. For all other zooplankton groups, except for salps and jellyfish, 𝑃𝑖 was fixed at 0.1. Salps and jellyfish are the two largest zooplankton groups, and 𝑃𝑖 was set at 0.01 to prevent these groups from dominating the biomass of the zooplankton community. 2.2.2 Running ZooMSS 485 Abundances of the zooplankton and fish groups were modelled with separate second order McKendrick-von Foerster equations, which we solved numerically using a second order semi- implicit upwind finite difference scheme (Press et al., 2007). For numerical implementation, we discretised the zooplankton and fish community size ranges into equal 0.1 log10 size intervals. The model is initialised with the same zooplankton community in each 5° grid cell 490 and integrated forward through time for 1000 years, with a half weekly time step. After experimenting with smaller and larger interval widths, we chose these values to discretise the weight and time intervals to ensure convergence in our numerical implementation, whilst minimising the time required to run the simulation. Large PPMRs and narrow feeding kernels have been shown to cause travelling waves across the spectrum through time (Datta et al., 495 2011; Heneghan et al., 2016). With the very high PPMRs of some of the zooplankton groups represented here, we found that our model settles into repeating cycles of travelling waves of abundance from small to large body sizes through time. Therefore, we average the last 500 22 years of abundance, biomass and growth rates to obtain the long-term average results presented here. 500 2.3 Model Assessment 2.3.1 Assessing the emergent global distribution of total zooplankton biomass We use two different estimates of the global distribution of zooplankton biomass to assess our model estimates. The first is from Strömberg et al. (2009) who created a global map of zooplankton carbon biomass based on a theoretical trophic transfer model driven by net 505 primary production and tuned to 4843 estimates of carbon biomass. The second global map is derived from all available zooplankton biomass estimates in the COPEPOD database (Supplementary Figure A3, n=196,707). As these estimates were measured and collected in different ways, we used a generalized additive model with zooplankton biomass as the response and a suite of predictors to standardize for measurement type and collection 510 methods, and included satellite chlorophyll a and sea surface temperature as environmental predictors (see Supplementary Information Section A2 for more details). Because this statistical model was driven by data we called this approach the Empirical Model. Both approaches have strengths and weaknesses: the Strömberg Model is based on theory and relatively few observations, whereas the Empirical Model is based on many more 515 observations but had to account for differences in sampling methods and types of biomass measurements. 2.3.2 Assessing emergent global distributions of zooplankton abundance The global distribution of each zooplankton functional group emerges from the model based on their relative fitness in the environmental conditions of each grid cell. To assess the 520 emergent distributions of the individual functional groups against global patterns of empirical abundance, we used general linear models to model the relationship between in situ sample data, and environmental variables and sample equipment. Similar to our comparison of the model’s biomass with the model from Strömberg et al., (2009), we compared the emergent distributions of abundance from ZooMSS with the empirical distributions from the 525 generalised linear models, using Pearson’s correlation coefficient and the root mean square error. 23 To build the generalised linear models of zooplankton abundance, we compiled zooplankton taxonomic data from multiple sources to create a database of 640,184 observations of zooplankton matching 7 of the taxonomic groups in ZooMSS (see Figure A3 in Supplementary 530 Information, for maps of data distribution). These groups were: euphausiids, omnivorous copepods, carnivorous copepods, larvaceans, chaetognaths, salps and jellyfish. The database was compiled from the COPEPOD database (O’Brien 2014), the Continuous Plankton Recorder (Richardson et al. 2006), the Integrated Marine Observing System (Eriksen et al. 2019, IMOS 2019), the Jellyfish Database Initiative (Condon et al. 2014) and the Warreen Data (Baird et al 535 2011). Data that did not match the 7 groups were discarded and remaining data were quality controlled to remove samples with missing or incorrect metadata (e.g. time, date, GPS location). GEBCO bathymetry and satellite-derived sea surface temperature and chlorophyll a climatologies from MODIS-Aqua was matched to each sample. Using this database of seven zooplankton groups we compiled, we generated a linear model 540 in R v3.5.3 (R Core Team 2019) for each zooplankton group to quantify global patterns and compare with ZooMSS. The same initial model was used for all groups. Predictors used in the initial model were environmental variables sea surface temperature (SST; °C), chlorophyll a (mg m-3), bathymetry (m), day of year (days), sample depth (m), hour of the day (h), mesh size (µm) and device type (Net or CPR). Natural splines using the splines package were fit to 545 bathymetry and SST, and a harmonic function was fitted to cyclical variables of day of year, and hour of the day. Day of year was standardised to the Northern Hemisphere, so the middle of summer was the same day of year in both hemispheres. Non-significant predictors were removed, and model residual plots were inspected. After viewing residual plots, we log- transformed the response for each zooplankton group (abundance) and chlorophyll a – to 550 improve assumptions of homogeneity of variance and normality. We tested a model using the Gamma family with a log-link, however the models did not converge. We do not have the original sample volumes for all the datasets, so we were unable to use counts as the response and volume as an offset in the model. Figures of the final statistical models for each of the seven functional groups can be found in the Supplementary Material A3 (Supplementary 555 Figures A6-A12). 24 2.3.3 Assessing modelled zooplankton growth rates Energy transfer through the zooplankton from phytoplankton to fish is driven by growth rates of the different zooplankton functional groups. To assess whether the model produced 560 realistic rates of energy transfer through the zooplankton groups, we compared modelled growth rates with empirical rates from Hirst et al. (2003) and Kiørboe and Hirst (2014). Evaluating growth rates in the model also allows us to assess whether a simple Type 1 functional response is a reasonable parameterisation of zooplankton feeding rates at this temporal and spatial scale. 565 2.4 Sensitivity analysis We conducted two sensitivity analyses. First, to assess the robustness of the model to our parameterisation of the abundances of the smallest zooplankton size classes, we conducted a sensitivity analysis of the total biomass of the zooplankton community when 𝑃𝑖 was varied by ±50% from chosen values. Second, to assess the sensitivity of fish community biomass to 570 the composition of the zooplankton, we calculated the change in total fish biomass when each of the zooplankton functional groups (excluding flagellates and ciliates) were removed individually. As part of this second analysis, we also assessed the sensitivity of total fish biomass when zooplankton were represented as a single functional group, with a fixed PPMR of 1000 (Kiørboe, 2008; Ward et al., 2012), a carbon content of 10% (Pauly and Christensen, 575 1995) and a feeding kernel width of 1.3 (Andersen et al., 2016b). 3 Results 3.1 Model assessment 25 3.1.1 Global zooplankton biomass 580 Figure 5 Comparison of the global zooplankton carbon biomass (mg C m −3) from the size- spectrum model, with maps derived from observations. a) ZooMSS total zooplankton carbon biomass map (mg C m−3); b) Global annual average (1998–2005) zooplankton carbon biomass map (mg C m−3), reproduced from Strömberg et al., (2009) (called the Strömberg Model); c) Correlation plot between a) and b) with Pearson’s correlation (r) and root mean square error 585 (RMSE) reported; d) Global zooplankton carbon biomass map based on an empirical statistical model using 196,707 observations from the COPEPOD database (called the Empirical Model, see Supplementary Information Section A2 for more information); e) Correlation plot between a) and d) with Pearson’s correlation (r) and root mean square error (RMSE) reported. In c) and e), the black dashed line is the one-to-one line, and the blue solid line is the line of 590 best fit. There is reasonable agreement in the spatial pattern between the global zooplankton biomass map from ZooMSS and that from the Strömberg and Empirical Models, with correlation coefficients of 0.85 and 0.9 respectively and RMSEs of 0.88 and 0.86 respectively (Figure 5). ZooMSS captures the global-scale patterns of total zooplankton biomass, showing 595 lowest levels of biomass in the oligotrophic ocean gyres, and highest levels in upwelling regions and coastal shelves. The distribution of global biomass in ZooMSS and the Strömberg and Empirical Models is similar to that of satellite-derived chlorophyll a at this scale. When comparing values from the ZooMSS model global biomass with the Strömberg model, two data groups are evident (Figure 5c): one above the line of best fit, which mainly contains 26 600 grid cells in ZooMSS in coastal and upwelling where zooplankton biomass is high; and one below the line of best fit, which has grid cells from oligotrophic regions where zooplankton biomass is low. At this scale, this variation between the biomass in ZooMSS and the Strömberg Model could be caused by the difference in underlying environmental drivers from the Strömberg study (which uses net primary production to drive zooplankton biomass) and our 605 study (which uses chlorophyll a concentration) or from differences in the structure of the two models. In the comparison between the biomass from ZooMSS and the Empirical Model, there were no similar groups obvious (Figure 5 e). In terms of absolute zooplankton biomass, the biomass ranges in ZooMSS (3–300 mg C m−3) is generally higher than from the Strömberg Model (0.5–30 mg C m−3, Figure 5 b). The 610 Strömberg zooplankton biomass is calculated from a simple trophic model, where the total zooplankton biomass depends on the trophic transfer efficiency between phytoplankton and zooplankton and tuned to a relatively small dataset on zooplankton biomass (n=4,843). It appears that the zooplankton biomass estimates from ZooMSS are shifted vertically from the Strömberg Model, evident in the relative position of the 1:1 line (Figure 5c). By contrast, the 615 biomass range in ZooMSS (3–300 mg C m−3) is fairly close to that from the Empirical Model (50–300 mg C m−3) based on a much larger global dataset of biomass measurements (n=196,707) than that of Strömberg et al. (2009). On the other hand, because the minimum values in ZooMSS are much lower (~3 mg C m−3) than from the Empirical Model (~70 mg C m−3), ZooMSS might underestimate zooplankton biomass when zooplankton biomass is low – 620 i.e. in low chlorophyll a areas. Overall, the large correlations and RMSEs between the pattern of global zooplankton biomass in ZooMSS and the Strömberg and Empirical Models, and the fact that ZooMSS absolute biomass range falls within estimates from the two models based on data, is evidence that the model produces reasonable estimates of the distribution of global zooplankton biomass. 625 3.1.2 Global zooplankton functional group abundances The numerical model shows good agreement with the spatial pattern of abundance from the statistical models, indicating that ZooMSS is capturing the broad patterns of zooplankton community composition (Figure 6). Correlations were positive for six of the seven groups. The 27 best performing groups were jellyfish, omnivorous copepods and euphausiids (r = 0.91, 0.82 630 and 0.79 respectively), followed by chaetognaths, carnivorous copepods and larvaceans (r = 0.39, 0.29 and 0.25 respectively). Salps were the only group that showed a negative correlation with ZooMSS. This could be due to limitations in the statistical model rather than ZooMSS, as we had fewer data than the other groups and the available data was clustered around Australia and the Southern Ocean (see Figure A3 in Supplementary Information). 635 For all groups, ZooMSS predicted that abundance increases with chlorophyll a concentration, with highest abundances in coastal and upwelling grid cells and lowest abundances in oligotrophic open ocean grid cells where chlorophyll a is lowest. The effect of temperature was secondary to chlorophyll a in ZooMSS, and so emergent abundance distributions agree well with empirical distributions when empirical abundance is strongly driven by chlorophyll 640 a, such as for omnivorous copepods, euphausiids and jellyfish. However, when a group’s empirical distribution was strongly driven by temperature the model’s distribution was poorly correlated, or even negatively correlated in the case of salps. Although the model was generally able to capture patterns in the empirical distributions of abundance, the reported absolute abundances from the model were much greater than those 645 from the statistical models. With the exception of jellyfish, the root mean square errors (RMSE) between abundances from ZooMSS and the statistical maps were all greater than 2, indicating that the average difference between the absolute abundances from ZooMSS and the statistical models was more than 2 orders of magnitude. For jellyfish, the RMSE was 0.53, indicating over a three-fold difference between the model jellyfish abundances and those 650 reported from the statistical model. 28 29 Figure 6 Assessing emergent abundance distributions from the size-spectrum model. The abundance distributions from the statistical models (left column, log10(ind m-3); see Supplementary Information Section A3 for more information), the emergent zooplankton 655 distributions from ZooMSS (centre column, log10(ind m-3)), and correlation plots (right column) between the two, with Pearson’s correlation (r) and root mean square error (RMSE) reported, for: a-c) Larvaceans, d-f) Omnivorous Copepods, g-i) Carnivorous Copepods, j-l) Euphausiids, m-o) Chaetognaths, p-r) Salps and s-u) Jellyfish. In each row, the dots in the correlation plot represent the abundance of that row’s zooplankton group, in a 5 grid square, 660 from the statistical model (x-axis) and ZooMSS (y-axis). The blue line in each correlation plot is the line of best fit, from a linear model of empirical versus ZooMSS abundances. 3.1.3 Zooplankton growth Figure 7 Zooplankton mass-specific growth rates (d−1) as a function of body size (g C) for a) Microzooplankton (heterotrophic flagellates and heterotrophic ciliates), b) Filter Feeders 665 (larvaceans and salps), c) Omnivores (omnivorous copepods and euphausiids) and d) Carnivores (jellyfish, carnivorous copepods and chaetognaths). For each group, dots are 30 empirical mass-specific growth rates (d−1) from Hirst et al. (2003) and Kiørboe and Hirst (2013), and lines represent modelled mass-specific growth rates. 670 The model’s mass-specific growth rates agree well with observed rates from Hirst et al. (1998) and Kiørboe and Hirst (2013) (Figure 7). For all groups, modelled growth rates fell within the range of observed growth rates, and, with the exception of large flagellates (Figure 7 a), model growth rates do not exceed the empirical data for any functional group. However, modelled growth rates for ciliates, larvaceans and jellyfish are slower than most of the empirical 675 observations for those groups (Figure 7 a, b, d). These results indicate that overall, our model is not overestimating growth rates and resulting energy flow through the zooplankton to higher trophic levels. There were some differences between the scaling of growth and body size between the model and empirical observations; with the exception of salps and larvaceans, modelled 680 growth rates decline less with body mass than observations do. Moreover, changes in modelled growth rates with body size are linear for the model zooplankton groups, with the exception of the carnivorous zooplankton, which are more variable (Figure 7 d). 3.2 Sensitivity analyses 685 3.2.1 Model sensitivity to boundary conditions for each group Table 4 Median and maximum changes (%) in total zooplankton biomass when the relative abundance of the smallest size class of each zooplankton group (𝑃𝑖 ) is varied by +50% or −50%, one at a time. −50% 𝑃𝑖 +50% 𝑃𝑖  Total Zooplankton Median Maximum Median Maximum Biomass Flagellates −5% +13% +4% +14% Ciliates −2% +3% +<1% +13% Larvaceans −4% +5% +3% +39% Omnivorous −5% +38% +12% +27% Copepods Carnivorous +5% +31% +3% +15% Copepods Chaetognaths −<1% +15% +<1% +5% Euphausiids +1% +9% +<1% +11% Salps +1% +15% +<1% +15% Jellyfish −<1% +<1% +<1% +9% 31 690 Zooplankton model biomass was robust to changes in the boundary conditions of all zooplankton groups. Total zooplankton biomass was altered by less than the change in 𝑃𝑖 for 695 all of the zooplankton groups, with total zooplankton biomass varying ≤15 % for flagellates, ciliates, chaetognaths, salps, euphausiids and jellyfish and by <40 % for all groups (Table 4). The median change in total zooplankton biomass was ≤ ±5 % when 𝑃𝑖 was changed by ±50 % for all groups, with the exception of omnivorous copepods, which saw a median increase of 12% in total zooplankton biomass when that group’s 𝑃𝑖 was increased by 50 %. 700 3.2.2 Sensitivity of fish biomass to different zooplankton groups Figure 8 Distribution of the change in fish biomass when a) Larvaceans, b) Omnivorous copepods, c) Carnivorous Copepods, d) Chaetognaths, e) Euphausiids, f) Salps and g) Jellyfish 705 are omitted from the model one at a time, compared to total fish biomass when all 32 zooplankton groups are present, across the 5° regions of the global ocean. The dashed red line shows the median change in fish biomass, and the solid black line at 0 % shows no change in total fish biomass. 710 Across the global ocean, total fish biomass was sensitive to the presence of different zooplankton groups (Figure 8). Removing omnivorous copepods or euphausiids caused declines of up to 80 % in total fish biomass, with median global declines of 59.3 % for omnivorous copepods and 46.5 % for euphausiids (Figure 8 b, e). When larvaceans were removed fish biomass declined by up to 50 %, however, the median global decline in the 715 absence of larvaceans was limited to 4.5 %, as the regions that had declines in fish biomass were balanced by the 50% of grid squares that had little change or a slight increase in fish biomass (Figure 8 a). Removing carnivorous copepods saw the greatest median increase in global fish biomass of 18.5%, with some grid squares recording increases in fish biomass of over 70% (Figure 8 c). The removal of chaetognaths, salps and jellyfish caused changes in total 720 fish biomass across the global ocean of over ±25 % (Figure 8 d, f, g). However, for each of these three groups there was little change in total fish biomass for at least 35 % of the global ocean when they were removed, and the global median change in fish biomass was at most ±1 %. 4 Discussion 725 Here, we described the Zooplankton Model of Size Spectra version 2 (ZooMSS), the first functional size-spectrum model of the global marine ecosystem that includes the body size, carbon content and size-based feeding characteristics of nine major zooplankton groups. The model requires only sea surface temperature and chlorophyll a concentration as inputs, and the spatial variation of the zooplankton community biomass and structure emerges based on 730 the functional traits of the zooplankton groups. Our model strongly suggests that zooplankton play a key role in regulating the biomass of fish in the ocean, and demonstrates how various zooplankton groups support energy transfer to higher trophic levels. The ability of ZooMSS to reproduce empirically-derived spatial patterns of global zooplankton biomass and abundance as well as growth rates gives confidence that the model provides a reasonable framework for 735 further exploration of the role of zooplankton in the ocean. 33 4.1 Importance of resolving zooplankton groups for higher trophic levels Accounting for functional diversity in marine models is important to better understand the unique role of zooplankton in mediating energy from phytoplankton to fish across the global 740 ocean (Mitra et al., 2014). In current global marine models focussed on higher trophic levels, the functional diversity of zooplankton is typically ignored, with zooplankton usually grouped together with phytoplankton as a single resource for small fish (e.g., Blanchard et al., 2009, 2012; Law et al., 2009; Datta et al., 2010), or taken as external inputs from nutrient- phytoplankton-zooplankton-detritus models (e.g., Maury, 2010; Christensen et al., 2015; 745 Petrik et al., 2019). In marine models focussed on lower trophic levels, more progress has been made in resolving the functional traits of different plankton groups (e.g., Ward et al., 2012, 2014; Prowe et al., 2018). However, these models do not resolve higher trophic levels such as fish, because they are not built to assess the unique and dynamic role of zooplankton in mediating energy from phytoplankton to higher trophic levels. By resolving the functional 750 traits of nine major zooplankton groups, our model has the unique ability to explore the impact of shifts in the phytoplankton and zooplankton community on higher trophic levels across global environmental gradients. Our results from the sensitivity analysis where we removed each zooplankton group one at a time indicate that not all zooplankton groups are equal, and that their role in mediating 755 energy from phytoplankton to higher trophic levels depends on their body size, carbon content and size-based feeding behaviour. When omnivorous zooplankton (omnivorous copepods and euphausiids) were removed from the model, it reduced fish biomass by up to 80%. Omnivorous zooplankton are thus important for promoting fish biomass, a likely consequence of their high PPMR, wide range of body sizes and high carbon content, which 760 means they efficiently transfer carbon to higher trophic levels. By contrast, when carnivorous zooplankton (carnivorous copepods, chaetognaths and jellyfish) were removed from the model, it increased fish biomass by up to 20%. Carnivorous zooplankton thus reduce fish biomass, a likely effect of their very low PPMR and diet of other zooplankton, which cause longer and more inefficient food chains where these groups dominate. When filter-feeding 765 zooplankton (larvaceans and salps) were removed from the model, there were small declines by up to 5% in fish biomass. 34 Larvaceans and salps share similar body size ranges with omnivorous copepods and euphausiids, but have much larger predator–prey mass ratios, meaning they can also access picophytoplankton for food. These groups are a direct energy pathway from 770 picophytoplankton to higher trophic levels, and so could be important for the productivity and efficiency of oligotrophic waters where picophytoplankton dominate (Brewin et al., 2010). On the other hand, these groups are gelatinous, meaning they are a less nutritious energy source for higher trophic levels and is part of the reason why median global fish biomass did not decline as much when these groups were removed from our model, 775 compared to omnivorous copepods and euphausiids. Increasing functional complexity in models is important for improving our understanding of key ecosystem processes in the marine environment. For instance, trait diversity in phytoplankton has a key role in ocean biogeochemistry (Fuhrman 2009), productivity (Chen et al., 2019) and carbon export (Guidi et al., 2012), and there are a growing number of models 780 that demonstrate the importance of trait diversity in structuring the phytoplankton community across environmental gradients (Acevedo-Trejos, 2015, 2018; Basu and Mackey, 2018; Dutkiewicz et al., 2019). Similarly, we have shown here that resolving zooplankton functional complexity is important to better understand how ecosystem efficiency is affected by shifts in phytoplankton community structure from oligotrophic ocean gyres to eutrophic 785 coastal and upwelling systems (Hansen et al., 1994; Boyce et al., 2015). 4.2 Model evaluation To evaluate the performance of ZooMSS, we focussed on comparing emergent properties of the model’s zooplankton, with empirical observations of total zooplankton biomass across the global ocean, as well as the distributions of abundance for seven of the zooplankton 790 groups, and empirical growth rates for all zooplankton groups. Across the three approaches the model performed well, reproducing the global distribution of zooplankton biomass and most of the zooplankton functional group abundances, as well as producing growth rates that were within the range of empirical observation for all zooplankton groups. The global distribution of total zooplankton biomass from ZooMSS was strongly correlated 795 with the empirical estimate from Strömberg et al., (2009), our estimate from the COPEPOD 35 database (see Supplementary Information) and also agrees with other model predictions of the distribution of zooplankton and fish biomass (Jennings et al., 2008; Ward et al., 2012; Harfoot et al., 2014; Petrik et al., 2019). What is more, total zooplankton biomass from ZooMSS falls within the range of the two empirical estimates. 800 The emergent distributions of abundance for six of the seven zooplankton groups that were assessed had a positive correlation with distributions that were derived from data. However, there was a difference in the ranges of absolute zooplankton biomass and abundance between the model and the empirical estimates, with ZooMSS projecting higher abundances than the empirical estimates. This would be explained first by the fact that parameters for 805 ZooMSS were selected from the literature, and were not tuned to estimates of abundance and biomass from the sample data. Moreover, there is a significant difference between sample data used to derive the empirical distributions, and how biomass and abundance are calculated from the ZooMSS. The empirical models were constructed using sample data taken with dozens of different gears and mesh sizes, over the past 60 years. However, there are 810 many mesh sizes and gear types not suitable to capture smaller zooplankton, biasing sampled abundances and biomass to larger zooplankton (Everett et al., 2017). At the same time, large zooplankton can actively avoid sample nets and avoid detection, thus leading to samples that underestimate the number of large zooplankton (Richardson et al., 2006). This means that the samples—and therefore the distributions predicted by the empirical models—are only 815 capturing a fraction of the actual total zooplankton. By contrast, biomass and abundance from ZooMSS was obtained by summing all zooplankton. This means we are counting 100% of the zooplankton from the size-spectrum model, while the sampled distributions from the empirical models only capture some fraction of total zooplankton. Mass-specific growth rates from ZooMSS agreed well with observed rates from Hirst et al. 820 (2003) and Kiørboe and Hirst (2013). However, there were two main discrepancies. The first is that for zooplankton groups other than salps and larvaceans, modelled growth rates decline less with body size compared to observations. The reason for this discrepancy could be in the difference between how growth rates are calculated for salps and larvaceans and the rest of the zooplankton groups; the prey size range for salps and larvaceans is fixed no matter the 825 predator body size, meaning that growth rates for these groups scale with body size based on 36 search volume only, which has a mass-specific scaling that is less than zero. By contrast, the other zooplankton groups have prey size ranges that shift with increasing predator size—prey get larger as predators get larger—and this could be offsetting the lower search rate for larger predators, as well as the lower abundance of larger prey. 830 The second discrepancy was that for carnivorous groups, the model had non-linear variation in growth rates with body size that was not evident in the empirical data. This variation in specific growth rates with body size for carnivores is caused by the variability in the density of their prey. Microzooplankton, filter feeders and omnivorous zooplankton all primarily eat phytoplankton, which is modelled as a static resource across time and body size, so their 835 growth rates are largely driven by predator search volume and size selectivity, not variability in their prey. However, carnivorous zooplankton only eat other zooplankton, which have varying abundance with time and body size, so variation in their growth rates is caused not only by search volume and size selectivity, but also by changes in their prey through time. 4.3 Caveats and next steps 840 4.3.1 Spatial and temporal dynamics The spatial and temporal resolution of ZooMSS means that the model is designed to explore steady-state conditions of the marine ecosystem across global environmental gradients. For this reason, we focused on traits such as body size, predator–prey mass ratio and carbon content because, at this spatial and temporal scale, these traits are important for how energy 845 moves through the marine food web (Hansen et al., 1994; Jennings et al., 2001; Woodward et al., 2005; Andersen et al., 2016a; McConville et al., 2017). However, this also means that we have made simplifying assumptions about dynamic processes such as reproduction and movement. This is in contrast to other fish-focussed global models, which although they resolve zooplankton simplistically, do resolve explicit reproduction (Blanchard et al., 2011; 850 Carozza et al., 2016; Petrik et al., 2019) or movement (Cheung et al., 2011; Maury, 2010). By fixing abundances of the smallest size classes for all functional groups, we assume constant recruitment, which means we are not explicitly resolving reproduction. This could have a stabilising effect on the dynamics of the model, and our sensitivity analysis of total zooplankton biomass to changes in the boundary conditions of each zooplankton group 855 indicates that our results are reasonably robust to the implicit parameterisation of 37 reproduction as constant recruitment. Zooplankton reproductive strategies are extremely complex, however a first step toward resolving reproduction in the zooplankton would be to follow the formulation of the community size-spectrum model framework (Andersen et al., 2016b), where reproductive complexity is ignored. For each functional group recruitment 860 could then be represented as a fraction of assimilated energy flowing from mature size classes to the smallest size class for that group. Finally, the spatial resolution of our model (5) is coarser than other global marine ecosystem models, which tend to use a 1 resolution (Tittensor et al., 2018). We use a 5 resolution because it allowed us to run the model quickly, while still being able to assess the emergent 865 patterns of zooplankton biomass and abundance across the global ocean. Similar to most global fish-focused models run with a 1 or less resolution (but see Maury, 2010), each 5 grid cell is run independently, which means we do not take into account the movement of plankton or fish between adjacent grid squares. This could introduce a possible bias due to movement of plankton by currents in some regions or apex predators by active transport, 870 which in the real world could cause discrepancies between plankton and fish abundance in certain areas. However, given the coarse resolution of our model output we expect this to have a minor impact, since most zooplankton are short-lived (less than a few weeks) and are unlikely to move outside of the large grid cells over their lifespan. Further, we have not included fish species in the model, but only used three broad, size-based groups to represent 875 all epipelagic fish, which precludes a good understanding of how fish might move among grid cells during large-scale movements and seasonal migrations. The 5 cells also mean that many fish would remain within a single grid cell during their life, particularly for the small and medium fish communities. In future models, passive movement by currents and active movement due to behaviour of fish could be implemented using existing advection and 880 diffusion algorithms (Maury, 2010; Castle et al., 2011; Watson et al., 2014). 4.3.2 Phytoplankton production and dynamics Similar to existing global marine ecosystem models (Blanchard et al., 2009; Maury, 2010; Christensen et al., 2015; Jennings and Collingridge, 2015; Carozza et al., 2016; Petrik et al., 2019), ZooMSS does not incorporate phytoplankton dynamics, instead representing primary 885 producers as a static abundance spectrum with slope, intercept and maximum size 38 determined by annual average chlorophyll a concentration. Satellite chlorophyll a represents in situ phytoplankton biomass, and so in a sense implicitly incorporates the processes of nutrient cycling and predation on the phytoplankton. Moreover, since ZooMSS was able to reproduce global zooplankton biomass and growth rates in the range of empirical estimates, 890 we believe a static resource for zooplankton that does not explicitly resolve phytoplankton dynamics and nutrient cycling in each grid cell is a reasonable compromise between realism and model complexity at this stage. Looking forward, resolving phytoplankton dynamics and the feedbacks between phytoplankton and zooplankton is critical to better resolving zooplankton in marine 895 ecosystem models. In any region, primary production sets the limits to growth for higher trophic levels. Limited primary production also induces competition amongst zooplankton competing for a finite resource. This competition would have implications for the stability of different zooplankton groups from oligotrophic to eutrophic regions that are not resolved in the current model. In any region primary production depends on nutrient supply. In coastal 900 and upwelling regions, a high supply of new nutrients from ocean mixing sustain productive systems, while in the open ocean nutrient recycling from zooplankton back to phytoplankton through the microbial loop plays a key role in nutrient supply for primary production (Azam et al., 1983). Finally, seasonal cycles of boom and bust in the phytoplankton is a major driver of variation of zooplankton reproduction, as many zooplankton groups time their 905 reproduction to coincide with phytoplankton blooms (Falkowski et al., 1988; Atkinson, 1996). These linkages and feedbacks between nutrient cycles, phytoplankton dynamics and zooplankton need to be considered if we wish to move beyond the current static representation of the plankton in this model. Processes of nutrient uptake, growth and mortality are strongly size-structured in the phytoplankton, and these size-based 910 relationships have been used to resolve the dynamics of the phytoplankton over large spatial scales (Follows et al., 2007; Fuchs and Franks, 2010; Ward et al., 2012, 2014; Cuesta et al., 2018). Future model developments could focus on resolving the size-structured functional traits of phytoplankton nutrient uptake, growth and mortality, or alternatively to shift from using unlimited phytoplankton biomass to limited primary production to fuel higher trophic 915 levels. 39 4.3.3 Fitting the model to data Emergent patterns of biomass and abundance from the model generally agree with empirical patterns, and the ZooMSS model estimate of absolute global zooplankton biomass falls within the range of two empirical estimates. However, ZooMSS’s absolute abundances are up to 920 several orders of magnitude higher than empirical estimates. Part of this discrepancy is because the sample data used to fit empirical abundance estimates are obtained with methods that are not able to capture all zooplankton (Richardson et al., 2006), whereas reported numbers from the model include all individuals across the entire size range of zooplankton. However, the magnitude of the difference between the model’s abundance 925 estimates and the empirical estimates highlights that the model’s parameters were not tuned to sample data, rather the functional traits incorporated into ZooMSS, were parameterized using values from the literature. This is in contrast to other global marine ecosystem models, which calibrate their parameter values so that their reported values are in the range of empirical estimates (Christensen et al., 2015; Carozza et al., 2017). For zooplankton, that kind 930 of calibration presents unique challenges given the homogeneity of sampling methods used to obtain empirical estimates of biomass and abundance. However, empirical maps of zooplankton abundance fit with sample data have been useful for constraining and assessing previous ecosystem model estimates (Everett et al., 2017), which means that the statistical models developed here could also be useful to constrain zooplankton abundance in future 935 iterations of ZooMSS. In the meantime, reported absolute abundances from this model should be read with these uncertainties in mind. 4.3.4 Temperature Effects of sea surface temperature were incorporated in the model as a multiplier on growth, mortality and diffusion terms, with the same temperature scaling (Q10) used for all functional 940 groups. We made this assumption because, despite studies elucidating temperature scaling for different zooplankton species (Hansen et al., 1997; Kiørboe and Hirst, 2014), we are unaware of any meta-analyses for temperature dependence of different processes for the broader taxonomic groups we use here. However, this means that changes in chlorophyll a, and the resulting shifts in the size structure of the phytoplankton community was the main 945 driver of biomass and abundance for all zooplankton groups, and so the model does not capture the unique temperature response of different functional groups. This is highlighted 40 in the contrast between the empirically derived abundance distributions, and the emergent distributions from ZooMSS; the emergent distributions of abundance from the model are driven by chlorophyll a, with more zooplankton abundance in high chlorophyll a waters and 950 less in oligotrophic waters. Emergent abundance distributions agree well with empirical distributions when empirical abundance is strongly driven by chlorophyll a. However, when a group’s empirical distribution was more strongly affected by temperature the model’s distribution was poorly correlated. This highlights the pressing need for greater synthesis of the extensive experimental zooplankton data on the rate of life processes at different 955 temperatures to inform modelling efforts. 4.4 Concluding remarks ZooMSS is the first functional size-spectrum model of the global marine ecosystem to resolve the functional traits of multiple zooplankton functional groups. Using body size ranges, size- based feeding behaviour and carbon content of nine of the most abundant zooplankton 960 groups, ZooMSS is capable of replicating global patterns in total zooplankton biomass and abundance, and resolving growth rates for different zooplankton groups. By showing how unique zooplankton groups have different impacts on fish biomass across the ocean, ZooMSS provides new insights into the role of different zooplankton in supporting higher trophic level biomass, and is a first step towards better resolution of zooplankton in marine ecosystem 965 models using the functional size-spectrum framework. In an era of environmental change, and given the pivotal role zooplankton play in food webs, improving the representation of zooplankton in ecosystem models will be critical to predict fisheries productivity and carbon cycling in the future ocean. 970 Code availability ZooMSS was written in R version 3.5. Scripts to run the model, as well as global forcings for sea surface temperature and chlorophyll a, are available for download at https://github.com/MathMarEcol/ZoopModelSizeSpectra. 975 Acknowledgements and Funding 41 RFH was funded by an Australian Government Research Training Program Scholarship the Spanish Ministry of Science, Innovation and Universities, through the Acciones de Programacion Conjunta Internacional (PCIN-2017-115). PS was supported by an Australian Government Research Training Program Scholarship. JDE was funded by Australian Research 980 Council Discovery Projects DP150102656 and DP190102293. Author Contribution Statement RFH, AJR, JDE and JLB defined the research question for the study. RFH built the model and ran the simulations, and analysed the results, with assistance from AJR, JDE, JLB and PS. SDB, ME, KT provided data for the model assessment. AJR and JDE assisted with model 985 assessment. RFH drafted the manuscript, with assistance and feedback from AJR, JDE, JLB, PS, IS, SDB, ME and KT. 42 References Acuña, J.L., López-Urrutia, Á., Colin, S., 2011. Faking giants: The evolution of high prey 990 clearance rates in jellyfishes. Science 333, 1627–1629. https://doi.org/10.1126/science.1205134 Agawin, N.S.R., Duarte, C.M., Agustí, S., 2000. Nutrient and temperature control of the contribution of picoplankton to phytoplankton biomass and production. Limnol. Oceanogr. 45, 591–600. https://doi.org/10.4319/lo.2000.45.3.0591 995 Andersen, K.H., Jacobsen, N.S., Farnsworth, K.D., 2016b. The theoretical foundations for size- spectrum models of fish communities. Can. J. Fish. Aquat. Sci. 73, 575–588. https://doi.org/10.1139/cjfas-2015-0230 Andersen, K.H., Prowe, A.E.F., Berge, T., Palacz, A., Olsson, K., Jacobsen, N.S., Andersen, K.H., Lindemann, C., Traving, S.J., Wadhwa, N., Visser, A.W., Gonçalves, R.J., Hylander, S., 1000 Martens, E.A., Neuheimer, A.B., Heuschele, J., Kiørboe, T., Sainmont, J., Hartvig, M., 2016a. Characteristic Sizes of Life in the Oceans, from Bacteria to Whales. Ann. Rev. Mar. Sci. 8, 217–241. https://doi.org/10.1146/annurev-marine-122414-034144 Acevedo, F. de, Dias, J.D., Braghin, L. de S.M., Bonecker, C.C., 2012. Length-weight regressions of the microcrustacean species from a tropical floodplain. Acta Limnol. Bras. 24, 01–11. 1005 https://doi.org/10.1590/s2179-975x2012005000021 Acevedo-Trejos, E., Brandt, G., Bruggeman, J., Merico, A., 2015. Mechanisms shaping size structure and functional diversity of phytoplankton communities in the ocean. Sci. Rep. 5, 17–20. https://doi.org/10.1038/srep08918 Acevedo-Trejos, E., Marañón, E., Merico, A., 2018. Phytoplankton size diversity and 1010 ecosystem function relationships across oceanic regions. Proc. R. Soc. B 285, 20180621. http://dx.dor.org/10.1098/rspb.2019.0621 Azam, F., Graf, J.S., 1983. The Ecological Role of Water-Column Microbes in the Sea. Mar. Ecol. Prog. Ser. 10, 257–263. Baird, M.E., Suthers, I.M., 2007. A size-resolved pelagic ecosystem model. Ecol. Modell. 203, 1015 185–203. https://doi.org/10.1016/j.ecolmodel.2006.11.025 Baird, M. E., J. D. Everett, and I. M. Suthers. 2011. Analysis of southeast Australian zooplankton observations of 1938–42 using synoptic oceanographic conditions. Deep Sea Research Part II: Topical Studies in Oceanography 58(5):699–711 Barnes, C., Irigoien, X., De Oliveira, J.A.A., Maxwell, D., Jennings, S., 2011. Predicting marine 1020 phytoplankton community size structure from empirical relationships with remotely sensed variables. J. Plankton Res. 33, 13–24. https://doi.org/10.1093/plankt/fbq088 Barton, A.D., Pershing, A.J., Litchman, E., Record, N.R., Edwards, K.F., Finkel, Z. V., Kiørboe, T., Ward, B.A., 2013. The biogeography of marine plankton traits. Ecol. Lett. 16, 522–534. https://doi.org/10.1111/ele.12063 1025 Basu, S., Mackey, K.R.M., 2018. Phytoplankton as Key Mediators of the Biological Carbon Pump : Their Responses to a Changing Climate. Sustainability. 10, 869, https://doi.org/10.3390/su10030869 Benedetti, F., Gasparini, S., Ayata, S.D., 2016. Identifying copepod functional groups from species functional traits. J. Plankton Res. 38, 159–166. 1030 https://doi.org/10.1093/plankt/fbv096 Blanchard, J.L., Andersen, K.H., Scott, F., Hintzen, N.T., Piet, G., Jennings, S., 2014. Evaluating targets and trade-offs among fisheries and conservation objectives using a multispecies 43 size-spectrum model. J. Appl. Ecol. 51, 612–622. https://doi.org/10.1111/1365- 2664.12238 1035 Blanchard, J.L., Heneghan, R.F., Everett, J.D., Trebilco, R., Richardson, A.J., 2017. From Bacteria to Whales: Using Functional Size Spectra to Model Marine Ecosystems. Trends Ecol. Evol. 32, 174–186. https://doi.org/10.1016/j.tree.2016.12.003 Blanchard, J.L., Jennings, S., Holmes, R., Harle, J., Merino, G., Allen, J.I., Holt, J., Dulvy, N.K., Barange, M., 2012. Potential consequences of climate change for primary production 1040 and fish production in large marine ecosystems. Philos. Trans. R. Soc. B Biol. Sci. 367, 2979–2989. https://doi.org/10.1098/rstb.2012.0231 Blanchard, J.L., Jennings, S., Law, R., Castle, M.D., McCloghrie, P., Rochet, M.-J., Benoît, E., 2009. How does abundance scale with body size in coupled size-structured food webs? J. Anim. Ecol. 78, 270–280. https://doi.org/10.1111/j.1365-2656.2008.01466.x 1045 Blanchard, J.L., Law, R., Castle, M.D., Jennings, S., 2011. Coupled energy pathways and the resilience of size-structured food webs. Theor. Ecol. 4, 289–300. https://doi.org/10.1007/s12080-010-0078-9 Bone, Q. (1991). The biology of chaetognaths. Oxford University Press, Oxford Bone, Q. (1997). The biology of pelagic tunicates. Oxford University Press, Oxford 1050 Bone, Q., Carre, C., Chang, P. 2003. Tunicate feeding filters. J. Mar. Biol. Ass. U.K. 83, 907-919. https://doi.org/10.1017/S002531540300804Xh Boudreau, P.R., and Dickie, L.M., 1992. Biomass spectra of aquatic ecosystems in relation to fisheries yield. Can. J. Fish. Aquat. Sci. 49, 1528-1538. https://doi.org/10.1139/f92-169 Boyce, D.G., Frank, K.T., Leggett, W.C., 2015. From mice to elephants: overturning the ‘one 1055 size fits all’ paradigm in marine plankton food chains. Ecol. Lett. 18, 504–515. https://doi.org/10.1111/ele.12434 Brewin, R.J.W., Sathyendranath, S., Hirata, T., Lavender, S.J., Barciela, R.M., Hardman- Mountford, N.J., 2010. A three-component model of phytoplankton size class for the Atlantic Ocean. Ecol. Modell. 221, 1472–1483. 1060 https://doi.org/10.1016/j.ecolmodel.2010.02.014 Brun, P., Payne, M.R., Kiørboe, T., 2016. Trait biogeography of marine copepods – an analysis across scales. Ecol. Lett. 19, 1403–1413. https://doi.org/10.1111/ele.12688 Castle, M.D., Blanchard, J. L., Jennings, S., 2011. Predicted Effects of Behavioural Movement and Passive Transport on Individual Growth and Community Size Structure in Marine 1065 Ecosystems. Advances in Ecological Research 45, 41–66. Carozza, D.A., Bianchi, D., Galbraith, E.D., 2016. The ecological module of BOATS-1.0: a bioenergetically constrained model of marine upper trophic levels suitable for studies of fisheries and ocean biogeochemistry. Geosci. Model Dev. 9, 1545–1565. https://doi.org/10.5194/gmd-9-1545-2016 1070 Carozza, D.A., Bianchi, D., Galbraith, E.D., 2017. Formulation , General Features and Global Calibration of a Bioenergetically-Constrained Fishery Model. PLos ONE. 12, e0169763 https://doi.org/10.1371/journal.pone.0169763 Chen, B., Smith, S.L., Wirtz, K.W., 2019. Effect of phytoplankton size diversity on primary productivity in the North Pacific: trait distributions under environmental variability. Ecol. 1075 Lett. 22, 56-66. https://doi.org/10.1111/ele.13167 Christensen, V., Coll, M., Buszowski, J., Cheung, W.W.L., Frölicher, T., Steenbeek, J., Stock, C.A., Watson, R.A., Walters, C.J., 2015. The global ocean is an ecosystem: Simulating marine life and fisheries. Glob. Ecol. Biogeogr. 24, 507–517. https://doi.org/10.1111/geb.12281 44 1080 Condon, R., Lucas, C., Duarte, C. M., Pitt, K. (2014) Jellyfish Database Initiative: Global records on gelatinous zooplankton for the past 200 years, collected from global sources and literature (Trophic BATS project). Biological and Chemical Oceanography Data Management Office (BCO-DMO). Dataset version 2014-08-28 doi:10.1575/1912/7191 Cuesta, J.A., Delius, G.W., Law, R., 2018. Sheldon spectrum and the plankton paradox: two 1085 sides of the same coin—a trait-based plankton size-spectrum model. J. Math. Biol. 76, 67–96. https://doi.org/10.1007/s00285-017-1132-7 Dam, H.G., Baumann, H., 2017. Climate Change, Zooplankton and Fisheries. Clim. Chang. Impacts Fish. Aquac. II, 851–874. https://doi.org/10.1002/9781119154051.ch25 Datta, S., Delius, G.W., Law, R., 2010. A jump-growth model for predator-prey dynamics: 1090 Derivation and application to marine ecosystems. Bull. Math. Biol. 72, 1361–1382. https://doi.org/10.1007/s11538-009-9496-5 Datta, S., Delius, G.W., Law, R., Plank, M.J., 2011. A stability analysis of the power-law steady state of marine size spectra. J. Math. Biol. 63, 779–799. https://doi.org/10.1007/s00285- 010-0387-z 1095 Daponte, M.C., Palmieri, M.A., Casareto, B.E., Esnal, G.B., 2013. Reproduction and population structure of the salp Iasis zonaria (Pallas, 1774) in the southwestern Atlantic Ocean (34°30′ to 39°30′S) during three successive winters (1999-2001). J. Plankton Res. 35, 813–830. https://doi.org/10.1093/plankt/fbt034 Deibel, D. (1998) Feeding and metabolism of Appendicularia, in: Bone, Q. (Ed.), The Biology 1100 of Pelagic Tunicates. Oxford University Press, New York, pp. 139–149. Dutkiewicz, S., Cermeno, P., Jahn, O., Follows, M.J., Hickman, A.E., Taniguchi, A.A., Ward, B.A., 2019. Dimensions of Marine Phytoplankton Diversity. Biogeosciences 1–46. https://doi.org/10.5194/bg-2019-311 Edwards, K.F., Litchman, E. Klausmeier, C.A., 2013. Functional traits explain phytoplankton 1105 community structure and seasonal dynamics in a marine ecosystem. Ecology Letters 16, 56–63. https://doi.org/10.1111/ele.12012 Eppley, R.W., 1972. Temperature and Phytoplankton Growth in the Sea Fish. Bull. 70 1063– 1085. Eriksen RS, Bonham P, Davies CH, Coman FE, Edgar S, McEnnulty FR, McLeod D, Miller MJ, 1110 Rochester W, Slotwinski A, Tonks ML, Uribe-Palomino J, Richardson AJ (2019) Australia’s Long-term Plankton Observations: The Integrated Marine Observing System National Reference Station Network. Frontiers in Marine Science 6: 161. 17 pp. Everett, J.D., Baird, M.E., Buchanan, P., Bulman, C., Davies, C., Downie, R., Griffiths, C., Heneghan, R., Kloser, R.J., Laiolo, L., Lara-Lopez, A., Lozano-Montes, H., Matear, R.J., 1115 McEnnulty, F., Robson, B., Rochester, W., Skerratt, J., Smith, J.A., Strzelecki, J., Suthers, I.M., Swadling, K.M., van Ruth, P., Richardson, A.J., 2017. Modeling what we sample and sampling what we model: Challenges for zooplankton model assessment. Front. Mar. Sci. 4, 1–19. https://doi.org/10.3389/fmars.2017.00077 Fautin, D.G., 2002. Reproduction of Cnidaria. Can. J. Zool. 80, 1735–1754. 1120 https://doi.org/10.1139/z02-133 Follows, M.J., Dutkiewicz, S., Grant, S., Chisholm, S.W., 2007. Emergent biogeography of microbial communities in a model ocean. Science 315, 1843–1846. https://doi.org/10.1126/science.1138544 Friedland, K.D., Stock, C., Drinkwater, K.F., Link, J.S., Leaf, R.T., Shank, B. V., Rose, J.M., 1125 Pilskaln, C.H., Fogarty, M.J., 2012. Pathways between primary production and fisheries 45 yields of large marine ecosystems. PLoS One 7, e28945. https://doi.org/10.1371/journal.pone.0028945 Finkel, Z. V., Beardall, J., Flynn, K.J., Quigg, A., Rees, T.A. V., Raven, J.A., 2010. Phytoplankton in a changing world: Cell size and elemental stoichiometry. J. Plankton Res. 32, 119–137. 1130 https://doi.org/10.1093/plankt/fbp098 Fuchs, H.L., Franks, P.J.S., 2010. Plankton community properties determined by nutrients and size-selective feeding. Mar. Ecol. Prog. Ser. 413, 1–15. https://doi.org/10.3354/meps08716 Fuhrman J.A., 2009. Microbial community structure and its functional implications. Nature 1135 459, 193–199. https://doi:10.1038/nature08058 Gascuel, D., Pauly, D., 2009. EcoTroph : Modelling marine ecosystem functioning and impact of fishing. Ecol. Modell. 220, 2885–2898. https://doi.org/10.1016/j.ecolmodel.2009.07.031 Gentleman, W.C., Neuheimer, A.B., 2008. Functional responses and ecosystem dynamics : 1140 how clearance rates explain the influence of satiation , food-limitation and acclimation J. Plankton Res. 30, 1215-1231. https://doi.org/10.1093/plankt/fbn078 Goldman, J.C., Carpenter, E.J., 1974. A kinetic approach to the effect of temperature on algal growth. Limnol. Oceanog. 19, 756-766. https://doi.org/10.4319/lo.1974.19.5.0756 Gusmão, L.F.M., McKinnon, A.D., 2009. Sex ratios, intersexuality and sex change in copepods. 1145 J. Plankton Res. 31, 1101–1117. https://doi.org/10.1093/plankt/fbp059 Guidi L, Stemmann L, Jackson GA, Ibanez F, Claustre H, Legendre, L., Picheral, M., and Gorsky, G., 2009. Effects of phytoplankton community on production, size and export of large aggregates: a world-ocean analysis. Limnol. Oceanogr. 54, 1951–63. https://doi.org/10.4319/lo.2009.54.6.1951 1150 Hall, S.J., Collie, J.S., Duplisea, D.E., Jennings, S., Bravington, M., Link, J., 2006. A length-based multispecies model for evaluating community responses to fishing. Can. J. Fish. Aquat. Sci. 63, 1344–1359. https://doi.org/10.1139/f06-039 Hansen, B., Bjørnsen, P.K., Hansen, P.J., 1994. The size ratio between planktonic predators and their prey. Limnol. Oceanogr. 39, 395–403. 1155 https://doi.org/10.4319/lo.1994.39.2.0395 Hansen, P.J., Bjørnsen, P.K., Hansen, B.W., 1997. Zooplankton grazing and growth: Scaling within the 2-2,000-μm body size range. Limnol. Oceanogr. 42, 687–704. https://doi.org/10.4319/lo.1997.42.4.0687 Harfoot, M.B.J., Newbold, T., Tittensor, D.P., Emmott, S., Hutton, J., Lyutsarev, V., Smith, M.J., 1160 Scharlemann, J.P.W., Purves, D.W., 2014. Emergent Global Patterns of Ecosystem Structure and Function from a Mechanistic General Ecosystem Model. PLoS Biol. 12. https://doi.org/10.1371/journal.pbio.1001841 Hartvig, M., Andersen, K.H., Beyer, J.E., 2011. Food web framework for size-structured populations. J. Theor. Biol. 272, 113–122. https://doi.org/10.1016/j.jtbi.2010.12.006 1165 Heneghan, R.F., Everett, J.D., Blanchard, J.L., Richardson, A.J., 2016. Zooplankton Are Not Fish: Improving Zooplankton Realism in Size-Spectrum Models Mediates Energy Transfer in Food Webs. Front. Mar. Sci. 3, 1–15. https://doi.org/10.3389/fmars.2016.00201 Henschke, N., Everett, J.D., Richardson, A.J., Suthers, I.M., 2016. Rethinking the Role of Salps in the Ocean. Trends Ecol. Evol. 31, 720–733. 1170 https://doi.org/10.1016/j.tree.2016.06.007 46 Heron, A., McWilliam, P., Dal Pont, G., 1988. Length-weight relation in the salp Thalia democratica and potential of salps as a source of food. Mar. Ecol. Prog. Ser. 42, 125–132. https://doi.org/10.3354/meps042125 Hirata, T., Hardman-Mountford, N.J., Brewin, R.J.W., Aiken, J., Barlow, R., Suzuki, K., Isada, T., 1175 Howell, E., Hashioka, T., Noguchi-Aita, M., Yamanaka, Y., 2011. Synoptic relationships between surface Chlorophyll-a and diagnostic pigments specific to phytoplankton functional types. Biogeosciences 8, 311–327. https://doi.org/10.5194/bg-8-311-2011 Hirst, A.G., Roff, J.C., Lampitt, R.S., 2003. A synthesis of growth rates in marine epipelagic invertebrate zooplankton. Adv. Mar. Biol. 44, 1–142. https://doi.org/10.1016/S0065- 1180 2881(03)44002-9 Hirst, A. G. and Bunker, A. J. 2003. Growth of marine planktonic copepods: Global rates and patterns in relation to chlorophyll a, temperature, and body weight. Limnol. Oceanogr., 48, 1988–2010. Hopcroft, R.R., Roff, J.C. and Bouman, H.A., 1998. Zooplankton growth rates: the larvaceans 1185 Appendiculara, Fritallaria and Oikopleura in trophic waters. Journal of Plankton Research. 20, 539–555. https://doi.org/10.1093/plankt/20.3.539 Huete-Ortega, M., Cermeño, P., Calvo-Díaz, A., Marañón, E., 2012. Isometric size-scaling of metabolic rate and the size abundance distribution of phytoplankton. Proc. R. Soc. B Biol. Sci. 279, 1815–1823. https://doi.org/10.1098/rspb.2011.2257 1190 IMOS (2019) IMOS National Reference Station (NRS) - Zooplankton Abundance. https://catalogue- imos.aodn.org.au/geonetwork/srv/eng/metadata.show?uuid=c13451a9-7cfe-091c- e044-00144f7bc0f4 Jennings, S., Collingridge, K., 2015. Predicting consumer biomass, size-structure, production, 1195 catch potential, responses to fishing and associated uncertainties in the world’s marine ecosystems. PLoS One 10, e0133794. https://doi.org/10.1371/journal.pone.0133794 Jennings, S., Pinnegar, J.K., Polunin, N.V.C., Boon, T.W., 2001. Weak cross-species relationships between body size and trophic level belie powerful size-based trophic structuring in fish communities. J. Anim. Ecol. 70, 934–944. 1200 https://doi.org/10.1046/j.0021-8790.2001.00552.x Jennings, S., Mélin, F., Blanchard, J.L., Forster, R.M., Dulvy, N.K., Wilson, R.W., 2008. Global- scale predictions of community and ecosystem properties from simple ecological theory. Proc. R. Soc. B Biol. Sci. 275, 1375–1383. https://doi.org/10.1098/rspb.2008.0192 Kawaguchi, S., Kurihara, H., King, R., Hale, L., Berli, T., Robinson, J.P., Ishida, A., Wakita, M., 1205 Virtue, P., Nicol, S., Ishimatsu, A., 2011. Will krill fare well under Southern Ocean acidification? Biol. Lett. 7, 288–291. https://doi.org/10.1098/rsbl.2010.0777 Kiørboe, T. (2008). A mechanistic approach to plankton ecology. Princeton University Press, Princeton Kiørboe, T., 2011. How zooplankton feed: Mechanisms, traits and trade-offs. Biol. Rev. 86, 1210 311–339. https://doi.org/10.1111/j.1469-185X.2010.00148.x Kiørboe, T., 2013. Zooplankton body composition. Limnol. Oceanogr. 58, 1843–1850. https://doi.org/10.4319/lo.2013.58.5.1843 Kiørboe T., 2016. Foraging mode and prey size spectra of suspension-feeding copepods and other zooplankton. Marine Ecology Progress Series 558, 15–20. 1215 https://doi.org/10.3354/meps11877 47 Kiørboe, T., Hirst, A.G., 2014. Shifts in Mass Scaling of Respiration, Feeding, and Growth Rates across Life-Form Transitions in Marine Pelagic Organisms. Am. Nat. 183, E118–E130. https://doi.org/10.1086/675241 Law, R., Plank, M.J., James, A., Blanchard, J.L., 2009. Size-spectra dynamics from stochastic 1220 predation and growth of individuals. Ecology 90, 802–811. https://doi.org/10.1890/07- 1900.1 Litchman, E., Ohman, M.D., Kiørboe, T., 2013. Trait-based approaches to zooplankton communities. J. Plankton Res. 35, 473–484. https://doi.org/10.1093/plankt/fbt019 López-Urrutia, Á., Harris, R.P., Smith, T., 2004. Predation by calanoid copepods on the 1225 appendicularian Oikopleura dioica. Limnol. Oceanogr. 49, 303–307. https://doi.org/10.4319/lo.2004.49.1.0303 Marañón, E., Cermeño, P., Rodríguez, J., Zubkov, M. V., Harris, R.P., 2007. Scaling of phytoplankton photosynthesis and cell size in the ocean. Limnol. Oceanogr. 52, 2190– 2198. https://doi.org/10.4319/lo.2007.52.5.2190 1230 Marañón, E., 2015. Cell Size as a Key Determinant of Phytoplankton Metabolism and Community Structure. Ann. Rev. Mar. Sci. 7, 241–264. https://doi.org/10.1146/annurev- marine-010814-015955 Maury, O., 2010. An overview of APECOSM, a spatialized mass balanced “Apex Predators ECOSystem Model” to study physiologically structured tuna population dynamics in their 1235 ecosystem. Prog. Oceanogr. 84, 113–117. https://doi.org/10.1016/j.pocean.2009.09.013 McGill, B.J., Enquist, B.J., Weiher, E., Westoby, M., 2006. Rebuilding community ecology from functional traits. Trends Ecol. Evol. 21, 178–185. https://doi.org/https://doi.org/10.1016/j.tree.2006.02.002 1240 Mcconville, K., Atkinson, A., Fileman, E.S., Spicer, J.I., Hirst, A.G., 2017. Disentangling the counteracting effects of water content and carbon mass on zooplankton growth. J. Plankton Res. 39, 246–256. https://doi.org/10.1093/plankt/fbw094 Menden-Deuer, S., Lessard, E.J., 2000. Carbon to volume relationships for dinoflagellates, diatoms, and other protist plankton. Limnol. Oceanogr. 45, 569–579. 1245 https://doi.org/10.4319/lo.2000.45.3.0569 Meyer, B., Teschke, M., 2016. Physiology of Euphausia superba, in: Volker, S., Biology and Ecology of Antarctic Krill. 1st ed, Springer International Publishing, pp. 145–175. Mitra, A., Castellani, C., Gentleman, W.C., Jónasdóttir, S.H., Flynn, K.J., Bode, A., Halsband, C., Kuhn, P., Licandro, P., Agersted, M.D., Calbet, A., Lindeque, P.K., Koppelmann, R., Møller, 1250 E.F., Gislason, A., Nielsen, T.G., St. John, M., 2014. Bridging the gap between marine biogeochemical and fisheries sciences; configuring the zooplankton link. Prog. Oceanogr. 129, 176–199. https://doi.org/https://doi.org/10.1016/j.pocean.2014.04.025 Mitra, A., Davis, C., 2010. Defining the “to” in end-to-end models. Prog. Oceanogr. 84, 39–42. https://doi.org/10.1016/j.pocean.2009.09.004 1255 Moreno-Ostos, E., Blanco, J.M., Agustí, S., Lubián, L.M., Rodríguez, V., Palomino, R.L., Llabrés, M., Rodríguez, J., 2015. Phytoplankton biovolume is independent from the slope of the size-spectrum in the oligotrophic Atlantic Ocean. J. Mar. Syst. 152, 42–50. https://doi.org/10.1016/j.jmarsys.2015.07.008 Neuheimer, A. B., Hartvig, M., Heuschele, J., Hylander, S., Kiørboe, T., Olsson, K. H., et al., 1260 2015. Adult and offspring size in the ocean over 17 orders of magnitude follows two life history strategies. Ecology 96, 3303–3311. https://doi.org/10.1890/14-2491.1 48 O’Brien, T.D. 2014. COPEPOD: The Global Plankton Database. An overview of the 2014 database contents, processing methods, and access interface. U.S. Dep. Commerce, NOAA Tech. Memo., NMFS-F/ST-37, 29 pp. 1265 Pauly, D., Christensen, V., 1995. Primary production required to sustain global fisheries. Nature 376, 255–257. Pearre, S., 1980. Feeding by Chaetognatha: The Relation of Prey Size to Predator Size in Several Species. Mar. Ecol. Prog. Ser. 3, 125–134. https://doi.org/10.3354/meps003125 1270 Peters, R.H., 1983. The Ecological Implications of Body Size. Cambridge University Press, Cambridge. Petrik, C.M., Stock, C.A., Andersen, K.H., van Denderen, P.D., Watson, J.R., 2019. Bottom-up drivers of global patterns of demersal, forage, and pelagic fishes. Prog. Oceanogr. 176, 102124. https://doi.org/https://doi.org/10.1016/j.pocean.2019.102124 1275 Polovina, J.J., Howell, E.A., Abecassis, M., 2008. Ocean’s least productive waters are expanding. Geophys. Res. Lett. 35, 2–6. https://doi.org/10.1029/2007GL031745 Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 2007. Numerical Recipes: The Art of Scientific Computing, Third Edition. Cambridge University Press, New York. Purcell J.E, and Arai M.N., 2001. Interactions of pelagic cnidarians and ctenophores with fish: 1280 a review. Hydrobiologia 451, 27–44. https://doi.org/10.1023/A:1011883905394 Prowe, A.E.F., Andersen, K.H., Kiørboe, T., Visser, A.W., Chiba, S., 2018. Biogeography of zooplankton feeding strategy. Limnol. Oceanogr. lno.11067. https://doi.org/10.1002/lno.11067 Richardson A.J., and Verheye, H.M., 1999. Growth rates of copepods in the southern Benguela 1285 upwelling system: the interplay between body size and food. Limnol. Oceanogr. 44, 382-392. Richardson, A. J., A. W. Walne, A. W. G. John, T. D. Jonas, J. A. Lindley, D. W. Sims, D. Stevens, and M. Witt. 2006. Using continuous plankton recorder data. Prog. Oceanog. 68, 27– 74. 1290 Rochet, M.J., Benoît, E., 2012. Fishing destabilizes the biomass flow in the marine size- spectrum. Proc. R. Soc. B Biol. Sci. 279, 284–292. https://doi.org/10.1098/rspb.2011.0893 Sato, R., Tanaka, Y., Ishimaru, T., 2003. Species-specific house productivity of appendicularians. Mar. Ecol. Prog. Ser. 259, 163–172. 1295 Sarmiento, J.L., Slater, R., Barber, R., Bopp, L., Doney, S.C., Hirst, A.C., Kleypas, J., Matear, R., Mikolajewicz, U., Monfray, P., Soldatov, V., Spall, S.A., Stouffer, R., 2004. Response of ocean ecosystems to climate warming. Global Biogeochem. Cycles 18. https://doi.org/10.1029/2003GB002134 Schmidt, K and Atkinson, A., 2016. Feeding And Food Processing Antarctic Krill (Euphausia 1300 Superba Dana), in: Volker, S., Biology and Ecology of Antarctic Krill. 1st ed, Springer International Publishing, pp. 175–224. Schnedler-Meyer, N.A., Mariani, P., Kiørboe, T., 2016. The global susceptibility of coastal forage fish to competition by large jellyfish. Proc. R. Soc. B Biol. Sci. 283. https://doi.org/10.1098/rspb.2016.1931 1305 Scott, F., Blanchard, J.L., Andersen, K.H., 2014. mizer: An R package for multispecies, trait- based and community size spectrum ecological modelling. Methods Ecol. Evol. 5, 1121– 1125. https://doi.org/10.1111/2041-210X.12256 49 Spitz, J., Mourocq, E., Schoen, V., Ridoux, V., 2010. Proximate composition and energy content of forage species from the Bay of Biscay: high- or low-quality food? ICES J. Mar. Sci. 67, 1310 909–915. https://doi.org/10.1093/icesjms/fsq008 Stock, C.A., Dunne, J.P., John, J.G., 2014. Drivers of trophic amplification of ocean productivity trends in a changing climate. Biogeosciences 11, 7125–7135. https://doi.org/10.5194/bg-11-7125-2014 Straile, D., 1997. Gross growth efficiencies of protozoan and metazoan zooplankton and their 1315 dependence on food concentration, predator-prey weight ratio, and taxonomic group. Limnol. Oceanogr. 42, 1375–1385. https://doi.org/10.4319/lo.1997.42.6.1375 Strömberg, K.H.P., Smyth, T.J., Allen, J.I., Pitois, S., O’Brien, T.D., 2009. Estimation of global zooplankton biomass from satellite ocean colour. J. Mar. Syst. 78, 18–27. https://doi.org/10.1016/j.jmarsys.2009.02.004 1320 Stuart-Smith, R.D., Bates, A.E., Lefcheck, J.S., Duffy, J.E., Baker, S.C., Thomson, R.J., Stuart- Smith, J.F., Hill, N.A., Kininmonth, S.J., Airoldi, L., Becerro, M.A., Campbell, S.J., Dawson, T.P., Navarrete, S.A., Soler, G.A., Strain, E.M.A., Willis, T.J., Edgar, G.J., 2013. Integrating abundance and functional traits reveals new global hotspots of fish diversity. Nature 501, 539–542. https://doi.org/10.1038/nature12529 1325 Taylor, W.D., 1978. Maximum growth rate, size and commonness in a community of bactivorous ciliates. Oecologia 36, 263–272. https://doi.org/10.1007/BF00348052 Terazaki, M., 2000. Feeding of Carnivorous Zooplankton, Chaetognaths in the Pacific, in: Dynamics and Characterization of marine Organic Matter. 1st ed. Terrapub, pp. 257–276. Tittensor, D.P. et al., 2018. A protocol for the intercomparison of marine fishery and 1330 ecosystem models: Fish-MIP v1.0. Geosci. Model Dev. 1421–1442. https://doi.org/https://www.geosci-model-dev-discuss.net/gmd-2017-209/ van Denderen, P.D., Lindegren, M., MacKenzie, B.R., Watson, R.A., Andersen, K.H., 2017. Global patterns in marine predatory fish. Nat. Ecol. Evol. 2. https://doi.org/10.1038/s41559-017-0388-z 1335 Watson, J. R., Stock, C. A., and Sarmiento, J. L., 2014. Exploring the role of movement in determining the global distribution of ma- rine biomass using a coupled hydrodynamic size-based ecosystem model, Prog. Oceanogr. 138, 521-532 doi:10.1016/j.pocean.2014.09.001, Ward, B.A., Dutkiewicz, S., Jahn, O., Follows, M.J., 2012. A size-structured food-web model 1340 for the global ocean. Limnol. Oceanogr. 57, 1877–1891. https://doi.org/10.4319/lo.2012.57.6.1877 Ward, B.E.N.A., Dutkiewicz, S., Follows, M.J., 2014. Modelling spatial and temporal patterns in size-structured marine plankton communities : top – down and bottom – up controls 36, 31–47. https://doi.org/10.1093/plankt/fbt097 1345 Wirtz, K.W., 2012. Who is eating whom? Morphology and feeding type determine the size relation between planktonic predators and their ideal prey. Mar. Ecol. Prog. Ser. 445, 1– 12. https://doi.org/10.3354/meps09502 Woodward, G., Ebenman, B., Emmerson, M., Montoya, J.M., Olesen, J.M., Valido, A., Warren, P.H., 2005. Body size in ecological networks. Trends Ecol. Evol. 20, 402–409. 1350 https://doi.org/10.1016/j.tree.2005.04.005 Woodworth-Jefcoats, P.A., Polovina, J.J., Dunne, J.P., Blanchard, J.L., 2013. Ecosystem size structure response to 21st century climate projection: Large fish abundance decreases in the central North Pacific and increases in the California Current. Glob. Chang. Biol. 19, 724–733. https://doi.org/10.1111/gcb.12076 50 1355 Zhang, C., Chen, Y., Ren, Y., 2015. Assessing uncertainty of a multispecies size-spectrum model resulting from process and observation errors. ICES J. Mar. Sci. 72, 2223–2233. https://doi.org/10.1093/icesjms/fsv086 Zhang, C., Chen, Y., Ren, Y., 2016. An evaluation of implementing long-term MSY in ecosystem-based fisheries management: Incorporating trophic interaction, bycatch and 1360 uncertainty. Fish. Res. 174, 179–189. https://doi.org/10.1016/j.fishres.2015.10.007 Zhang, L., Thygesen, U.H., Knudsen, K., Andersen, K.H., 2013. Trait diversity promotes stability of community dynamics. Theor. Ecol. 6, 57–69. https://doi.org/10.1007/s12080-012- 0160-6 Zhou, M., Carlotti, F., Zhu, Y., 2010. A size-spectrum zooplankton closure model for ecosystem 1365 modelling. J. Plankton Res. 32, 1147–1165. https://doi.org/10.1093/plankt/fbq054 51