Timing of Calanus finmarchicus diapause in stochastic environments

In environments with strong seasonality, many herbivorous zooplankton remain active only during the productive season and undergo a period of inactivity and suppressed development termed ‘diapause ’ during the unproductive season. The ability to time the diapause entry and exit in response to the seasonality of the environment is thus essential for their survival. However, timing of diapause may become challenging when environmental conditions vary stochastically across shorter and longer timescales, and particularly when zooplankton lack external cues to predict these variations. In this study, we used a novel individual-based model to study the emerging patterns of diapause timing of the high-latitude marine herbivorous copepod Calanus finmarchicus under shorter-(6-h) and longer-term (interannual) environmental stochasticity. The model simulated growth, development, survival and reproduction (income breeding) of a C. finmarchicus population over multiple calendar years and traced the emergence of behavioral responses and life history strategies. The emergent timing of diapause entry and exit were robust to shorter-term environmental stochasticity, which was manifested through morphological (i.e., body and energy reserve sizes) and behavioral plasticity (i


Introduction
Characterizing variability in natural environments and understanding the adaptations of organisms to these changes are fundamental aspects of behavioral and evolutionary ecology.Environmental changes occur both across time (temporal environmental heterogeneity) and space (spatial environmental heterogeneity) (Pigliucci, 2001).A part of this heterogeneity follows well-established cyclic patterns, such as the diel and seasonal variability of solar irradiance; or trends, such as the gradual change of climate across latitude and altitude.Organisms generally respond to environmental heterogeneity in three main ways: (i).When environmental changes are long-lasting and the selection pressures are persistent, standing genetic variation together with recombination and mutation will produce locally adapted genotypes through adaptive evolution ('adaptive tracking') (Byers, 2005).(ii).When selection pressures fluctuate across shorter spatio-temporal scales, a strategy termed 'phenotypic plasticity' becomes potent.This involves the production of different phenotypes by a single genotype that characterize various physiological, morphological, behavioral and life history responses (Ghalambor et al., 2007).(iii).In environments with highly unpredictable variability, a strategy termed 'bet hedging' reduces fitness fluctuations of genotypes between generations (Slatkin, 1974).
This may occur through the production of phenotypes that perform reasonably well under 'good' and 'bad' environmental scenarios (conservative bet hedging) or by spreading the risk of certain behavioral and life history decisions among an array of different phenotypes (diversified bet hedging) (Simons, 2011).
Zooplankton in high-latitude marine environments are particularly exposed to spatial and temporal environmental heterogeneity.In terms of time, the daylength (photoperiod) may transit from near-continuous darkness (polar night, ~0L:24D) to near-continuous daylight (midnight sun, ~24L:0D) within a year.This seasonality results in a productive season with well-demarcated windows of high primary production (e.g., spring-summer) and an unproductive season (e.g., autumn-winter).Although these irradiance-driven environmental variations are largely predictable, notable stochastic oscillations occur atop these diel and seasonal patterns.For example, shorter-term variations of cloud cover can attenuate subsurface irradiance and momentarily influence the efficiency of visually orientating pelagic predators (Eiane et al., 1997;Ryer and Olla, 1999).In addition, the upper-pelagic ambient temperature and primary production may also vary in the shorter-term depending on the fluctuations of irradiance, extent of convective vertical mixing, horizontal advection and grazing pressure (Brainerd and Gregg, 1995;Cushing, 1990;Wang et al., 2005;Wood and Corcoran, 1966).In the longer-term, timing and the duration of the productive season may vary between years depending on the fluctuations of ambient temperature, water column stratification, ice-retreat timing and nutrient supply (Uitz et al., 2010;Wassmann et al., 2006).In terms of space, zooplankton encounter environmental gradients (e.g., temperature, salinity, density, dissolved oxygen, food availability, predation risk) during their routine vertical excursions through the water column.Zooplankton also continuously drift with water currents, which can either be of shorter spatial extent (e.g., cross-shelf exchange processes: Torgersen and Huse, 2005) or of broader trans-latitude extent (e. g., Ji et al., 2012).Although vertical and latitudinal environmental gradients follow generic trends and are predictable to some extent, frequent changes in vertical mixing and mesoscale phenomena, such as eddies, hydrographical fronts and varying bottom depths can induce stochastic environmental variations in space (Mann and Lazier, 2006;McGillicuddy, 2016).
Adaptations of high-latitude zooplankton to predictable cyclic environmental variations are relatively well-understood.These adaptations are most pronounced among the herbivore community (Hagen and Auel, 2001;Varpe, 2012).The food (phytoplankton) availability for herbivorous zooplankton is seasonally limited and restricted to the upper pelagial (photic zone).The earliest available food emerges in the spring when the upper pelagial is typically colder and perhaps under ice cover.This may reduce grazing, assimilation and growth rates of herbivorous zooplankton (Huntley and Lopez, 1992;Romare et al., 2005).Although the primary production generally extends towards summer, the near-constant summertime illumination elevates the light-dependent (visual) predation risk in the upper pelagial.Zooplankton tend to minimize visual predation risk by descending to deeper, darker layers to take refuge during daytime while ascending to the near-surface layers to feed during the nighta behavior known as 'diel vertical migration' (DVM) (reviewed in Bandara et al., 2021;Brierley, 2014).Although DVM is effective against visual predation risk, it typically shrinks the daily foraging window, which leads to slower growth and development rates (Bandara et al., 2018;Loose and Dawidowicz, 1994).Because of higher predation risk or fading food concentrations towards autumn, most herbivorous zooplankton descend to deeper waters and enter a state of hibernation termed 'diapause' for 6-8 months of the year (reviewed in Bandara et al., 2021;Baumgartner and Tarrant, 2017).
Investigations of the adaptations of high-latitude zooplankton to stochastic environmental variability are mostly conducted in the shorter-term and focus on the plasticity of zooplankton behavior to stochastic environmental oscillations.For example, Eiane and Parisi (2001) and Record and Young (2006) used empirical data and simulation models to study the changes of zooplankton DVM in response to variations in cloud cover.Their findings suggest that zooplankton can instantaneously modify their DVM behavior to increase grazing time in near-surface layers even at midday when the skies are relatively cloudy and darka classic example of the plasticity of DVM behavior (see Bandara et al., 2021 for a review).Similar cases of DVM plasticity have been observed and predicted in response to sudden and often unpredictable subsurface irradiance attenuation episodes caused by solar and lunar eclipses (e.g., Strömberg et al., 2002;Tarling et al., 1999), smoke from wildfire (e.g., Urmy et al., 2016), suspended matter and algal blooms (e.g., Fiksen and Carlotti, 1998;Williamson et al., 2020).
Investigations of the influence of longer-term environmental stochasticity on diapause timing of high-latitude zooplankton are rare, and the present knowledgebase is largely theory-driven.For example, Ji (2011) viewed the model-predicted variations of diapause entry timing of the high-latitude copepod Calanus finmarchicus in the light of phenotypic plasticity, which allows the animals to cope with spatial and temporal (interannual) environmental stochasticity expected in the western north Atlantic.Further, a model for C. finmarchicus by Fiksen (2000) demonstrated prioritized reserve accumulation and early diapause entry in stochastic settings (see also Kvile et al., 2018;Varpe and Ejsmond, 2018).Since diapause entry usually negates the attainment of sexual maturity and production of an additional generation or two within the same productive season (see also Kaartvedt, 2000;Varpe and Fiksen, 2010), such strategies are viewed as acts of conservative bet hedging, which can be essential for survival when the growth potential (i.e., food availability and temperature) varies stochastically between years.Empirical evidence in this regard exists in freshwater literature, where diapausing stages (i.e., resting eggs) of rotifers are produced earlier when the length of the growing season varies stochastically over time (e.g., Franch-Gras et al., 2019;Tarazona et al., 2017).In contrast, the role of bet hedging in diapause exit is poorly understood, except in the case of the hatching of resting eggs from egg banks (Evans and Dennehy, 2005).In addition, there is no broad agreement about the mechanism(s) by which high-latitude herbivorous zooplankton time their diapause exit in relation to the interannually variable timing of pelagic primary production.The timing of diapause exit becomes even more enigmatic when diapausing populations are advected into different environments by moving water masses (Espinasse et al., 2016;Rullyanto et al., 2015).This is because timing and the duration of primary production vary significantly across space (Daase et al., 2013;Falk-Petersen et al., 2009) and diapausing populations seem to lack perceivable cues to predict the timing of pelagic primary production, particularly when they occupy habitats several thousand meters below the productive near-surface waters (Østvedt, 1955).
Our research focuses on investigating the adaptations of a predominantly herbivorous high-latitude copepod Calanus finmarchicus to environmental stochasticity.C. finmarchicus is a species that encounters pronounced environmental variability in their natural habitats since its geographical range spans across ca.40 • of latitudes across the North Atlantic towards the high-Arctic (Conover, 1988;Fleminger and Hulsemann, 1977;Melle et al., 2014).In this study, we present an individual-based model that can simulate the life cycle of C. finmarchicus in great spatial, temporal and biological resolution and use it to investigate its adaptations to environmental stochasticity.In particular, we invoke the potential ecological roles of adaptive tracking, phenotypic plasticity and bet hedging as strategies that aids C. finmarchicus to time its diapause in environments where the irradiance, temperature, food availability and predation risk vary stochastically on diurnal, seasonal and interannual basis.

Model overview
The present model is an extension of our previous work Bandara et al. (2018) and(2019), which were aimed towards investigating the adaptations of high-latitude herbivorous copepods (Calanus spp.) to predictable cyclic diel and seasonal environmental heterogeneity (deterministic settings).In this study, we added a level of stochastic environmental variations atop these predictable cyclic patterns.Further, we upgraded the Genetic Algorithm based strategy-oriented construct of the predecessor models to a more versatile individual based construct Fig. 1.Simplified conceptualized construct of the model.The model focuses on a herbivorous copepod occupying a high-latitude seasonal environment (A).Model copepods comprise a 'genome' of five 'genes' that represent various attributes of their behavior and life history.The 'genome' information of adult male and female copepods recombines and mutates to produce new 'genotypes' (B).The growth and development, survival and reproduction of copepods carrying these 'genotypes' (= 'phenotypes') are simulated in artificial seasonal environments, which may vary stochastically in the shorter-and longer-term.These simulations are performed for 100 calendar years and depending on the model environment, their 'genes' may eventually be fixed in the population (C).and made notable changes and upgrades to the growth and development, survival and reproductive submodels (see Appendix A1 in Supplementary material for a detailed comparison of the present model and its predecessors).
The model comprises two entities: model copepods and the model environment.Model copepods are simplified representations of Calanus finmarchicus in terms of morphology (morphometry), behavior and life history (see Section 2.2.1.1 for details).The behavioral and life history strategies of model copepods are determined by five evolvable attributes ('genes') that collectively form their 'genome'.The 'genome' of each copepod contains information on energy allocation patterns, body size, and the timing of DVM, diapause and seasonal vertical migration (SVM) (Fig. 1A, B, Table 1).The model environment is a 500 m-deep unidimensional seasonal setting resolved to 1 m bins (see Section 2.2.1.2for details).
This model follows an open-ended life cycle simulation approach.Accordingly, starting from an artificially seeded batch of eggs, the growth and development, survival and reproduction of model copepods are simulated in model environments at 6 h temporal and 1 m vertical spatial resolution (see Section 2.2.2 for details).During reproduction, the 'genome' information of male and female copepods is recombined and mutated to generate new 'genotypes' that may produce 'phenotypes' with different behavioral and life history strategies (i.e., each mating male and female pair can produce a range of 'phenotypes': Fig. 1B).The model environments may vary stochastically in the shorter term (between 6-h intervals) or longer-term (interannually) during the simulation.Simulations are performed for 100 calendar years (i.e., several hundred generations, assuming a ≤ 1 year generation time for C. finmarchicus:).Depending on the simulated environmental dynamics, model copepods with favorable strategies may attain higher fitness and their 'genes' may eventually be fixed in the population over time (Fig. 1C).

Model description
2.2.1.Entities 2.2.1.1.Model copepods.C. finmarchicus is a predominantly herbivorous copepod species that inhabits the North Atlantic and Arctic oceans (Conover, 1988).Viable C. finmarchicus populations can range from the Gulf of Maine (~43 • N) (Fish, 1936) to north of the Svalbard archipelago ( > 80 • N) (Daase and Eiane, 2007).This broad trans-latitude geographical distribution is due to their association with Atlantic water masses, which transport C. finmarchicus populations from their center of distribution in the North Atlantic towards the Arctic Ocean.This exposes them to notable variations in temperature, irradiance, food and predation environments.
C. finmarchicus has a life cycle with 13 developmental stages (egg, six naupliar stages, five copepodite stages and adults).In higher latitudes, it typically follows an annual life cycle (one generation per year), whereas in lower latitudes, it may complete 2-3 generations each year (Melle et al., 2014).Towards the end of the feeding season, late-juvenile stages (copepodite stages IV and V) of C. finmarchicus store lipids and descend to deeper waters for diapause (Falk-Petersen et al., 2009).Diapause duration varies between 4 and 6 months, during which the stored lipids are metabolized at a considerably lower rate than for metabolism during the parts of the year when they are active (Hirche, 1996a;Maps et al., 2013) (Fig. 1A).C. finmarchicus has many similarities with other herbivorous copepods such as Calanoides acutus in the Southern Ocean and Calanus glacialis in the northern hemisphere.

Model environment.
The model environment characterizes three variables: irradiance, temperature and food concentration.Their formulation is originally based on Cottier et al. (2010) and are described in detail in Bandara et al. (2018) and ( 2019).In the most basic form, the Active metabolic rate Eq. ( 14) Basal metabolic rate (sizedependent) Eq. ( 13) Basal metabolic rate (temp.dependent) Eq. ( 14) Parameter for satiation food conc.

Growth and development.
The growth and development of model copepods follow the formulations developed for C. finmarchicus by Maps et al. (2012a) and Bandara et al. (2019).All rates below are given as hourly estimates, which are adjusted to the temporal resolution of the model (6 h) in calculations.The body mass (µg C) is the sole proxy of body size used in this model.The somatic growth rate of individual i at timestep t in depth bin z (G, µg C ind − 1 h − 1 ) is estimated in carbon units as a balance between the assimilation (A, µg C ind − 1 h − 1 ) and metabolic rates (B, µg C ind − 1 h − 1 ) as, The assimilation rate is a product of the ingestion rate (I, µg C ind − 1 h − 1 ) and the assimilation coefficient (a: Table 1).The growth equation can thus be rearranged into, At the reference temperature of − 2 • C, the ingestion rate relates with the structural mass of the copepods (W c, µg C) as, Here, the terms b and m represent mass coefficient and exponent of ingestion (Table 1).The ingestion rate relates exponentially with the ambient temperature as, Here, I ´is the maximum temperature-dependent ingestion rate, and c and n are temperature coefficient and exponent of ingestion (Table 1).Warmer years in reddish hue and colder years in blueish hue-synonymous with color representation in Fig. 6.
The ambient food concentration (F, µg C l − 1 ) maps the temperaturedependent ingestion rate between 0 and 1 as, Here, the ingestion rate becomes solely temperature-dependent above a satiation food concentration, which is given by the sizedependence of the parameter d as, The satiation food concentration thus increases with the body mass of the copepods and falls in the range of 75-125 µg C l − 1 (cf.Campbell et al., 2001;Huntley and Boyd, 1984).
The growth formulations in Eqs.
(1)-6 do not apply to eggs and the first two nauplii stages (NI and NII), which do not feed.Instead, they encounter a negative growth as the energy reserves are catabolized to meet the energetic demands (Marshall and Orr, 1972).The development of non-feeding stages is thus solely temperature-dependent (Corkett et al., 1986) and occurs following a Bělehrádek function as, Here, D is the estimated development time in 6 h time intervals and thus the scalar (4) on the right-hand side of the Eq. ( 7).The stagespecific (1 ≤ j ≤ 3) values for the parameter q are based on the estimates of Campbell et al. (2001) (Table 1).The development of feeding stages (NIII-adult: 3 < j ≤ 13) occurs as the structural mass of copepods (W c ) exceed a stage-specific critical molting mass (W j ).We estimated the W j values following the growth formulations of Maps et al. (2012a) for C. finmarchicus.At each environmental setting, we produced the minimum and maximum estimates for W j (i.e., W j min and W j max ) by running the Maps et al. (2012a) model for annual maximum and minimum temperatures at non-limiting food concentrations.Consequently, W j min and W j max estimates fluctuate between years depending on the interannual ambient temperature and food concentration variations.The 'gene' α (0-1) maintains the variability of body size within the copepod population.The environment-specific critical molting mass for any developmental stage (NIII onwards) is estimated as, Here, based on the value of the 'gene' α, the ontogenetic body mass trajectories of copepods tend to occupy a fixed fraction of the estimated minima and maxima.

Predation risk.
Similar to the study by Bandara et al. (2019), the predation risk in this model possesses two components: the light-and size-dependent visual predation risk and the non-visual (tactile) predation risk.The light-dependency of visual predation risk is given as, Here, L t,0 and L t,z are the estimated irradiance (µmol m − 2 s − 1 ) on the modelled sea surface and at depth z (500 m ≤ z > 0 m).The constant ψ represents the water column light attenuation coefficient (Table 1).We scaled L t,z to obtain a probability metric (L ′ t,z ) that ranges between 0.1 − 0.9, which offers non-zero probability for survival and death respectively at the highest and lowest levels of irradiance (Bandara et al., 2018(Bandara et al., , 2019)).The visual predation risk (M v , i.e., the probability of death by visual predation) is thus given as, Here, K (10 − 7 − 0.25) is a scalar for visual predation risk that adjusts the visual predation risk to the temporal resolution of the model.The asymptotic exponential function at the right-hand side of Eq. ( 10) represents the size-dependence of visual predation risk (Fig. 3A), where W cs (µg C) is the total mass of the copepod (i.e., the sum of the structural mass, W c and the energy reserve mass, W s ).
Unlike the light-and size-dependent visual predation risk, the modelled non-visual predation risk (M n , i.e., the probability of death by non-visual predation: range = 10 − 7 − 10 − 3 ) is kept constant over time and depth for simplicity (Eiane and Parisi, 2001).

Diel vertical migration (DVM).
In this model, the shorter-term vertical behavior of the modelled copepods is driven by their photoreactive behavior (Ringelberg, 2010).Thereby, the model copepods react to absolute light intensity and tend to avoid irradiance levels above an individual-specific threshold, determined by the 'gene' β (0 − 1).If the maximum irradiance of the current environmental setting is L max µmol m − 2 s − 1 , the threshold irradiance of a given individual at a given time (L threshold , µmol m − 2 s − 1 ) is estimated as, The asymptotic exponential function at the right-hand side of the Eq. ( 11) highlights the size-dependency of the threshold irradiance, which decreases (and light sensitivity and risk averseness increases) over a copepod's lifespan (Fig. 3B).Accordingly, at any given time, the model copepods tend to remain at a depth, which provides the maximum growth potential (Eq. ( 1)) below the estimated irradiance threshold (Bandara et al., 2019).Given the diel periodicity of modelled irradiance (Fig. 2A, B), the above photoreactive behavior leads to a classic DVM pattern, i.e., the occupation of shallower waters during nighttime and retreat to deeper waters during daytime.For each copepod at any given time, we predicted the depth with maximum growth potential below the threshold irradiance deterministically and assumed that the neutrally buoyant copepods reach this depth via cruising (van Someren Gréve et al., 2017) where individuals swim at a constant velocity (Fig. 3C) as, 2.2.2.2.3.Energy storage.The entire surplus assimilation (Eq.( 1)) of younger developmental stages (NIII-CIII) is allocated to somatic growth (indicated by the structural mass, W c , µg C).Late-juvenile CIV and CV stages can allocate an individual-specific fraction of the surplus assimilation to build up an energy reserve (indicated by the energy reserve mass, W s , µg C), which may occupy up to the entirety of the individual's structural mass.Here, 'gene' γ (0 − 1) defines the pattern of energy allocation, where the entire surplus assimilation is allocated to structural growth if γ = 0 and energy reserves if γ = 1.The adult stages of this model do not build up energy reserves but inherit energy reserves accumulated in their pre-adulthood.

Diapause and seasonal vertical migration (SVM).
Each model copepod possesses an equal probability to either (i).perform a seasonal descent into deeper waters, undergo diapause, ascend back to nearsurface waters (SVM) and then develop into adults or (ii).directly develop into adulthood without diapause and SVM.As in our predecessor models, we used the state of the energy reserve as a proxy of timing of diapause entry and exit, a reasonable assumption for the primarily income breeding C. finmarchicus (Varpe and Ejsmond, 2018).Here, copepods descend to deeper waters when their stores reach an individual-specific fraction of the structural mass (either at CIV or CV stage), which is determined by the 'gene' δ (0 − 1).The selection of diapause habitat is simple and similar to that of Carlotti and Wolf (1998) and Bandara et al. (2019) where individuals are randomly placed at depth below the mixed layer (300− 500 m).During diapause, the vertical position of each copepod passively varies by a random amount at each time step (± 25 m).The model copepods utilize stored reserves to meet the energy requirements during diapause, which occurs at a lower rate compared to the basal metabolic rate.Diapause terminates when an individual-specific fraction of the stored reserves is exhausted, which is determined by the 'gene' ε (0 − 1).As C. finmarchicus is not known to possess > 1-year life cycles (reviewed in Bandara, 2014;Falk-Petersen et al., 2009), we did not model a diapause strategy among adult stages, and thus the maximum generation time simulated in this model is 1-year.The copepods that develop directly into adulthood may produce one or more generations within a calendar year.
2.2.2.2.5.Metabolism.The metabolic rate is the sum of the basal metabolic rate (B b , µg C ind − 1 h − 1 ) and the active metabolic rate (B a , µg C ind − 1 h − 1 ).These hourly estimates are adjusted to the temporal resolution of the model (6 h) during calculations.At the reference temperature − 2 • C, the relationship of B b with the body mass is given as, Here, f and o are mass coefficient and exponent of respiration (Table 1).The relationship between the basal metabolic rate and the ambient temperature is given as, Here, B ′ b is the temperature-dependent basal metabolic rate, and g and p are temperature coefficient and exponent of metabolism, respectively (Table 1).The active metabolism (B a ), which is assumed to account for 1.5 B b , is added when the individuals are swimming via cruising (Bandara et al., 2019).During diapause, the B a is 0 and B b occurs at 75% reduced rate (Maps et al., 2013).At each time interval, the metabolic costs are deducted from the gross assimilation (Eq.( 1)).When metabolic demand is larger than assimilated energy, energy reserves are mobilized.

Starvation risk.
Starvation sets in when the metabolic demands exceed energy reserves.When a model copepod starves, the structural mass is catabolized to meet the metabolic demands.This loss of structural mass increases the starvation risk (M s , i.e., the probability of death by starvation) as, Here, W cx is the catabolized structural mass expressed as a proportion to the maximum structural mass attained prior to structural catabolization.The above function approaches an upper asymptote (= 1.00, Fig. 3D) as W cx reaches 0.50 (the Chossat's rule: Chossat, 1843).

Population ceiling.
The maximum population size (P max ) allowed in this model is 10 6 , which is solely dependent upon the available computational resources for model simulation-not on natural resource (food) limitation due to the lack of two-way coupling between the feeding and phytoplankton dynamics in this model.To avoid abruptly cutting off the population at P max , when the population size at a given time (P t ) approaches P max , a computational-resource-dependent artificial mortality (M a ) set in as, Therefore, M a is negligible at smaller population sizes and begins to increase exponentially as P t approaches P max (Fig. 3E).M a operates irrespective of copepod's internal states, such as the body mass, energy reserve mass and the developmental stage.

Fig. 3. A:
The relationship between visual predation risk (M v ), total mass (W cs ) of the copepod and the scalar K at highest level of ambient irradiance (i.e., L' = 0.9: Eq. ( 10)).B: Relationship between irradiance sensitivity (presented as a threshold irradiance, I threshold ) of the copepod with the total mass (W cs ) and the attribute value of 'gene' β (Eq.( 11)).C: Relationship between the cruising velocity (U) and the structural mass (W c ) (Eq. ( 12)).D: Relationship between the modelled starvation risk (M s ) and the proportion of catabolized structural mass (W cx ) (Eq. ( 15)).E: The emergence of computational-resource-dependent artificial mortality (M a ) as the simulated population size (P t ) approaches the ceiling P max (Eq.( 16)).

Reproduction
CV stages possess an equal probability of molting into either an adult male or a female (sex ratio = 1).We assumed that adult stages reproduce immediately following the final molt.The possibility of sex switching at the adult stage (sexual dimorphism: Svensen and Tande, 1999) is disregarded for simplicity.The adult males in this model do not feed (Mauchline, 1998) but use energy reserves to meet their metabolic demands, and thus are generally short-lived.Females mate only once in their lifetime (Hirche, 1996b;Marshall and Orr, 1972;Titelman et al., 2007) with a random male.However, modelled males can mate with more than one female (Nicholls, 1933).When mating, a copy of the male 'genome' is transferred to the female.We assumed that adult females follow a pure income breeding strategy (cf.Varpe et al., 2009), where the egg production is entirely dependent on the food intake (Hirche, 1996b).The egg production (E i,t ) is estimated as, Here, W R is the matter allocated to the egg production and W E is the unit egg mass (Table 1).During egg production, the male and female 'genes' recombine at a probability of 0.70 per 'gene' following a heuristic crossover method (Haupt and Haupt, 2004;Michalewicz, 1996).Recombined genetic information undergoes mutation at a probability of 0.20 per 'gene' using random replacement (Eiben and Smith, 2003).Females can produce a maximum of 1000 eggs (Carlotti and Hirche, 1997) and are 'killed' afterwards, hence assuming strict semelparity (Varpe and Ejsmond, 2018) to allow space in the simulated population.
As the maximum allowed population size (P max ) in this model is limited to 10 6 , the total egg production of the population at a given time may exceed P max .This excess egg production is systematically cut off using a fecundity-proportional selection method as, Here, E i,t is the estimated egg production of a given mated female at a given time and P t is the estimated population size at that time.

Model operation
The model initializes at time t (at 0000 h of 1 January of year 1) with the generation of a model environment for the entire calendar year (up to 1800 h on 31 December).In stochastic simulations (see below), the model environment may vary between time intervals and/or between calendar years.The simulation begins with the seeding of 100 eggs with random 'gene' combinations to random depths (< 100 m) of the water column at each time interval throughout the first calendar year of the simulation (i.e., a total of 146,000 eggs).At each time interval, the somatic growth, developmental progression, reserve build-up (in latejuvenile stages) and the vertical position of each copepod is simulated and are updated into five state variables.Thereafter, each copepod is assessed for survival as, Here, S is a binary state variable (survivorship), which takes the default value 1 ('alive') and is set to 0 ('dead') if the estimated survival probability is smaller than the individual-and time-specific random number (ω it ) drawn from a uniform probability distribution.When a copepod is presumed 'dead' (i.e., S i,t = 0), all the corresponding state variable values are erased to allow space for a new individual in the simulated population.At each time interval, when the survival assessment of all copepods is complete, the population size (P t ) is estimated and updated into a state variable (population size) as, The model employs two types of tracking (bookkeeping) variables.High-resolution trackers are updated at each time interval and record the total and stage-specific population sizes, the number of diapauseentries, number of diapause-exits and the number of directdevelopments (i.e., individuals that develop directly into adulthood without diapause), along with their structural and energy reserve masses into a set of unidimensional arrays.Further, the vertical distribution of the population along with the distribution of temperature, food concentration and irradiance throughout the water column are recorded into a set of two-dimensional arrays at each time interval.In stochastic environments, high-resolution trackers are written to the disk at the end of each calendar year of the simulation.In the deterministic environment, high-resolution trackers are only recorded and written to the disk at the final calendar year of the simulation.Low-resolution trackers keep an annual mean of the 'genome' information of the population and are written to disk at the end of the final calendar year of the simulation.The model simulations continue over multiple years and terminate as t approaches a termination time, which is typically 1800 h of 31 December of the year 100.

Simulation experiments
Model simulations were performed under four main scenarios.(i).A basic model simulation was performed in a deterministic environment (Fig. 2A-C).(ii).Model simulations were performed in an environment with shorter-term (6-h) stochasticity (Fig. 2D-F).In this environment, the mean ambient temperature, food concentration, predation risk (K = M n = 10 − 5 ) and the timing and duration of pelagic primary production remained constant between years.(iii).Model simulations were performed in an environment with both shorter-and longer-term stochasticity in ambient irradiance, temperature and food availability (Fig. 2G).(iv).The model was run in the same environment in the above scenario iii, but with the addition of uniformly random interannual variations to visual (K, range = 10 − 7 − 0.25: Fig. 2H) and non-visual predation risks (M n , range = 10 − 7 − 10 − 3 : Fig. 2I).

Model development, analysis and archiving
The model was developed in FORTRAN 95 (ISO/IEC 1539-1:1997) with the support of OpenMP application program interface modules for parallelization.The model simulations were performed in a custom-built gaming rig with a liquid-cooled Intel® Core™ i9-7920X processor with 24 nodes (workers) running at an overclocked turbo frequency of ca.4.4 GHz.The model outputs were analyzed using R™ version 3.4.1 (R Core Team, 2017) and RKWard™ version 0.7.0 (Rödiger et al., 2012) in Ubuntu version 18.04.The data presented in both the deterministic and shorter-term stochastic simulation experiments are outputs from the final calendar year of simulations.In contrast, in the simulation experiments with shorter-and longer-term environmental stochasticity, we used the mean values of environmental variables and emerging behavioral and life history attributes from each calendar year in data presentation.To interpret the environmental correlates of emerging strategies in stochastic environments, the above values were standardized (i.e., centered by mean and divided by standard deviation) and used in a Principal Component Analysis (PCA), which was run in a correlation matrix.Raw data of these simulations along with the source code of the model can be downloaded at https://git.io/JclNW.

Deterministic model environment
The abundance of the simulated C. finmarchicus population peaked in spring (April-May) and then in autumn (July-September) (Fig. 4A).The first peak was dominated by eggs and younger developmental stages (Fig. 4B).The proportion of late-juvenile stages (copepodite stages IV and V) began to increase towards mid-May.By late June, a fraction of CIVs and CVs descended to deeper waters and entered diapause (Fig. 4C,  H).The size of the energy reserve of these diapause-entries was relatively low (mean ≈ 11 µg C: Fig. 4D), which was metabolized over a relatively short duration (ca.1-3 months).These individuals, therefore, ascended to the surface waters again relatively quickly, after what seems to be a short summertime diapause.Consequently, the number of diapause-exits with nearly spent energy reserves increased from July to August (1000-1600 exits d − 1 : Fig. 4E, F).
The proportion of adult males and females gradually increased from July onwards (Fig. 4B).Approximately 34% of these adult stages emerged following the final molt of summertime diapause-exits.The rest of the adults matured from late-juvenile stages and completed their final molt without undergoing diapause (Fig. 4G).A late summer egg production peak (although less prominent compared to that in the spring) coincided with the emergence of adults in July and contributed to a second generation (Fig. 4B).The proportion of second-generation CIV and CV stages began to increase towards autumn, and ca.60% of them migrated to deeper waters for diapause with considerably higher energy reserves compared to the summertime diapause-entries (mean W s ≈ 80 µg C) (Fig. 4B-D, H).Some second-generation CIVs and CVs did not undergo diapause and developed directly to adults and produced a third generation in late August and early September (Fig. 4B).While the thirdgeneration copepods developed under diminishing food conditions, some of their late-juvenile stages (predominantly CIVs) descended to deeper waters with relatively low energy reserves and entered diapause by mid-September (mean W s ≈ 21 µg C) (Figs. 2C and 4C, D, H).Individuals that did not reach an overwintering stage or those that failed to build sufficient energy reserves starved to death in the upper pelagial as the modelled food supply faded into the winter.
Diapause exit continued during the autumn-to-winter transition (October-December) in smaller numbers (< 200 exits d − 1 , mostly CIV stages).The energy reserves of these diapaus-exits were spent, and, Fig. 5. Comparison of emergent life history attributes and their timing in the deterministic and shorter-term stochastic model environments.The color shading indicates the extent (on a percentage scale) to which the attributes had deviated under shorter-term environmental stochasticity (cf.Fig. 4).
K. Bandara et al. consequently, they died in the upper pelagial due to starvation (Fig. 4E,  H).In contrast, the diapause-exits that emerged between late winter and early spring (February-May, ca.400 exits d − 1 ) were predominantly large CVs with higher energy reserves (Fig. 4F).These had post-diapause energy reserves that prevented starvation and aided the elevation of their structural mass to attain sexual maturity and mate, so that spawning commenced as soon as food levels increased in early March.

Stochastic model environments
3.1.2.1.Shorter-term stochasticity.Under shorter-term environmental stochasticity (Fig. 2D-F), the fraction of late-juvenile (CIV and CV) stages that entered diapause increased slightly (2%) relative to the simulation in the deterministic environmental setting (Fig. 4K, O and diapause index in Fig. 5).Although the timing of diapause entry and exit Fig. 6.Principal Component Analysis (PCA) biplots indicating the relationships between environmental variables and the emergent life history attributes (annual means) in simulations performed under shorter-and longer-term environmental stochasticity.Each point signifies a year in the simulated timeseries, indicated by subscript year numbers.A: Simulation with interannually consistent, modest predation risk.B: Simulation with variable visual predation risk.C: Simulation with variable non-visual predation risk.Emergent population dynamics of selected years in each of these environmental scenarios are presented in Figs.7-9.
were not markedly different, the structural mass and energy reserve levels of diapause-entries and diapause-exits differed significantly from the deterministic simulation (Figs.4K-N and 5).Particularly, the 2nd and 3rd quarters of diapause-entries, which endured the longest diapause duration were on average 20%-40% smaller in structural mass and carried ca.40% less lipid reserves (Figs.4L and 5).As a result, when they emerged from diapause between January and late June (i.e., 1st quarter of diapause-exits), they had on average < 9 µg C of reserves remaining, which was < 50% compared to those in the deterministic environment (Fig. 5).In contrast, the mean structural and energy reserve masses of the 1st quarter of diapause-entries that endured the shortest diapause duration were ca.20% higher than in the deterministic setting.Many of them emerged in the late summer and early autumn (2nd quarter of diapause-exits) with higher energy reserves (Fig. 5).These energy reserves allowed for reduced starvation, especially since feeding opportunities gradually deteriorated towards autumn (Fig. 2F) and these larger individuals had to frequently evacuate the productive near-surface waters as they performed diel vertical migrations to reduce predation risk (Fig. 4P).The mean size at sexual maturity of the simulated C. finmarchicus population decreased by ca.12% and the timing of egg production occurred ca. 12 d later in the stochastic model environment (Fig. 5).

Longer-term stochasticity.
In all model environments where longer-term (interannual) stochasticity was added to the shorter-term stochasticity, food concentration and the duration of the productive season were positively correlated to ambient temperature (altogether, the growth potential: eigenvectors I-III: Fig. 6, cf.Fig. 2G-I).Consequently, in warmer years, the timing of the pelagic bloom occurred earlier compared to colder years (eigenvector IV: Fig. 6).
When both visual and non-visual predation risks were modest and constant between years, the mean timing of diapause entry and the structural and energy reserve masses of diapause-entries were negatively correlated with the growth potential (eigenvectors IX, X and XI in Fig. 6A).This shows that in (warmer) years with higher growth potential, the simulated C. finmarchicus population entered diapause earlier in the year at smaller body size and with less reserves (Fig. 7C, D & K, L).Similarly, size at sexual maturity correlated negatively with the growth Fig. 7. Predicted dynamics of the simulated C. finmarchicus population in two years with different growth potential but with constant modest visual and non-visual predation risks subjectively selected from the simulated time series (6-h estimates).This is synonymous with the data presented in Figs.2G and 6A.DEN: diapause entry, DEX: diapause exit, DD: direct development (without diapause).Q1-Q4 are times at which 25%, 50%, 75% and 100% of the population enters or exits diapause and their structural (W c ) and energy reserve masses (W s ).potential (eigenvector XV: Fig. 6A).
The mean timing of diapause exit correlated with structural and energy reserve masses of diapause-exits, but not with the growth potential (eigenvectors XII, XIII, XIV: Fig. 6A).This is because the early diapause-exits, which ascended to near-surface waters during the springsummer transition were large developmental stages (mainly CVs) with high energy reserves, while those emerged during the autumn-winter transition were smaller CIVs and CVs with nearly depleted reserves (Fig. 7E, F & M, N).The diapause index correlated negatively with the growth potential (eigenvector VIII: Fig. 6A) and suggests that a higher proportion of late-juvenile stages remained in diapause until feeding opportunities emerged in the following productive season in (colder) years with lower growth potential (Fig. 7E, M).In (warmer) years with higher growth potential, most late-juvenile stages directly developed to adulthood and produced additional generation(s) (Fig. 7G & O).Although the timing of egg production (eigenvector XVI) correlated positively with the timing of diapause exit (Fig. 6A), the lack of a strong positive correlation (e.g., alike that between temperature and food concentration) suggests that non-diapausing late-juvenile stages that directly develop into adults contributed to a substantial fraction of the total egg production irrespective of the modelled interannual variations of growth potential.
Interannual stochastic variability in both visual or non-visual predation risks (eigenvectors VI and VII: Fig. 6) had little effect on the mean timing or structural and energy reserve masses at diapaus entry (Fig. 6B,  C, and panels C, D & E, K of Figs. 8 & 9).However, visual predation risk correlated positively with the mean timing of diapause exit (Fig. 6B), indicating that diapause of the simulated C. finmarchicus population terminated later in years with higher visual predation risk and vice versa (Fig. 8E, M).This predicted shift in the timing of diapause exit primarily accounts for the fraction of C. finmarchicus population that exits diapause during late summer and autumn because the mean timing of diapause exit always fell between late-summer and early-autumn, when the largest number of copepods exited diapause in all model simulations (Fig 8E , M). Non-visual predation risk correlated negatively with the timing of diapause exit (Fig. 6C).Therefore, in years with higher nonvisual predation risk, diapause exit of the simulated C. finmarchicus population occurred relatively earlier than in years with lower risk Fig. 8. Predicted dynamics of the simulated C. finmarchicus population in two years with different visual predation risk but with average growth potential and modest non-visual predation risk subjectively selected from the simulated time series (6-h estimates).This is synonymous with the data presented in Figs.2H and 6B.DEN: diapause entry, DEX: diapause exit, DD: direct development (without diapause).Q1-Q4 are times at which 25%, 50%, 75% and 100% of the population enters or exits diapause and their structural (W c ) and energy reserve masses (W s ).

Discussion
Two main diapause strategies emerged in the simulated C. finmarchicus population irrespective of the level of environmental variability: (i). a mid-to late spring diapause entry as CIVs or smaller CVs with limited energy reserves and a late summer to autumn diapause exit, which resembles an oversummering strategy (Wang et al., 2003), and (ii).a summer-autumn diapause entry, predominantly as CVs with larger reserve loads and a diapause exit during the winter-spring transition, which resembles a classic overwintering strategy (Sømme, 1934).As a result of both strategies being present in the population, diapause entry occurred throughout the productive season and diapause exit occurred asynchronously throughout the year in varying numbers.Additionally, and as a third strategy, a part of the population did not undergo diapause at all and instead developed directly to adults in the attempt to reproduce within the same productive season.The proportions of the diapausing and non-diapausing populations varied between years despite the equal probabilities forced at the diapausing and non-diapausing decision.This occurred due to the differential survival rates associated with the diapausing and non-diapausing strategies given the interannual stochastic variations of the model environment.

The curious case of oversummering
In high-latitude herbivorous copepods, overwintering is a state of dormancy that increases survival during unproductive winter periods where food is scarce (Sømme, 1934).Overwintering has been studied extensively since the pioneering work of Gran (1902).Oversummering is comparatively less well studied.Oversummering C. sinicus in the Yellow Sea is probably the most well-known among marine copepods (Pu et al., 2004;Wang et al., 2005), where CVs occupying the Yellow Sea Cold Water Mass enters diapause as the upper pelagial warms up and feeding opportunities deteriorate towards summer (Li et al., 2004).Recent studies point to an oversummering strategy among some C. finmarchicus populations in the Western North Atlantic.For example, an individual based modeling study of C. finmarchicus in the Gulf of Maine by Maps  et al. (2012b) predicted that a part of the simulated population that enters diapause in mid-summer returns to the upper pelagial in autumn while using a late pelagic algal bloom for feeding and reproduction.Their model-predicted oversummering timing and durations align well with the field data (Durbin et al., 1997;Saumweber and Durbin, 2006).However, the environmental dynamics (e.g., ambient temperature and the timing and duration of the productive season) of these Northwest Atlantic study locations considerably differ from those of the present model and direct comparisons to field data should be made with caution.In the Northeast Atlantic, Russell (1926) noted a mid-summer seasonal descent of C. finmarchicus populations of the English Channel and described it as a strategy to avoid harmful levels of summertime irradiance.Further, Kaartvedt (2000) highlighted that C. finmarchicus populations in the Norwegian Sea may enter diapause in summer primarily due to the elevated predation risk imposed by seasonally-migrating planktivorous fish.However, as none of these studies continued towards autumn, we do not know if these summer diapause entries had ascended back to the upper pelagial later in the same productive season.From the classic of Østvedt (1955) to more recent studies, such as Strand et al. (2020), we see CIV and CV stages remaining in the near-surface layers of the Norwegian Sea towards autumn.However, it is unclear whether these late developmental stages are a part of a non-diapausing multigenerational development or remnants of a brief summer diapause.
A main problem of entering diapause earlier in the year (late spring and early summer) is that these diapause-entries require substantial energy reserves to remain dormant until the following productive season.For example, a copepod that enters diapause in late May or early June at the deterministic model environment has to remain dormant for ca. 10 months before the pelagic primary production commences in the following year (Fig. 2C).This is unlikely for most copepods given the modelled diapause metabolic rate at 0 • C (Fig. 2B, Eq. ( 14)) and their lower energy reserve sizes (Fig. 4D).Instead, energy reserve exhaustion drove large numbers of these early diapause-entries to the upper pelagial between late-summer and autumn (Fig. 4E, F, M, N), where they successfully produced another generation since food was still available at the time (Figs.2C, F and 4B, J).In the model, this evolved into an adaptive strategy in both deterministic and stochastic settings, since the summertime escape from the upper pelagial reduced the mortality risk imposed by visual predation and there was sufficient food for the growth and development of their offspring until autumn.This somewhat resembles the summer diapause of freshwater zooplankton in response to vertebrate (Pijanowska and Stolpe, 1996) and invertebrate (Strickler and Twombly, 1975) predation, where resting eggs are produced prior to the elevation of seasonal predation risk-a well-studied conservative bet hedging strategy (García-Roger et al., 2017).Despite this interesting model prediction that highlights underlying trade-offs in the annual routine of these copepods, we are uncertain about how likely it is that such briefer summer descents are to be expected for Northeast Atlantic C. finmarchicus.This is because contrast to our model, summer-autumn food conditions in nature are far from guaranteed (due to grazing pressure and nutrient limitation), and consequently, most early diapause-entries that exit diapause in late summer and autumn with exhausted energy reserves usually starve to death without contributing to another generation (Pepin and Head, 2009).

Shorter-term environmental stochasticity
The shorter-term stochastic oscillations of ambient temperature and food concentration collectively induced growth potential variations between the simulated 6-h timesteps Fig. 2E, F and Eqs. ( 1)-( 6).However, diapause timing of the simulated C. finmarchicus population was robust to these shorter-term growth potential variations Fig. 5).This robustness was manifested through the plasticity of behavioral and morphological (morphometric) attributes.First, copepods abandoned their routine DVM behavior (that emerged in the deterministic environment: Fig. 4H) and spent more time in near-surface waters, especially during time intervals with higher simulated cloud cover to elevate foraging effort (Fig. 4P).Nonetheless, this did not fully compensate for the growth potential variability, which is indicated by the lower mean structural and energy reserve masses of late-juvenile CIV and CV stages compared to the deterministic setting (Figs.4D, L and 5).Since it takes lesser time to attain a lower body mass at a prescribed growth potential in this model (Eqs.( 1)-( 8), the lower structural and energy reserve masses of late-juvenile stages compensated for potential delays of diapause entry caused by environmental variability.However, this body size plasticity came at a significant cost in the form of reduced survival and reproductive output.At their springtime ascent, these smaller postdiapause stages had little energy reserves (Figs. 4F,N and 5).This reduced the survivorship among the earliest diapause-exits that ascended to near-surface waters preceding the onset of the pelagic primary production.Further, as the attainment of sexual maturity generally requires a significant elevation of the structural mass from late-juvenile copepodite stages (Eq.( 8)), the relatively smaller diapause-exits took a longer time to develop towards sexual maturity and caused the 12d delay in the mean timing of egg production (Fig. 5).However, a further delay in egg production did not occur as the C. finmarchicus reached sexual maturity at ca. 12% lesser size (Fig. 5).This caused ca.20% drop in egg production since the fecundity of the modelled copepods is proportional to their body size (Eq.( 17)).

Longer-term environmental stochasticity
4.2.2.1.Diapause entry.The mean timing of diapause entry was highly plastic to the modelled interannual variations of growth potential and occurred earlier in warmer years with earlier onset and longer duration of pelagic primary production and vice versa (Figs.6A, cf.Fig. 7).In a ca.two-decade timeseries study, Espinasse et al. (2018) reported a similar phenological shift in the seasonal appearance of CVs among two distinct C. finmarchicus populations inhabiting northern Icelandic and western Spitsbergen shelves.Further, a five-year study on a co-occurring copepod (C.glacialis) reported ca.3-week delay in the appearance of lipid-accumulated late-juvenile stages during abnormally colder years when the pelagic primary production occurred relatively late in the White Sea (Pertsova and Kosobokova, 2010).Growth potential implications on the diapause entry of C. finmarchicus are also found along spatial gradients.For example, Melle et al. (2014) reported a notable (30-60 d) delay in C. finmarchicus diapause entry in the northwest Atlantic compared to northeast Atlantic.The earlier onset of primary production in the warmer upper pelagial of the northeast Atlantic generally warrant faster growth rates and a large fraction of C. finmarchicus population enters diapause between late May and July.In contrast, the pelagic primary production of the northwest Atlantic occurs relatively late and slower growth rates sustained in the colder upper pelagial generally delays their diapause entry.Similarly, in colder seasonally ice-covered high-Arctic fjords where the period of net positive pelagic primary production commences in mid-June and is limited in duration to 2-3 months, the timing of C. finmarchicus diapause entry can shift to August-September (Arnkvaern et al., 2005;Bandara et al., 2016).
It has been hypothesized that elevated predation risk imposed by planktivorous fish drive parts of the C. finmarchicus populations in the Norwegian Sea to diapause depths by mid-summer, thus negating the possibility to produce an additional generation or two in the autumn Kaartvedt, 2000;Varpe and Fiksen, 2010).In contrast, our predecessor model (Bandara et al., 2019), which assumed an annual life cycle for C. finmarchicus at ~70 • N, highlighted the possibility of elevated visual predation risk to cause a relative delay in the timing of diapause entry due to the reduced growth rates induced by diel vertical migration.Consequently, the modelled copepods required more time to develop into CIVs or CVs with sufficient energy reserves for overwintering.However, in the present model, the visual predation risk did not significantly influence (i). the predicted mean timing of diapause entry, (ii).the mean structural and energy reserve masses of diapause-entries, (iii).generation time or (iv).the extent of summer diapause (Fig. 6B, cf.Fig. 8).In this model, visual predation risk varied interannually following a uniformly random probability distribution (Fig. 2H).In the studies of Varpe and Fiksen (2010) and Bandara et al. (2019), a higher predation risk was sustained or assumed across multiple generations.For example, Bandara et al. ( 2019) maintained elevated visual predation potential for ca.400 generations, thus maintaining a constant top-down selection pressure over an adequate time period for life histories of the simulated copepods to evolve.Since the top-down selection pressure of the stochastic versions of the present model varied randomly between years, it is unlikely that life history responses, such as smaller body masses and reduced energy stores are evolved in higher predation risk environments.This adaptive evolution (adaptive tracking) of body size contrasts the body size plasticity predicted under uniformly random interannual fluctuations of bottom-up selection pressures (i.e., temperature and food availability: Fig. 7), where the response of copepods' body size (and hence the reserve size) is near-instantaneous due to the strict food and temperature dependence of the growth and developmental formulations (Eqs.( 1)-( 8).

Diapause exit.
The proximate drivers of the diapause exit and seasonal ascent of high-latitude copepods are not well-understood and an array of internal and external cues are hypothesized (reviewed in Bandara et al., 2021;Baumgartner and Tarrant, 2017;Miller et al., 1991).According to recent studies, lipid-modulated endogenous timers play a crucial role in diapause exit of C. finmarchicus (Häfker et al., 2018;Skottene et al., 2019).Similarly, the modelled proximate driver of diapause exit in this study is the exhaustion of energy reserves (represented by the 'gene' ε).This, in association with the lipid-based diapause entry attribute, (i.e., 'gene' δ, Table 1) created more or less a continuum of diapause exit strategies that spanned asynchronously throughout the year in both deterministic and stochastic model environments (e.g., Fig. 4E, M).Observations on asynchronous diapause exit and seasonal ascent of Calanus spp.have been noted in several high-latitude year-round zooplankton investigations, particularly fjord environments (e.g., Bandara et al., 2016;Darnis and Fortier, 2014;Dezutter et al., 2019;Wallace et al., 2010).Further, wintertime upper pelagial occupation of Calanus spp. is documented across numerous high-latitude locations, such as western and northern Norwegian fjords (Balinö and Aksnes, 1993;Falkenhaug et al., 1997), open waters of the Norwegian and Barents Seas (Bathmann et al., 1990;Heath and Jónasdóttir, 1999;Pedersen et al., 1995), Greenland Sea (Hirche, 1991) and in fjords on the Svalbard archipelago (Berge et al., 2009;Darnis et al., 2017).It is unclear why Calanus spp.ascend to the upper pelagial several months prior to the commencement of pelagic primary production.Some studies are insightful of the ability of winter-active predators to drive diapausing populations out of the deeper pelagial (Błachowiak-Samołyk et al., 2015;Sims et al., 2003).Results of the simulations performed under variable non-visual predation risk also points towards the same direction, where the diapausing C. finmarchicus population tended to ascend earlier at higher levels of non-visual predation risk (Fig. 6C, cf.Fig. 9).However, a major contrast between the simulated Calanus population and those in the field is that when driven to the upper pelagial, the latter may survive the unproductive winter given their omnivorous feeding strategy (Cleary et al., 2017;Ohman et al., 1998;Søreide et al., 2008).
In this model, most individuals that exited diapause and ascended to the upper pelagial during the unproductive season (October-February) died from starvation as they were modelled as strict herbivores .It may be questioned how such a 'wasteful' strategy that causes mass wintertime mortality becomes prominent in a model of artificial evolution.When the mutation probability of 'gene' ε was reduced to 0.02 (1/10th of the normal value) in a trial simulation performed in a deterministic model environment, the year-round diapause exit largely diminished and was limited to the spring and mid-summer (data not presented).However, when the same modification was attempted in a trial simulation performed in an environment with shorter-and longerterm stochasticity, the model collapsed within a few years due to the inability of copepods to maintain a viable population over the simulated timeseries.This suggests that the year-round diapause exit was crucial to the year-to-year survival of the simulated C. finmarchicus population, especially since the copepods in diapause cannot predict the time at which the pelagic primary production starts in the following year.
The asynchronous year-round diapause exit predicted in this model represents a diversified bet hedging strategy, where parents (parental genotypes), through random mating, recombination and mutation, produce a series of offspring with diverse diapause exit strategies (phenotypes), across which the mortality (starvation) risk of ascending to the upper pelagial at the 'wrong' food-deprived time of the year is spread (Philippi and Seger, 1989;Slatkin, 1974).Similar diversified bet-hedging strategies are common among many freshwater zooplankton (e.g., rotifers and cladocerans) that undergo diapause as resting eggs.Here, parents produce resting eggs through sexual reproduction, which sinks to the sediment and remain viable for many years and act as 'germ banks' from which the offspring emerge asynchronously in time (Brendonck and De Meester, 2003).The population-level consequences of such bet hedging strategies are paramount, particularly because synchronous hatching of resting eggs can be detrimental to the population when there are pronounced stochastic variations in the environment and if there are no widespread proximate cues to predict those variations (Hairston Jr, 1996).Therefore, the asynchronous diapause exit of marine zooplankton predicted by our model and the asynchronous hatching of resting eggs from egg banks seems to be analogous strategies selected for by environmental variability.As more detailed in-situ data are about to become increasingly available via autonomous vehicles (gliders and floats equipped with imaging devices, e.g., Lombard et al., 2019), this predicted outcome could soon be tested in various ecosystems around the world.

Conclusions
Most marine zooplankton live in a near-continuous flow and faces uncertainty in foraging and predator encounter probabilities in the shorter-term (Seuront et al., 2004;Visser et al., 2008).Along with longer-term interannual variations of temperature, food and predation environments, they occupy a world of highly variable bottom-up and top-down selection pressures across time and space (Lurling and de Senerpont Domis, 2013).Understanding adaptations of zooplankton to environmental heterogeneity is not an easy task.However, mechanistic simulation models provide a cost-effective supplement to field studies and laboratory experiments in this regard.The present model extends the capabilities of our predecessor models and sheds new light on the timing and body size dynamics of the high-latitude copepod C. finmarchicus against stochastic environmental variations that occur atop predictable cyclic diel and seasonal patterns (see also Fiksen, 2000).Although spatial environmental heterogeneity is not explicitly formulated in this model, its computationally efficient simulation framework opens future opportunities to couple this unidimensional individual based biological model to a three-dimensional ocean circulation model towards generating predictions that can directly be validated with field data.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could potentially influence the research reported in this paper.
model runs in a deterministic environment (Fig.2A-C).To simulate shorter-term environmental stochasticity, we modified the deterministic model environment by introducing uniform random variations in the range of ±25% to temperature and food concentration at each 6-h time interval (Fig.2E, F).Similarly, uniform random variations were introduced to irradiance in the range of -90% and 0% to simulate the attenuation of incident irradiance due to cloud cover (Fig.2D).To simulate the longer-term (interannual) environmental stochasticity, we introduced uniform random variations to the maximum sea surface temperature in the range of ± 3 • C, which in turn drives the modelled maximum food concentration in the range of ± 45 µg C l − 1 and the timing and duration of thermal stratification and pelagic primary production in the range of ± 75 d (Fig.2G-I).The above data ranges were defined based on (i).environmental data collected by autonomous surface and underwater vehicles and FerryBoxes operating in the northern Norwegian sea (68 • -71 • N) between 2018 and 2019 and (ii).ERA-Interim reanalysis archives of European centre for Medium-Range

Fig. 2 .
Fig. 2. Modelled diurnal and seasonal variation of irradiance, temperature and food concentration in the deterministic (A-C) and shorter-term (6-h) stochastic (D-F) environments.Panels G-I show the modelled longer-term (interannual) stochasticity in temperature, food and visual (scalar K, Eq. (10)) and non-visual predation (M n ) environments used in three simulation experiments.Temperature, food concentration and duration of the productive season are shown as mean deviations.Warmer years in reddish hue and colder years in blueish hue-synonymous with color representation in Fig. 6.

Fig. 4 .
Fig. 4. Predicted annual dynamics of estimated population size (A, I), developmental stage composition (B, J), number of diapause entries (DENs: C, K), number of diapause exits (DEXs: E, M) and their structural (W c : D, L) and energy reserve masses (W s : F, N), number of direct-developing (non-diapausing) individuals (DDs: G, O) and the vertical distribution of the population (H, P) in the deterministic (left panels) and shorter-term stochastic environments (right panels) (6-h estimates).gray-shaded regions shows the time at which 25%, 50%, 75% and 100% (Q1-Q4) of the population enters and exits diapause and their structural and energy reserve masses.Relative variations of food concentration and temperature (excluding shorter-term variability) at the surface (z = 0 m) are shown on the top (cf.Fig. 2).

Fig. 9 .
Fig. 9.Predicted dynamics of the simulated C. finmarchicus population in two years with different non-visual predation risk but with average growth potential and modest visual predation risk subjectively selected from the simulated time series (6-h estimates).This is synonymous with the data presented in Figs.2I and 6C.DEN: diapause entry, DEX: diapause exit, DD: direct development (without diapause).Q1-Q4 are times at which 25%, 50%, 75% and 100% of the population enters or exits diapause and their structural (W c ) and energy reserve masses (W s ).

Table 1
List of definitions, values and units of the terms used in the model.