Should we account for mesozooplankton reproduction and ontogenetic growth in biogeochemical modeling?

Mesozooplankton play a key role in marine ecosystems as they modulate the transfer of energy from phytoplankton to large marine organisms. In addition, they directly influence the oceanic cycles of carbon and nutrients through vertical migrations, fecal pellet production, respiration, and excretion. Mesozooplankton are mainly made up of metazoans, which undergo important size changes during their life cycle, resulting in significant variations in metabolic rates. However, most marine biogeochemical models represent mesozooplankton as protists-like organisms. Here, we study the potential caveats of this simplistic representation by using a chemostat-like zero-dimensional model with four different Nutrient-Phytoplankton-Zooplankton configurations in which the description of mesozooplankton ranges from protist-type organisms to using a size-based formulation including explicit reproduction and ontogenetic growth. We show that the size-based formulation strongly impacts mesozooplankton. First, it generates a delay of a few months in the response to an increase in food availability. Second, the increase in mesozooplankton biomass displays much larger temporal variations, in the form of successive cohorts, because of the dependency of the ingestion rate to body size. However, the size-based formulation does not affect smaller plankton or nutrient concentrations. A proper assessment of these top-down effects would require implementing our size-resolved approach in a 3-dimensional biogeochemical model. Furthermore, the bottom-up effects on higher trophic levels resulting from the significant changes in the temporal dynamics of mesozooplankton could be estimated in an end-to-end model coupling low and high trophic levels.


Introduction
Zooplankton plays a pivotal role in the marine ecosystems because they are at the interface between primary producers and higher trophic levels. They are as such responsible for a large fraction of the energy transfer to fish (Beaugrand et al. 2010;Carlotti and Poggiale 2010;Mitra et al. 2014). They are also key players in marine biogeochemical cycles as they convert a large fraction of the organic matter they consume to dissolved inorganic carbon and nutrient pools. Furthermore, they modulate the amount of organic matter that is exported to the interior of the ocean either through the sinking of fecal pellets and dead organisms (Henschke et al. 2016;Steinberg and Landry 2017) and through diel and seasonal vertical migration (Emerson et al. 1997;Packard and Gomez 2013;Aumont et al. 2018).
Assessment of trophic flows through zooplankton is difficult to perform due to insufficient data and methodological difficulties (Buitenhuis et al. 2006;Everett et al. 2017). Nevertheless, several estimates on the global scale have been proposed for the magnitude of feeding by different zooplankton size classes on primary autotrophic producers, i.e., phytoplankton. Microzooplankton feeding is found to consume from 50% up to almost 80% of primary production (PP) and exerts thus a major top-down control on phytoplankton (Calbet and Landry 2004;Schmoker et al. 2013). Mesozooplankton respiration has been estimated as 10% to about 30% of global PP, both phytoplankton and microzooplankton contributing importantly to the diet of mesozooplankton (Calbet 2001;Hernandez-Len and Ikeda 2005).
A first separation that should be made within zooplankton is the distinction between protozoan unicellular organisms which can be referred to as microzooplankton and metazoan pluricellular organisms which are in general larger animals included in mesozooplankton and macrozooplankton (Carlotti and Poggiale 2010;Mitra et al. 2014). Despite the fact that protozoan zooplankton can have quite complex lifecycles, they share many similarities with phytoplankton except for their heterotrophy. Their size changes are limited to about doubling or halving of their cell volume, which involves small changes in their metabolic rates during their life history. On the other hand, metazoan zooplankton undergo large changes in their size, up to several orders of magnitude, over their life cycle. As a consequence, their metabolic rates change strongly over their life cycle history. Furthermore, intraguild predation (cannibalism) is common in mesozooplankton, for instance among copepods and chaetognaths (Brodeur and Terazaki 1999;Basedow and Tande 2006;Ohman et al. 2008).
Mathematical models of marine biogeochemistry and plankton dynamics have gained a considerable interest over the last decades and are now a central plank in ecological, biogeochemical, and climate change studies. Most of the biogeochemical models used today find their roots in the NPZ formalism popularized by Fasham et al. (1990). Their complexity has strongly increased by adding multiple limiting nutrients and multiple phytoplankton and zooplankton functional groups or size classes, e.g., Moore et al. (2002); Le Queré et al. (2005); Follows et al. (2007); Ward et al. (2012); Aumont et al. (2015). However, most of the effort has focused on the description of the nutrient-phytoplankton link. In comparison, the representation of zooplankton has received much less attention and is very often restricted to employing a few size-classes, generally only 2 corresponding to micro-and mesozooplankton, both considered as protist-like organisms (Sailley et al. 2013). Indeed, for the vast majority of past and current biogeochemical models, zooplankton act as a closure term for nutrient and carbon fluxes.
Physiologically Structured Population Models (PSPMs) (Metz and Diekmann 1986;De Roos 1997) is a class of models that are designed to investigate the dynamics of population controlled by size-dependent processes. PSPMs are based on a description of relevant processes at the individual level such as feeding, development, reproduction and mortality which depend on the environment (for instance, the amount of food) and on the state of the organism itself. As size controls most, if not all, of these processes, including feeding interactions, PSMPs are often constructed as size-resolved models, (e.g., Persson et al. 1998;van Kooten et al. 2007). Indeed, this subclass of PSPMs has been used to improve our understanding of the effects of ontogenetic growth at the population level (De Roos 2018) as well as at the community level (De Roos 2020). In systems where size-structured predators feed on unstructured preys, cyclic cohort cycles may appear with periods depending on the competitive advantage between juveniles and adults De Roos et al. 2008). This type of model may also exhibit alternative stable states, at least when adult and juvenile stages feed on different resources (Guill 2009).
The complexity of size-structured PSPMs makes their application at the community level challenging, because they require a large set of species-specific parameters. As a more tractable alternative, aggregated stage-based approaches have been derived from PSPMs (De Roos et al. 2008;Guill 2009;Soudijn and De Roos 2017). However, even if they proved to be able to capture some of the properties of the PSPMs, this is not always the case and their predictions depend on the number of stages they represent (Guill 2009). Purely size-based models have been proposed to represent communities which reduce the complexity of PSPMs by assigning a single trait to differentiate between species. This trait can be size at maturation (Hartvig et al. 2011) or maximum size (Maury and Poggiale 2013). The structuring role of size at the community level is supported by marine observations that often show continuous and quasi-linear size spectra in logarithmic space (Sheldon et al. 1972;Quinones et al. 2003;Zhou et al. 2009;Thompson et al. 2013). The dynamics of the abundance distribution is generally described by the well-known McKendrick von Foerster equation that combines an advection term along the size dimension corresponding to growth and some sink terms due to predation and other mortality processes (Guiet et al. 2016). Using these different approaches, numerous models at the population or community levels have been developed to study zooplankton and their response to both biotic and abiotic changes in their environment (e.g., Carlotti et al. 2000;Flynn and Irigoien 2009;Zhou et al. 2010;Pinceel et al. 2016). Some of these models have been integrated in biogeochemical models (e.g., Fennel 2001;Baird and Suthers 2007;Zhou et al. 2010;Serra-Pompei et al. 2020). However, to our knowledge, the consequences of explicitly describing ontogenetic growth and reproduction of mesozooplankton have not been extensively explored, especially concerning biogeochemistry and plankton dynamics. Furthermore, no spatially resolved biogeochemical modeling study includes such a detailed description of mesozooplankton community.
In this paper, we propose to investigate the consequences of an explicit description of individual growth and reproduction of a mesozooplankton community on the temporal dynamics of a simple NPZ-type plankton model. To achieve this goal, we have developed a suite of four different model configurations in which the description of mesozooplankton goes from the dominant representation in biogeochemical models, as protist-like organisms, to a consistent size-based formulation with reproduction and ontogenetic growth. These different model configurations are used in a zero-dimensional chemostat setting in which the incoming nutrient levels are varied with time in order to assess the impact on the temporal dynamics of nutrient and plankton biomass.

