Intra-specific diversity and adaptation modify regime shifts dynamics under environmental change

1.   Introduction

Ecosystem management is a major concern due to climatic and anthropogenic pressures on the environment. To study such pressures, ecologists developed tools based on early warning signals (EWs, ^[1]) to anticipate changes in ecosystems, such as regime shifts. They depict modifications from one state (e.g., a forest), to another (e.g., barren lands). For instance, such tool was applied to fishery management in the Caribbean ^[2]. Fisheries target the parrot-fish Sparisoma viride, a major grazer of algae. Overfishing of the parrot-fish is suspected to induce an increase in algae, and a loss of corals ^[2]. Such shift is irreversible, so that recovering the initial coral cover requires a greater reduction in fishing pressure than the increase in fishing needed to shift. Thus, regime shifts depict catastrophic ecosystem shifts, and restoring the system requires strong recovery efforts ^[3].

However, species can adapt to changes in the community. For instance, in planktonic systems, predation pressure by zooplankton selects for resistant phenotypes in phytoplankton ^[4]. Furthermore, selection can alter population dynamics. In a chemostat experiment, ^[4] showed that selection of predation resistance in phytoplankton leads to cycles of predator and prey densities. Zooplankton preyed on vulnerable phytoplankton and selected resistant prey. Then, zooplankton density decreased due to a lack of vulnerable prey. As zooplankton were low in density, vulnerable phytoplankton grew again, and cycles were driven by the evolution of predation resistance in phytoplankton. Another example is cuckoos and their interactions with bird hosts. Cuckoos must select vulnerable nests to lay eggs. However, hosts can learn to recognize cuckoo eggs, and reject them. Unmasked cuckoos then need to find other nests, resulting in co-evolution of cuckoos and their hosts ^[5]. Interactions between ecology and evolution brings further complexity to community dynamics. Selection alters life-cycle of species ^[6], and alters community structure. Thus, evolutionary changes in community structure can alter regime shifts. However, research on interactions between evolution and regime shifts is at its infancy. The researchers in ^[7] argued that regime shifts depend on evolution, and can lead to both facilitated or delayed shifts. In extreme cases, evolution can cause regime shifts ^[8].

Evolution relies on phenotype selection and heritable traits diversity. However, the size of diversity is also affected by extreme events. For instance, storms shape diversity by eliminating rare phenotypes and selecting more resistant ones. The environment thus eliminates large deviations from the optimal phenotype, leading to a decrease in genetic diversity ^[9]. However, high genetic diversity favors the survival of populations under extreme events. Reference ^[10] showed that such diversity increases the survival of populations of Zoster marina under heat waves. Thus, trait diversity in populations was argued to increase the capacity to absorb environmental change and maintain the population's state. This idea is related to population survival and resilience ^[10,11]. Resilience should increase with the probability of presence of resistant phenotypes, and the size of diversity. On the other hand, more complex interactions between resilience and diversity may exist. Diversity represents an investment, as some diverging phenotypes are more vulnerable to current environmental conditions ^[12]. Furthermore, environmental changes can select non-dominant phenotypes, and cause regime shifts ^[13].

Diversity depends on environmental change ^[14], and may play a central role in regime shift dynamics. Fast environmental change erodes diversity and adaptation, and induces higher mortality rates. Such increase may in turn lead to populations extinction. Thus, more knowledge between heritable trait diversity, evolution, and regime shifts is crucial to better anticipate shifts and establish management strategies.

The link between heritable trait diversity in populations and regime shifts is thus unclear, and is still a subject of debate in ecology. Here, we used different approaches to explore the influence of evolution and diversity on regime shifts dynamics. First, we used the adaptive dynamics ^[15] to explore the role of evolution. This approach assumes that evolutionary and ecological time scales are distinct. However, here we relaxed this assumption and allowed evolution to occur on ecological time scales. Second, we also used a partial differential equations model to simulate trait diversity dynamics in populations during regime shifts. In this second approach, evolution occurs through selection in the available trait diversity. Thus, we used simulations of heritable trait diversity to validate assumptions made on eco-evolutionary time scales. The different approaches, and their models, are described in the "Materials and methods" section.

Here, we used three numerical models to explore interactions between regime shifts and evolution. First, we investigated how evolution affects regime shifts thresholds, especially in a scenario of ecosystem recovery. Second, using a stochastic model we asked if environmental variability modifies the probability of regime shifts. Then, we discussed interactions between stochasticity, population resilience, and evolution. Third, using simulations of heritable trait diversity and a partial differential equations model, we asked if diversity can delay regime shifts through the selection of the "fittest" trait in the available genetic stock. Finally, findings are summarized in the "Conclusions" section.

2.   Materials and method

2.1.   Regime shifts

