Life-history evolution in response to foraging risk, modelled for Northeast Arctic cod ( Gadus morhua )

Foraging behaviour is known to be a key element in ecology and evolution. Increased foraging intensity increases energy intake, which is useful for growth and reproduction but comes at the cost of higher mortality risk due to increased exposure to predators. Here, we investigate these trade-offs through an individual-based, mechanistic modelling framework adapted to the Northeast Arctic Cod. The model incorporates a series of life-history traits, survival trade-offs, and heritability, which allow evolution to occur and optimal strategies to emerge due to individual trait combinations and their fitness consequences. By altering the relationship between foraging intensity and mortality risk, we find that increased risk causes evolution towards lower foraging effort leading to lower growth and in turn, earlier maturation and a faster pace of life. These results build on previous studies by demonstrating behavioural evolution without direct anthropogenic stressors. Natural mortality among fish is poorly understood, and these results highlight an interesting point of further research that could help future modelling approaches make more accurate assumptions about natural mortality and its components.


The ecology of foraging
All animals require energy to survive and reproduce, as energy is used in a multitude of processes ranging from basal metabolism to tissue repair and synthesis.While some animals can passively acquire energy even when stationary, most animals will need to actively seek out food, an activity known as foraging.Foraging provides an interesting behavioural framework to consider the cost/benefit of performing an activity (Emlen, 1966;MacArthur and Pianka, 1966).
On one hand, the more time spent foraging, the more energy an individual could feasibly obtain -energy that could be used for faster growth, stores, or reproduction, likely leading to more offspring.This makes the benefit of foraging clear, and had there been no downsides, one would expect individuals to forage continuously, barring other limitations such as digestive capacity (Fall and Fiksen, 2020).
However, there are costs associated with foraging, which take on two primary forms: first off, foraging is active and requires energy.The energy gained must exceed the energy cost to be worth foraging.Secondly, foraging individuals are more exposed to predation since they must venture away from potential hiding spots to find food.This leads to an increase in mortality with an increase in foraging effort (Toms et al., 2010).
Individual differences in willingness to engage in risky behaviour have been found in a wide variety of species (Harcourt et al., 2009;Ward et al., 2004;Wilson et al., 1994).This risk-taking also includes a willingness to accept greater risks when foraging (Dammhahn and Almeling, 2012), which can remain consistent within individuals, suggesting that this is an intrinsic trait.Furthermore, studies have highlighted interactions between these behavioural patterns and environmental factors such as temperature (Biro et al., 2010;Killen et al., 2013), making understanding these behaviours increasingly relevant in our warming climate.
Such traits have been linked to personality, often categorised as being either 'shy' or 'bold' (Toms et al., 2010), with risk-averse individuals generally categorised as 'shy' and risk-accepting individuals being categorised as 'bold.'Alternative frameworks have also been suggested, such as risk-taking being characterised by a fear/hunger trade-off (Budaev et al., 2018).Regardless of terminology, personality differences are well documented within fish, and there is good reason to suspect at least partial heritability of personality traits (Dochtermann et al., 2015;Wilson et al., 1994).