Model description
Model structure This study is based on the simple biogeochemical model developed by Mayersohn (2018). The state variables of this model represent a resource or a group of organisms whose temporal evolution is governed by an ordinary differential equation (ODE). These state variables are expressed in terms of their phosphorus content ( molP L −1 ). The original version is a classical NPZ model with three major types of compartments: nutrient, phytoplankton, and zooplankton. The version developed for this study includes a subdivision of the zooplankton compartment to explicitly represent mesozooplankton size classes, and thus sexual reproduction and growth. Figure 1 summarizes the different (sub)compartments of the model and their interactions.
The model contains 4 main compartments (dotted round squares in Fig. 1), expressed in terms of biomass or concentration, and each comprising one or more subcompartments (simple round squares in Fig. 1). Phosphate is the only resource considered here. Phytoplankton P is divided into three size classes P i where i ∈ [0, 2] . Zooplankton is divided into two main compartments according to size: microzooplankton Z and mesozooplankton M . Microzooplankton is further divided into strictly heterotrophic protists U (i.e., strict microzooplankton) and N S 2 size classes representing juvenile multi-cellular zooplankton J i with i ∈ [0, The maximum size mesozooplankton compartment is also referred to as A max . The symbols are summarized in Table 1.
Our system represents a chemostat experiment in a zerodimensional (0-D) setting. This system can be considered similar to an idealized two-layer ocean model with instantaneous remineralization of organic matter under the mixed layer (Evans et al. 1985;Thingstad et al. 1996;Ward et al. 2012). Inorganic phosphorus (phosphate) is supplied to the system at a specified dilution rate (t) and concentration C(t), which can be either constant or vary with time, to mimic vertical mixing events with nutrient-rich deep waters. An equal outflow dilutes inorganic phosphorus and the organic living compartments in the system. Available Photosynthetically Active Radiation (PAR) I(t) is specified as an external environmental variable and is supposed homogeneous within the system. The external environmental conditions are represented by bold arrows in Fig. 1.

General system dynamics
The following set of equations governs the temporal evolution of the inorganic resource and phytoplankton groups: Fig. 1 Model structure. The dotted round squares are the major functional compartments. Plain round squares are the sub-compartments (i.e., Resources, Phytoplankton and Zooplankton). The thin arrows represent the flows between compartments or subcompartments. Bold arrows represent the external forcings. The large round white dotted square is the modeled system, which is similar to a chemostat. The small round blue dotted square is the external medium with which our system is mixed. "Repro." claims for "Reproduction" Growth rate of phytoplankton r P i depends on light and on nutrient availability ( Table 1). The only loss term for phytoplankton is predation by the different zooplankton groups. The inorganic resource R is taken up by phytoplankton and is resupplied by mortality and excretion from zooplankton.
The temporal evolution of the different zooplankton groups is computed according to the following set of equations: The formulation for the different terms used in this set of equations is detailed in Table 1 (Parts IV and V). All zooplankton groups feed on all three phytoplankton size classes.
In addition, mesozooplankton feed on strict microzooplankton and on juveniles. Grazing is a function of the amount of food available ( T Z for strict microzooplankton and juveniles, T M for mesozooplankton). It is parameterized according to a Holling type III function that improves the stability of the model (Gentleman and Neuheimer 2008). Feeding on

Amount of food available for mesozooplankton
Mean maximum ingestion rate of microzooplankton Mean maximum ingestion rate of mesozooplankton r z = g z S Z Ingestion rate of microzooplankton compartment z ∈ Z the different available preys is described by a proportional formulation which implies no-switching (Gentleman et al. 2003).
For each mature mesozooplankton A s , part of the assimilated food w is allocated to reproduction and is transferred to the juvenile sub-compartment J s . The remainder of the assimilated food is used for growth, resulting in a transfer between adjacent size classes at a rate v. The value of this parameter depends on the number of size classes and the assumed size distribution within each size class (see Supplementary Material Section 1). For the largest size class of mature mesozooplankton A max , no size growth is possible. Thus, the amount of assimilated food that would lead to an increase in size is channeled to the inorganic resource (we refer to this term as growth closure). Furthermore, we assume that the deepwater mass contains only mature life stages for mesozooplankton.
The overall dynamics of mesozooplankton M (Eq. 2) and microzooplankton Z (Eq. 1) compartments, on which we focus in this study, is obtained by summing the terms of the different sub-compartments: Size-based parameterization Phytoplankton is divided into three classes of increasing size. The nutrient half-saturation constant increases with size (e.g., Eppley et al. 1969;Edwards et al. 2012). Maximum growth rate is kept constant with size, despite observational studies that suggest a unimodal size scaling of maximum phytoplankton growth rate (e.g., Wirtz 2011;Maran et al. 2013). Zooplankton preference for phytoplankton is also set to vary with size: microzooplankton groups are assumed to preferentially feed on small phytoplankton, while mesozooplankton compartments are more likely to feed on large phytoplankton. Food preference is constant in each major zooplankton compartment (microzooplankton and mesozooplankton).
We consider a subdivision of the two metazoan zooplankton compartments into size classes of equal width in logarithmic space. The center of each size class is defined as in logarithmic space and is therefore constant. We define the length factor of the size class X s as 2s+1 2N s which goes from 0 to 1.
The maximum ingestion rate of the different zooplankton classes is set according to an allometric relationship proposed by Hansen et al. (1997). The half-saturation constant used in the grazing parameterization is supposed constant as observations suggest no significant variations with size (Hansen et al. 1997). The transition rate v between the different size-classes was computed by assuming that the slope of the biomass size-spectrum within each size-class is constant in a log-log space. It is set to -3 following the seminal study of Sheldon et al. (1972), which corresponds to an approximate constant biomass in logarithmically equal size intervals. The expressions for the transition rate and for the maximum ingestion rate are shown in Table 2 and their computation is detailed in the supplementary material (Sections 1 and 2). The size-dependent formulations used in our standard model configuration are listed in Table 2.

Alternative models
The model described above is hereafter referred to as the "standard model". It relies on an explicit description of the size distribution of metazoan zooplankton, including ontogenetic growth and reproduction. Furthermore, maximum ingestion rates vary with size following the allometric relationship proposed by Hansen et al. (1997). To investigate the effects of this representation of mesozooplankton ontogenetic growth and reproduction on the model behavior, three alternative model formulations were developed for comparison. The specificities of these models are presented in Fig. 2. Detailed equations are available in the supplementary materials (Section 3). The first alternative model is a "constant-rate" model, in which maximum ingestion rates are constant across all size classes of each major zooplankton compartment. The geometric mean of the maximum ingestion rates of microzooplankton and mesozooplankton, respectively, is assigned to the juvenile and mature size classes, i.e., g Z and g M (Table 2). By comparison with the standard model, this model configuration is designed to evaluate the effects of allometric scaling of maximum ingestion rates on the temporal dynamics of the system. The second alternative model is a "two-life-stage" model ( Fig. 2b) in which there is no size discretization of juveniles and mature organisms. Consequently, the representation of metazoan zooplankton is limited to two life stages: juveniles and mature organisms. In this model, adult mesozooplankton allocates part of their assimilated food to reproduction, which is represented as a flux of organic matter from mature mesozooplankton to the juvenile life stage. As a result of ontogenetic growth, juveniles are transferred to adult mesozooplankton at a rate v computed according to the formulation presented in Table 1. The value of v is here less than 1 in contrast to the first two model configurations that have a better resolved size distribution. Thus, food assimilation leads to an increase in the biomass of each metazoan compartments. As mature mesozooplankton is represented by a single compartment, the specific treatment of the largest size-class as applied in the standard and constant-rate models is not used here. Comparing the standard model with the two-life-stage model helps to investigate the effect of an increased size resolution on the temporal dynamics of the system. The last alternative model is a "no-life-stage" model ( Fig. 2a), in which metazoans are represented by a single compartment corresponding to mesozooplankton. Thus juvenile and mature organisms are supposed to have the same metabolic rates and the same predation behavior. There is no internal predation within the mesozooplankton compartment. This model is very similar to the classical approach adopted by most of the PFT-type biogeochemical models. In this case, the representation of both microzooplankton and mesozooplankton is similar and corresponds to a formalism used for protists whose typical mode of reproduction is based on cell division. Length factor for phytoplankton P i Minimal metazoan zooplankton body length l max Maximal metazoan zooplankton body length Transition rate between the mesozooplankton size-classes Mesozooplankton preference for phytoplankton P i g M Geometric mean of the maximum mesozooplankton ingestion rate g Z Geometric mean of the maximum microzooplankton ingestion rate Length factor for generic microzooplankton U Maximum ingestion rate of the zooplankton size-class X s

Numerical experiments
Each model configuration was run for 4 years using a simple Euler forward integration scheme with a time step of 1 hour. Light intensity, organism concentrations in the incoming waters and the dilution rate are kept constant throughout the simulations. Light intensity is set to a non-limiting value of 60 W m −2 . Dilution rate is maintained at 0.1 d −1 . Organisms concentrations in the inflow waters are set to a very small value to avoid any extinction of communities. The temporal dynamics of the different model configurations is investigated by periodically varying the incoming resource concentration according to a step function with a period T of 2 years. Incoming nutrients are varied between a null value and a non-limiting concentration of 1 molP L −1 .
The model parameters and their values are listed in Table 3. In the standard and constant-rate models, metazoans are split into 20 size-classes, equally distributed between juveniles and mature organisms. The reproduction parameter is set to 0.3 (Kooijman 2013). The size parameterization defines metazoans in the model as a community of individuals ranging from 10 to 4000 μ m with an mean egg to adult ratio of 1/20. All the other parameter values have been assigned to typical values found in the literature. Reference for each parameter is provided in Table 3. Some parameters of the standard model (denoted by a ⋆ in Table 3) have been varied by ±50% around their standard value to produce a series of sensitivity simulations discussed in the last section of the results. "two-life-stage" model configuration. Left red arrow represents maturation, right red arrow represents reproduction. In a. and b., light grey items denote the unchanged parts of the model presented in Fig. 1. c. Parameteriza-tions of the maximum ingestion rate of mesozooplankton. Plain lines with dots represent the standard model parameterization (allometric), and dashed lines with squares represent a configuration in which the maximum ingestion rate is constant within each main compartment. Vertical grey lines are delimiting the 20 size-classes An analytical study of the size spectra at equilibrium allows a comparison of the theoretical continuous spectra with the simulated discrete spectra. This comparison presented in the supplementary materials (Section 4) confirms that 1) the discretization into 20 size classes is sufficient enough to derive a good representation of the size growth dynamics according to the McKendrick Von Foerster equation and, 2) the formulation of the transition rate v is robust to a deviation from the assumption that size spectrum has a slope of -3 in log-log space.