Two types of environmental change are defined here. First, populations can suffer a disturbance, a sudden and temporary change. Disturbances are spontaneous events in time, such as a sudden increase in mortality due to storms. Second, if the change persists in time, it is called a perturbation. A linearly increasing fishing pressure over years is such an example. Large perturbations can thus lead to regime shifts (see Figure 1(b)). Under low amplitudes of changes, populations can persist, and maintain the community's current structure ^[3,16]. This absorption capacity of populations is described using the concepts of stability, and resilience ^[16]. First, stability defines the propensity to recover a population's equilibrium after environmental changes. The higher the stability, the faster the recovery. Second, resilience defines the amplitude of change populations can absorb before changing state, or shifting. The more resilient a population, the more change it can absorb. Resilient systems can suffer large disturbances and maintain their structure. In Figure 1, high resilience translates to very large basins of attraction. For instance, at low environmental parameter values, shifting from the upper basin of attraction to the lower basin requires a substantial decrease in population density (see Figure 1(b)). If the limit between the basins is crossed, populations are not resilient enough to recover the high density equilibrium. Thus, disturbances that exceed resilience thresholds lead to a shift to the other state ^[17]. Shifts in population density were modeled using the hysteresis, as in Figure 1. The environmental change thresholds before regime shifts are also displayed as tipping points. Thus, populations can shift due to a disturbance in density (a change of basin of attraction) or due to a perturbation (the system crosses a tipping point, see Figure 1(b)).

Figure 1.  (a) Depiction of the influence of evolution and stochasticity on regime shift dynamics, based on ^[7] and ^[17]. The hysteresis depicts two stable states for the same range of an environmental parameter. Here, the states are high population density, and low density. The system can shift from one state to the other by crossing a tipping point. These tipping points are thresholds of parameter values, beyond which the population shifts. Evolution (colored lines) modifies these tipping points, and thus thresholds for shifts. Shifts can also occur if the population changes its basin of attraction. The upper basin of attraction is pictured on the upper part. While inside this basin (density is higher than the dotted line), the population will remain attracted to the high density state. On the contrary, if the density decreases and enters the lower basin, the population becomes attracted to the low density state. In (b), examples for regime shifts due to perturbation and disturbance are shown. Following a perturbation, the system can cross a tipping point and shift. In the case of a disturbance, as a sudden increase in mortality, the system can shift due to entering the lower basin of attraction.

DownLoad: Full-Size Img PowerPoint

These tipping points can be delayed, or be crossed earlier due to evolution ^[7] (see Figure 1). Thus, regime shifts depend on environmental change and ongoing evolution. Furthermore, evolution can also modify basins of attraction, and thus resilience ^[18]. For instance, adaptation of prey to an increase in predation pressure can then delay shifts to low prey densities. Here, we explore relations between evolution and regime shifts using an adaptive dynamics model of a prey population (adaptive dynamics theory is presented in the next section).

2.2.   Eco-evolutionary model

2.2.1.   Population dynamics model

The model used in this work is inspired from ^[19] and we use a very similar form. It includes prey intrinsic growth rate and resource competition, such as in ^[20]. Predation is represented through a Holling Ⅲ functional response.

Prey population at any time $t$, abundance $N(t)$, in (2.1), depends on the intrinsic growth rate $r$, on intra-specific competition and on predation.

The parameter $B$ is similar to the parameter $\beta$ in Eq (3) of ^[19], and the parameter $a$ corresponds to $\dfrac{1}{\alpha^2}$ in ^[19]. $a$ is actually proportional to the attack rate, and $B$ represents the maximal predation intensity per unit of time. In our study, predation pressure $B$ is regarded, from the prey population point of view we adopted, as an external environmental factor, that can vary continuously through different levels. "Environmental parameter" will always refer to this predation pressure $B$.

2.2.2.   Prey adaptive trait $s$

We considered prey evolution through adaptation of one of its phenotypic trait, that we denote $s$, that takes continuous real values. In the following, we will develop formulation choices based on the hypothesis that this trait corresponds to the mean body size of prey individuals. Yet, it may be replaced by any specific trait that would fit to our dependencies hypotheses.

This phenotypic trait naturally impacts life history traits in prey demography model (2.1) in a way that translates evolutionary trade-off or side effects. In this study, we use the adaptive dynamics framework ^[15] to develop and analyse model of prey evolution. The first step relies on the explicit formulation of the dependencies of population dynamics model parameters on the adaptive trait we focus on. Most parameters are thus replaced form then by functions, without changing their names for commodity:

Intrinsic growth rate:

We hypothesize that the intrinsic growth rate decreases with body size. This reflects the cost of growing big or can be seen as a very common allometric relationship between population growth rate and body size. Furthermore, we make the assumption that competition gives advantage to individuals of greater body size. This is natural to assume that it comes with a cost that negatively impacts growth rate. All this considered, We chose a simple, general formulation:

This function decreases from a maximal value $r_{max}$, with a speed determined by $D_r$ and tends to $0^+$ from very large body sizes.

Competition

In the next section, we explain how adaptive dynamics introduces diversity in the prey population, in order to model natural selection. Equation (2.2) will be replaced by a system of ODE to asses the dynamics of subpopulations of prey that present a diversity in their trait $s$. We hypothesize that in such a setup, competition $\alpha$ will depend on the trait difference $s_1-s_2$ between an individual of type $s_2$ that exerts competition and an individual of type $s_1$ that undergoes the competitive effect and it is asymmetric in favor of the greater. Within individual of a same subpopulation (or equivalently when the population is monomorphic in $s$), individuals are of identic type and the competition intensity simplifies to $\alpha(0)$ following this logic. Otherwise, following ^[20], it is given by the function:

where $\delta_s = s_1-s_2$ stands for the trait difference between interacting individuals. This is a decreasing sigmoid function that varies from $c_{\alpha}$ to $0$ for very large negative of positive trait difference. Parameter $k$ tunes the steepness of the curve (intensity of the competitive asymmetry) while $V$ slides the inflection point and therefore the convexity in $\delta_s = 0$.

Predation rate

We hypothesize that individuals of high body size benefit from some advantage with regards to reduced predation, whereas smaller individuals are consumed at higher rate (see ^[21,22]). More precisely, $a(s)$ in taken as:

$1/K_c^2$ corresponds to the minimum predation rate that undergo large individuals and $D_c$ stands as a "four times minimum predation rate" parameter.

2.2.3.   Invasion fitness function and canonical equation of adaptive dynamics

Adaptive dynamics theory ^[15] models evolution dynamics through a sequence of substitution steps of an in-place resident type by a successful mutant type, that then becomes the new resident type. The success, or failing, of a random mutant lineage is assessed in a setup for natural selection in which resident and mutant types compete, through the mutant-resident system:

that explicitly quantifies intra- and inter-type interactions, with $s$ and $s_{mut}$ respectively the resident and mutant adaptive trait value, and $N$ and $N_{mut}$ the resident and mutant population size.

From there, we infer the invasion fitness function that quantifies the initial growth rate of the mutant lineage when it is initially rare and appears in the resident population at its demographic equilibrium $\hat N(s)$:

A mutant type with a negative invasion fitness will certainly fails whereas a positive invasion fitness provides it a positive invasion probability. The sequence of randomly appearing mutant lineages and their selection against the former resident type (when their fitness is positive) gives the evolutionary dynamics of the trait $s$.

Fitness gradient gives the direction of evolution of the adaptive trait:

Lastly, canonical equation of adaptive dynamics provides a smooth approximation of the dynamics of the expectation of the adaptive trait evolution on the long-term evolutionary timescale:

2.2.4.   Population size equilibria for fixed $s$

In the reference "no evolution" scenari, $s$ has a fixed value. Beside its trivial equilibrium $N = 0$, Eq (2.2) for monomorphic population exhibits one or three real positive equilibria, depending on parameters and trait $s$ values. The value of this equilibria can be computed analytically using Caradano's formula (see Supplementary 1). $N = 0$ is always unstable. In cases where there is only one other equilibrium, $\hat N(s)$, it is always stable. In cases where there are 3 other equilibria, the smaller and the higher ones are stable whereas the middle one is unstable and delimits the two basins of attraction. In Eqs (2.7)–(2.9), $\hat N(s)$ stands either for one or for the other of these stable equilibria, depending on the "demographic state" of the population. This possibility of bistability allows for a hystereris scheme (see Figure 1).

2.2.5.   Coupled eco-evolutionary dynamics

Adaptive dynamics theory setup relies on the assumption that demography is fast with regards to evolution. In this study, relative timescales between evolutionary adaptation, environment change and demography are a key factor that lead to the variety of system behaviors and scenario we obtained. In order to allow ourselves to intermingle these timescales, we stepped aside from the strict application of adaptive dynamics hypothesis and loosen the strict timescale separation between demography and evolution, by setting up a coupled system of ODE combining population dynamics (2.2) and canonical equation for evolution (2.9):

Notice that invasion fitness $f(s_{mut}, s, N)$ (and thereby fitness gradient) is now a function of $N(t)$ too, since in this setting $N$ remains a state variable of the system and is no longer fixed at its (or one of its) equilibrium value $\hat N(s)$.

From (2.10), several numerical approaches were derived. First, we performed direct integrations of ODE. Second, we added random effects on the prey populations. and third, we simulated heritable trait diversity during regime shifts (see the protocol in Supplementary 2). All simulations were carried out under Python version 3.10.9.

2.2.6.   Prey individual based model (IBM)

Disturbances, such as storms, impact prey by inducing sudden peaks in mortality. Hence, small populations are expected to be most sensitive to environmental variability. An IBM was thus constructed using a Gillespie algorithm, to explore the interactions between stochasticity and regime shifts.

Birth and death processes occur at each time step $\Delta t$ following a multinomial law. Each prey follows one of 8 outcomes : die, survive without offspring, or survive and have up to 6 offspring. A number of 6 offspring was an extremely rare case, and was considered sufficient for simulations. Here, we write $\lambda$ and $\mu$, the rates of reproduction (or survival), and of mortality. They were deduced from the growth (intrinsic growth rate) and mortality (competition and predation) terms in $\frac{dN}{dt}$ (2.10).

System (2.11) describes the probabilities used in the Gillespie algorithm. The probability for survival is $a$, and $b$ is for mortality.

Using probabilities in (2.11), we computed the probability of having a number $k$ of offspring for a trait value $s$, given by $P_s^k$. Offspring were computed as $P_s^k = (1-a).(1-b).b^k$. Thus, $P_s^k$ gives the number of offspring for each different trait value, and at each time step.

