Multiple‐batch spawning as a bet‐hedging strategy in highly stochastic environments : an exploratory analysis of Atlantic cod

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Multiple‐batch spawning as a bet‐hedging strategy in highly stochastic environments : an exploratory analysis of Atlantic cod © 2021 the Authors Published version Hočevar, Sara; Hutchings, Jeffrey A.; Kuparinen, Anna Hočevar, S., Hutchings, J. A., & Kuparinen, A. (2021). Multiple‐batch spawning as a bet‐hedging strategy in highly stochastic environments : an exploratory analysis of Atlantic cod. Evolutionary Applications, 14(8), 1980-1992. https://doi.org/10.1111/eva.13251 2021 1980 | Evolutionary Applications. 2021;14:1980–1992.wileyonlinelibrary.com/journal/eva 1 | INTRODUCTION Natural environmental conditions change continuously across multiple spatial and temporal scales. While some of the fluctuations can be predictable and of low magnitude, others are uncertain, occurring with varying degrees of intensity, periodicity and stochasticity. Stochasticity can be integral to the shaping of genotypes, phenotypes and populations, either concomitantly or not Received: 18 October 2020 | Revised: 30 April 2021 | Accepted: 3 May 2021 DOI: 10.1111/eva.13251 ORIGINAL ARTICLE Multiplebatch spawning as a bethedging strategy in highly stochastic environments: An exploratory analysis of Atlantic cod Sara Hočevar1 | Jeffrey A. Hutchings1,2,3,4 | Anna Kuparinen1 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2021 The Authors. Evolutionary Applications published by John Wiley & Sons Ltd 1Department of Biological and Environmental Science, University of Jyväskylä, Jyväskylä, Finland 2Department of Biology, Dalhousie University, Halifax, NS, Canada 3Institute of Marine Research, Flødevigen Marine Research Station, His, Norway 4Department of Natural Sciences, University of Agder, Kristiansand, Norway Correspondence Sara Hočevar, Department of Biological and Environmental Science, University of Jyväskylä, Jyväskylä, Finland. Email: sara.s.hocev[email protected] Funding information Academy of Finland, Grant/Award Number: 317495; Natural Sciences and Engineering Research Council of Canada; European Research Council, Grant/Award Number: COMPLEXFISH 770884 Abstract Stochastic environments shape lifehistory traits and can promote selection for riskspreading strategies, such as bethedging. Although the strategy has often been hypothesized to exist for various species, empirical tests providing firm evidence have been rare, mainly due to the challenge in tracking fitness across generations. Here, we take a ‘proof of principle’ approach to explore whether the reproductive strategy of multiplebatch spawning constitutes a bethedging. We used Atlantic cod (Gadus morhua) as the study species and parameterized an ecoevolutionary model, using empirical data on sizerelated reproductive and survival traits. To evaluate the fitness benefits of multiplebatch spawning (within a single breeding period), the mechanistic model separately simulated multiplebatch and singlebatch spawning populations under temporally varying environments. We followed the arithmetic and geometric mean fitness associated with both strategies and quantified the mean changes in fitness under several environmental stochasticity levels. We found that, by spreading the environmental risk among batches, multiplebatch spawning increases fitness under fluctuating environmental conditions. The multiplebatch spawning trait is, thus, advantageous and acts as a bethedging strategy when the environment is exceptionally unpredictable. Our research identifies an analytically flexible, stochastic, lifehistory modelling approach to explore the fitness consequences of a riskspreading strategy and elucidates the importance of evolutionary applications to lifehistory diversity. KEYWORDS Atlantic cod, bethedging, environmental stochasticity, fitness, multiplebatch spawning, riskspreading | 1981 HOČEVAR Et Al. (for a detailed review, see Lenormand et al., 2009). Environmental stochasticity influences the ecological processes of a population, determines the rate and direction of its evolutionary change (Frank & Slatkin, 1990; May, 1973) and can even lead to its extinction (Lande, 1993). Hence, to mitigate challenges arising from prevailing environmental uncertainty, organisms have evolved a diversity of lifehistory strategies (Kussell & Leibler, 2005; Maynard Smith, 2010; Meyers & Bull, 2002; Moran, 1992; Stearns, 1976; Tufto, 2015). One of these is bethedging (Gillespie, 1973, 1974, 1975; Slatkin, 1974). Bethedging is a costly genotypic strategy that maximizes longrun or geometric mean fitness across generations by trading off the arithmetic mean in reproductive output and its variance (Cohen, 1966; Gillespie, 1974; Lewontin & Cohen, 1969; Seger & Brockmann, 1987; Simons, 2002; Yoshimura & Clark, 1991). In other words, bethedging can act as a ‘portfolio effect’ (Markowitz, 1952) through which the diversification of assets, here partitioning of offspring among batches, reduces the risk and stabilizes the returns, that is geometric mean fitness of a genotype. Organisms can spread risk among their offspring on a temporal or spatial scale, in a conservative or diversified way, or even as a complex combination of all the above (Haaland et al., 2019; Scheiner, 2014). While conservative bethedging maximizes fitness by reducing the variance in fitness at the individual level, diversifying bethedging does so by reducing the correlation in expected fitness among individuals in the same population (Starrfelt & Kokko, 2012). Examples of bethedging strategies appear in a wide range of systems and forms (Childs et al., 2010; Philippi & Seger, 1989) such as iteroparity (Cole, 1954; Ranta et al., 2002), seed dormancy (Cohen, 1966; Simons, 2009), seed dispersal (Beckman et al., 2018; Snyder, 2011), flowering schedule (Simons & Johnston, 2003), timing of sexual reproduction (Tarazona et al., 2017), embryonic diapause (Furness et al., 2015) and hatching asynchrony (Laaksonen, 2004). Simons (2011) extensively reviewed over 100 studies on bethedging, categorizing them based on the strength of the empirical evidence. Although bethedging life histories have been reported for a variety of species and hypothesized for even more, from bacteria (Beaumont et al., 2009) to vertebrates (Lips, 2001; Mahony & Thumm, 2002), the strength of the evidence for most has been limited or, as Simons (2011) put it, the evidence has been elusive. He proposed six, ranked evidence conditions that need to be met: (I) recognize a bethedging trait; (II) monitor the unpredictable environment; (III) observe differences in the trait among populations; (IV) demonstrate differences in fitness dynamics; (V) validate whether the trait is favoured under relevant varying environments; and (VI) test the optimality of the trait under a range of conditions of fluctuating selection (Simons, 2011). Few studies possess sufficient empirical support to fulfil the highest three and most datademanding categories of evidence for bethedging (i.e. categories IV– VI; Simons, 2011). The majority of those that do fulfil these conditions are on plants (Childs et al., 2010; Simons & Johnston, 2003). This general lack of evidence can be attributed to the very considerable challenges of recognizing the adaptive significance of a trait that is bethedged and the difficulty of tracking acrossgenerational fitness in a stochastically fluctuating environment. In our study, we attempt to overcome these challenges and aspire to provide support for or against the fifth evidence category on bethedging significance of multiplebatch spawning strategy. Multiplebatch spawning is a reproductive strategy common among marine fishes, such as gadoids and flounders, for example haddock (Melanogrammus aeglefinus), pollock (Pollachius virens), whiting (Merlangus merlangus), halibut (Hippoglossus hippoglossus) and dab (Limanda limanda) (Murua & SabridoRey, 2003). Yet, the fitness benefits of the multiplebatch strategy have not been comprehensively explored. To sustain population resilience in a stochastic environment, bethedging could be crucial for multiplebatch spawning fish populations. A strategy of broadcast spawning on multiple spawning grounds, multiple times (Kjesbu, 1989) and over prolonged periods (Hutchings & Myers, 1994; Kjesbu et al., 1996) might act as a portfolio effect by reducing the risk of complete reproductive failure. The production of multiple egg batches within a spawning season, the number of which increases with female weight and body size (Kjesbu et al., 1996; Roney et al., 2018), could enable a batch spawner to spread the environmental risk among its offspring and mitigate the fitness consequences of environmental fluctuations. As a tradeoff in diversification, the variance in reproductive output of a multiplebatch spawner could be lower, boosting the acrossgenerational geometric mean fitness, at the expense of producing a lower average number of offspring. To tackle the question of whether multiplebatch spawning yields the predicted fitness benefits of a bethedging strategy, we used Atlantic cod, G. morhua (Linnaeus, 1758), as a focal species in this study. Atlantic cod, one of the most studied batch spawning fish species, has been speculated to be a conservative bethedger (e.g. Hutchings & Rangeley, 2011), but never in fact tested for it. Here, we test this hypothesis by expanding an ecoevolutionary model parameterized for cod (Kuparinen et al., 2012). Our primary objectives are to (i) observe how multiplebatch spawning affects populational dynamics; (ii) evaluate the fitness consequences of multiplebatch spawning within a spawning season, under different levels of environmental stochasticity; (iii) inspect the variance in reproductive output within generations; and (iv) analyse the proportion of successful spawning seasons. 2 | MATERIALS AND METHODS There can be several risk distribution strategies acting on different stages or processes in a species at any given time. This complication has potential to obscure the fitness consequences of any one component of the bethedging strategy (Simons, 2011). Thus, we focused solely on the component of multiplebatch spawning. 1982 | HOČEVAR Et Al. 2.1 | Multiplebatch spawning and environmental stochasticity We examined whether multiplebatch spawning constitutes a bethedging strategy by exploring its ecoevolutionary impacts on fitness dynamics under varying levels of environmental stochasticity affecting batch survival. We did so by implementing an individualbased mechanistic model developed by Kuparinen et al. (2012) which characterizes the ecoevolutionary dynamics and demographic processes of Atlantic cod. The main evolving trait of the model is body size, which fits our research design since batch production, fecundity and spawning duration are sizerelated traits. Given that the model's configuration and parameterization have been thoroughly described elsewhere (Kuparinen et al., 2012, 2014), we outline below only the main features (Table 1). Here, we focus on a detailed description of newly implemented batch spawning strategies and components of generated environmental stochasticity in batch survival. We simulated the fecundity of every mature female at the start of each spawning season through juvenile production and survival. The production of eggs was a positively dependent function of the female's weight derived from the empirically based length (L) – w e i g h t ( W) relationship W t=3.52 ⋅10 −6 ⋅L 3.19 t . The egg number was calculated based on the empirical relationship as eggs = ( 0.48 ⋅ ( Wt+0.37 ) 1.45 +0.12 ) ⋅10− 6 (Hutchings, 2005; Kuparinen et al., 2012). We portrayed the lifehistory strategy of multiplebatch spawning cod in the model as the occurrence of multiple reproductive events within a single spawning season. To highlight how the strategy feeds back on the population dynamics of cod, we compared it to a separately simulated hypothetical population with no such riskspreading strategy, that is singlebatch spawners. While individuals in a multiplebatch spawning population distribute their eggs among several batches, individuals in a singlebatch spawning population only have one reproductive event within each spawning season and, thus, practise the tactic of ‘placing all your eggs in one basket’ (Figure 1). The number of eggs that a multiplebatch spawning female produced in one spawning season was distributed evenly among batches. The number of batches depended on body size. At each annual step, the number of batches for each mature female was derived from an empirically based nonlinear regression fit between the fork length Lfork( t ) and batch number Nbatches , where N batches = 21.156 1+exp ( 55.014 −Lfork(t) 10.141 ) (parameterized based on an empirical data set collected from Risør fjord in coastal Skagerrak by Roney et al., 2018; Figure 2). Before the nonlinear regression was fitted, the lower value of 0 batches at size 25 cm and higher value of 21 batches at 100 cm were added to the experimentally gathered data to account for somatic constraints. We based these constraints on data of maximum observed batch number in captive Norwegian coastal cod (Kjesbu et al., 1996) and limited the function by setting the maximum available number of batches to 21 to prevent the continuous increase of produced batches with a female's size. TABLE 1 Summarizing the main underlying empirically derived variables of the ecoevolutionary model Variable Description Equation Value Unit Source LLength calculated each year for every individual following the Von Bertalanffy growth curve L(t) =L∞− ( L∞−L 0) ⋅e− kt log (k)=−0.609 −0.013 ⋅ L∞ L(t) cm von Bertalanffy (1938) Kuparinen et al. (2012) WWeight calculated each year for every individual from the lengthweight relationship Wt =3.52 ⋅10 −6 ⋅L 3.19 t W(t) kg Kuparinen et al. (2012) Lmat Length at maturity Lmat = 0.66 ⋅ L∞ Lmat cm Jensen (1997) Nbatches Number of spawned batches calculated each year for every mature female N batches = 21.156 1+exp ( 55.014 −Lfork(t) 10.141 ) Nbatches(t) Roney et al. (2018) Neggs Number of spawned eggs, calculated each year for every mature female N eggs = ( 0.48 ⋅ ( Wt+0.37 ) 1.45 +0.12 ) ⋅10− 6 Neggs(t) Hutchings (2005) CostsMBS Costs of multiplebatch spawning strategy, calculated each year for every mature female spawning more than one batch Costs MBS = N batches ∑ batch =1 1−0.00523 ⋅(batch −1 ) CostsMBS(t) Roney et al. (2018) | 1983 HOČEVAR Et Al. To make multiplebatch spawning a costly trait, we introduced costs to batch survival (rightsided bar plot on Figure 2). We set these costs based on empirical findings, following the spawning dynamics of 73 wildcaught Norwegian coastal cod in Skagerrak and their offspring quality (Roney et al., 2018). Larval length and yolksac volume of offspring spawned in an experimental spawning basin at the Institute of Marine Research Flødevigen exhibited declining trends during the spawning period. Given that later spawned batches produced smaller larvae at hatch, we added an assumption that shedding of the first or single batch had no associated costs, while for every consecutive batch, shed within the same season, the batch spawning costs increased. Because larvae length tends to correlate with survival probability, we applied the tradeoff in costs of batch production and riskspreading potential in a gradually decreasing survival probability of each batch from 1.00 to 0.89 for the first to 21st produced batch, respectively (Roney et al., 2018; following the mortality function described by Pepin, 1991). Therefore, while singlebatch spawners experienced no such batch spawning costs, as they could shed only one batch per spawning season, the average costs of multiplebatch spawners varied depending on how many batches an individual has shed in a spawning season. We introduced environmental stochasticity to batch survival as the environmental pressure (0.05, 0.10, 0.15, 0.20, 0.25) and environmental fluctuations (0, 0.0001, 0.001, 0.01) and varied them in a full factorial manner (Figure 1). We define environmental pressure as a change in the mean batch mortality rate that could be driven, for example by any combination of ecological, environmental or anthropogenic factors. On the other hand, environmental fluctuations are defined as a change in the variance of batch mortality rate and can be generated by the stochasticity about the ecological, environmental or anthropogenic factors. Hence, we separately exposed multiplebatch and singlebatch spawning populations to each of 20 simulated environmental scenarios and compared their fitness performances. Higher rates of environmental pressure and fluctuations were dismissed upon trial testing, since the stochasticity became overwhelming and drove the populations simulated in our study towards extinction. When the added environmental fluctuations were preset to zero, the scenario was considered nonstochastic. Under these circumstances, we applied a constant environmental pressure to a batch or batches as a success probability in a Bernoulli trial to determine the survival of an entire batch or group of batches. While the environmental pressure under such nonstochastic scenarios remained constant, the outcome of successful survival of a batch or batches could still vary among seasons and individuals as it was newly drawn in every spawning season for every mature female. By contrast, the stochastic scenario was characterized by the presence of a variance in the form of continuous environmental fluctuations around mean rates of environmental pressure (Figure 1). The final environmental stochastic rate applied to batches was preadjusted, based on the rate of environmental pressure and environmental fluctuation, using a beta distribution parameterized by the α and β shape parameters (see Tables S1 and S2), and was drawn for every batch in every spawning season. To derive the final survival outcome of each batch, we applied the final environmental stochastic rate as a success probability in a Bernoulli trial to each batch separately, meaning that, for each batch, we drew a random number (0 or 1) as to whether the batch either dies or survives (i.e. a predator or environmental disaster destroyed a whole batch). Subsequently, to determine the final number FIGURE 1 The schematic diagram demonstrates the multiplebatch and singlebatch spawning cod populations simulated under a stochastic environment with varying rate of environmental pressure (0.05– 0.25) and environmental fluctuations (0– 0.01). The first graph is illustrating the probability density function for the beta distribution where mean environmental pressure applied to batch survival equals 0.20 and fluctuates depending on the environmental fluctuation rate (0, 0.0001, 0.001, 0.01). Correspondingly, the second graph is demonstrating the random, betagenerated batch mortality rates, drawn from each probability density function 1984 | HOČEVAR Et Al. of offspring, we summed the number of eggs from survived batches for every multiplebatch spawning female and multiplied the sum with a natural survival rate from the egg stage to 3yearold recruit estimated to be 1.13 ⋅10−6 for northern cod (Hutchings, 2005). The same process was adopted for the singlebatch spawning population. Thus, we created a combination of 40 different scenarios, comprising multiplebatch spawning populations under 20 distinct conditions of environmental stochasticity and singlebatch spawning populations under 20 distinct conditions of environmental stochasticity (Figure 1). 2.2 | Mechanistic model of Atlantic cod The mechanistic model (Kuparinen et al., 2012) follows the life stages of each individual fish in a population at annual time steps and combines genetic and optimization approaches through the use of heritable growth trajectories. The trajectories were derived from least square fits of empirically gathered 258 cod growth trajectories (Kuparinen et al., 2012), using the von Bertalanffy growth model L( t ) =L ∞ − ( L ∞ −L 0) ⋅e −kt , where L( t ) is the length of a fish at age t , L∞ is the asymptotic body length, L0 is the length at t=0 and k (year − 1) is the growth coefficient which describes the rate at which L∞ is reached (von Bertalanffy, 1938). Two observed associations underpin the model: (i) the observed negative correlation between L∞ and k , where log (k)=−0.609 −0.013 ⋅ L∞ , and (ii) the ratio of the length at maturity Lmat and L∞ , where Lmat =0.66 ⋅ L∞ (Jensen, 1997) when 30 cm ≤ L∞ ≤ 120 cm (Kuparinen et al., 2012). Each individual carried a genotype of 10 unlinked, diploid loci with 2 alleles (0 and 1) that followed classical Mendelian inheritance. The sum of these 10 loci, that could range from 0 to 20, coded for the genotypic value of L∞ and, thus, allowed for evolution of growth to occur. Ten loci were sufficient in describing the trait distribution smoothly; adding additional loci did not affect the simulations. Final phenotypic value, generated as an environmental variation (s.d. = 3.5) around the genotypic trait value, coded for the phenotypic L∞ value that provided a basis for the estimation of other relevant sizebased traits. To initiate the external fertilization process, a mature male was randomly assigned to a mature female, and the sex of offspring was determined by a 50/50 Bernoulli trial (Kuparinen et al., 2012). In addition to the demographic processes of reproduction and survival described in the batch spawning component of the model, densitydependent growth and natural mortality were also simulated at each annual time step on the individual level, and the state of the population was tracked accordingly. Growth of each individual was defined by its von Bertalanffy parameters ( L∞ and k ) and additionally altered by densitydependent population dynamics. Time available for growth at each annual step was bounded between 0 and 1. If the population reached or exceeded the carrying capacity, the individual's growth time and progress along its growth trajectory was reduced in accordance with a logistic equation e 15 −17.6 ⋅c ( 1+e15 −17.6 ⋅c )−1 , where c is the ratio between the biomass of the population and its carrying capacity. Populations of multiplebatch and singlebatch spawners had the same preset carrying capacity in all scenarios. At low population density, the individual's time spent on growth within one year was close to 1, allowing an individual almost a full annual growth increment along its von Bertalanffy curve (see Figure S1). Therefore, the population density affected fecundity by regulating the individual's growth, which impacts (i) the time the individual needs to reach 66% of its asymptotic length and mature and (ii) the age when reproduction starts. An instantaneous rate of natural mortality rate of 0.15, which was not applied until individuals had reached 3 years of age (see above), was assumed to be equal for all individuals of age 3 years or older (Kuparinen et al., 2012). If the individual was mature, the mortality rate was additionally increased by 0.10 to account for the survival cost of reproduction (following Kuparinen et al., 2012), resulting in an instantaneous rate of 0.25, which corresponds with the estimated natural mortality of many cod populations (Beverton et al., 1994). Using the binomial distribution, the model simulated the survival of every individual at each annual step with the maximum lifespan set to 25 years. 2.3 | Simulation design To achieve reproducibility of the code and results, we initialized a pseudorandom number generating sequence in a repeatable manner FIGURE 2 Relationship between mature female size, batch number and batch survival probability. The empirically gathered data (black) on the abundance of batches were plotted against the fork length of mature female cod and added constraints (red). A sigmoid curve given by a solid line was fitted to the data set and used in the simulation process. Vertical blueshaded area indicates the female size smaller than 250 mm, which were not considered mature in our model. Greyshaded bar chart on the right side of the figure illustrates the batch survival probability assigned to every batch according to the order at which female has shed it in the spawning period (values based on Roney et al., 2018). As a result, the batch survival probability is decreasing with increasing order number at which batch is shed (yaxis) | 1985 HOČEVAR Et Al. before each scenario run in R software (R Core Team, 2019) and allowed the loop to iterate and produce 50 replica simulations. As bethedging is predicted to have the greatest benefits in stochastic environments (Simons, 2011), we produced several combinations of runs to generate quantitative data (Figure 1). In the run of each scenario, we initialized the simulations with a preadapted cod population to separately simulate a population consisting of only multiplebatch spawners and a population consisting of only singlebatch spawners under each of the environmentally stochastic scenarios for 5000 years. This time interval was sufficient as it was beyond the time needed for the populations to reach their dynamic ecoevolutionary equilibriums (Figure S2). To investigate the potential underlying feedbacks of the bethedging strategy on life histories, we looked into the dynamics of the last 2500 years of the population in each run by recording the population variables and lifehistory traits at each annual step. The primary recorded output data were total number of recruits produced in one year, and annual population averages of L∞ and k , abundance, biomass and mortality rate. These outputs of 50 replica simulations per run were then summarized across replicates by recording the mean, coefficient of variance and standard deviation of each variable in every year. For the last 300 years of each simulation process, we tracked the individual fitness of every mature fish and recorded its total number of successfully shed batches, batch size and realized reproductive output over the individual's lifespan. The total reproductive output of every mature individual was recorded as the sum of the realized number of 3yearold recruits produced in the lifetime of the original individual. To explore whether multiplebatch spawning is a bethedging strategy of Atlantic cod, we looked into acrossgenerational fitness elements of multiplebatch spawning populations under simulated environmentally stochastic scenarios and compared them to a singlebatch spawning population. Acrossgenerational fitness elements included the arithmetic mean fitness WAM , variance in arithmetic mean fitness among generations or cohorts, and geometric mean fitness WGM . We measured the arithmetic mean fitness WAM and its variance across generations for each scenario as an average realized lifetime reproductive output of a generation and the acrossgenerational geometric mean in fitness WGM as the nth root of the product of average realized lifetime reproductive output of every generation following WGM = ( W AM1 ⋅W AM2 ⋅…⋅W AM n )1∕n (Seger & Brockmann, 1987), where n is a number of generations or cohorts. Fitness outputs were pooled together per run, and the mean, variance and coefficient of variation of each variable were recorded. Statistical analyses of relationships and trends were done using Welch's twosample t test (Welch, 1947), nonparametric Kruskal– Wallis test (Kruskal & Wallis, 1952) and simple linear regression (Kenney & Keeping, 1962). All simulations and analyses were performed in the opensource software R (R Core Team, 2019), and figures were produced using the package tidyverse (Wickham et al., 2019). 3 | RESULTS 3.1 | Population dynamics of each of the spawning strategists Separately simulated multiplebatch and singlebatch spawning cod populations exposed to 20 variations of environmentally stochastic scenarios (Figure 1) were analysed to investigate whether the costly spawning of multiple batches might be considered a bethedging strategy in Atlantic cod. In all 2000 simulations (2 spawning strategies * 20 different environmental scenarios * 50 replica simulations of each spawning population), the populations adapted from initial conditions to the new specific environments in fewer than 2500 years, maintaining a stable dynamic thereafter (Figure S2). Mean population size for all scenarios was 5311 individuals (s.d. 623 individuals across scenarios) and exhibited a declining trend (up to a 33.10% decrease) as the environmental fluctuations and pressure on batch survival or batch mortality rate increased. Coefficient of variation in population size was greater overall for singlebatch spawning populations (Table S3) and significantly increased with increasing environmental pressure on batch survival: from 2.28 (CV) under the least pressured environmental scenarios to 3.20 (CV) under the most pressured environmental scenarios in multiplebatch spawning populations, and from 2.37 (CV) under the least pressured environmental scenarios to 3.52 (CV) under the most pressured environmental scenarios in singlebatch spawning populations. Environmental fluctuations had no significant effect on either population size or CV in population size of singlebatch spawners but did affect multiplebatch spawners (Table S4), causing a significant difference between scenarios with and without any applied environmental fluctuations. In particular, population size was more variable under scenarios with no environmental fluctuations as it varied from s.d. 165 individuals to s.d. 134 individuals under scenarios with environmental fluctuations. However, there was no significant variation among scenarios with different environmental fluctuation rates. The same dynamics were observed for population biomass and its CV (Tables S3 and S4). Multiplebatch spawning cod populations had significantly lower realized average mortality rate compared to singlebatch spawning populations (Figure S3a). While the realized mortality rate did not differ among different environmental fluctuating rates in singlebatch spawning populations, it increased significantly in multiplebatch spawning populations when environmental fluctuations were applied to each batch separately (Table S4). 3.2 | Fitness components Under increasing environmental fluctuations, multiplebatch spawners experienced a significant increasing trend in longrun geometric mean fitness WGM , resulting in a higher longrun WGM under most 1986 | HOČEVAR Et Al. unpredictable and uncertain environmental conditions (Figure 3a). Both strategies had an increasing trend in WGM with increasing environmental pressure (Figure 4a), but the relationship was significant, albeit weak, only in multiplebatch spawning populations (Table S6). The difference in WGM across all environmentally stochastic scenarios differed significantly between multiplebatch and singlebatch spawning populations by being lower in the multiplebatch spawning population (Table S5) as their WGM values were lower under conditions when the environmental fluctuations were absent or low, and/ or the environmental pressure was weak. The variance in WAM acrossgenerations was higher overall in multiplebatch spawning populations (Table S5) and had a significant decreasing trend with increasing environmental fluctuations (Figure 3c; Table S6). This significance and decreasing trend were due to higher variance in WAM among generations under nonstochastic conditions, while the variance among singlebatch spawning cod generations increased in the presence of environmental fluctuations. In contrast, environmental pressure had no significant effect on the variance in WAM among generations of either multipleor singlebatch spawners. Multiplebatch spawners had overall significantly higher arithmetic mean in fitness WAM compared to singlebatch spawners (Table S5) due to greater realized reproductive output in the absence of environmental fluctuations (Figure 3b). While the environmental pressure had no significant effect on WAM for either strategy, the presence of environmental fluctuations significantly decreased the WAM of multiplebatch spawners (Table S6). This effect resulted in a lower WAM for multiplebatch spawners when exposed to the highest environmental pressure and environmental fluctuations, hence, experiencing elevated environmental uncertainty and mortality (Figure 3b and Figure 4b). The environmental scenarios where the three fitness components: (i) high longrun geometric mean fitness WGM , (ii) low arithmetic mean in fitness WAM and (iii) low acrossgenerational variance in WAM overlap illustrate that multiplebatch spawning is a bethedging strategy (Figure 4 greyshaded area). 3.3 | Variance in reproductive output within generations Multiplebatch spawning populations had a lower variance in withingenerational reproductive output than singlebatch spawning populations (Figure S3 and Table S3). The two components of environmental stochasticity— the environmental fluctuations and environmental pressure applied to batch mortality— had a significant effect on the variance in withingenerational reproductive output of multiplebatch spawning populations (Table S4). FIGURE 3 Trends of three fitness components: acrossgenerational geometric mean in fitness ( WGM ), acrossgenerational average reproductive output or arithmetic mean in fitness ( WAM ) and acrossgenerational variance in arithmetic mean fitness (variance in WAM ) plotted against increasing rate of environmental fluctuations (xaxis), and grouped by increasing environmental pressure applied to batch survival (colours). Linear regression trends, that were based on the average value of each observed variable across all generations per every run of each spawning type, are illustrated by a solid line for multiplebatch spawners and by a dashed line for singlebatch spawners. Greyshaded area indicates the environmental conditions under which a costly multiplebatch spawning evolutionary outperforms the singlebatch spawning strategy | 1987 HOČEVAR Et Al. However, post hoc Dunn's test revealed that the effects were significant only due to the increased variance under the most pressured environments (0.25) and predictable environmental conditions (0). The scenarios with environmental fluctuations (0.0001, 0.001, 0.01) did not differ significantly in the withingenerational variance of reproductive output. 3.4 | Spawning success The proportion of successful spawning seasons in a population, when an individual produced at least one successfully surviving offspring, reaching age 3 (recruitment age when juveniles become catchable by fishing), per season, differed between populations of multiplebatch and singlebatch spawners (Table S5). The average frequency of such occasions was significantly higher and more consistent in multiplebatch spawning populations when environmental conditions were less predictable and more stressful (Figure 5a,b). In contrast, the scenarios with no environmental fluctuations (0) and low pressure applied to batch survival (0.05) were more favourable to singlebatch spawning populations which exhibited higher spawning success under such conditions (Figure 5a,b). The success probability of mature fish was on average higher and more predictable in the presence of the multiplebatch spawning strategy (Figure 5c). 4 | DISCUSSION It has been hypothesized that multiplebatch spawning in fishes might comprise a bethedging strategy and yield high fitness returns (e.g. Hutchings & Rangeley, 2011). In the present study, we used Atlantic cod as a model species and extended the ecoevolutionary mechanistic model of Kuparinen et al. (2012) to theoretically and empirically explore the hypothesis by evaluating the fitness consequences of such a riskspreading trait. The most interesting finding to emerge from the simulations of our empirically parameterized ecoevolutionary model is that the costly multiplebatch spawning strategy can constitute a bethedging trait under sufficiently uncertain natural environments. The multiplebatch spawning strategy of individuals exposed to fluctuating environmental conditions served to reduce the variance in arithmetic mean fitness across generations, reflecting the decreasing variance in offspring output within generations. The fitness of spawners under stochastic environments is governed by the geometric mean in their reproductive success (Gillespie, 1974) rather than the average mean, which fails to account for environmental variability (Lewontin & Cohen, 1969). We followed the acrossgenerational fitness and found that the multiplebatch spawning strategy maximizes geometric mean fitness by lowering the acrossgenerational variance in arithmetic FIGURE 4 Trends of three fitness components: acrossgenerational geometric mean in fitness ( WGM ), acrossgenerational average reproductive output or arithmetic mean in fitness ( WAM ) and acrossgenerational variance in arithmetic mean fitness (variance in WAM ) plotted against increasing rate of environmental pressure (xaxis), and grouped by increasing environmental fluctuations applied to batch survival (colours). Linear regression trends, that were based on the average value of each observed variable across all generations per every run of each spawning type, are illustrated by a solid line for multiplebatch spawners and by a dashed line for singlebatch spawners. Greyshaded area indicates the environmental conditions under which multiplebatch spawning strategy constitutes a bethedging, as the WGM increases at the cost of reduced WAM