Results
For all model simulations, we only analyze the last two years of the 4-yr simulations. Hereafter, we refer to the period following the switch from No-Nutrient Mode to High-Nutrient Mode as HNM and the period following the switch from High-Nutrient Mode to No-Nutrient Mode as 0NM.
In the No-Nutrient mode (0NM), the temporal dynamics across all model versions is similar (Fig. 3), with a sharp decrease for all compartments (resource and biomass) towards low values after ∼ 50 days. This response is dictated by the mixing term with inflow waters that contain no nutrients and virtually no phytoplankton and zooplankton. The delay in the response varies depending on the compartment, i.e., more rapid for nutrient concentrations and less so for phytoplankton and zooplankton biomasses.
Hereafter, we focus our analysis on the response to the High-Nutrient mode (HNM) and present the results of the standard model first, and then a comparison to the alternative models.

Standard model
The standard model output is represented in orange on Fig. 3. In HNM, about 70% (i.e. 0.7 mol P L −1 ) of the total phosphorus content of the system is not converted into biomass at equilibrium, so that phytoplankton growth is not limited by nutrient. Contrary to 0NM, the temporal dynamics varies a lot across compartments, with equilibrium values reached much earlier for nutrients, phytoplankton, and microzooplankton biomasses than for mesozooplankton. A short phytoplankton bloom is observed during the first 10 days. During this phase, up to 65% of the total phosphorus (P) content of the system is converted into phytoplankton biomass. Following the phytoplankton bloom, microzooplankton biomass also increases sharply reaching 55% of the total P content by the end of the phytoplankton bloom. The high grazing pressure exerted by microzooplankton leads to a rapid collapse of phytoplankton at t = 10 days. Microzooplankton biomass then decreases, being limited by prey abundance, and reaches an equilibrium ∼ 1 month after the switch to the high-nutrient mode. To summarize, both microzooplankton and phytoplankton experience a rapid and intense early bloom phase that is then followed by a rapid decline. Apart from the mesozooplankton compartment, a quasi-equilibrium is reached by the end of the first month in HNM.
For mesozooplankton, the temporal dynamics is markedly different from that of phyto-or microzooplankton. No outburst is observed and the temporal evolution is indeed much slower (Fig. 3). Consequently, equilibrium takes several months ( > 1 year) to be reached.
To get some insights into the mesozooplankton response to HNM, we explore the temporal dynamics of the 20 different size-classes of metazoans, for both juveniles (in the microzooplankton compartment) and mature individuals (in the mesozooplankton compartment). Figure 4 shows the evolution of the relative abundance of biomass in the different size-classes over time (i.e., the size-spectrum time evolution). During the first 9 days of HNM, biomass is concentrated in the mesozooplankton size classes, as the inflow waters contain only mature organisms.
Reproduction of these mature organisms produces juveniles and induces the propagation of a cohort across the size spectrum as individuals keep progressing in their sizeclasses. This first cohort reaches the maximum size-class after about two months in HNM. When this cohort enters the mature size-classes domain (red points on Fig. 4), it triggers reproduction that leads to the production of new juveniles and thus, generates a second cohort. During that second cohort, the relative biomass of the mature sizeclasses reaches lower values than the first cohort because of the progressive build-up of juveniles (Fig. 3). The start of a third and a fourth cohort is observed before the system reaches equilibrium.
The temporal evolution of mesozooplankton in HNM can be split into 5 phases (Fig. 5). Phase I is defined as the early response (fast growth) and is followed by a moderate growth phase (Phase II), when the first cohort reaches mesozooplankton compartment. Between t ∼ 30 and t ∼ 90 days , the first cohort gradually fades away before a second cohort starts, resulting in a mean mesozooplankton growth close to zero between t ∼ 50 and t ∼ 80 days. The time lapse between the two first cohorts defines Phase III. We characterize this phase by three indices: its duration DP3, the day of its onset TP3 and its amplitude AP3 (ie the change in growth during the phase). Further details about the calculation of those parameters are provided in supplementary materials (Section 6). For the standard model, DP3 equals 62 days, TP3 equals 61 days and AP3 equals 0.03 molP L −1 d −1 . Phase III is then followed by a second moderate growth phase (Phase IV) during which the growth speed gradually decreases until the mesozooplankton biomass reaches an equilibrium at t eq > 1 year (Phase V) (not shown on Fig. 5, discussed in the last section of the results). Figure 3 shows that all the model versions (standard and alternative model structures) seem to tend towards the same equilibrium state in HNM. Furthermore, the temporal evolution of nutrient concentrations, phytoplankton and microzooplankton biomasses is similar across all model versions. With the exception of mesozooplankton, on which we focus hereafter, differences between models are weak suggesting that the introduction of mesozooplankton reproduction and size variations in a simple 0-D biogeochemical model such as ours has no significant influence on either lower trophic levels or nutrient availability.