The number of individuals simulated depends on a scaling factor, as $N \cdot \text{Scaling factor}$. The scaling factor modifies the influence of stochasticity on prey densities. The greater the scaling factor, the lower the influence of stochasticity. To maintain homogeneity of units between density and parameters, mortality and reproduction probabilities were computed using $\frac{N}{\text{Scaling factor}}$.

2.2.7.   Prey population model with heritable trait diversity

Interactions between heritable trait diversity and population resilience is still unclear. First, high diversity can increase population resilience. Resilience may increase with phenotypic diversity, as a greater array of responses to environmental changes are available ^[10,11]. Other researchers argue that the link between resilience and genetic diversity is unclear ^[18,23]. The fitness of a population depends on non-linear interactions between genetic diversity and the environment.

Our last approach uses heritable trait diversity in the prey population, using partial differential equations ^[24]. This partial differential equations (PDE) model uses the same functions as in Eq (2.10), and integrated a parameter of stochastic fluctuations of population density.

Thus, in (2.12), the dynamics of a prey $\dfrac{\partial}{\partial t}N(t, s)$ depends on its own trait $s$, and on the entire diversity of prey with a different trait $x$. Evolution then acts through selections from the diversity of available trait values. Mutation events are assumed to provide new diversity by slightly modifying the trait value from parents to offspring. The diffusion term $\Delta N(t, s)$ represents incoming and outgoing mutations for the trait $s$, weighted by the mutation rate $\mu$. Environmental variability is introduced as an additional mortality $\bar{n}$ (2.12). Mortality increases if $N(t, s)$ is of the same order as $\bar{n}$, such that rare traits are eliminated by environmental variability. This model was implemented as a numerical scheme under Python version 3.10.9.

The trait distribution, ranging from 0 to $s_{max}$, is divided into boxes of equal length. Each box $i$ represents a range of traits of index $i$, and contains a prey density $N_i$. Functions use the median trait within boxes. Thus small and numerous boxes allow for better precision of the trait distribution.

Mutations appear in the trait distribution, and are computed as diffusion between neighboring boxes. Diffusion is modeled as the second derivative of population density $\frac{\partial N}{\partial t}$ weighted by the mutation rate $\mu$, such as a flux of offspring between boxes. To close the distribution boundaries, the two boxes at the borders of the distribution can diffuse only toward the center.

Environmental variability was introduced as an additional mortality $\bar{n}$ (2.12). Parameter $\bar{n}$ is applied to every box such that mortality increases if density $N_i$ is too low (e.g., of the same order of $\bar{n}$). Such mortality parameter allows to observe the influence of environmental variability on heritable trait diversity and density in a deterministic model.

For better numerical stability during simulations, Eq (2.13) are defined by the difference between growth and mortality rates in (2.14) (\(R_i^k - \sum_{ j } K_{ i j }\)) at every time step $\Delta_t$. If growth is higher, equation (+) is used to compute density $N_i^k$ at the next time step (using $N^k_i$). Otherwise, equation (-) is used (using $N^{ k + 1 }_i$).

3.   Results and discussion

3.1.   Evolution direction and steady state

As argued by the researchers in ^[7], evolution could both delay, or lead to early regime shifts. Here, we speculate on the effects of evolution under an increase of predation pressure (Figure 2). An increase in predation pressure decreases the population density of prey. Such change in abundance influences evolution by modifying the pressures of mortality and birth on fitness ^[25]. Here, a reduction in population density reduces encounters, and, thus, competition pressure between prey ^[26,27,28].

Figure 2.  Hypotheses on resilience changes through evolution of the prey. First, if evolution is driven by competition (upper part), an increase in predation pressure decreases competition pressure on prey. The prey can invest more energy in reproduction, and the high population density state becomes more resilient. Second, if adaptation is driven by resistance to predation (lower part), an increase in predation pressure induces a selection of more resistant prey. The population reproduction rate decreases, leading to a more resilient low density state.

DownLoad: Full-Size Img PowerPoint

An increase of predation pressure increase could thus have one of two effects. First, the increase in predation pressure leads to a decrease of competition in prey. Less competition pressure allows prey to invest more energy in reproduction. In such cases, the access to resources is improved for prey, and favors the resilience of high population densities (Figure 2). In the second case, an increase of predation pressure induces selection of resistant prey ^[4]. Such evolution costs a decrease of the reproduction rate, and favors early shifts to low population densities.

Moreover, ^[29] observed that a reduction in size at maturity, paired with increased fecundity, led to increased resilience of fish populations under fishing pressure. Thus, qualitative knowledge on ecological and evolutionary interactions should help predict population regime shifts on eco-evolutionary time scales.

3.2.   Population shifts and recovery under evolution

Intense management efforts can restore the population's initial state after a regime shift. For instance, in the Caribbean, overfishing of the parrot-fish (Sparisoma viride) is expected to lead to a dominance of algae and a loss of corals ^[2]. Efforts to reduce fishing pressure, and restore parrot-fish populations, reduce algal dominance and preserve coral reefs.