Anticipating evolution
Fish evolve in response to stressors, and the evolution of various lifehistory traits has the potential to not only alter demographic parameters such as growth (Crozier and Hutchings, 2014;Enberg et al., 2012;Holt and Jørgensen, 2014) and age-at-maturity (Dieckmann and Heino, 2007;Heino et al., 2002b), but these changes may, in turn, affect the relationship between recruitment and indices such as spawning stock biomass (Enberg et al., 2010), which in turn can lead to mismanagement if not accounted for.It becomes apparent then that proper consideration of evolution should be included in management and conservation efforts (Jørgensen et al., 2007).
Given that individuals will differ in willingness to forage under varying degrees of risk and that these differences are heritable, we have reason to suspect that time spent foraging (and subsequent energy acquisition) would evolve in response to varying foraging risk (Stephens, 2008; Fig. 1), which might vary due to changing community structure or physical environment.Building on this, changes in energy intake might in turn lead to shifts in optimal life history.
Mechanistic models are tools that aim to simulate the natural processes and trade-offs that impact survival and subsequent lifetime fecundity.When properly calibrated, this allows us to calculate optimal strategies for any scenario we wish to test by assuming that strategies maximising lifetime fecundity will be favoured by evolution.
When changing input parameters (such as temperature or foraging risk), new optimal strategies might emerge, indicating the direction we can expect evolution to drive the population.
Among mechanistic models, Individual Based Models (IBMs, sometimes referred to as Agent-Based Models) are particularly suited to studying evolution, as they can simulate evolution by including interindividual variability in inherited traits that influence fitness, allowing individuals to reproduce and die according to their life histories (e.g., Enberg et al., 2009).However, to our knowledge, no evolutionary IBMs for fish have been created that include realistic physiology with oxygen limitation.Such physiological frameworks have been included in other types of models, for instance, dynamic optimization (Holt and Jørgensen, 2014), but including it in an IBM will allow us to more Increased foraging effort increases both energy intake and mortality, but since energy intake experiences diminishing returns various risk scenarios will have different optimal foraging efforts, here assumed to be where the energy gained in relation to mortality experienced is largest.realistically simulate actual population responses to stressors rather than describing optimal behaviour.Additionally, few evolutionary models have included costs associated with foraging, leading to a lack of mechanistic understanding of how risk might impact optimal life history.

Aim of the study
This study aims to build on previous evolutionary modelling efforts using Atlantic cod (Gadus morhua), specifically the Northeast Arctic (NEA) cod stock, as the focal population.The NEA cod stock is not only commercially important, accounting for 796,000 tonnes of catches annually between 2010 and 2020 (ICES, 2021), but is also well-researched, with a large body of literature available on maturation schedule (Heino et al., 2002b), temperature-growth response (Björnsson et al., 2007;Brander, 1995), behaviour (Freitas et al., 2015) as well as expected evolution in response to both climate (Holt and Jørgensen, 2014) and fisheries (Heino et al., 2002b;Jørgensen et al., 2007Jørgensen et al., , 2009) ) to name a few.
Through the lens of NEA cod, we explore the relationship between foraging risk and the evolution of foraging effort and maturation schedule, as well as the emerging life histories.In doing so, we develop an IBM that includes explicit energy budgeting sensitive to temperature, suitable for studying evolutionary responses to fisheries and climate.As such, the model does include complexity that is not strictly necessary to study the effects of foraging risk but is useful for increased versatility in future model applications.

Model description
We introduce an individual-based model based on the framework previously published by Enberg et al. (2009), adding explicit energetics and oxygen budget for NEA cod within a stochastic food environment, drawn mainly from Holt and Jørgensen (2014).The inclusion of explicit energetics and oxygen budget allows for better accounting of what effect temperature can be expected to have on the population for use in future applications of the model.Each individual carries traits that impact their life history through maturation, growth, reproductive investment, and mortality, while population dynamics and fitness are allowed to emerge from the interactions of all these individuals (Fig. 2).Each of these sections are described in more detail below, and a full list of parameters can be found in Table 1.

Purpose
The model aims to simulate the evolution of the maturation scheme and somatic energy allocation (life history) by including genetic traits that influence the phenotypic expression, with a measure of stochasticity to mimic environmental variability, which impacts survival.The model also introduces the inherited trait "Appetite," which takes the form of desired energy intake, in turn governing foraging activity.As such, individuals with high appetites are comparable to bold individuals on the bold/shy spectrum (Toms et al., 2010) or hungry individuals in the fear/hunger framework (Budaev et al., 2018), while the inverse is true for individuals with low appetites.
Surviving individuals reproduce, passing on their genetic traits to the next generation.This method allows for evolutionary patterns to emerge over time, as individuals with favourable traits will survive and produce more offspring over their lifetimes, simulating natural selection.
Using this model, we will test evolutionary responses to varying degrees of foraging risk, represented by the variable 'foraging risk exponent,' e 2 , as explained below.

Process overview
The model simulates 2000 years in annual time steps.Every year, individuals are first faced with the possibility of maturing, after which energetics are calculated.Based on trade-offs in these sections, mortality/survival probability is determined.Recruitment (and inheritance) is then performed based on gonadal energy allocation, followed by model upkeep, removing dead individuals and ageing the remaining ones.Individuals are allowed to live for 20 years, after which they are removed, and recruits are added as 1-year-olds.