Comparison between the standard and alternative models
The temporal evolution of the mesozooplankton biomass is displayed on Fig. 5 for all model versions. Overall, the response is slower in the twenty-size-class models (standard and constant-rate models) and in the two-life-stage model than in the no-life-stage model. Consequently, the implementation of reproduction seems to delay the response of mesozooplankton to a sudden increase in food abundance resulting from a phytoplankton bloom and microzooplankton outburst. Using the phases defined for the standard model in the first section of the results, we summarize the differences in the temporal evolution of mesozooplankton between the different model versions (Fig. 5). First, growth is faster in Phase I (initial response) for the constant-rate than for the two-life-stage and standard models. In the no-life stage model, there is no Phase I. Second, Phase III is not simulated in the no-lifestage and two-life-stage models, so that Phases II and IV are combined into a single pre-equilibrium phase (moderate growth phase). In the constant-rate model, Phase III Fig. 6 Time evolution of the mesozooplankton sources terms, sinks terms and mean size. All sink and source terms are normalized by the current mesozooplankton biomass M . Refer to Eq. 2 for the definition of each term cannot be visually identified. Using the indices defined above, we identify however that its temporality is close to the standard model (TP3 = 55 days, DP3 = 55 days) but its amplitude is much lower (AP3 < 0.001 molP L −1 d −1 ). Third, Phase V (equilibrium ) is not reached before the end of HNM for the standard, constant-rate and two-life-stage models ( t eq > 1 year). In the no-life-stage model, Phase V is reached much faster ( t eq ∼ 100 days).