However, evolution is suspected to affect populations or ecosystem recovery ^[7]. After a regime shift, populations adapt to new conditions, and the path of recovery changes. Adaptation to new conditions modify the tipping points toward the initial state.

The adaptive dynamics model (2.10) simulated scenarios of regime shifts and recovery. Simulations started from the eco-evolutionary equilibrium of the prey on the low density state (or high density). Then, a constant negative (or positive) predation pressure change rate of $\pm 10^{-3}$ was applied to induce shifts. Different maximum predation changes were tested to cross past the tipping point and mimic a delay before initiating recovery. Then, recovery began by reversing the sign of the predation pressure change rate, and recovering the initial $B$ value. These different recovery experiments were tested under different rates of prey evolution. Then, relative changes of prey density between initial eco-evolutionary equilibria and final ecological equilibria after recovery were computed as $\frac{N_{final} - N_{initial}}{N_{initial}}$ (see Figure 3).

Figure 3.  Experiments of predation change and recovery after a regime shift. On the bottom, the white shows the bistability region in the predation and prey trait space. The grey areas show the regions with only one density state (low or high). After reaching thresholds in change of predation pressure (see $\bigstar$), the change rate's sign is inverted to start recovery. Dots ($\bullet$) indicate the starting point of shift experiments from the eco-evolutionary equilibrium. For panel (a), the starting point of all simulations is $B = 1$ and $s = 4.5$. For panel (b), starting points are $B = 5.7$ and $s = 33.2$. Arrows show examples of shift then recovery without evolution (in black), or with evolution (in red). In (a) and (b), different evolution rates are shown against different changes in predation pressure. Relative changes of prey density between the initial eco-evolutionary equilibria and the final ecological equilibria are reported in the color plots. The parameter values are $r_{max} = 2$, $D_r = 0.001$, $c_{\alpha} = 0.2008$, $k = 0.001$, $V = 0.64$, $K_c = 4.06$, $D_c = 40.6$, and a change of predation pressure of $\pm 10^{-3}$.

DownLoad: Full-Size Img PowerPoint

In the case without evolution, there was no difference between prey density at the initial state and density after recovery ^[30]. However, enabling prey evolution modified the tipping points ^[7], and induced differences of density after recovery. In Figure 3, tipping points depended on both predation pressure and on the prey trait $s$. For low $s$ values, the bistability region encompassed a large range of predation pressure, while tipping points were closer as $s$ increased. Thus, the position of tipping points depended on the ratio between the rate of predation pressure change, and the rate of evolution. Tipping points were further modified as evolution was faster. Thus, evolution modified the path of recovery in the trait-predation space as prey continued to evolve after regime shifts ^[7].

Modifications of the path of recovery depended on the regime shift considered. First, let us consider the shift from high prey densities (low predation pressure) to low densities (high predation pressure, see case Figure 3(a)). On high densities evolution was much faster due to large population sizes and higher mutation counts. On the contrary, while on low densities, evolution speed decreased. In Figure 3(a), further predation change after shifting to low densities had little effect on the path of recovery. Relative changes of density were mostly influenced by fast evolution rates while on high densities. After recovery, initial eco-evolutionary equilibria could be recovered given sufficient time and could result in little changes of density. Thus, delays before initiating recovery from shifts to low densities had little effect on changes in final density.

In the opposite case in Figure 3(b), starting from low densities, evolution was much faster after shifting to high densities. Thus, the path of recovery mostly depended on the time spent after the shift. In this case, both delays before initiating recovery, and rates of evolution and predation change influenced density equilibria after recovery. Furthermore, retrieving initial eco-evolutionary equilibria after recovery to low densities was difficult due to slow evolution. Thus, modifications of tipping points and population density brought by delays of recovery were long-lasting in this case.

As argued by the researchers in ^[30], regime shifts were sensitive to the ratio between the rate of evolution and the rate of environmental change. Fast rates of environmental change reduced modifications brought by evolution. Recovery strategies should thus consider the initial state of the system (and ecological interactions), the rate of environmental change, and generation time or the rate of evolution. For instance, in Figure 3 the two types of regime shifts showed different dynamics of evolution and of recovery. These differences mostly depended on the position of high densities, if these were the initial condition, or resulted from shifts. Thus, rates of evolution and of environmental change, as well as interactions between density and evolution speed, alter a population's density after restoration.

Here, we discussed the effects of heritable trait diversity. However, phenotypic diversity in populations may also result from phenotypic plasticity, which unfolds on shorter time scales than heritable trait selection. Thus, changes in tipping points are expected to happen on ecological time scales, and to be faster in the case of plasticity. Thus, phenotypic plasticity would reduce prey density changes after recovery. However, regime shifts to low population densities would also induce a loss of genetic diversity, and thus alter the extent of plasticity ^[18]. We speculate that interactions between phenotypic plasticity and regime shifts would require to study genetic diversity as well, but such questions are beyond the scope of this study.

3.3.   Limit cycles maintained by evolution