Initializing the model
To initialize the model, we first ran it for 5000 years to let the population stabilize using a foraging risk exponent e 2 = 1.8 (explained in more detail below).5000 years is used to make sure that the initial population is evolutionarily stable and does not carry trends into the simulation experiments.The risk exponent value (e 2 = 1.8) was chosen because it best replicated observed growth patterns and proportion mature-at-age when comparing the model against observed data (ICES, 2021; Fig. 3) from the Barents Sea and Lofoten between 2004 and 2018.While the fit is not perfect, it should be noted that the NEA cod has been fished for many generations, likely explaining this small discrepancy, as the proportion of mature individuals better fit with the data from 1946, before fisheries were intensified (Fig. 3, B).The resulting population (305,470 individuals) was used as a starting point for all future runs.For every value of e 2 tested, the model was run for 2000 years, and each run was done 20 times with different random seeds to act as replicates.Two thousand years was chosen because that was the earliest point in which inherited traits were stabilised, and evolutionary trends were visible.

Input
The model includes several sources of stochasticity: phenotypic expression of genetic traits, variability in the food environment, temperature and recruitment, and offspring variability compared to midparental values.For simplicity, all of these are assumed to be normally distributed around their respective values with a standard deviation of 5%, except for offspring variability (described below), as sensitivity analysis showed that none of them were driving results.Temperature is centred around 4 • C, which is the average annual temperature in the Barents Sea as measured along the Kola hydrographical transect (0 -200 m, 1900(0 -200 m, -2020(0 -200 m, , Boitzov et al. 2012)).

Maturation
The maturation scheme is based on a Probabilistic Maturation Reaction Norm (PMRN) Heino et al., 2002b), shown in Eqs.(1)-((3).The probability of maturation for individual i, p mat (i) within year y is dependent on the intercept and slope of the maturation reaction norm, where L y is the length of individual i in year y, L p50 is the length at 50% maturity (Eq.( 3)), and θ is based on reaction norm width as follows, where w is the PMRN width, p u and p l being the upper and lower probability bonds, respectively (generally set to 0.75 and 0.25).The length at 50% maturity, L p50 , is calculated as, where L I and L S are the intercept and slope of the reaction norm, respectively, and Age y− 1 is the age prior to year y.While Enberg et al. ( 2009) allowed both intercept (L I ) and slope (L S ) to evolve, they also found that the evolution of slope added little explanatory value.Hence, in this study, only the L 1 is an inherited trait.Slope (L S ) is set to 2 cm year − 1 , and width (w)¸is set to 40 cm, as this resulted in a reasonable maturation schedule when compared to available data (ICES, 2021).When interpreting results, higher L I means late maturation at a larger size, while low L I means earlier maturation at a smaller size.The slope parameter L S influences how dependant the maturation process is on length and age relatively: if slope is 0, it is only the size that impacts maturation probability, whereas increasingly negative slope values mean that age start playing increasingly important role in the maturation probability.For more information on PMRNs, please refer to Heino et al. (2002b).

Energetics
All individuals (i) forage until they reach their predetermined energy intake, ϕ, determined as follows, where ζ is the appetite described earlier for individual i, W is weight, and β is an allometric exponent.This model then includes a modified energetics framework based on the Wisconsin Bioenergetics Framework Hewett and Johnson, 1992), previously parametrized for NEA cod by Holt and Jørgensen (2014), Eqs. ( 5)-((12) (Eq.( 6) modified to account for density dependence).
For every year, a food environment ω y is set, where χ 1 is a stochastic function normally distributed around 1 with a standard deviation of 0.05, and Σ(W) is the total population biomass to simulate competition for food (density dependence).The total energy intake for individual i in any year, ϕ y (i), is then determined by, where f int is the foraging intensity, which individuals scale to reach the desired energy intake set by ζ, and c F1 and c F2 are constants.For ease of interpretation, f int is given in units of standard metabolic rate, B SMR .
Based on the food intake we can construct an energy budget for individual i, where N y is the net available energy for growth in year y, B SDA is the energy lost by specific dynamic action (digestion, etc.), B SMR is the standard metabolic rate based on weight and temperature, B ϕ is the energetic cost associated with foraging, and c R is the efficiency of converting energy to new tissue.B SDA , B SMR and B ϕ are calculated as follows, where c SDA , c SMR and c ϕ are constants (Table 1) and f(T) is an Arrhenius function of temperature, c T1 , c T2 and c T3 are constants (Table 1).
The available energy, N y , is then divided based on the maturity status of the fish; for immature fish, all the energy is allocated to somatic growth G s , where ρ s is the energy density of somatic tissue.
For mature individuals, a somatic growth allocation function, p t , is used to describe the decreasing allocation to somatic growth as the fish ages (Quince et al., 2008), resulting in increased gonadal allocation, where c D1 is the initial investment to somatic growth, c D2 is the rate of somatic growth allocation decay, A is the current age and A m is the age at maturation.Initial investment, c D1 , is an evolving trait in the model.Energy is then divided between somatic growth, G s , and gonad growth, G g , as in Holt and Jørgensen (2014), where B M is the energy used for spawning migration and ρ g is the energy density of gonadal tissue.Mature individuals of NEA cod undertake annual spawning migrations of 780 km each way to spawn near the Lofoten islands (Opdal et al., 2011).B M is modelled based on the swimming speed bioenergetics published by Ware et al. (1978), where optimal swimming speed is a function of length (body lengths per second), U opt (i) = c u * L b2 (i), where c u and b 2 are constants.The cost of transport is then calculated as a function of length and optimal swimming speed, c M (i) = c COT * L b3 (i) * U opt b 4 (i), where c COT , b 3 and b 4 are constants (Table 1).The total energetic cost of spawning migration then becomes, where D M is the migratory distance, undertaken twice for the round-trip.
A minimum amount of energy allocated to gonads, N y * (1− p t ), is required for the fish to undertake the spawning migration, here calculated as having enough energy to reach a gonadosomatic index of at least 0.1 post-migration.This value was chosen in order to achieve a realistic number of spawning-skippers, based on personal experience.If an individual does not meet this energy threshold, it instead allocates all energy to somatic growth and skips spawning for the year.Growth in length is derived from the growth in weight by allometric scaling, where k is a constant (Table 1).
The maximum oxygen uptake V max of individual i is calculated as, where V 1− 4 are all constants, and W som is the somatic weight (excluding gonad weight).Oxygen consumption V y is defined as the sum of the metabolic processes, where B growth is metabolic work associated with converting energy, c R ,

Mortality
Since many of the activities accounted for in the energetics section above impact mortality, the trade-offs must be included in the model.As such, mortality in the model is split into several compartments, based on previous modelling efforts by Holt and Jørgensen (2014), Eqs. ( 21)-( 26), where M fixed is a fixed rate of size-independent natural mortality, and the remaining factors are the mortalities associated with predation, foraging, reproduction and respiration, respectively.For any given individual i, the probability of surviving, S(i), in a given year is, Predation mortality is the primary size-dependant mortality, in the form of the increased risk of being preyed upon as smaller individuals.
where c predation and e 1 are constants (Table 1), and L is the length of the fish in cm.Foraging mortality comes primarily in the form of increased risk of being preyed upon when increasing foraging intensity (f int ), and is therefore related directly to M predation , where c ϕ and e 2 are constants.The risk associated with foraging is governed by e 2 , the foraging risk exponent, which is the focus of this paper.By changing e2 we alter the relationship between foraging effort and mortality as seen in Fig. 4. Reproductive mortality relates to decreased swimming capacity due to the change in form factor and additional risk-taking behaviour during courtship and mating.
where GSIref is the GSI at which Mreproduction = Mpredation, and e3 is a constant.Respiration mortality is the result of limited respiratory capacity in relation to respiratory demand (essentially exhaustion), where c respiration and e 4 are constants (Table 1).

Recruitment
The number of recruits in year t is calculated using Beverton-Holt recruitment Beverton and Holt, 1993) based on total fecundity, ΣQ(i, t), adapted for Atlantic cod by Enberg et al. (2009), Eqs. ( 27) & ((28), as follows, where c R1 and c R2 are constants (Table 1) determining survival at low fecundity and strength of density dependence, respectively, and e λR (t)  describes inter-annual environmental variability.In practice larger individuals are expected to dedicate more energy to gonads, making individual fecundity weight-dependant.Inheritance of traits is simplified in the model, as it is in Enberg et al., 2009, due to the often highly polygenic nature of life-history traits (Roff, 1993;Conner and Hartl 2004).For every recruit, two parents are  24)), in order to make the figure general.chosen from mature individuals by random sampling, weighted by gonad growth, G g (i).This makes a fish more likely to be chosen as a parent the more it contributed to the total fecundity of the population.Evolving traits are then calculated based on mid-parental values as shown here for appetite, ζ where χ 2 is a stochastic function normally distributed around 1 with a standard deviation of 0.14.This stochastic function aims to simulate the effects of mutation, segregation and recombination, generating heterogeneity in offspring for selection to act upon.The standard deviation of 0.14 was chosen as it yielded emergent heritability of approximately 0.2 calculated as the linear correlation between mid-parental values and offspring values for both length-and age at maturation, which is within the range typically seen for life history traits (Gjedrem 1983;Law 2000;Carlson and Seamons 2008).Note that the model does not separate between males and females, but only has a single sex.This is considered to be an acceptable simplification when male and female life-histories and demography are similar (Dunlop et al., 2009) -which we believe is the case for NEA cod.