Drivers of mesozooplankton dynamics
To explain the differences between the model configurations, we investigate the different terms of the mesozooplankton equation (Fig. 6, Eq. 2) and demonstrate the role of explicit reproduction and of explicit variations in the mean size of mesozooplankton.
First, we detail the processes involved in the first 15 days (Phase I) in HNM. For all models including reproduction (two-life-stage, constant rate, standard), the input from maturation increases sharply to reach 0.5 d −1 during the second week, and then rapidly decreases to zero at t = 15 days. This rapid decrease marks the shift from Phase I to Phase II and happens just after the drop of microzooplankton biomass, which is explained by the sharp decline in food availability after all phytoplankton have been grazed by zooplankton. In the size-resolved models (standard and constant-rate), during the first week of HNM, net growth dominates mesozooplankton dynamics. Then Phase I actually starts at t = 7 days, and is mostly explained by the maturation input from the first wave of juveniles. There is no Phase I in the no-life stage model because it does not include reproduction.
Second, we detail the processes involved after t = 15 days in HNM and before the equilibrium is reached. The main difference between the different models is the occurrence of a Phase III in the size-resolved models which is completely absent in the two-life-stage model, and of very low amplitude in the constant-rate model. Phase III is explained by a cohort dynamics characterized by the successive propagation of waves of individuals along the size spectrum. Phase III is actually the period between the first two of these cohorts.
In the standard model, Phase III occurs between t ∼ 30 days and t ∼ 90 days and its onset (TP3) corresponds to a net growth rate minimum. Large size variations in the mean size of mature organisms during Phase III are shown on the bottom panel of Fig. 6 but also in Fig. 4. A higher mean size induces a lower mean maximum growth rate in the standard model as a consequence of the allometric scaling of growth rate. Thus, the cohorts materialize in the standard model dynamics by a strong variation of the net growth rate related to their ontogenetic growth. In the constant-rate-model, a slight reduction of the maturation term coincides with Phase III which occurs slightly earlier than in the standard model (TP3 = 55 days vs 61 days in the standard model), and has a slightly lower duration (DP3 = 55 days vs 62 days in the standard model). A low departure from zero of the growth closure term coincides with the end of phase III. As growth rate does not vary in the constant-rate model, net growth rate variations are lower than in the standard model. Hence, maturation and growth closure terms dominate the constantrate model dynamics during phase III, but the values of the dominant rates are much lower than in the standard model, which explains its lower amplitude. Thus, in the constantrate model, the cohorts materialize through subtle variations in maturation and closure terms.
Finally, we analyze the timing of equilibrium (Phase V) for each model version. Mesozooplankton equilibrium is reached once all terms of Eq. 2 are at equilibrium. Differences in this timing are related to the fact that the term reaching equilibrium last differs across models. In the nolife-stage model, Phase V is reached once mortality does not evolve anymore (at t ∼ 100 days). However, the key process that controls the temporal dynamics of this model is net growth rate, which reaches equilibrium very quickly. Mortality is just following the temporal evolution of the mesozooplankton biomass as it is a density-dependent closure term. In the two-life-stage and constant-rate models, net growth is reaching equilibrium as quickly as in the nolife-stage model. However, maturation continues to evolve until the very end of HNM which explains their very long equilibrium time.
As a conclusion on the global behavior of the 4 model configurations, slower temporal dynamics for all life-stage models (two-life-stage, constant-rate and standard) than for the no-life-stage model is induced by the introduction of reproduction, which converts part of the energy assimilated by mesozooplankton into juveniles. This results in lower net growth offset by input from maturation, which introduces a delay. The appearance of cohort dynamics in models with explicit size representation is more marked in the case of the standard model as a consequence of the allometric scaling of growth rate and explains the faster dynamics in this model during the first two months.