Cycles of prey density were observed under no change in predation pressure. Prey evolution was unable to reach a stable eco-evolutionary equilibrium, and thus evolution drove eco-evolutionary cycles of the trait and of population density. Cycles of density also influenced evolution speed, as evolution was faster on high densities ^[31]. Furthermore, evolutionary constraints differed between high, and low prey densities ^[32]. For instance, high densities implied high rates of encounter between prey, and predators. More encounters resulted in more events of predation and competition. To cope with competition and predation, prey evolved towards higher trait values. However, high $s$ values led to a reduction in reproduction. Eventually, such cost was too high, and the population shifted to low densities. For low densities, interactions between prey and predators weakened due to fewer encounters ^[4,32]. prey suffered less predation and competition, and the trait evolved back to lower values. Thus, evolution eventually led prey to shift back to high densities, and the cycle began anew.

As emphasized by these eco-evolutionary cycles, evolution did not necessarily stabilize the system ^[8]. In Figure 4, evolution could not satisfy the trade-off between competition, resistance to predation, and reproduction. As in ^[33], eco-evolutionary cycles depended on the dependence of intra-specific competition on the trait. The trait needed to satisfy a compromise. Fitness depended on the investment in competition and resistance to predation, or in reproduction. Thus, the trait's response depended on prey density, and on rates of encounter. Changes in encounters modified ecological interactions between prey, and predators. Such changes altered evolutionary pressures, and led to eco-evolutionary cycles of prey densities.

Figure 4.  In (a), limit cycle (red line) in the predation and prey trait space. The bistability region is white, and regions with only one density state (low densities or high densities) are grey. Adaptation maintains the limit cycle (in red), without environmental change (constant predation). The limit cycle is also shown in the density-trait space. The prey's adaptation maintains subsequent regime shifts at different speeds depending on population density. High densities are in red, and low densities are in blue. First, at high densities, evolution is fast. Large populations induce evolution towards high trait values. Second, at low densities evolution is much slower. Small populations induce evolution towards low trait values. Population density (b), trait values (c), and predation pressure (d) are shown at the beginning of the experiment. Density shifted to high values, inducing an increase of trait values at constant predation pressure. The parameter values are $r_{max} = 2.575$, $D_r = 0.25$, $c_{\alpha} = 0.001$, $k = 2.5$, $V = 0.55$, $K_c = 10$, $D_c = 1$, and $B = 8.6$.

DownLoad: Full-Size Img PowerPoint

Phenotypic plasticity may reduce the amplitudes and length of cycles in density as trait values change faster ^[34]. Plasticity induces fast adaptation, and such cycles could happen on ecological time scales. In such cases, plasticity could dampen changes in density ^[34], which would stabilize predator-prey interactions. However, plasticity may not completely prevent shifts and oscillations in population density. Decreases in population density leads to a loss of genetic variance and may in turn slow down the effects of plasticity ^[18].

In extreme cases, shifts to low densities led prey to extinction ^[33]. Extinction risks should be even more likely in cases of sudden, stochastic disturbances. Environmental variability, such as storms or disease outbreaks, could affect population survival. In such cases, a shift to low densities driven by evolution, followed by a storm, could lead to extinction ^[35]. Furthermore, low population densities favor inbreeding and reduction of genetic diversity ^[18]. Risks of extinction should be even higher, as low genetic diversity reduces chances of survival against disturbances, such as diseases.

3.4.   Stochasticity and early regime shifts

As time series of population densities are highly variable, we ask if such variations would set the prey in one or the other of the attraction basins ^[36]. To test this idea, prey density was monitored in IBM simulations (see Eq (2.11)) at constant predation pressure $B$ and constant trait $s$. Evolution and changes of predation were not allowed in this experiment. Simulations of prey density for different values of predation were then monitored (Figure 5).

Figure 5.  Frequency distributions of prey densities during regime shifts experiments using the IBM model, with no evolution. Colors go from blue, for low population densities, to red for high densities. For each experiment, the trait is set to $s = 2$, and the scaling factor to $10$. The predation pressure $B$ is set to constant values of $22.0$, $22.5$, or $23.1$. IBM simulations were run until $t_{final} = 10^5$ (2000 samples). These three experiments explore the flickering of densities depending on the relative width of attraction basins, controlled by $B$. Flickering appears as the population flips to a different state depending on the basins of attraction and stochasticity ^[36]. The parameter values are $r_{max} = 1$, $D_r = 0.01$, $c_{\alpha} = 0.026$, $k = 0.1$, $V = 0.75$, $K_c = 14$, $D_c = 0.01$, and $B = 8.6$.

DownLoad: Full-Size Img PowerPoint

Under stochastic simulations, prey densities oscillated from one attraction basin to the other, thus shifting populations between density states (Figure 5). These stochastic fluctuations depended on population densities. Large densities (or a large scaling factor) decreased the amplitude of fluctuations, and thus the probability to shift.

This probability also depended on predation pressure $B$, as predation modified the width of basins of attraction. For instance, low $B$ values increased the upper basin width in Figure 5. Thus, the resilience of high population densities increased with the size of the upper basin. Such changes in population resilience modified the probability of early regime shifts due to disturbances.