Results
We found that an increase in foraging risk leads to evolution towards lower appetite and, thereby lower growth (Fig. 5, A).Comparing the runs over time (2000 years), appetite rapidly (within 250 years) diverged and stabilised around new means, indicating strong selection.The highest mean appetites were found for the lowest risks (~ 16,100 kJ kg − 1 ) while increasing risk decreased mean appetite up to 12% (~ 14,200 kJ kg − 1 ) in the highest risk scenario.A consistent trend for all results was that values of foraging risk (e 2 ) between 1.0 and 1.8 were more similar than the highest values (e 2 = 2.2 and 2.6).This is likely because foraging risk increases exponentially with increasing e 2 .Increasing the foraging risk even higher leads to the population crashing due to too high mortality in the early years (results not shown).
In extension of the reduced appetite, growth was consistently lower in the higher-risk scenarios, as seen by the lower size-at-age (Fig. 5, B).Individuals in the lowest risk scenario grow to reach 100 cm in length by age 10, while individuals in the highest risk scenario reached 9 cm.There were no clear differences between the three lowest risk scenarios, implying that growth is not significantly impacted until risk reaches a certain threshold.
Mortality increased notably in the two high-risk scenarios, a trend that is particularly apparent for the younger fish (Fig. 5, C), likely due to mortality being scaled by length-dependant predation mortality (Eq.( 23)).For the one-year-old fish, the rate of mortality was 0.28 y − 1 in the lowest risk scenario, but 0.61 y − 1 in the highest risk -more than a twofold increase.While this difference diminishes as the fish ages it doesn't disappear until the fish reaches 14 years of age, at which point the fish have already matured and reproduced (Fig. 6).This also translates to a smaller population size, measured both as total biomass and as the number of individuals (Appendix C).
The changes in energy acquisition and mortality lead to changes in life history, specifically by favouring a reduction in PMRN intercept after stabilization (Fig. 6, A).In this case, the clustering of low-risk scenarios compared to high-risk scenarios is even more pronounced than for appetite.The average intercept for low-risk scenarios ranges from 113 to 107 cm, while the highest-risk scenario leads to an intercept of 95 cm.
Unlike the other evolving traits (appetite and PMRN intercept), the initial gonadal allocation c D1 did not seem to change significantly with increasing risk (e 2 ) (Fig. 6, B).While the higher e 2 values did at first lead to an increase in the initial allocation, allocation started decreasing again once PMRN intercept and appetite had stabilised (around year 4-500) (Appendix B), and even when the difference between risks is largest, average allocation only differs by 1 percentage point.After stabilization, differences appear negligible.The initial increase in gonadal allocation in high-risk scenarios is likely driven by the sudden increase in mortality, making investing in reproduction over growth more beneficial.Once appetite stabilises, this pressure is relaxed, and allocation once again decreases.
To see how these changes in PMRN intercept translated into changes in the actual life histories, we considered age and length at maturation, which follows the same trend (Fig. 6, C&D).As risk increases, mean length and age at maturation decrease, with this effect being more pronounced for higher risks.However, while the overall trend is similar, the variability of length and age at maturation is notably higher than for the PMRN intercept, implying that other processes also influence the maturation schedule.This can happen because the PMRN intercept is directly inherited, while the emergent traits (such as length/age at maturation) result from more processes and traits also involving stochasticity.This is the case for e.g.age-at-maturation, which is not directly controlled by the PMRN intercept but is also influenced by appetite that controls energy acquisition, which drives growth and in turn, maturation.