Sensitivity to the parameters
The biological parameters of our model configurations have been set to constants and chosen from previous studies (Table 3). Note that these parameter values are not well constrained by measurements and observations, either in situ or in the lab, and that they show large variations within and between species and groups.The modelling results presented above may depend on the choice of these parameter values. To explore the sensitivity of our results to the choice of parameter values, we performed a sensitivity analysis by systematically varying selected parameters by ±50%. The dynamics of juveniles and mature metazoan organisms show the greatest differences between our four model configuration simulations. Therefore, we focus here on key zooplankton parameters: gross growth efficiency ( ), maximum ingestion rates ( g Z and g M ), quadratic mortality rate of mature organisms (m) and the fraction of food allocated to reproduction (w). Figure 7 shows the changes induced by varying parameter values on five key metrics (mesozooplankton biomass at equilibrium and time to reach equilibrium, and the three metrics used to define Phase III, DP3, TP3, AP3).
The mesozooplankton biomass at equilibrium (EQ) in the 4 model configurations is sensitive to all parameters except that of reproduction. For all parameter values, EQ is similar across the three alternative model versions (constantly rate, two-life-stage and no-life-stage). In the standard model, EQ is slightly larger than in the other model versions for almost all parameter values. In the supplementary materials (Section 5), we demonstrate this result analytically by manipulating the set of equations of the different model setups.
The time to reach equilibrium (Q90, defined as the number of days to reach 90% of EQ) is also sensitive to all parameters except quadratic mortality. An increase in the maximum ingestion rate or in the growth efficiency leads to a decrease in Q90 for all model configurations because both result in an increase in the mesozooplankton growth rate. When the lowest value of these two parameters is prescribed, mesozooplankton go extinct in the standard, the constant-rate and in the two-lifestage model configurations as the achieved net growth rate remains lower than the dilution rate. But despite the sensitivity to parameter values, the sensitivity analysis confirms one of our main results: the no-life-stage model always simulates the fastest adjustment to equilibrium, while the standard model shows longer Q90 than any other model version, except when using the highest reproduction parameter value. Fig. 7 Sensitivity of the mesozooplankton dynamics to the four parameters that have been varied by ±50% . Four indicators of mesozooplankton dynamics are represented : EQ the biomass at equilib-rium, Q90 the number of days to attain 90% of the mesozooplankton biomass at equilibrium, DP3 the duration of the plateau (in phase III) and TP3 the day of its onset The cohort dynamics, i.e., the occurrence of Phase III, is simulated in the constant-rate and the standard model configurations for all parameter values, except for the lowest values of the maximum ingestion rate, growth efficiency, and quadratic mortality (for the constant-rate model). In these three cases, a Phase III can not be properly defined, either because mesozooplankton go extinct (the constant-rate model) or because its biomass remains too small (the standard model) to properly identify cohorts. The amplitude AP3 of the Phase III is at least 5 times lower in the constant-rate model than in the standard model for all parameter values. In the standard model, AP3 significantly increases as growth efficiency and microzooplankton maximum growth rate increase. This can be related to the fact that we are considering absolute amplitude, thus directly related to the speed of population growth, which increases as those two parameters increase. The difference in the time of onset TP3 and in the duration DP3 between the constant-rate and the standard model does not exceed 25% , which suggests that the first two cohorts present a similar temporality in both models. DP3 and TP3 systematically decreases when growth efficiency and microzooplankton maximum growth rate increases. Indeed, higher net growth accelerates the dynamics of the cohorts, which are therefore closer together (DP3 decreases) and shorter (TP3 decreases). TP3 and DP3 in the standard model are relatively insensitive to reproduction and quadratic mortality.