Thus, environmental changes modified the resilience of both low and high density states ^[3]. Populations with low resilience were sensitive to sudden disturbances, such as storms, that induced early shifts ^[17,36]. Implementing evolution added even more complexity. Evolution not only altered tipping points ^[7], but also the basins of attraction. Population resilience to disturbances was thus also sensitive to evolution (see Figures 1 and 4). As argued by the researchers in ^[17], anticipating changes in resilience could be a strategy for sustainable ecosystems management. Thus, management in the long term should consider evolutionary processes.

3.5.   Diversity and environmental variability

Here, we asked if regime shifts thresholds depend on prey phenotypic diversity. Thus, used the PDE model to simulate heritable trait diversity 2.13. The model integrated evolution and an additional mortality rate due to environmental variability. Different values of evolution rate and environmental variability were tested. For each one of these values, the PDE model simulated shifts at constant change rate of predation pressure and monitored the predation pressure $B$ at regime shifts.

The initial conditions of PDE simulations were estimated as the eco-evolutionary equilibria of the adaptive dynamics model 2.10. Trait equilibria were used as the mean value of the trait distribution, and density equilibria gave total population densities. Simulations then waited for the prey population to reach its diversity equilibrium before initiating the change of predation pressure (see the detailed protocol in Supplementary 2).

Shift thresholds depended on both the mutation rate, or trait diversity, and environmental variability. In Figure 6, higher environmental variability induced regime shifts at lower predation pressure. High diversity further amplified this change in shift thresholds. Shifting to high densities required stronger decreases of predation pressure, and recovering high densities required further effort. Such a delay for the shift to high densities was not due to evolution, as the mean trait did not change significantly before shifts (see Figure 6(b), (c)) On the contrary, shifts to low densities occurred at lower predation pressure values. Thus, high trait diversity and high environmental variability seemed to favor shifts to low densities (Figure 6(a)).

Figure 6.  Monitoring experiments of predation pressure $B$ at regime shifts using heritable trait diversity. The time step used is $\Delta_t = 0.01$, with a constant change rate of predation pressure of $10^{-2}$. (a) shows tipping points measured, either by decreasing $B$ (shift to high densities) or increasing $B$ (shift to low densities). Mutation rates and environmental variability $\bar{n}$ both modify tipping points. (b) and (c) show results for one simulation using $\bar{n} = 0.1$, and at different diversity levels (black for no mutations, and red for $\mu = 10^{-4}$). In (b), the population mean trait does not change significantly until the shift is reached. The parameters are $r_{max} = 1$, $D_r = 0.01$, $c_{\alpha} = 0.026$, $k = 0.1$, $V = 0.75$, $K_c = 14$, $D_c = 0.01$, and $\Delta s = 0.01$.

DownLoad: Full-Size Img PowerPoint

Different phenotypes have different sensitivities to environmental variability. The mean phenotype, or mean trait, is considered as the most optimized value. A different phenotype thus suffers stronger mortality. Under a highly variable environment, vulnerable phenotypes with the lowest densities suffer even stronger mortality. For instance, rare and vulnerable phenotypes are even more likely to go extinct following a severe storm. Therefore, high trait diversity, paired with strong environmental variability, increased the overall population mortality by eliminating the more vulnerable phenotypes. This increase of mortality translated to early shifts to the low densities. Thus, interactions between trait diversity and environmental variability modified shift thresholds.

Large population densities were more complex to describe due to fast evolution. Shift thresholds were modified by both evolution ^[7], and environmental variability. The two either had a similar effect, such as both favoring shifts to low densities, or had opposite effects. Thus, selection could lessen this effect and increase the resilience of high densities (see Figure 2). Such interaction depended on the choice of parameters. Further exploration would be needed to disentangle the effects of selection, increased mortality due to diversity, and regime shift dynamics.

Under environmental change, phenotypic and genetic diversity in populations favor resilience ^{[10,11,37,38,39,40]}. High diversity increased the pool of available phenotypes, and thus of available responses to environmental change. However, increased diversity had a cost through interactions with the environment (Figure 6). The most vulnerable phenotypes, in a given context, risked extinction under high environmental variability (e.g., storms, diseases) ^[17]. Such mortality events increased the probability of shifts to low densities.

Furthermore, interactions between diversity and environmental variability were not constant through time, as vulnerable phenotypes were eliminated during selection events ^[41]. Environmental variability and the recovery speed of intra-specific diversity are essential in determining modifications of shift thresholds. Too slow diversification leads to maladaptation, and ultimately higher risks of extinction ^[39]. Thus, there would be an equilibrium between diversification and population resilience. Rates of diversification could be of major importance to anticipate population shifts under environmental change. First selection events are fast, as they rely on the genetic standing stock ^[42], and would result in decreased genetic diversity ^[9]. Then, new mutations and recombination events would contribute to new genetic diversity. New mutations provide few genetic backgrounds, which would reduce genetic diversity under selection ^[42]. Genetic diversity may also increase due to recombination events, by producing different genotypes and combinations of alleles, at the risk of breaking advantageous combinations ^[43]. Thus, rates of recombination and mutations should depend on a balance between decreasing the risks of maladaptation, while not decreasing population resilience.