Discussion
We found that with increasing foraging risk, evolution selected for lower appetite and earlier maturation, both leading to lower growth.This is expected, as higher mortality associated with energy acquisition should select for individuals with lower appetite, which forage less.Further, the fitness benefit associated with late maturation is offset by lower overall survival due to higher foraging risk, causing a selection towards earlier maturation.
While it seems intuitive that an increase in foraging risk and a decrease in energy intake would lead to reduced growth and higher mortality, respectively, these points still warrant further consideration.Since we are not enforcing a lower limit for appetite, the fish can  theoretically maintain lower mortalities even in high-risk scenarios by foraging less and accepting slower growth.Similarly, they could have accepted higher mortality risks to maintain growth, but it seems that the optimal strategy is a combination of the two.
These results highlight a critical life-history trade-off related to optimal foraging behaviour.Given that overall mortality is higher for smaller individuals, a selection towards fast growth early in life (McGurk, 1986;Shepherd and Cushing, 1980;Ware, 1975) motivates foraging.However, foraging is a risky activity, so high foraging effort will lead to higher mortality risk -this dissuades foraging (Stephens, 2008).This is further complicated by larger individuals having higher relative fecundity (Hixon et al., 2014;Morita et al., 1999), further motivating foraging.However, fisheries may nullify this factor by targeting larger individuals (Jørgensen et al., 2009).
A similar trade-off has been considered by Claireaux et al. (2018), who used a dynamic optimization to simulate evolutionary endpoints when fisheries target certain behaviours on the bold/shy spectrum.They found that any increase in fishing mortality, even unselective, caused evolution towards reduced growth and age/length-at-maturity.However, when explicitly targeting bold/hungry individuals, the selection became stronger, and also selected for reduced foraging effort.
This is an interesting comparison to our study.While the sources of mortality (predators vs fisheries) and the modelling framework (IBM vs dynamic optimization) differ, the experiments are quite analogousincreasing mortality disproportionately for foraging individuals.In both studies, this leads to lower foraging rates and a faster pace of life.As such, our study expands on Claireaux et al. (2018) by using an IBM to simulate actual life histories and by showing that this sensitivity to foraging risk is not limited to fisheries but is fundamental to optimal foraging theory.
It is a commonly found trend that higher mortality leads to earlier maturation (e.g., Heino et al., 2002a;Law, 2000Law, , 2007;;Jørgensen et al., 2009;Claireaux et al., 2018;Andersen et al., 2018).When life expectancy is reduced, the associated benefits of large body size and late maturation diminish (Heino et al., 2015).This aligns with our findings that high-risk scenarios and increased mortality select for earlier maturation life histories.
Adaptation to an increase in external mortality, like fishing mortality, has been shown to lead to increased natural mortality (Jørgensen and Fiksen 2010;Jørgensen and Holt 2013), and even though in our study, the increase in natural mortality is due to adaptation to higher foraging risk, which could be caused by for example increased predator population, the bottom line is that while a population is adapting to increased external mortality, it can lead to an increase in total natural mortality rates.
Our results (along with those of Claireaux et al., 2018) highlight the need to consider the coupling between energy acquisition and mortality.It becomes clear that when mortality is coupled with foraging, attempts to decrease mortality risk associated with being small comes at the cost of increased mortality risk from foraging.Other studies have found that increased selection pressure for small or intermediate sizes (e.g., gillnet fishing) can favour later maturation since late-maturing fish will more quickly grow out of the vulnerable size range (Huse, 2000;Jørgensen et al., 2009;Stepputtis et al., 2016).This trend, however, is only found at moderate fishing pressure, allowing fish to grow past the vulnerable size range.
In the past couple of decades, fisheries selection for behavioural patterns has been getting more attention, as different gear types have been shown to select for bold/shy behaviour (Cooke et al., 2007;Uusi-Heikkilä et al., 2008;Biro and Stamps, 2008;Diaz Pauli and Sih, 2017;Arlinghaus et al., 2017;Claireaux et al., 2018).As such, changes in behavioural patterns and the effects they have on life history evolution are becoming known.Our results further build on this by showing that the evolution of behaviour (e.g., foraging activity) is not limited to direct anthropogenic stressors such as fisheries but can also result from other changes to the ecosystem.For instance, foraging risk might change due to new species altering community structure, or changes in the physical environment that either reduce opportunities for seeking shelter or make foraging more time-consuming.Knowledge of what behavioural phenotype is favoured by natural selection would allow managers to make better-informed decisions about which gears to use.
The current model is an amalgamation of two earlier models: an ecoevolutionary individual-based model with rich population dynamics including density dependence, several evolving life history traits, but no explicitly physiological relationship climate (Enberg et al., 2009), and an optimization model with detailed physiological mechanisms including oxygen budget (Holt and Jørgensen 2014).Even though not fully taken advantage of in the current study, this model is a powerful tool for modelling the concurrent contemporary evolution and eco-evolutionary dynamics of the NEA cod in relation to fishing and climate warming.Future studies should include investigations of the relative importance of these anthropogenic drivers and which harvest strategies might be most suitable and least harmful for the long-term sustainability of the fished population.
In summary, this study shows that increasing foraging risk leads to a decrease in appetite, foraging activity and energy acquisition.These changes in turn, lead to an increase in mortality and subsequent selection for faster life histories.This builds on previous studies showing that foraging behaviour is subject to evolution in response to stressors and expands the field by demonstrating that even without direct anthropogenic influence (such as fisheries), ecological factors may select for bold or shy foraging strategies.For fisheries and their managers, this avenue of potential evolution should be addressed, as knowing which behavioural phenotype is most favoured locally allow for more informed decision-making.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Conceptual illustration of the survival trade-off when foraging.Increased foraging effort increases both energy intake and mortality, but since energy intake experiences diminishing returns various risk scenarios will have different optimal foraging efforts, here assumed to be where the energy gained in relation to mortality experienced is largest.