Limitations and prospects for improvement
Our study is based on a series of different model configurations that are used to explore the role of reproduction and ontogenetic growth on the mean state and the temporal dynamics of a plankton system. We employ a very simple 0-D physical framework similar to a chemostat in which the nutrient content of the inflow water is periodically switched between high and low concentrations. This framework has been used in many previous studies either to study the equilibrium states of a modeled system (eg., Thingstad et al. 1996;Ward et al. 2011;Prowe et al. 2012;Grigoratou et al. 2019) or its temporal dynamics, (e.g., Eisenhauer et al. 2009;Banas 2011;Kloosterman et al. 2016). Obviously, chemostat analogues are highly idealized configurations that cannot reproduce the complex and very diverse environmental conditions of the ocean. Furthermore, results are sensitive to the prescribed dilution rate and to the concentrations of the different organisms in the inflow waters that, in particular, controls the viability of the different modeled species and their response to imposed changes (Sommer 1986;Chatterjee and Pal 2013). However, owing to its simplicity, this physical framework represents a useful tool to analyze and understand the complex behavior of nonlinear systems developed to represent plankton ecosystems.
Our size-resolved description of metazoan organisms does not take into account interspecific and intraspecific differences. The actual body size is the sole structuring trait. Consequently, organisms of the same size are assumed to have the same metabolic and physiological characteristics as well as the same feeding behaviors. Different approaches have been proposed in the literature to introduce diversity into size-spectrum models. Some studies have adopted a mixed formalism between groupor species-based models and size-spectrum models (Blanchard et al. 2009;Maury 2010;Scott et al. 2014). The community is divided into several groups, each being modeled using a size-spectrum approach. A weakness of this mixed formalism is that within each community, body size remains the sole trait and therefore, inter-and intraspecific differences between organisms belonging to the same group remain disregarded. As an alternative to that approach, recent developments have focused on the introduction of a trait-based description of diversity. The most common additional traits are maximum size (Maury and Poggiale 2013) and size at maturity (Hartvig et al. 2011;Castellani et al. 2013). Intraspecific variability is generally represented by the addition of a diffusion term in the McKendrick von Foerster equation (Benoit and Rochet 2004;Datta et al. 2010).
The introduction of a size-resolved representation of mesozooplankton population produces fluctuations of the total biomass of metazoans which correspond to the propagation of successive cohorts through the whole size spectrum (Fig. 4). These fluctuations have been evidenced and extensively studied previously, often using PSPMs, (e.g., McCauley and Murdoch 1987;Persson et al. 1998). These studies showed that these cycles occur when juveniles have a competitive advantage over adults, for instance as a consequence of a higher mass specific ingestion rate (Persson et al. 1998;De Roos et al. 2008;Persson and De Roos 2013). Our size-resolved configurations do assume such an advantage for juveniles as the size scaling factor of ingestion rate is negative. Furthermore, when an allometric relationship is applied to each size class, cyclic cohort dynamics is stronger than when a uniform value is prescribed for juveniles and mature organisms. At equilibrium, cycles disappear despite the competitive advantage of juveniles which seems in contradiction with the previously quoted studies. Our results are difficult to compare to these previous studies which explore the properties of size-structured population models as there are significant differences. For instance, the ressource on which adults and juveniles are able to feed is very often assumed to be the same in many theoretical studies which is not the case in our model experiment: mature organisms feed both on microzooplankton and phytoplankton whereas juveniles only feed on phytoplankton. Furthermore, since our model represents a community, organisms can reproduce over a large size range defined by our mesozooplankton compartment. Finally, a quadratic mortality rate is applied which differs between the juvenile and the mature classes. Previous studies have shown that the way mortality is represented may significantly impact the behavior of the model, in particular its stability (Murdoch et al. 2003;van Kooten et al. 2007).
The temporal response of the community size-spectrum to short time-scales or seasonal perturbations has been investigated with single size-spectrum approaches that are similar to ours. They show qualitatively similar results, i.e., the propagation of successive waves of biomass from small to large organisms (Pope et al. 1994;Maury et al. 2007;Zhou et al. 2010). The inclusion of some representation of both intra-and interspecific diversity has been shown to increase the intrinsic stability of size-spectrum models (Datta et al. 2010;Zhang et al. 2013). The only study to our knowledge that has investigated the impact of seasonality with a sizespectrum model that includes a description of interspecific biodiversity is the one of Datta and Blanchard (2016). They showed that the simulated annual mean community size spectra are identical in seasonal and non seasonal setups. In addition, each modeled species exhibits propagating waves from smaller to larger individuals, similarly to what is simulated for the community by single size-spectrum approaches. Unfortunately, they did not explore the impacts of biodiversity on the community's response to seasonal fluctuations. Thus, the impact of both intra-and interspecific diversity on the response of an ecosystem to changes in the environmental conditions remains to be studied, at least at intra-seasonal to seasonal time scales. Therefore, we do not know how the representation of diversity would alter our results. This could be investigated in a future study.
The seasonal variability of the environment has resulted in adaptations in marine organisms, including zooplankton, in particular at high latitudes, such that their feeding activity, reproduction and maintenance may display strong seasonal rhythms (McNamara and Houston 2008). Concerning reproduction, our model formulation assumes that zooplankton are income breeders and thus, allocate part of the available resources directly to reproduction. This choice was dictated by simplicity. However, organisms may adopt an alternative reproduction strategy which is termed capital breeding. During the feeding season, an individual may allocate energy to reserves that are used later in the year to reproduce (Jönsson 1997;Varpe et al. 2009). Capital breeding has been suggested to be a better strategy in highly seasonal environments with a short feeding season, such as the polar oceans (Sainmont et al. 2014). Obviously, changing the breeding mode in our model configuration would significantly change the simulated temporal dynamics of juveniles and mature organisms as reproduction would become temporally decoupled from feeding.