Thus, the diversity-variability link presented here supported the idea of complex interactions between population resilience and diversity ^[23,44]. A population's fitness not only depends on its genetic diversity, but also on the interactions between genetic diversity and the environment ^[18]. Adding more heritable traits or more species to the model, and further diversity, should lead to more complex dynamics ^[33]. Dynamics of species would then depend on environmental variability, eco-evolutionary feedbacks, and predators as well as competitor evolution ^[44].

4.   Conclusions

Here, we showed that evolution influenced population dynamics under environmental change. Evolution and ecology act together to shape communities and should not be considered separately ^[4]. Furthermore, eco-evolutionary feedbacks depend on trait diversity for the selection process, and congruent time-scales between ecology and evolution, so that diversification is fast enough ^[45]. Such interactions can occur on ecological time scales due to fast generations, such as in planktonic systems. Thus, monitoring mutation and diversification rates seem mandatory to quantify such eco-evolutionary dynamics. Moreover, the extent of such interactions also depends on the rate of environmental change ^[30]. Assessing eco-evolutionary dynamics is the next step towards anticipating population shifts in the contemporary context of environmental change.

Evolution thus modifies ecological interactions, and the position of tipping points ^[7]. The researchers in ^[46] argued that evolution can also lead to cross tipping points and cause shifts. Thus, evolution does not necessarily increase system resilience and the link between evolution and population resilience seems complex. Selection can negatively impact the population and lead to the population's suicide ^[35]. Furthermore, evolution modified population resilience, and thus the probability to shift following a disturbance. The researchers in ^[36] argued that stochastic events, such as storms, could lead to early regime shifts. Thus, evolution could increase risks of extinction due to stochasticity. Assessing evolution seems crucial to anticipate changes in populations and risks of evolutionary suicide in variable environments ^[23,35].

Under environmental change, first selection events are fast as they rely on genetic standing stocks ^[42]. New genetic diversity is thus needed for subsequent selections, as low diversity is expected to induce maladaptation and increase extinction risks ^[39]. Increases of genetic diversity are created by gene flows between populations, although gene flows decrease the speed of natural selection ^[47], or by new mutations and recombination events. However, mutations can be deleterious, and recombination can also break advantageous genotypes ^[43]. Too high mutation and recombination rates should thus produce vulnerable phenotypes and reduce population resilience.

Heritable trait diversity had two major outcomes on the system. First, diversity enhanced the probability for species to recover following a disturbance. A greater phenotypic diversity increases the probability for an optimal phenotype being present. Second, diversity had a cost. Too much diversity negatively impacted populations. Diversity is an investment, and the presence of vulnerable phenotypes increases the population mortality rate. heritable trait diversity thus alters population resilience under disturbances and environmental change.

Thus, diversity not only depends on the mutation rate, but also on environmental variability ^[12]. For instance, high variability can promote a large trait diversity that can cope with different situations. Populations can adopt a mixture of phenotypes adapted to storms, and phenotypes adapted to more stable conditions. The proportion of phenotypes depends on the frequency of storms ^[12]. A change in the frequency of extreme events should influence phenotypic diversity, population resilience, and regime shift dynamics. Thus, we speculate that genetic diversity follows a balance between the rate of diversification and population resilience.

On top of trait diversity in prey, predators and competitors also evolve and have their own dynamics. Species evolution depends on ecological interactions and in turn impacts community and ecosystem dynamics ^[44]. Here, we assumed that gaining an advantage in competition also reduces predation pressure. However, both processes can also go opposite evolutionary directions. Predators may target the "fittest" trait with the most abundance, which may counteract the advantages of competition ^[48]. Furthermore, interactions such as cooperation between different species of prey, or even within populations, may hamper competition and predation pressures ^[49]. Ecological interactions may thus have a stabilizing effect on populations through the feedback on evolutionary pressures ^[50]. Thus, modifications of ecological interactions between species would modify evolutionary pressures ^[44,51], modify the dynamics of tipping points, and of basins of attraction (and resilience). As an added level of complexity, species may also change niche during development. Such ontogenetic shifts may affect population resilience, such as observed in fish populations ^[29], and even alter community and ecosystem dynamics ^[52]. Replacing population eco-evolutionary dynamics in an ecosystem context would require further study, which is beyond the scope of this study.

Hysteresis is a useful tool for stock management ^[53], and considering evolution together with ecological interactions is necessary to better anticipate regime shifts. Furthermore, trait diversity needs to be considered, especially in variable environments. Trait diversity (or genetic diversity) induced non-linear interactions with population resilience. Resilience depends on the fitness of the different traits in a given environment ^[18], and deviations from the "fittest" trait result in greater mortality, a decrease in population resilience, and higher risks of extinction. As the link between diversity and population resilience is non-linear, more knowledge on the feedback between diversity and ecology is needed. Thus, evolutionary dynamics, such as presented here, could provide a new tool to anticipate population regime shifts. Thus, we need to better assess links between diversity and regime shift dynamics to manage ecosystems. The present framework may represent a small step towards integrating complex eco-evolutionary dynamics and genetic diversity in a context of population regime shifts under environmental change.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was supported by an internship at Aix-Marseille Université.

Conflict of interest

The authors declare there is no conflicts of interest.