Fig. 2 .
Fig. 2. Conceptual schematic highlighting critical processes in the model relating external factors and heritable traits to fitness outcomes through energetics.Solid arrows indicate direct influence, and dashed arrows indicate the energy flow.The shaded '%' area represents the proportion of energy dedicated to somatic/ gonadal growth.

Fig. 3 .
Fig. 3. Initial population compared to ICES data (squares).(A) Average weight-at-age ±SD of the final population compared to ICES averages from 2014 to 2021.(B) Average proportion mature-at-age ±SD of the last 200 years of stabilisation (Appendix B), compared to both ICES averages from 2014 to 2021 (black squares) and ICES values from 1946 (hollow squares).1946 values were not available for weight-at-age.

Fig. 4 .
Fig. 4. Changes in base foraging mortality in response to increasing foraging effort, shown for several values of foraging risk.Note that the mortality values on the yaxis does not yet include the scaling by predation mortality, M predation (see Eq. (24)), in order to make the figure general.

Fig. 5 .
Fig. 5. Changes in energy acquisition and resulting growth/mortality differences.(A) Average genetic appetite ±SD over the model's runtime.(B) Average length ±SD as a function of age at the end of the 2000-year runtime.(C) Total mortality ±SE (n = 20) at the end of the 2000-year runtime.For plots B and C, values for foraging risks of 1.4 and 2.2 are omitted, as they didn't fall outside visible trends, but the dense clustering made the plots difficult to read.

Fig. 6 .
Fig. 6.Changes in maturation schedule in response to varying foraging risk.Values shown are all stabilised averages at the end of the 20 model runs.(A) Average PMRN intercept ±SD.(B) Average genetic initial gonadic allocation ±SD.(C) Average length at maturation ±SE.(D) Average age at maturation ±SE.Time trajectories of changes in these traits are shown in Appendix B.

Table 1
Model parameters.