Implications
Mesozooplankton occupy a key position in most ocean food chains. As predators of smaller plankton, they can influence the dynamics of microzooplankton and phytoplankton. As producers of fecal pellets (Steinberg and Landry 2017), and since they are able to migrate vertically in the water column (Emerson et al. 1997), they also control the transfer of matter from the ocean surface to depth, and as such may influence the marine carbon and nutrient cycles. As preys for larger organisms, mesozooplankton also control the transfer of energy to the higher trophic levels (Beaugrand et al. 2010). By eating and being eaten, mesozooplankton both exert a top-down control (on marine biogeochemistry and lowertrophic levels) and a bottom-up control (on higher trophic levels) (Kiørboe 1997). Bottom-up effects induced by the size structuration of a population with size-dependent lifehistory traits have been summarized by De Roos (2020). In particular, a predator feeding on a size structured prey may suffer from an Allee effect (i.e., reduced reproductive success at low population densities) with the emergence of a stable state where the predator population collapses, a behavior that is not observed in the absence of size structuring (De Roos and Persson 2002).
Here, we show that the explicit representation of meatozoans reproduction and size-growth in a idealized 0-D biogeochemical model strongly affect the temporal dynamics of mesozooplankton. However, we do not highlight any effects on the temporal dynamics of phytoplankton and nutrients when including reproduction and growth in size. Previous studies have on the contrary demonstrated potential top-down effects. It has been shown that early life stages of copepods, such as nauplii, can represent a significant part of microzooplankton (up to 30%) (Quevedo and Anadón 2000;Safi et al. 2007) and hence induce a top down control of mesozooplankton abundance (in this case with the species M. leidyi) on microzooplankton (McNamara et al. 2013). Also, mesozooplankton have been shown to be responsible for a significant part of the vertical export of carbon due to fecal pellets (Stamieszkin et al. 2015). Exploring those latter processes implies a vertically resolved representation of biogeochemical processes, which is absent in our model. This could explain why we do not observe any effect of our different representations of mesozooplankton on nutrients and lower-trophic levels.
As an explicit representation of individual size growth and reproduction of mesozooplankton in our model induces a delay in the mesozooplankton response to a phytoplankton bloom, we could expect that in a model coupled with upper trophic levels, the response of mesozooplankton predators would also be delayed. Indeed, numerous examples of bottom-up effect of mesozooplankton on upper trophic levels have been described which calls to the developments of robust zooplankton submodels in end-to-end marine ecosystem models (Rose et al. 2010;Mitra et al. 2014). As copepods dominate mesozooplankton in many systems (Gallienne 2001), studies on those species may give a representative overview on the bottom-up control of mesozooplankton on upper trophic levels. Bi et al. (2011) focused on the correlation between water-type affinities of copepod communities and salmon type predominancy along the Washington and Oregon coast and have showed that the predominant salmons species shifts followed switchs between cold-and warm-water copepod species, which suggest that salmon population dynamics is controlled by prey availability. In the Mediterranean Sea, Saraux et al. (2019) have shown that variations in plankton community around 2008, especially marked by a decrease in copepod abundance, have had a direct impact on fish biomass and on the mean size and weight of sardine and anchovy, highlighting again an important bottom-up effect.

Conclusion
The explicit representation of metazoan life-cycle in our simple biogeochemical model does not affect the temporal dynamics of non-mesozooplankton compartments. However, mesozooplankton temporal dynamics is changed, mostly through a few months delay in the response to an increase in food availability when reproduction is taken into account. In model versions with reproduction, we also identify an early acceleration phase, in which rates are particularly high, but which happens before mesozooplankton biomass reaches a significant value. Lastly, accounting for ontogenetic growth induces cohort dynamics in metazoans. The comparison of our 4 model versions enables to get some insights into the main factors inducing these changes.
First, as a consequence of size structuration, the temporal dynamics is affected by the size-based parameterization of the maximal growth rates. At the start of HNM, when zooplankton biomass is limited, metazoan individuals are mostly juveniles, i.e., they are part of the microzooplankton and therefore benefit from high growth rates. Those fast growing individuals are then matured into adult animals, inducing a faster initial increase of mesozooplankton to nutrient input as compared to a model without reproduction. In addition, temporal changes in the size distribution of mesozooplankton lead to the appearance of cohort dynamics, totally absent in model configurations without discretization of metazoan size.
The second factor influencing the temporal dynamics is the introduction of non-predation terms between microzooplankton and mesozooplankton size-classes, i.e., reproduction and maturation. Differences in the temporal dynamics of those two opposite terms induce a delay in the mesozooplankton response to equilibrium: reproduction, a loss term for mature mesozooplankton, reaches its equilibrium faster than the gain term, maturation. This results in a net loss term and reduces mesozooplankton growth before equilibrium is reached, and hence retards this equilibrium when compared to a model without any explicit representation of reproduction.
Mesozooplankton play a key role in marine ecosystems because of their intermediate position in trophic food webs. Explicit representation of mesozooplankton reproduction and ontogenetic growth should be considered when modeling the upper trophic levels as the dynamics of upper trophic levels feeding on mesozooplankton could be affected, particularly in regions where the seasonality is marked. Numerous top down effect on lower trophic levels (particularly microzooplankton) and biogeochemistry (particularly on the vertical carbon export), not testable in the context of our 0-D study, would be expected in a vertically resolved model. Size spectrum approaches are interesting avenues that could provide all the tools necessary to improve mesozooplankton community representation in biogeochemical models.

Supplementary information
The online version of this article (https:// doi. org/ 10. 1007/ s12080-021-00519-5) contains supplementary material, which is available to authorized users. the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.