<figure><img src="image_1.png" /></figure>

HAL
open science

# Evolution between two competing macrophyte populations along a resource gradient leads to collapse in a bistable lake ecosystem

Sirine Boucenna, Gael Raoul, Vasilis Dakos

► To cite this version:

Sirine Boucenna, Gael Raoul, Vasilis Dakos. Evolution between two competing macrophyte populations along a resource gradient leads to collapse in a bistable lake ecosystem. Theoretical Population Biology, 2024, 164, pp.23-36. 10.32942/X2532Z. hal-04794239

HAL Id: hal-04794239

https://hal.science/hal-04794239v1

Submitted on 22 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Highlights

Evolutionary outcomes arising from bistability in ecosystem dynamics

Sirine Boucenna, Gael Raoul, Vasilis Dakos

— Rapid evolution of the growth depth of submerged macrophytes promotes diversifica-
tion and coexistence with competitive species.

— Rapid evolution of the growth depth of submerged macrophytes in a bistable ecosystem doesn't prevent collapse and leads to evolutionary oscillations and evolutionary suicide.

— Co-evolution of the growth depth of floating and submerged macrophytes plays a sta-
bilizing role and delays the evolutionary suicide of the submerged macrophytes.
Evolutionary outcomes arising from bistability in ecosystem dynamics

Sirine Boucenna<sup>a</sup>, Gael Raoul<sup>b</sup>, Vasilis Dakos<sup>a</sup>

<sup>a</sup>Institut des Sciences de l'Evolution de Montpellier (ISEM), Université de Montpellier, CNRS, IRD, EPHE, Montpellier, France

<sup>b</sup>Centre de Mathématiques Appliquées, École Polytechnique, Palaiseau, 91120, France

**Abstract**

While it is known that shallow lakes ecosystems may experience abrupt shifts (ie tipping points) from one state to a contrasting degraded alternative state as a result of gradual environmental changes, the role of evolutionary processes and the impact of trait variation in this context remain largely unexplored. It is crucial to elucidate how eco-evolutionary feedbacks affect abrupt ecological transitions in shallow lakes. These feedbacks can significantly alter the dynamics of aquatic plants competition, community structure, and species diversity, potentially affecting the existence of alternative states or either delay or expedite the thresholds at which these ecological shifts occur. In this paper, we explore the eco-evolutionary dynamics of submerged and floating macrophytes in a shallow lake ecosystem under asymmetric competition for nutrients and light. We use adaptive dynamics and a structured population model to analyze the evolution of the growth depth of the submerged and floating macrophytes population, which influences their competitive ability for the two resources. We show how rapid trait evolution can result in complex dynamics including evolutionary oscillations, extensive diversification and evolutionary suicide. Furthermore, we find that the co-evolution of the two competitive species can play a stabilizing role, while not significantly affecting the overall evolutionary dynamics. Overall, this study shows that evolution can have strong effects in the ecological dynamics of bistable ecosystems.

*Keywords:* evolution, adaptive dynamics, structured population model, evolutionary suicide, evolutionary cycles, co-evolution.
# 1. Introduction

Extensive theoretical work and empirical observations have revealed the potential for ecosystems to undergo abrupt and dramatic transitions, shifting from one stable state to a contrasting degraded alternative state [23, 36, 48]. Some of the best-known examples of ecosystems with alternative stable states include lakes transitioning from a clear to a turbid water state [44], drylands transitioning from vegetation to desert [29], or coral reefs shifting from a healthy coral dominance to being overgrown by macro-algae [18, 24]. These transitions occur when external conditions exceed critical thresholds known as "tipping points" [34, 51]. Tipping points come in various forms, such as bifurcations, noise-induced shifts, and rate-dependent transitions [3]. In particular, bifurcation tipping points typically correspond to pitchfork, subcritical Hopf, or saddle-node (fold) bifurcations [16] [40]. It is important to understand tipping point responses, because tipping points not only lead to abrupt changes but these changes can be hard to reverse. That is, if the conditions change in the opposite direction, the system will remain in the alternative state until it reaches another tipping point. The difference between the two tipping points is known as hysteresis and the existence of multiple alternative states is referred to as bistability [43]. As a result, there is mounting concern regarding the crossing of tipping points in environmental conditions, which can trigger catastrophic and potentially irreversible shifts in ecosystems.

In shallow lakes, alternative stable states refer to different ecological states that a lake system can attain and persist in. These states can be distinguished based on their physical, chemical, and biological characteristics, such as water clarity, phytoplankton biomass, and dominance of different species of aquatic plants and animals. The two most common alternative stable states in shallow lakes are related to a clear-water and a turbid-water state [47]. A clear-water state is characterized by a low nutrient concentration (typically phosphorus and nitrogen), which results in minimal algal growth, a dominant vegetation cover of submerged
macrophytes, such as pond weeds or water milfoil, which stabilizes the sediment, provides shelter for phytoplankton consumers and altogether enhances water clarity. A turbid-water state is characterized by a high nutrient concentration that leads to excessive algal growth that shades out submerged macrophytes and reduces their growth. Loss of submerged macrophytes can lead to increased sediment suspension that further reduces water clarity and further inhibits the growth of submerged macrophytes [4]. This chain of events create a strong positive feedback that results in the emergence of alternative stable states in shallow lakes where shifts between clear and turbid water states occur in response to changes in external drivers, such as nutrient input or extreme weather conditions including droughts [45, 48] [17].

Various models of alternative stable states have been proposed in shallow lakes. The most common case describes the competition between phytoplankton and macrophytes in order to determine the level of nutrient loading at which a transition may occur, showed a highly non-linear response and hysteresis [26]. An aquatic Lotka-Volterra type food web model in shallow lakes was also used by [32] and showed the occurrence of alternative stable states. During eutrophication, the macrophyte-dominated clear-water state marked by a low level of chlorophyll-a disintegrates abruptly when the critical phosphorus loading is reached, shifting the system to a phytoplankton-dominated state. Another nutrient-phytoplankton-zooplankton model was used in [53] where discontinuous Hopf bifurcations were found and therefore bistable phenomenon occurred with the presence of a stable equilibrium and a stable limit cycle. In another model of competition between blue-green filamentous cyanobacteria and algae used in [42, 50], alternative stable states occurred in response to turbidity where a cyanobacteria dominance state was found in shallow lakes. Additionally, a predator-prey model of plankton dynamics describing algae and zooplankton densities [46] showed that switches from one regime to another occur abruptly at a critical fish density. Bistability was also described in a model of competition for nutrients and light between free-floating macrophytes and submerged macrophytes in a shallow lake [49]. Bistability in this ecosystem
is a consequence of a strong positive feedback between free-floating macrophytes and sub-
merged macrophytes, where floating macrophytes take precedence over light on the surface
of the water. Dark anoxic conditions beneath dense layers of floating vegetation offer little
habitat opportunity for flora and fauna and can have significant adverse effects on fishing
and shipping in tropical lakes [49]. Conversely, submerged macrophytes can thrive at lower
nutrient levels, thereby reducing nutrient availability and diffusion to the surface [25]. Thus,
in response to nutrient loading, a shallow lake can shift from a coexistence state between
submerged and floating macrophytes to a floating macrophytes dominance state.

But while ecological tipping points between alternative states have received extensive atten-
tion especially in the case of shallow lakes, there has been a significant neglect of evolutionary
dynamics and trait changes, which also play a significant role in shaping the dynamics and
stability of ecosystems [22] [35]. Large-scale ecological shifts triggered by tipping points can
have profound impacts on species interactions, resource availability, and overall ecosystem
functioning. These shifts can create selective pressures that influence the adaptive evolution
of species traits. At the same time, trait dynamics, including genetic variation, mutation,
and selection, can influence the capacity of populations to respond to ecological shifts [12].

Only little theoretical work exists in studying how trait evolution and trait variation af-
fects of shapes the bistable dynamics of ecosystems. For instance, fast environmental change
mediated by high trait variation and eco-evolutionary feedbacks could drive regime shifts
in ecosystems even before tipping points are crossed [8]. The introduction of evolutionary
process can qualitatively affect the response of ecosystems with alternative stable states to
increasing environmental stress and makes these ecological systems also sensitive to the rate
at which environmental changes occur. Another study about the influence of evolution on the
collapse and recovery pathways in bistable ecosystems using quantitative genetics in shallow
lakes found that adaptive evolution has the capacity to either advance or delay thresholds at
which a tipping points occur [7]. These alterations of trajectories of ecosystem collapse and recovery were proposed as a potential consequence of trait change in a bistable ecosystem [12]. Additionally, [13] found that including evolution of a trait that relates to competitive performance in a system where multiple ecological equilibria coexist can lead to adaptive reversals that drive cyclic alternation between alternative stable states. Work on arid ecosystems shows that the evolution of local facilitation does not lead to critically low population densities that would put the population at risk of extinction from random processes, but populations evolve towards trait values that are not viable [28]. Interestingly, these populations consistently exhibit alternative stable states. This finding aligns with the results of [21], where, in the case of a metapopulation model, if the transition to extinction is continuous, the extinction boundary acts as a repellent for evolutionary dynamics. In such cases, there is a very low population density near the extinction boundary, allowing only mutants that are viable without the resident population to successfully invade. This mechanism keeps the strategy within a viable range, preventing extinction due to adaptation. Thus, a sudden transition to extinction, represented by alternative stable states, is a necessary condition for the evolution towards extinction, although it is not sufficient on its own.

Moreover, other studies have shed light on selection-induced regime shifts and the role of evolution in the occurrence of bistability. We know for example that a small ecological change can set the stage for a future tipping point where significant reduction in mortality rates initiated the evolution of body size, leading to a population shift over the long term [9]. Furthermore, studies indicated that the results of evolutionary rescue models can be influenced by the choice of fitness function [39, 31]. They showed that there exist evolutionary tipping points that, when crossed by increasing the rate of environmental change, abruptly lead the population to extinction. This feature (among many others) was also characterized analytically in [19] that investigated more broadly the influence of the selection functions on the adaptation of sexually and asexually reproduction populations to a changing envi-
ronment. Some other complex frameworks have been shown to lead to tipping points. For
example, including an age structure to the population [10] [11] allows for feedback loops bet-
ween the dynamics of the demography and the trait dynamics to create multiple co-existing
equilibria that promote evolutionary tipping points. Another recent work [2] identified the
necessary ecological mechanisms for the evolutionary emergence of alternative stable states
using adaptive dynamics to study the evolution and diversification of a population of sub-
merged macrophytes in a shallow lake. It showed that the emergence of different phenotypes
and alternative stable states required a trade-off between nutrient availability and light avai-
lability, asymmetry in competition among individuals, and the presence of a priority effect.
But aside from these three mechanisms, the occurrence of a tipping point also depends on
the relative position of the different phenotypes.

In this paper, we integrate evolutionary dynamics into ecosystem models exhibiting ecologi-
cal bistability to gain a deeper understanding of the influence of evolutionary processes on
bistable ecosystems. We address two research questions. First, what are the evolutionary out-
comes in an ecosystem that exhibits bistability? Particularly, does evolution drive the system
towards or away from regions with alternative stable states, and if so, how does evolution
affect the ecological dynamics of the ecosystem? Second, what role does co-evolution play in
shaping these evolutionary outcomes, and how does it influence the ecological dynamics of
a bistable ecosystem? We investigate these two questions in the context of bistable shallow
lakes. We consider a mathematical model of competition for light and nutrients between
two types of aquatic plants : floating and submerged macrophytes. We then use adaptive
dynamics and a structured population model to study the evolution of a phenotypic trait of
the two plants that characterizes their competitive ability, and to investigate the evolutio-
nary and co-evolutionary behavior of the two types of macrophytes and their effects on the
ecological dynamics of a shallow lake ecosystem.
## 2. Materials and methods

### 2.1. Model description

We describe the eco-evolutionary dynamics of the populations of two dominant plant types: submerged macrophytes **S** and floating macrophytes **F** in a shallow lake that compete for nutrients and light. Both species are characterized by a phenotypic trait **z** that represents their growth depth and underlies their competitive ability for the two resources. We assumed in the first part of our study that (by definition) the floating macrophytes are fixed at the surface of the lake and have a phenotypic trait $z_f = 0$, while that the submerged macrophytes can grow along the water column with a phenotypic trait $z_s \in [0, z_b]$ where $z_b$ represents the depth of the lake.

We used an ecological model based on [49, 2], defined by the following ordinary differential equations (ODEs):

$$ 
\begin{cases}
\frac{dS(t,z)}{dt} = rS \frac{n(z,S,F)}{n(z,S,F)+h_s} \frac{1}{1+a(z,z)S+a(z,0)F+W(z)} - L_s(z)S \\
\frac{dF(t)}{dt} = rF \frac{n(0,S,F)}{n(0,S,F)+h_f} \frac{1}{1+a(0,0)F} - L_f(0)F
\end{cases}
\quad (2.1) $$

Where $S(t, z) \in \mathbb{R}^+$ denotes the population size of the submerged macrophytes at time $t \ge 0$ with a trait $z \in [0, z_b]$ and $F(t) \in \mathbb{R}^+$ is the population biomass of the floating macrophytes. $r$ is the maximum growth rate of the two populations, and the function $n$ represents the nutrient limitation with half saturation $h_s$ (respectively $h_f$) for **S** (respectively for **F**) and mortality functions $L_s(z)$ and $L_f(z)$. Light limitation depends on water turbidity $W(z)$ at the level $z$ and the shading effect induced by conspecifics is represented by the function $a$.

floating macrophytes are found on the water surface with their roots not attached to the sediment of the lake. However, submerged macrophytes completely grow underwater with roots attached to the bottom. This difference makes submerged macrophytes have a priority for access to nutrients, which reduces their availability in shallower depths, while floating
macrophytes have an advantage on the access to light due to their position and have a sha-
ding effect on deeper submerged macrophytes, which reduces their access to light [49] [2].
Also, we assumed that the submerged macrophytes do not exert a shading effect on the
floating macrophytes. Consequently, the competition for light in the population of floating
macrophytes is solely intrapecific, meaning it occurs among individuals of the same species.

Based on the above, attenuation of light followed the Beer-Lambert law as described in *W* 2.1, where *w* determines the strength of light attenuation with depth *z*. Nutrients stored at the sediment of the lake diffuse to the surface with nutrient diffusion strength determined by *u*. Submerged macrophytes' priority effect on nutrients that negatively impacts floating macrophytes is represented by the function *q* which increases with the growth depth *z* according to the value *p*. In turn, the floating macrophytes limit submerged macrophytes growth because of the shading effect, whose steepness is determined by the parameter *b* which represents the strength of asymmetry in competition for light.

Lastly, we defined the functions Ls and Lf so that the submerged macrophytes do not grow
near the surface and to prevent the floating macrophytes from reaching deep levels to keep
the two plant types distinct. We set the same inflection point for the two step functions
at the trait threshold zthres below/above which the population has a high mortality and
cannot survive. In order to tally with previous studies [49, 2], we assumed that the mortality
functions had the same minimum value of 0.05.
<table><caption>TABLE 1 – Definitions of the functions of the trait-dependent mechanisms in the model</caption><thead><tr><th>Function</th><th>Mechanism</th></tr></thead><tbody><tr><td>n(z, S, F) = <sup>Nu(z)</sup>&frasl;<sub>1+q(z)S+q(0)F</sub></td><td>Nutrient limitation</td></tr><tr><td>a(z<sub>i</sub>, z<sub>j</sub>) = <sup>2a<sub>0</sub> exp(b(z<sub>j</sub>-z<sub>i</sub>))</sup>&frasl;<sub>1+exp(b(z<sub>j</sub>-z<sub>i</sub>))</sub></td><td>Shading effect</td></tr><tr><td>Nu(z) = <sup>Nu<sub>0</sub></sup>&frasl;<sub>u&radic;2&pi;</sub> exp(<sup>-(z<sub>b</sub>-z)<sup>2</sup></sup>&frasl;<sub>u<sup>2</sup></sub>)</td><td>Nutrient concentration available</td></tr><tr><td>q(z) = q<sub>0</sub> exp(pz)</td><td>Rate of nutrient absorption</td></tr><tr><td>W(z) = W<sub>0</sub> exp(wz) - W<sub>0</sub></td><td>Light limitation</td></tr><tr><td>L<sub>s</sub>(z) = l<sub>max<sub>s</sub></sub> - (l<sub>max<sub>s</sub></sub> - l<sub>min<sub>s</sub></sub>) <sup>-z<sub>thres</sub></sup>&frasl;<sub>z<sub>thres</sub><sup>i</sup>+z<sup>i</sup></sub></td><td>S mortality function</td></tr><tr><td>L<sub>f</sub>(z) = l<sub>max<sub>f</sub></sub> - (l<sub>max<sub>f</sub></sub> - l<sub>min<sub>f</sub></sub>) <sup>-z<sub>thres</sub></sup>&frasl;<sub>z<sub>thres</sub><sup>i</sup>+z<sup>i</sup></sub></td><td>F mortality function</td></tr></tbody></table>

<table><caption>TABLE 2 – Model variables and parameters, their definition and default values</caption><thead><tr><th>Variable/Parameter</th><th>Definition</th><th>Default value</th></tr></thead><tbody><tr><td>z</td><td>Growth depth</td><td>∈ [0, z<sub>b</sub>]</td></tr><tr><td>z<sub>b</sub></td><td>Bottom of the lake</td><td>10</td></tr><tr><td>r</td><td>Maximum macrophyte growth rate</td><td>0.5</td></tr><tr><td>l<sub>max<sub>s</sub></sub></td><td>Maximum of the mortality function of S</td><td>0.1</td></tr><tr><td>l<sub>min<sub>s</sub></sub></td><td>Minimum of the mortality function of S</td><td>0.05</td></tr><tr><td>l<sub>max<sub>f</sub></sub></td><td>Maximum of the mortality function of F</td><td>0.05</td></tr><tr><td>l<sub>min<sub>f</sub></sub></td><td>Minimum of the mortality function of F</td><td>0.1</td></tr><tr><td>h<sub>s</sub></td><td>Submerged macrophyte dependency on water nutrient</td><td>0.1</td></tr><tr><td>h<sub>f</sub></td><td>floating macrophytes dependency on water nutrient</td><td>0.2</td></tr><tr><td>Nu<sub>0</sub></td><td>Total nutrient concentration in the water in the absence of plants</td><td>50</td></tr><tr><td>q<sub>0</sub></td><td>Surface macrophyte impact on water nutrient content</td><td>0.005</td></tr><tr><td>W<sub>0</sub></td><td>Baseline light attenuation and light limitation</td><td>1</td></tr><tr><td>a<sub>0</sub></td><td>Intra-specific light competition coefficient</td><td>0.01</td></tr><tr><td>p</td><td>Strength of priority effect</td><td>0.5</td></tr><tr><td>z<sub>thres</sub></td><td>Threshold trait and inflection point of the mortality functions</td><td>2</td></tr><tr><td>i</td><td>Strength of slope of the mortality functions</td><td>10</td></tr><tr><td>w</td><td>Strength of light attenuation</td><td>0.1</td></tr><tr><td>b</td><td>Strength of asymmetry in intra-specific competition for light</td><td>10</td></tr><tr><td>&sigma;<sup>2</sup></td><td>Initial phenotypic variance of S</td><td></td></tr><tr><td>&epsilon;</td><td>Seeding value</td><td>0.001</td></tr><tr><td>u</td><td>Nutrient diffusion value</td><td>∈ [0, 10]</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>FIGURE 2.1 – Mechanism functions : (a) : Water turbidity $W$ that increases with the depth. Impact of the plants of nutrient content and priority effect of the submerged macrophytes $q$. Nutrient availability in the water column $Nu$ when the nutrient diffusion $u = 5$. (b) : Shading effect for light competition according to the distance between conspecifics. (c) : Mortality functions of $S$ and $F$ that prevent the submerged macrophytes to grow above the inflection point threshold $z_{thre}$ and floating macrophytes to grow below it.</figcaption></figure>

### 2.1.1. *Adaptive Dynamics - Evolution of Submerged macrophytes only*

We used Adaptive Dynamics (AD) to investigate the evolution of the phenotypic trait $z$ of the submerged macrophyte population [14] [20] [38] [15]. We recall that the floating macrophytes do not evolve and are fixed at the surface of the lake. Initially, the populations $S$ and $F$ are at their ecological equilibria determined by the external conditions. We define the fitness function of a mutant $S_m$ with trait $z_m$ as its long-time growth rate in an environment set by the resident populations $(S^*, F^*)$:
$$s(z_m, z) = r \frac{n(z_m, S^*, F^*, S_m)}{n(z_m, S, F, S_m) + h} \frac{1}{1 + a(z_m, z)S^* + a(z_m, 0)F^* + W(z_m)} - L_s(z_m) \quad (2.2)$$

where $S^*$ and $F^*$ are the ecological equilibria of the populations of the submerged and floating macrophytes.

The direction of evolution is determined by the sign of the selection gradient defined as :

$$\frac{\partial}{\partial z_m} s(z, z) = \frac{\partial}{\partial z_m} \left[ s(z_m, z) = r \frac{n(z_m, S^*, F^*, S_m)}{n(z_m, S, F, S_m) + h} \frac{1}{1 + a(z_m, z)S^* + a(z_m, 0)F^* + W(z_m)} - L_s(z_m) \right]_{z=z_m} \quad (2.3)$$

A Singular point $\bar{z}$ is a trait value at which the selection gradient vanishes

$$\frac{\partial}{\partial z_m} s(\bar{z}, \bar{z}) = 0. \qquad (2.4)$$

In the context of AD, a singular strategy refers to an evolutionary stable strategy (ESS) that cannot be invaded by any alternative strategy in a population. In our model, this means that if the majority of submerged macrophytes have a trait corresponding to the ESS, no rare mutant with a different trait can grow.

In order for $\bar{z}$ to be a local ESS, $\bar{z}$ should be a singular strategy meaning that it should satisfy the condition [15] :

$$\frac{\partial^2}{\partial z_m^2} s(\bar{z}, \bar{z}) < 0 \qquad (2.5)$$

The last condition is related to the stability of $\bar{z}$ for monomorphic populations. A singular strategy is also Convergence Stable Strategy (CSS) if the resident monomorphic population S has a strategy $\tilde{z}$ close to $\bar{z}$ such that successive beneficiary mutations bring the population's trait towards $\bar{z}$ [15]. In order to have a CSS, the selection gradient should be positive for $z < \bar{z}$ and negative for $z > \bar{z}$. Thus, a strategy $\bar{z}$ is a CSS if it is a singular strategy and

$$0 > \left[ \frac{\partial}{\partial z} \left( \frac{\partial}{\partial z_m} s(z_m, z) \right) \right]_{z=z_m} = \left[ \frac{\partial^2 s}{\partial z \partial z_m} + \frac{\partial^2 s}{\partial z_m^2} \right]_{z=z_m=\bar{z}}. \quad (2.6)$$
An interesting case occurs when mutants can invade the resident population both for larger and smaller traits. This mutual invasibility case is possible when $\bar{z}$ is CSS but not ESS, resulting in a branching point (BP) where from monomorphic population more than one phenotypes emerge. For a BP to occur, the singular strategy must apply [15] :

$$ \frac{\partial^2}{\partial z^2} s(\bar{z}, \bar{z}) > \frac{\partial^2}{\partial z_m^2} s(\bar{z}, \bar{z}) > 0 \qquad (2.7) $$

### 2.1.2. Structured population model of phenotypic emergence - Evolution of Submerged macrophytes only

The emergence of multiple phenotypes and their specific characteristics during branching points are not captured by AD. Instead AD only allows to estimate long-term evolutionary attractors and their classification [15]. To understand the nature of emerged phenotypes, we use a quantitative genetics approach, initially developed by [30] which considers a continuous trait distribution and frequent mutations (generally assumed Gaussian with a fixed variance). In our study we used a structured population model with selection-mutation obtained as a limit in large populations of a probabilistic individual-based model describing the evolution of the adaptive trait [6]. In that sense, the structured population refers to the emergence and potential coexistence of multiple phenotypes that all belong to the submerged macrophyte type. To simplify the analysis and in the context of rapid evolution, we neglected the effect of mutations by seeding submerged macrophytes across the water column. In other words, we initially consider a normal distribution of the submerged macrophytes with additional very small phenotypic input uniformly distributed along the trait set. This means that the initial condition resembles a Gaussian distribution with non decaying tails. Thus, model dynamics were generated only by competition and selection by a partial differential equation for $S$ at time $t \ge 0$ with trait $z \in [0, z_b]$ which is a trait continuous version of the previous model :
$$ \partial_t S(t, z) = \left( r \frac{Nu(z)}{Nu(z) + h(1 + \int q(x)S(t, x)dx + q0F)} \frac{1}{1 + \int a_0S(t, x)dx + a(z, 0)F + W(z)} - L_s(z) \right) S \quad (2.8) $$

We assumed that the floating macrophytes still follow the same dynamics as in 2.1:

$$ \frac{dF(t)}{dt} = rF \frac{n(0, \bar{S}, F)}{n(0, \bar{S}, F) + h_f} \frac{1}{1 + a(0, 0)F + a(0, \bar{z})\bar{S} + W(0)} - L_f(0)F $$

where $\bar{S}(t) = \int_0^{z_b} S(t, x)dx$ is the total size of the population of the submerged macrophytes and $\bar{x} = \frac{\int_0^{z_b} xS(t, x)dx}{\int_0^{z_b} S(t, x)dx}$ is the mean trait of the population.

Initially the population S follows a Gaussian distribution centered around the dominant mean trait $z_0$ with sample seeding along the trait space $[0, z_b]$ introduced through $\varepsilon$:

$$ S(0, z) = S^* \Gamma_{\sigma^2}(z - z_0) + \varepsilon $$

where $S^*$ is the initial ecological steady state (resident population size), $\Gamma_{\sigma^2}(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp -\frac{x^2}{2\sigma^2}$ is the Gaussian Kernel and $z_0$ is the initial mean trait.

## 2.2. Co-evolution of Submerged and Floating macrophytes

We also investigated a co-evolutionary scenario involving both floating and submerged macrophytes. In this co-evolutionary scenario, there was a constraint imposed on the growth limits of the two plant types. Specifically, floating macrophytes were restricted from growing at deep levels, while submerged macrophytes were restricted from growing at shallow levels. We implemented this constraint by using step mortality functions 2.1.

Following AD, we made the assumption that the two populations of submerged and floating macrophytes, characterized by their initial traits $z_s$ and $z_f$ respectively, are at their ecological
equilibrium, denoted as $S^*$ and $F^*$. A mutant submerged plant, denoted as $S_m$, emerges with a trait value of $z_{s_m}$, and a mutant floating plant, denoted as $F_m$, appears with a trait value of $z_{f_m}$. The new system is described as follows :

$$ (S') \left\{ \begin{aligned} \frac{dS(t, z_s)}{dt} &= rS \frac{n(z_s, S, S_m, F, F_m)}{n(z_s, S, S_m, F, F_m) + h_s} \frac{1}{1 + a_0 S + a(z_s, z_{s_m}) S_m + a(z_s, z_f) F + a(z_s, z_{f_m}) F_m + W(z_s)} - L_s(z_s)S \\ \frac{dS_m(t, z_{s_m})}{dt} &= rS_m \frac{n(z_{s_m}, S, S_m, F, F_m)}{n(z_{s_m}, S, S_m, F, F_m) + h_s} \frac{1}{1 + a(z_{s_m}, z_s) S + a_0 S_m + a(z_{s_m}, z_f) F + a(z_{s_m}, z_{f_m}) F_m + W(z_{s_m})} - L_s(z_{s_m})S_m \\ \frac{dF(t, z_f)}{dt} &= rF \frac{n(z_f, S, S_m, F, F_m)}{n(z_f, S, S_m, F, F_m) + h_f} \frac{1}{1 + a_0 F + a(z_f, z_{f_m}) F_m + a(z_f, z_s) S + a(z_f, z_{s_m}) S_m + W(z_f)} - L_f(z_f)F \\ \frac{dF_m(t, z_{f_m})}{dt} &= rF_m \frac{n(z_{f_m}, S, S_m, F, F_m)}{n(z_{f_m}, S, S_m, F, F_m) + h_f} \frac{1}{1 + a_0 F_m + a(z_{f_m}, z_f) F a(z_{f_m}, z_s) S + a(z_{f_m}, z_{s_m}) S_m + W(z_{f_m})} - L_f(z_{f_m})F_m \end{aligned} \right. \quad (2.9) $$

We defined the respective fitness functions of the mutants $S_m$ and $F_m$ as :

$$ \left\{ \begin{aligned} s_s(z_{s_m}, z_s) &= r \frac{n(z_{s_m}, S^*, 0, F^*, 0)}{n(z_{s_m}, S^*, 0, F^*, 0) + h_s} \frac{1}{1 + a(z_{s_m}, z_s)S^* + a(z_{s_m}, z_f)F^* + W(z_{s_m})} - L_s(z_{s_m}) \\ s_f(z_{f_m}, z_f) &= r \frac{n(z_{f_m}, S^*, 0, F^*, 0)}{n(z_{f_m}, S^*, 0, F^*, 0) + h_f} \frac{1}{1 + a(z_{f_m}, z_f)F^* + a(z_{f_m}, z_s)S^* + W(z_{f_m})} - L_f(z_{f_m}) \end{aligned} \right. $$

Where the respective fitness Gradients are defined by :

$$ \left\{ \begin{aligned} \frac{\partial s_s}{\partial z_{s_m}}(z_s, z_s)|_{z_s=z_{s_m}} \\ \frac{\partial s_f}{\partial z_{f_m}}(z_f, z_f)|_{z_f=z_{f_m}} \end{aligned} \right. \quad (2.10) $$

To assess the final outcome of the co-evolutionary process, we used a structured population model to describe the dynamics of the submerged macrophytes and the floating macrophytes populations :

$$ \left\{ \begin{aligned} \partial_t S(t, z_s) &= \left( r \frac{Nu(z_s)}{Nu(z_s) + h(1 + \int q(x)(S(t,x) + F(t,x))dx)} \frac{1}{1 + \int a(z_s, x)S(t,x)dx + \int a(z_s, x)F(t,x)dx + W(z_s)} - L_s(z_s) \right) S(t, z_s) \\ \partial_t F(t, z_f) &= \left( r \frac{Nu(z_f)}{Nu(z_f) + h(1 + \int q(x)(S(t,x) + F(t,x))dx)} \frac{1}{1 + \int a(z_f, x)F(t,x)dx + \int a(z_f, x)S(t,x)dx + W(z_f)} - L_f(z_f) \right) F(t, z_f) \end{aligned} \right. $$
Similarly to the analysis of the submerged macrophytes' only evolution (section 2.1.2), we assumed that initially both the submerged macrophytes and the floating macrophytes populations were normally distributed around their respective resident mean traits. Additionally, we assumed a uniform seeding of these populations along the lake. This seeding process ensured that individuals with varying traits were present throughout the lake, allowing for potential evolutionary changes in both populations.

## 2.3. Numerical methods

In this section, we explain the numerical methods we followed, as analytical expressions for the eco-evolutionary dynamics of our shallow lake ecosystem were not possible to derive.

The analysis took place in two steps: first, we solved the ODE equations 2.1 in order to find the different equilibria of the system and determined the ecological steady states $(S^*, F^*)$ for each nutrient diffusion value $u \in [0, 10]$ and trait $z \in [0, z_b]$ using the python function **odeint** from the `scipy` package. This yields the two-dimensional bifurcation analysis where we identified regions of qualitatively different behaviors corresponding to different stable equilibria of the system for each $(u, z)$ combination. The second step used AD for finding the singular strategies (i.e., the roots of the gradient fitness 2.4) for different environmental conditions set by the parameter coefficients and the ecological equilibria $(S^*, F^*)$ using the Newton-Raphson method. In the presence of alternative stable states, the ecological equilibria used in the fitness equation 2.2 are the states with a positive population size of submerged macrophytes. We did not consider the alternative state of the floating macrophytes dominance since the submerged macrophytes are the evolving species and a resident population is needed. Subsequently, we determined the nature of the distinct singular strategies by assessing the stability and convergence conditions through the evaluation of second derivatives (see equations 2.5, 2.6, and 2.7). We completed the bifurcation analysis by incorporating the different CSS, ESS, and branching points for each nutrient diffusion value $u$.
For the structured population model simulations 2.8, we used the explicit Euler method with Neumann boundary condition that specifies the normal derivative at the boundary to be 0 with a numerical discretization of the finite trait space $[0.1, z_b = 10]$ using a step $\Delta z = 0.1$. This analysis allows us to gain insights into the dynamics and trends exhibited by the population over an extended period, providing a holistic perspective on the evolutionary processes and their ultimate consequences. We used arbitrary long simulation times of $N = 30000$, $N = 50000$ and $N=70000$ for different examples of nutrient diffusion values $u$ and the initial conditions $S(0, z) = S^*\Gamma_{\sigma^2}(z - z_0) + \varepsilon$ and $F(0) = F^*$ where $\varepsilon$ represents the seeding of the submerged macrophytes uniformly distributed across the entire lake. We assumed $\varepsilon \ll 1$. This assumption is made because we neglected mutations in our analysis, and our objective was to enable the growth of distinct traits in the submerged macrophyte population.

We used the total population size $\bar{S}(t)$, the emergence of different phenotypes, and the density distribution $S$ to describe the evolutionary outcome of the population of submerged macrophytes. For each of the 100 nutrient diffusion values $u \in [0.1, 10]$, we ran simulations for a long arbitrary number of iterations $N = 30000$ and considered the final population size and the trait distribution to determine the asymptotic behavior of the population $S$. If the total population size was below a certain threshold ($< \varepsilon$), then it was considered extinct. The presence of an evolutionary cycle was determined by analyzing the difference between the maximum and minimum values of the total mass of submerged macrophytes over the final time steps. If this difference remained relatively large or displayed recurring patterns of oscillation, it indicated the existence of cyclic behavior. Lastly, if the population size stabilized, then we assumed that the submerged macrophytes have evolved towards a monomorphic or polymorphic state according to the final trait distribution of the population.
3. Results

3.1. *Ecological dynamics*

<figure><img src="image_4.png" /><figcaption>FIGURE 3.1 – **Top** : Bifurcation diagram in a (u,z)-plane, where the trait of submerged macrophytes and the nutrient diffusion z, u ∈ [0, 10], with z = 0 being the surface and z_b = 10 the bottom of the lake. The different colors represent the different stable equilibria of the system 2.1. For low nutrient diffusion, either none of the two populations grows (I) or only the submerged macrophytes are able to invade (III). When u increases, a coexistence state appears (V), a region where the floating macrophytes outcompete the submerged (II), and two different bistability zones emerge. In the region (VI), the two alternative stable states are coexistence or an only **F** state and in the region (IV), the alternative stable states are a floating macrophytes dominance or a submerged macrophytes’ dominance. **Bottom** : Equilibrium biomass of the submerged and floating macrophytes (S,F) for specific values of nutrient diffusion u according to the growth depth of the submerged macrophytes z. (a) u = 3 : steady-states for low nutrient diffusion, where only the submerged macrophytes S can grow starting from a certain depth α. (b) u = 7 : we observe the different equilibria according to the trait z : when the submerged macrophytes are near the surface, only the floating macrophytes **F** grow while they get outcompeted when **S** are near the bottom. A coexistence stable state is possible when the submerged macrophytes are in intermediate levels, and a small bistable zone (VI) occurs between two tipping points z ∈ [δ, η] with only **F** and coexistence as alternative stable states. (c) u = 8 : Another bistability emerges (IV) of either submerged macrophytes dominance or a floating macrophytes’ dominance when z < γ.</figcaption></figure>
Our study focused on a specific scenario characterized by high nutrient concentration ($Nu_0$ = 50), low light attenuation ($w$=0.1), and a strong asymmetry in competition for light ($b$=10). This parameters choice was deliberate because it allowed us to identify a range of ecological equilibrium states and wide regions of bistability (different scenarios are illustrated in SI A. Fig. S1 highlighting less rich ecological equilibria.). These environmental parameter settings give rise to a bifurcation diagram in the (u,z)-plane, where the ecological dynamics of shallow lakes are strongly influenced by nutrient diffusion rate $u$ (Fig. 3.1). Figure 3.1 illustrates six distinct ecological steady states that differ qualitatively. When nutrient diffusion rate is low, either none of the two populations grows for the non-availability of the resources (I in Fig. 3.1) or only the submerged macrophytes invade (III in Fig. 3.1), outcompeting the floating macrophytes because of their priority effect on the nutrients. Increasing the level of nutrient diffusion rate ($u$ > 5) allows the floating macrophytes to grow (V in Fig. 3.1)and even outcompete the submerged macrophytes that are at higher depths (II in Fig. 3.1). A coexistence state also emerges for certain ranges of (u,z) (V in Fig. 3.1). Higher nutrient diffusion rate generates two additional regions with bistability. In these bistable regions (IV and VI in Fig. 3.1), different alternative stable states are present. In the first region (IV), the two alternative states possible are either only floating or only submerged macrophytes presence, whereas in the second bistable region (VI), a coexistence state is possible along with a floating macrophytes dominance alternative state. These alternative stable states result from different mechanisms; a trade-off between nutrients and light, where the parameter of nutrient diffusion rate $u$ determines the distribution of the nutrients in the water column and the decrease of luminosity in deeper levels is determined by light attenuation $w$ [49], a priority effect of the submerged macrophytes for nutrient access, determined by $p$ and asymmetrical competition for light [2].

Bifurcation plots of the submerged and floating macrophytes biomass for specific values of nutrient diffusion rate give more insights on the ecological equilibria of the system. In the
case of nutrient diffusion rate $u = 3$ (Fig. 3.1 b1, b2), the floating macrophytes never grow and the submerged macrophytes start growing at trait $z = \alpha$. When the nutrient diffusion rate $u = 7$ (Fig. 3.1 c1, c2), a coexistence state and a floating dominance alternative state are present (bistable region (VI)) where the submerged macrophytes switch abruptly from a coexistence state to their extinction when their depth trait $z$ approaches the lake surface and crosses the tipping point $z = \eta$. When the trait $z$ of the submerged macrophytes change in the opposite direction and increases from 0, the population $S$ shifts from extinction to coexistence at the point $z = \sigma$. When the nutrient diffusion rate increases to $u = 8$ (Fig. 3.1 d1,d2), a floating macrophytes dominance is the unique stable state when the trait $z$ is low, but when $z \in [\theta, \gamma]$, a coexistence alternative stable state emerges (bistable region VI). Above that threshold, when $z > \gamma$, the floating macrophytes go extinct and the remaining alternative states are floating macrophytes dominance or submerged macrophytes dominance (bistable region IV).

## 3.2. Evolutionary asymptotic analysis

<figure><img src="image_5.png" /><figcaption>FIGURE 3.2 – Bifurcation diagram with singular strategies obtained as roots of the gradient fitness $2.4$ from the AD framework for each nutrient diffusion rate $u$. The red line represents the final mean trait of the long time distribution of the submerged macrophytes population obtained from the structured population model simulations. The red hatches represent the range of the present phenotypes in the limit distribution of the population of submerged macrophytes $S$, which indicates the trait variation of the submerged macrophytes. The line and the hatches stop when the population $S$ doesn’t persist and goes extinct.</figcaption></figure>
### 3.2.1. Singular evolutionary strategies from Adaptive Dynamics

For the evolutionary analysis, starting only with the evolution of submerged macrophytes, we show how selection operates in the presence of floating macrophytes fixed to the surface when nutrient diffusion rate $u$ increases and the submerged macrophytes population $S$ is precluded from growing near the surface. Figure 3.2 shows the singular evolutionary strategies of the evolving submerged macrophytes on top of the underlying ecological states in the $(u,z)$-plane. When nutrient diffusion rate $u$ is less than 4.5, floating macrophytes are absent, resulting in a situation akin to the model proposed by [2] for the evolution of submerged macrophytes in the absence of floating macrophytes. The submerged macrophytes population evolves towards a CSS characterized by a monomorphic population when $u$ is very low $u < 1$. Alternatively, when $1 < u < 4.5$, the population reaches a branching point, leading to the emergence of a polymorphic population with multiple phenotypes (III in Fig. 3.2). When $u$ exceeds 4.5, different ecological scenarios are possible, and the evolutionary strategies appear to lie precisely on the bifurcation points between feasible regions of the submerged macrophytes (V, VI in Fig. 3.2) and the extinction state of $S$ (II in Fig. 3.2). In other words, the population of the submerged macrophytes will evolve towards the surface, crossing either the coexisting region (zone V) when $u < 7.5$ or the bistability region characterized by the alternative states of coexistence or only floating macrophytes (zone VI) when $u > 7.5$.
<figure><img src="image_6.png" /><figcaption>FIGURE 3.3 – **top :** (a) : Final total population size of submerged **S** and floating macrophytes **F** from the structured population model simulations as a function of nutrient diffusion rate **u**, where we observe different evolutionary outcomes :</figcaption></figure>

$u \in [0, U_d = 0.5[$ : Evolution of the submerged macrophytes **S** towards a monomorphic steady-state in the absence of the floating macrophytes. $u \in [U_d, U_{co} = 5.8[$ : Diversification of the submerged macrophytes **S** towards a polymorphic steady-state in the absence of the floating macrophytes. $u \in [U_{co}, U_c = 6.78[$ : Evolution of **S** towards a stable equilibrium of coexistence with **F**. $u \in [U_c, U_s = 7.45[$ : Beginning of oscillations where is represented the minimum and maximum of the population size during a cycle (amplitude of oscillation). $u \in [U_s, z_b = 10]$ : Beginning of evolutionary suicide, where **S** experiences evolutionary suicide.

**bottom :** Total population size of submerged $\int S(t, z)dz$ and floating macrophytes $F(t)$ through time using the structured population model (b1-f1). Emergence of different phenotypes (trait $z$) through time quantified in intervals of 0.5 of the trait as $\int_{i}^{i+0.5} S(t, z)dz$ using a threshold of detection $\epsilon$ for different nutrient diffusion rate values $u$ (b2-f2). This leads to different evolutionary outcomes : b1, b2) $u = 3$ : Evolution towards an equilibrium of diversification of the population **S** without the presence of the floating macrophytes. c1, c2) $u = 6$ : Evolution towards a steady-state of coexistence between the submerged and the floating macrophytes. $d1, d2$) $u = 7$ and $e1, e2$) $u = 7.4$ : Emergence of evolutionary oscillations of different amplitude and $f1, f2$) $u = 7.45$ : Evolution towards the extinction of the submerged macrophytes.
### 3.2.2. Evolutionary oscillations and evolutionary suicide from the structured population model

While AD allowed us to determine the expected phenotypic trait towards which the population of submerged macrophytes will evolve, simulations of the structured population model provided us with the long-term trait distribution of the submerged macrophytes population in the cases when the singular strategy was a BP. We plotted the final mean trait (red line in Fig. 3.1) and the range of all present traits (red hatched area) of the asymptotic distribution of the simulated submerged macrophytes phenotypes in Fig. 3.1 (a similar result is found if we consider the median trait instead of the mean). Interestingly, we found that the adaptive dynamics singular strategies and the evolutionary outcomes from the structured population model approaches are similar only to a certain extent. Specifically, the diversification of the submerged macrophytes population in the structured population model coincides with the predicted branching points from AD (Fig. 3.2). However, in the structured population model, the mean trait of the submerged macrophytes population does not match the trait values predicted by AD except for very low nutrient diffusion rates $u < 1$, where the CSS predictions of AD coincided with the single final dominant trait observed in structured population model. Nevertheless, the diversification initiated at lower nutrient diffusion rates compared to the initial branching points predicted by the AD framework. This discrepancy is due to two distinct factors : firstly, the AD approach foresees, for the majority of nutrient diffusion rates, the emergence of branching points, but it does not precisely specify the ultimate positions of the emerged traits [15]. Conversely, AD considers small, rare mutations that we did not incorporate into our structured population model. Instead, we have neglected the effects of mutations but allowed the simultaneous emergence of different traits within the population (see section 2.1.2). This assumption can be likened to the scenario of very frequent and substantial mutations (We refer to Fig. S4 in the SI. D for the final mean trait of the population with frequent mutations), or to the presence of a seed bank in the lake [27] [37], which allowed the rapid emergence of different phenotypes leading to polymorphisms
with potential dominant traits. Consequently, these emergent phenotypes influenced fitness through competition, potentially explaining the disparity with AD predictions.

But a closer look at the eco-evolutionary dynamics from the structured population model simulations revealed that the trajectories of the emerged phenotypes were not always stationary but oscillating (panel a in Fig. 3.3). In fact, we found that selection towards increased nutrient diffusion rate led to 4 distinct evolutionary behaviors before eventually reaching evolutionary suicide (Fig. 3.3). Firstly, there is the convergence towards a stable equilibrium characterized by either a monomorphic or polymorphic population of submerged macrophytes dominance. Secondly, the submerged macrophytes can evolve to a stable equilibrium of coexistence with the floating macrophytes. Thirdly, evolution can manifest as cyclic behavior with the presence of both the submerged and the floating macrophytes oscillating at bigger amplitude as *u* becomes higher. Finally, the submerged macrophyte population evolves to extinction (evolutionary suicide or Darwinian extinction [52]).

At low nutrient diffusion rate ($u < U_{co} = 5.8$) (Fig. 3.3 a), the submerged macrophytes evolve into a submerged macrophytes dominance state. When $u \in [U_{co}, U_c = 6.78[$, the evolution of **S** leads to their coexistence with the floating macrophytes **F**. In panels (Fig. 3.3 b1,b2) when the nutrient diffusion rate $u = 3$ and the initial population of submerged macrophytes has a mean trait $z_0 = 8$, the evolution leads to stable equilibrium in the absence of the floating macrophytes (Fig. 3.3 b1). The initial seeding allows submerged macrophytes with different traits to grow simultaneously and rapidly, resulting in a large trait variation in the population **S** which leads to a polymorphism with different traits present from 4 to 8 (Fig. 3.3 b2).

When the nutrient diffusion rate exceeds $U_{co}$, the floating macrophytes are able to grow and coexist with the submerged macrophytes. In panel (Fig. 3.3 c1,c2) when $u = 6$, the population size of **F** increases and the submerged macrophytes evolve towards a coexistence steady-state
(Fig.3.3 c1). Different phenotypes start growing towards the surface within the population **S**, until reaching the inflection point of the mortality function $L_S$ 2.1 $z_{thres} = 2$ level at which the population concentrates with a wide phenotypic variance (Fig. 3.3 c2).

At nutrient diffusion rate ($u \ge U_c$), evolutionary oscillations emerge with larger and larger amplitudes as $u$ increases. In panel (Fig. 3.3) d1,d2 when $u = 7$, starting from a resident population **S** with initial mean trait $z_0 = 8$ in the absence of the floating macrophytes, these latter grow rapidly and the two population sizes of **S** and **F** start oscillating (Fig. 3.3 d1). We observe the rapid emergence of new traits at higher depths up to the turning threshold feature $z_{thres} = 2$ defined by the step mortality function, in which the majority of the population settles down as it has the best access to light, while the rest of the population stays at deeper depths where it grows, increasing the total population size of **S**, and then disappears due to the intrapecific competition for light and shading effect (Fig. 3.3 d2).

At nutrient diffusion rate $u = 7.4$ in panel (Fig. 3.3 e1,e2), a similar scenario occurs, with evolutionary oscillations resulting from dynamics driven by competition and selection. At this $u$ rate, however, the population of floating macrophytes reaches high total biomass due to higher nutrient diffusion rate, resulting in a strong shadowing effect on the submerged macrophytes which almost die out before being rescued by the emergence of deeper phenotypes, resulting in larger fluctuations in one limit cycle. When the nutrient diffusion rate $u$ ramps up a little more to $U_s = 7.45$ (Fig. 3.3 f1,f2), it traverses a homoclinic bifurcation where the limit cycle grows such that the population of submerged macrophytes eventually transitions to the alternative stable state of extinction due to the bistability present. Then all trajectories belong to the attraction basin of the extinction state of the submerged macrophytes. In other words, selection leads to the evolutionary suicide of the submerged macrophytes population. It is worth noting that in the absence of the floating macrophytes, the submerged macrophytes persist under the same conditions. SI. C (Fig.S3) summarizes the distributions
of the submerged macrophytes population at different time steps for the examples of nutrient
diffusion rates studied above.

3.3. Co-evolution of submerged and floating macrophytes

In the previous section we presented the case where only the submerged macrophytes popu-
lation evolved. In this section, we explored the role of co-evolution by repeating the previous
analysis while allowing for the evolution of the floating macrophytes as well. Note that we
imposed a step mortality function for both macrophytes types with an arbitrary inflection
point at z_{thres} = 2 to prevent the submerged macrophytes from growing above the threshold
of z_{thres} and to prevent the floating macrophytes from growing below it. We did this because
in the absence of such constrain, we would observe each macrophytes type to evolve into its
competing type; a situation of reversal that ecologically would be meaningless.
<figure><img src="image_7.png" /></figure>

FIGURE 3.4 – Bifurcation diagrams for ecological equilibria of **F** and **S** with traits $z_f$ and $z_s$, respectively, when the nutrient diffusion rate $u = 7$ in (a) and $u = 7.6$ in (b). The yellow arrows indicate the direction of evolution of the trait $z_f$ of **F** and trait $z_s$ of **S** in each dotted region, according to the sign of their respective fitness gradients. We used the structured population model under the same set of parameter to simulate the co-evolution of **S** and **F** starting from initial phenotypes $z_{f0} = 0$ and $z_{s0} = 4$ (red dots). The results showed that when $u = 7$ in (a), the floating macrophytes population remained at the surface with a final trait of $z_f = 0$, while the submerged macrophytes evolved to greater depths, resulting in diversification with final mean trait of $z_s = 2.56$ (red star) in (a). However, in (b) when $u = 7.6$ the submerged macrophytes evolved to similar depths but didn't persist as this position fell within a bistable region where the extinction of **S** was an alternative stable state (red star).

We used another examples of initial positions $z_{f0} = 1$ and $z_{s0} = 4$ (blue dot) when $u = 7$ where the floating macrophytes evolved towards the surface, reaching a final trait $z_f = 0$ while the submerged macrophytes grew at greater depths, later moving deeper again with a final mean trait $z_s = 3.1$ where they coexist with **F** (blue star) in (a). In the case of $u = 7.6$ in (b), we studied the scenario with initial positions $z_{f0} = 1$ and $z_{s0} = 5$ where both populations took the same trajectory towards a final position of the floating macrophytes $z_f = 0$ which outcompeted the submerged macrophytes which evolved towards the surface and then went back to deeper depths in the bistability region (blue star) in (b).

In Figure 3.4, we first mapped the ecological equilibria in the ($z_f, z_s$)-plane for two nutrient diffusion rates $u = 7$ and $u = 7.6$ (panels a and b). In both cases, floating macrophytes cannot grow towards the bottom of the lake due to a high mortality rate (II). However, when both traits $z_f$ and $z_s$ are sufficiently low, the floating macrophytes can either outcompete the submerged macrophytes (I) or coexist with them (III), depending on the depth of the submerged macrophytes $z_s$. Relative bistability regions appear (IV) where the two alternative stable states are a coexistence or a floating macrophytes dominance state. Interestingly, the two populations never occupy the same depth. When $z_f$ is set to 0, we revert to the
scenario where the floating macrophytes were fixed at the surface, and a bistability zone
is observed. In the first case of $u = 7$ (Fig. 3.3 a1), this zone widens as the depth of the
floating macrophytes increases, until reaching approximately $z_f \sim 2$. In the second scenario
$u = 7.6$ (Fig. 3.3 b1), we observe a similar result with a bigger bistable region which inversely
becomes narrower as $z_f$ increases.

In these specific examples, we did not find any evolutionary equilibria or singular strategies
using AD. This implies that there is no combination of traits that results in both the floating
and submerged macrophytes having zero fitness gradients simultaneously (equations 2.10).
Instead, we only used the direction of selection based on the signs of the fitness gradients,
which indicate the trajectory towards which floating and submerged macrophytes are likely
to co-evolve (yellow arrows in Fig. 3.4). To find the actual co-evolutionary trajectory, we used
the structured population model. For example, by initially setting the traits of the submerged
and floating macrophytes to ($z_{s0} = 4, z_{f0} = 0$) (red dots in Fig. 3.4, we estimated the
evolutionary simulated trajectory (red dashed line) showing that the floating macrophytes
stayed at the surface while the submerged macrophytes diversified towards a mean trait of
$z_s = 2.56$ when $u = 7$ (red cross in a in Fig. 3.4) while in panels b1,b2 when $u = 7.6$,
the submerged macrophytes evolved towards a mean trait $z_s = 3.01$ which lies within the
bistable region (IV) where the population **S** went extinct. In another example when the
nutrient diffusion rate $u = 7$, we set the initial populations at initial mean traits $z_{s0} = 4$ and
$z_{f0} = 1$ (blue dot in panels a1,a2), the floating macrophytes evolved towards the surface and
the submerged macrophytes grew at higher depths until reaching a mean trait of $z_s \sim 2.5$
before evolving back to deeper levels with a final mean trait of $z_s = 3$ (blue star in a).
Similarly, in the case of $u = 7.6$, we studied the evolution of **S** and **F** initially set at $z_{s0} = 5$
and $z_{f0} = 1$ (blue dot in panels b in Fig. 3.4) where the two populations adopted a similar
trajectory but in this case leading to the suicide of the submerged macrophytes (blue star
in b).
<figure><img src="image_8.png" /><figcaption>FIGURE 3.5 – **Top**: (a) : Final total population size, when both **S** and **F** co-evolve, as a function of the nutrient diffusion rate *u*. The oscillations zone disappears and is replaced by a stable equilibrium of coexistence. Moreover, the floating plants start growing at a smaller nutrient diffusion rate *u* = *U*<sub>co</sub> = 5.7 evolutionary suicide of the submerged macrophytes is delayed and occurs for a higher *u* = *U*<sub>s</sub> = 7.6. **Bottom** : Total population size of the submerged **S** and floating macrophytes **F** through time and presence of different phenotypes in the submerged and floating macrophytes populations for different nutrient diffusion rates : *u* = 7 : b1) : The two population evolve towards a coexistence steady-state. b2) the submerged macrophytes population asymptotically concentrate around the inflection point of the mortality function *z*<sub>thres</sub> = 2 with the presence of other phenotypes bellow. b3) While we observe within the floating macrophytes population the emergence of low density population bellow the surface between 0 and 2 during the first time steps but which rapidly disappear leaving only the phenotypes at the surface *z*<sub>*f*</sub> = 0. *u* = 7.6 : c1) The floating macrophytes evolve towards a high density population outcompeting the submerged macrophytes that undergo evolutionary suicide. c2) The submerged macrophytes diversify with the emergence of different phenotypes at high levels before they get outcompeted by the floating macrophytes and go extinct. c3) The floating macrophytes stay at the surface after a few deeper phenotypes appeared and then vanished afterwards.</figcaption></figure>
Using the structured population model, we simulated the co-evolution of the submerged macrophytes and the floating macrophytes considering the final population sizes of **S** and **F** ($\int S(t_f, z)dz$ and $\int F(t_f, z)dz$ respectively) Fig. 3.5. On first sight, the co-evolutionary outcomes were similar to case when only submerged macrophytes evolved (panel a in Fig. 3.3 vs panel a in Fig. 3.5). This means that submerged macrophytes transitioned from monomorphic or polymorphic populations to a state of coexistence with floating macrophytes to extinction (panel a in Fig. 3.5). Nevertheless, the growth of the floating macrophytes occurred at a slightly smaller nutrient diffusion value compared to when only the submerged macrophytes evolve ($u = U_{c0} = 5.7$ Fig. 3.5 vs $u = U_{c0} = 5.8$ Fig. 3.3). In the co-evolution scenario when the nutrient diffusion rate was set at $uU_{c0} = 5.7$, the floating macrophytes evolved below the surface and diversified into a polymorphic population (Fig. S5 b3,b4) reaching a low but positive total population size (Fig. S5 b1), whereas they remained extinct when fixed to the surface and only the submerge macrophytes evolve (Fig. S5 a1). Strikingly, however, we no longer observed the oscillatory dynamics we observed in the case of only the submerged macrophytes evolution (panel a in Fig. 3.5). In addition, the extinction of submerged macrophytes is delayed and occurs for a higher nutrient diffusion rate as both plants co-evolve. To understand the lack of oscillations and the delayed evolutionary suicide of submerged macrophytes, we considered some examples in detail.

When the nutrient diffusion rate $u = 7$ explored in (b1,b2,b3 in Fig. 3.5), where evolutionary cycles were observed for when only the submerged macrophytes evolved (panels d2,d2 in Fig. 3.3), we extend the analysis to include the evolution of floating macrophytes. Starting from the resident traits ($z_{s0} = 8, z_{f0} = 0$), the co-evolutionary dynamics result in stabilized populations without oscillations (panel b1 in Fig. 3.3). The submerged plants evolve towards the surface, diversify and stabilize with a concentration around the trait threshold set by the inflection point of the mortality function $z_{thres} = 2$. We observe that individuals of the floating macrophytes population start growing at lower levels than the surface during
the initial time steps $\int_{0.5}^{2} F(t, z)dz > threshold$ (panel b2 in Fig. 3.3), as they have positive fitness as can be seen in SI. J Fig. S10). Nonetheless, they do not reach high densities before disappearing. Although the difference is not substantial compared to the scenario where only S evolved. The same observation is made for the example of nutrient diffusion rate $u = 6.8$ in SI. G Fig.S7 vs SI. H Fig.S8 where the oscillations disappeared in the co-evolution scenario. These observations suggests that the inclusion of the evolution of floating macrophytes may contribute to the non-replication of the same result without their evolution, which in our study led to evolutionary oscillations. In the co-evolutionary context, floating macrophytes are not confined to the surface alone and can influence the fitness of submerged macrophytes through their shading effect. A similar behavior unfolds within the floating macrophytes population when the nutrient diffusion rate is set to $u = 7.6$ (panels c1,c2,c3 in Fig. 3.3). However, in this case, the submerged macrophytes found themselves outcompeted by the floating macrophytes that attained a notably high biomass (panel c1 in Fig. 3.3). It's worth noting that the phenomenon of evolutionary suicide among the submerged macrophytes, when considering co-evolution with floating macrophytes, occurred at a higher nutrient diffusion rate. Furthermore, this occurred without the emergence of cyclic behavior, which contrasts with the situation when the submerged macrophytes evolved in isolation.

# 4. Discussion

In this study, we explored the effects of evolution and co-evolution in a shallow lake eco-system exhibiting ecological bistability. The model developed in this study goes beyond the ones proposed by [49] and [2] as it integrates the competition between floating and submerged macrophytes and their characterisation by a phenotypic trait in order to describe the eco-evolutionary dynamics of the ecosystem. Bistability arises due to the positive feedback resulting from competition between floating and submerged macrophytes for nutrients and light [43, 49, 4]. Consistent with previous results, we find that ecological alternative stable states occur for high rates of nutrient diffusion, causing shifts from coexistence or exclusive presence of submerged macrophytes to a floating macrophytes dominance. We find that
the evolution of a phenotypic trait representing the growth depth of macrophytes in our shallow lake ecosystem leads to evolutionary oscillations and evolutionary suicide. Instead, co-evolution of both macrophyte types stabilises the dynamics but still leads to evolutionary suicide of the submerged macrophytes. Below we discuss our results in detail.

4.1. *Evolution of the submerged macrophytes can lead to evolutionary branching, oscillations and evolutionary suicide*

Our analysis shows that submerged macrophytes' growth depth mostly evolves towards the water surface for varying levels of nutrient diffusion rate : basically evolving towards a higher availability of light as light is the limiting resource for the submerged macrophytes. For a range of very low nutrient diffusion rates (0<u<1.5 Fig. 3.2 red line) when the submerged macrophytes evolve towards the bottom of the lake, we observe an initial decrease in the total population size of the submerged macrophytes (Fig. 3.3). This arises because, when the nutrient diffusion rate is very low (near zero), nutrient concentration accumulates at the bottom of the lake, providing an advantage to the deep submerged macrophytes. As the diffusion rate slightly increases (e.g. from 0.1 to 0.2), nutrients disperse to shallower depths, reducing their availability at the bottom. Within this range of nutrient diffusion rates, as the submerged macrophytes evolve toward the bottom, higher nutrient diffusion leads to reduced nutrient access in deep levels, explaining the decrease in the total biomass of the submerged macrophytes population (SI. B Fig. S2).

For a range of intermediate nutrient diffusion rates when the floating macrophytes are ab-
sent (1<u<4.5, Fig. 3.2), submerged macrophytes evolve towards branching points, leading
to diversification of the monomorphic population into a polymorphic population as found in
[2]. Diversification can be attributed to the advantageous priority effect held by submerged
macrophytes in terms of nutrient access. As nutrient diffusion rate increases,the depth of the
singular strategy decreases (goes towards the surface) since higher diffusion allows shallower
depths to experience a significant enrichment in nutrients. However, the evolutionary out-
come is not a stable strategy anymore but a branching point. As nutrient availability widens at different depths of the lake, it enables the emergence and almost stable coexistence of different phenotypes within a specific depth range (Fig. 3.3).

Interestingly, when the nutrient diffusion rate exceeds a certain value ($u>4.5$, Fig. 3.1), the evolutionary equilibria of the submerged macrophytes based on AD (BP in Fig. 3.2) coincide with the critical depth traits that represent bifurcation points between bistability (VI Fig 3.2) and floating macrophytes dominance (II in Fig. 3.2). Simulations of the structured population model, though, show that submerged macrophytes evolve to depths where they can either dominate or coexist with the floating macrophytes either in stable coexistence or in oscillations. Especially, for diffusion rates $4.5<u<7.5$ Fig. 3.2, despite the existence of a wide range of phenotypes in which submerged macrophytes could potentially dominate the lake on their own, they still evolve towards depths that enable their competing floating macrophytes to invade and coexist with them.

But what does drive the oscillating evolutionary patterns when nutrient diffusion rate exceeds a specific threshold $u \ge U_c$ (Fig. 3.3)? At this threshold of diffusion rate, limitation due to nutrient competition diminishes. Consequently, competition for light becomes the primary driver influencing population dynamics. As we assumed continuous seeding throughout the water column, the emergence of different phenotypes of submerged macrophytes at various depths is possible. Submerged macrophytes evolve to the surface but they cannot be maintained there in high abundance because of the shading effect from the floating macrophytes. Competition for light reduces the population size of submerged macrophytes, and simultaneously increases the population of floating macrophytes as they get more access to nutrients. At the same time though, the reduced size of submerged macrophytes also reduces intraspecific shading of submerged macrophytes and enables the deeper submerged phenotypes to grow again. As a result phenotypes at different depths emerge and then disappear
triggering the oscillating patterns between submerged and floating macrophytes. Similar cycles between high and low abundance has been observed in Antarctic krill populations ([41]). In that case, intraspecific food competition limits the growth of krill within the same species when krill populations reach a critical abundance, but when krill abundance is low, growth relies more on environmental conditions [41].

When the nutrient diffusion rate increases the amplitude of the oscillations becomes larger until a threshold ($u \ge U_s$, Fig. 3.3), where evolutionary oscillations become so large that the submerged macrophytes cross to the extinction state (Fig. 3.3 f1,f2). More in general, at high diffusion rates, the evolutionary trajectory tends towards depth regions where submerged macrophytes face extinction rather than to regions where they can be rescued (or more properly persist). This is due to the fact that at such high range of diffusion rate, nutrient limitation is not anymore present. Instead, the only way to increase fitness is by growing to shallower depths where light is more available, but where the only feasible equilibrium is the one of a dominant floating macrophytes population and an extinct submerged ones. This is a classical case of evolutionary suicide. Evolution to extinction was put forward by [1] in which a prey population went extinct when decreasing foraging activity. Evolutionary suicide through the existence of a catastrophic bifurcation was highlighted by [21] in structured metapopulation models with evolving dispersal. Other examples showed that evolutionary suicide could also occur through non-catastrophic bifurcations, as depicted in a model of pathogen evolution with frequency-dependent transmission in [5]. In our simulations of the structured population model, the submerged macrophytes experience a catastrophic collapse as evolutionary oscillations become so large that they cross a homoclinic saddle bifurcation leading to the extinction of the submerged macrophytes. Similar dynamics have been shown in a predator prey model where a sudden ecological shift from a coexistence limit cycle to predator extinction can occur through catastrophic homoclinic bifurcations [33].
4.2. Co-evolution of the floating macrophytes prevents evolutionary cycles and delays evolu-
tionary suicide of submerged macrophytes

Strikingly, co-evolution does not affect the overall qualitative behaviour of the system com-
pared to the evolution of the submerged macrophytes only. One might have expected that,
under the influence of co-evolution, the floating macrophytes would start growing at much
lower nutrient diffusion values (u << Uco = 5.8 Fig. 3.3), allowing them to establish them-
selves at deeper depths where the floating macrophytes can grow at levels below the surface.
However, the floating macrophytes fail to grow at these low nutrient diffusion rates, because
under low nutrient diffusion, floating macrophytes are able to grow at deeper depths only if
submerged macrophytes would occupy the surface. This reversal is unachievable given that
the submerged macrophytes are restricted below a certain threshold (i.e. they reside in the
deeper depths of the lake, SI. I Fig. S9 b,c,d,e,f). Despite similar evolutionary outcomes,
co-evolution triggers floating macrophytes growth at a slightly lower nutrient diffusion value
where they diversify at deeper depths (SI. E Fig. S5), compared to scenarios where they are
fixed at the surface. The floating macrophytes continue evolving towards phenotypes below
the surface for a limited range of nutrient diffusion rates leading to higher total biomass
compared to scenarios where floating macrophytes are fixed at the surface (SI F. Fig. S6 a1).
However, at sufficiently high nutrient diffusion rates deep traits vanish after a brief period,
leaving only surface traits (SI I. Fig. S9). This is due to the fact that nutrients are no longer
a limiting resource and staying at the surface becomes the optimal strategy for the floating
macrophytes.

Yet, co-evolution of the floating and submerged macrophytes has two more subtle important
effects. First, co-evolution prevents the occurrence of evolutionary oscillations, but leads to
the establishment of a stable coexistence equilibrium. Comparing dynamics in the nutrient
diffusion ranges where evolutionary oscillations disappear (Fig. 3.3 d vs (Fig. 3.4 b), in the
case of co-evolution the floating macrophytes population size is smaller than when floating
macrophytes do not evolve. Co-evolution allows the emergence of deeper phenotypes within the floating macrophytes population leading to nutrient competition on the floating macrophytes at the surface. This competition decreases the overall biomass of floating macrophytes and, thus, diminishes their shading impact on submerged macrophytes. In its turn, this results in a higher biomass of submerged macrophytes at shallow depth ($z \sim 2$), which reduces the abundance of deeper phenotypes due to increased shading. As a result the recurrent emergence and decline of submerged macrophytes are no longer possible when the floating macrophytes co-evolve.

Interestingly, this stabilization persists even though there is just a brief emergence of deeper phenotypes within the floating macrophytes populations and asymptotically the floating macrophytes are fixed at the surface. One might question why, in this case, evolutionary oscillations do not occur at a higher nutrient diffusion rate where the floating macrophytes biomass is high enough to trigger oscillatory dynamics as described in the no co-evolution case? This is because deeper phenotypes within the floating macrophytes population appear with lower density and for shorter time as the nutrient diffusion rate increases (SI. H Fig. S8 c and Fig. 3.5 b3, c3). By the time their impact fades away, the nutrient diffusion rate becomes that high that the evolutionary suicide of the submerged macrophytes occurs. This can be seen in the non-persistent oscillations just before the extinction of submerged macrophytes outcompeted by the floating macrophytes (Fig. 3.5 c1).

The reduction in floating macrophytes' biomass driven by the emergence of deeper phenotypes is also responsible for the second difference observed in the co-evolution case: the evolutionary suicide of the submerged macrophytes is delayed to a higher nutrient diffusion rate ($u = U_s = 7.45$ Fig. 3.3 vs $u = U_s=7.6$ Fig. 3.5). The emergence of deeper floating macrophytes phenotypes due to co-evolution suppresses oscillatory dynamics, leads to the growth of the submerged macrophytes to the shallowest possible depth, and allows the sub-
merged macrophytes to persist longer even if only for a narrow additional range of nutrient diffusion rates.

### 4.3. Conclusion

Evolution of the growth depth of submerged macrophytes in a bistable shallow lake ecosystem destabilised the stable coexistence between submerged and floating macrophytes. We found an emergence of evolutionary oscillations that eventually lead to the evolutionary suicide of the submerged macrophytes as nutrient diffusion increased and nutrient limitation became weak. Co-evolution of both floating and submerged macrophytes was stabilising by preventing the occurrence of evolutionary oscillations, and delayed, but it did not prevent, the evolutionary suicide of the submerged macrophytes. The fact that evolution ultimately resulted in the submerged macrophytes extinction despite the availability of growth depths where submerged and floating macrophytes could coexist suggests that adaptive evolution does not necessarily prevent collapse in bistable ecosystems.

## 5. Data Availability

The codes to reproduce the figures of this article are available at https://github.com/SirineBoucenna/Evolutionary-model-shallow-lakes

### Acknowledgments

We thank Nicolas Loeuille and Sepideh Mirrahimi for valuable comments and suggestions. This research was supported in part by the International Centre for Theoretical Sciences (ICTS) for the program "Tipping Points in Complex Systems" (code: ICTS/tipc2022/9).

### Références

[1] Peter A ABRAMS et Hiroyuki MATSUDA. "The evolution of traits that determine ability in competitive contests". In: *Evolutionary Ecology* 8 (1994), p. 667-686.
[2] Alice Nadia ARDICHVILI, Nicolas LOEUILLE et Vasilis DAKOS. "Evolutionary emergence of alternative stable states in shallow lakes". In : *Ecology Letters* 26.5 (2023), p. 692-705.

[3] Peter ASHWIN et al. "Tipping points in open systems : bifurcation, noise-induced and rate-dependent examples in the climate system". In : *Philosophical Transactions of the Royal Society A : Mathematical, Physical and Engineering Sciences* 370.1962 (2012), p. 1166-1184.

[4] Irmgard BLINDOW, Anders HARGEBY et Gunnar ANDERSSON. "Alternative stable states in shallow lakes : what causes a shift ?" In : *The structuring role of submerged macrophytes in lakes*. Springer, 1998, p. 353-360.

[5] Barbara BOLDIN et Eva KISDI. "Evolutionary suicide through a non-catastrophic bifurcation : adaptive dynamics of pathogens with frequency-dependent transmission". In : *Journal of Mathematical Biology* 72 (2016), p. 1101-1124.

[6] Nicolas CHAMPAGNAT, Régis FERRIÈRE et Sylvie MÉLÉARD. "Unifying evolutionary dynamics : from individual stochastic processes to macroscopic models". In : *Theoretical Population Biology* 69.3 (2006), p. 297-321.

[7] P Catalina CHAPARRO PEDRAZA et al. "Adaptive evolution can both prevent eco-system collapse and delay ecosystem recovery". In : *The American Naturalist* 198.6 (2021), E185-E197.

[8] P Catalina CHAPARRO-PEDRAZA. "Fast environmental change and eco-evolutionary feedbacks can drive regime shifts in ecosystems before tipping points are crossed". In : *Proceedings of the Royal Society B* 288.1955 (2021), p. 20211192.

[9] P Catalina CHAPARRO-PEDRAZA et André M de ROOS. "Ecological changes with minor effect initiate evolution to delayed regime shifts". In : *Nature Ecology & Evolution* 4.3 (2020), p. 412-418.
[10] Olivier COTTO et Ophélie RONCE. "Maladaptation as a source of senesence in habitats variables in space and time." In : *Evolution* 68.9 (sept. 2014), p. 2481-2493.

[11] Olivier COTTO et al. "Maladaptive shifts in life history in a changing environment". In : *The American Naturalist* 194.4 (2019), p. 558-573.

[12] Vasilis DAKOS et al. "Ecosystem tipping points in an evolving world". In : *Nature ecology & evolution* 3.3 (2019), p. 355-362.

[13] Fabio DERCOLE, Régis FERRIÈRE et Sergio RINALDI. "Ecological bistability and evolutionary reversals under asymmetrical competition". In : *Evolution* 56.6 (2002), p. 1081-1090.

[14] Ulf DIECKMANN et Richard LAW. "The dynamical theory of coevolution : a derivation from stochastic ecological processes". In : *Journal of Mathematical Biology* 34 (1996), p. 579-612.

[15] Odo DIEKMANN et al. "The dynamics of adaptation : an illuminating example and a Hamilton-Jacobi approach". In : *Theoretical Population Biology* 67.4 (2005), p. 257-271.

[16] Alexandra ERBACH, Frithjof LUTSCHER et Gunog SEO. "Bistability and limit cycles in generalist predator-prey dynamics". In : *Ecological Complexity* 14 (2013), p. 48-55.

[17] Bjørn A FAAFENG et Marit MJELDE. "Clear and turbid water in shallow Norwegian lakes related to submerged vegetation". In : *The structuring role of submerged macrophytes in lakes*. Springer, 1998, p. 361-368.

[18] Carl FOLKE et al. "Regime shifts, resilience, and biodiversity in ecosystem management". In : *Annu. Rev. Ecol. Evol. Syst.* 35 (2004), p. 557-581.

[19] Jimmy GARNIER et al. "Adaptation of a quantitative trait to a changing environment : New analytical insights on the asexual and infinitesimal sexual models". In : *Theoretical Population Biology* 152 (2023), p. 1-22.
[20] Stefan AH GERITZ et al. "Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree". In: *Evolutionary Ecology* 12 (1998), p. 35-57.

[21] Mats GYLLENBERG et Kalle PARVINEN. "Necessary and sufficient conditions for evolutionary suicide". In: *Bulletin of Mathematical Biology* 63 (2001), p. 981-993.

[22] Nelson G HAIRSTON JR et al. "Rapid evolution and the convergence of ecological and evolutionary time". In: *Ecology Letters* 8.10 (2005), p. 1114-1127.

[23] Crawford S HOLLING. "Resilience and stability of ecological systems". In: *Annual review of ecology and systematics* 4.1 (1973), p. 1-23.

[24] Terry P HUGHES et al. "Spatial and temporal patterns of mass bleaching of corals in the Anthropocene". In: *Science* 359.6371 (2018), p. 80-83.

[25] G Evelyn HUTCHINSON. "Chapter 32 : The algal benthos". In: *A treatise in limnology III Limnological Botany* (1975), p. 509-571.

[26] Jan H JANSE. "A model of nutrient dynamics in shallow lakes in relation to multiple stable states". In: *Shallow Lakes' 95 : Trophic Cascades in Shallow Freshwater and Brackish Lakes* (1997), p. 1-8.

[27] K KAPLAN et Th MUER. "Beobachtungen zum diasporenreservoir im Bereich ehemaliger Heideweiher". In: *Floristische Rundbriefe* 24 (1990), p. 38-45.

[28] Sonia KÉFI et al. "Evolution of local facilitation in arid ecosystems". In: *The American Naturalist* 172.1 (2008), E1-E17.

[29] Sonia KÉFI et al. "Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems". In: *Nature* 449.7159 (2007), p. 213-217.

[30] Motoo KIMURA. "A stochastic model concerning the maintenance of genetic variability in quantitative characters." In: *Proceedings of the National Academy of Sciences* 54.3 (1965), p. 731-736.
[31] Christopher A KLAUSMEIER et al. "Ecological limits to evolutionary rescue". In: *Philosophical Transactions of the Royal Society B* 375.1814 (2020), p. 20190453.

[32] Jan J KUIPER et al. "Food-web stability signals critical transitions in temperate shallow lakes". In: *Nature Communications* 6.1 (2015), p. 7727.

[33] Sami O LEHTINEN. "Ecological and evolutionary consequences of predator-prey role reversal: Allee effect and catastrophic predator extinction". In: *Journal of Theoretical Biology* 510 (2021), p. 110542.

[34] Timothy M LENTON et al. "Tipping elements in the Earth's climate system". In: *Proceedings of the national Academy of Sciences* 105 (2008).

[35] Nicolas LOEUILLE. "Influence of evolution on the stability of ecological communities". In: *Ecology Letters* 13.12 (2010), p. 1536-1545.

[36] Robert M MAY. "Thresholds and breakpoints in ecosystems with a multiplicity of stable states". In: *Nature* 269.5628 (1977), p. 471-477.

[37] DG MCFARLAND et SJ ROGERS. "The aquatic macrophyte seed bank in Lake Onalaska, Wisconsin". In: *Journal of Aquatic Plant Management* 36.JAN. (1998), p. 33-39.

[38] Johan AJ METZ, Roger M NISBET et Stefan AH GERITZ. "How should we define 'fitness' for general ecological scenarios?" In: *Trends in Ecology & Evolution* 7.6 (1992), p. 198-202.

[39] Matthew M OSMOND et Christopher A KLAUSMEIER. "An evolutionary tipping point in a changing environment". In: *Evolution* 71.12 (2017), p. 2930-2941.

[40] Lucia RUSSO et Konstantinos G SPILIOTIS. "Bifurcation analysis of a forest-grassland ecosystem". In: 1738.1 (2016).

[41] Alexey B RYABOV et al. "Competition-induced starvation drives large-scale population cycles in Antarctic krill". In: *Nature Ecology & Evolution* 1 (2017).
[42] M SCHEFFER, S RINALDI et L MUR. "On the dominance of filamentous blue-green
algae in shallow lakes". In : (1994).

[43] Marten SCHEFFER et al. *Ecology of shallow lakes*. T. 1. Springer, 1998.

[44] Marten SCHEFFER et Stephen R CARPENTER. "Catastrophic regime shifts in ecosystems : linking theory to observation". In : *Trends in Ecology & Evolution* 18.12 (2003), p. 648-656.

[45] Marten SCHEFFER et Erik JEPPESEN. "Regime shifts in shallow lakes". In : *Ecosystems* 10.1 (2007), p. 1-3.

[46] Marten SCHEFFER, Sergio RINALDI et Yuri A KUZNETSOV. "Effects of fish on plankton dynamics : a theoretical analysis". In : *Canadian Journal of Fisheries and Aquatic Sciences* 57.6 (2000), p. 1208-1219.

[47] Marten SCHEFFER et al. "Alternative equilibria in shallow lakes". In : *Trends in Ecology & Evolution* 8.8 (1993), p. 275-279.

[48] Marten SCHEFFER et al. "Catastrophic shifts in ecosystems". In : *Nature* 413.6856 (2001), p. 591-596.

[49] Marten SCHEFFER et al. "Floating plant dominance as a stable state". In : *Proceedings of the National Academy of Sciences* 100.7 (2003), p. 4040-4045.

[50] Marten SCHEFFER et al. "On the dominance of filamentous cyanobacteria in shallow, turbid lakes". In : *Ecology* 78.1 (1997), p. 272-282.

[51] Egbert H VAN NES et al. "What do you mean, 'tipping point' ?" In : *Trends in Ecology & Evolution* 31.12 (2016), p. 902-904.

[52] Colleen WEBB. "A complete classification of Darwinian extinction in ecological interactions". In : *The American Naturalist* 161.2 (2003), p. 181-205.
[53] Tianran ZHANG et Wendi WANG. "Hopf bifurcation and bistability of a nutrient-phytoplankton-zooplankton model". In: *Applied Mathematical Modelling* 36.12 (2012), p. 6225-6235.

**6. Supplementary information**

The following subsections gather supplementary figures of the main text.

**6.1. SI A. Bifurcation diagrams under different parameter conditions**

<figure><img src="image_9.png" /><figcaption>FIG. S1 – Bifurcation diagrams for different parameter conditions than used in our study. **a**: Higher light attenuation *w* = 0.2. **b**: Smaller nutrient concentration $Nu_0$ = 10. **c**: Smaller nutrient concentration $Nu_0$ = 10 and higher light attenuation *w* = 0.3.</figcaption></figure>
## 6.2. SI B. Steady states of the submerged macrophytes for very low nutrient diffusion rates

<figure><img src="image_10.png" /><figcaption>FIG. S2 – Ecological equilibria representing the biomass of the submerged macrophytes population for three small nutrient diffusion rates. We can see that the maximum of the population size decreases as *u* increases because the higher the diffusion the less the nutrients concentration in the bottom.</figcaption></figure>

## 6.3. SI C. Distribution of the submerged macrophytes population from the structured population model

<figure><img src="image_11.png" /><figcaption>FIG. S3 – Density distributions of the submerged macrophytes populations along the trait set of the examples studied in Fig. 3.3 at different time steps. In **a**), we observe the population evolving from a monomorphic population at the depth *z* = 8 to a polymorohism with several trait present. In **b**), **c**) and **d**) where we observed evolutionary oscillations, we illustrate the distributions of **S** at the exterma in a cycle. In **d**) the population experienced extinction afterwards.</figcaption></figure>
6.4. SI D. *Comparison of the final mean traits from the structured population model with and without mutations*

<figure><img src="image_12.png" /><figcaption>FIG. S4 – Final mean traits obtained from the structured population model in the case of neglecting the mutations with initial seeding (from our study in Fig. 3.2) and in case of frequent Gaussian mutations with a probability of mutation 0.1 and a variance 0.05 where in this scenario we did not observe oscillations and the final distributions of the submerged macrophytes population were a unique Gaussian with a certain phenotypic variance.</figcaption></figure>
6.5. SI E. An example of nutrient diffusion rate $u = 5.7 : Only S evolution vs co-evolution

<figure><img src="image_13.png" /><figcaption>FIG. S5 – Comparison of the scenario where only the submerged macrophytes evolve **left** and where the floating and submerged macrophytes co-evolve **right** when the nutrient diffusion rate $u = 5.7$. We observe in panel (a1) that floating macrophytes remain at 0 and don’t grow, while in the case of co-evolution the total population size of the floating macrophytes increases in panel (b1) where different phenotypes of the floating macrophytes population emerge $0 < z_f < 2$ in panel (b3) and then the ones closer to the surface disappear. Panels (a4) and (b4) represent the initial and final distributions of the submerged macrophytes population and in case of co-evolution of the floating macrophytes population as well.</figcaption></figure>
6.6. SI F. An example of nutrient diffusion rate $u = 5.9 : Only S evolution vs co-evolution

<figure><img src="image_14.png" /><figcaption>FIG. S6 – Comparison of the scenario where only the submerged macrophytes evolve **left** and where the floating and submerged macrophytes co-evolve **right** when the nutrient diffusion rate $u = 5.9$. We observe in panel (a1) that floating macrophytes that are fixed at the surface grow and reach low population size, while in the case of co-evolution the total population size of the floating macrophytes increases and reaches higher values in panel (b1). Different phenotypes of the floating macrophytes population emerge $0.5 < z_f < 2.5$ in panel (b3) whereas they are fixed to the surface when only the submerged macrophytes evolve in panel (a3). Panels (a4) and (b4) represent the initial and final distributions of the submerged macrophytes population and in case of co-evolution of the floating macrophytes population as well.</figcaption></figure>
6.7. SI G. An example of nutrient diffusion rate $u = 6.8$ at the beginning of the cycle region

<figure><img src="image_15.png" /><figcaption>FIG. S7 – (a) : Total population size of the submerged macrophytes and the floating macrophytes when only S evolve for $u = 6.8$ and (b) : the presence of different traits z of the submerged macrophytes, which start from initial depth $z_0 = 8$, with the emergence of small amplitude evolutionary cycles.</figcaption></figure>
6.8. SI H. The example of nutrient diffusion rate $u = 6.8$ in case of co-evolution

<figure><img src="image_16.png" /><figcaption>FIG. S8 – (a) : Total population size of the submerged macrophytes and the floating macrophytes when **S** and **F** co-evolve for $u = 6.8$. Emergence of different phenotypes of the submerged macrophytes population $z_s$ in (b) and of the floating macrophytes population $z_f$ in (c) through time.</figcaption></figure>
6.9. SI I. Bifurcation diagrams in (z<sub>s</sub>, z<sub>f</sub>)-plane for different nutrient diffusion rates

<figure><img src="image_17.png" /><figcaption>FIG. S9 – Bifurcation diagrams in (z<sub>f</sub>, z<sub>s</sub>)-plane with ecological steady-states of the floating and submerged macrophytes when their respective growth depths z<sub>f</sub> ∈ [0.5] and z<sub>s</sub> ∈ [0.10] for different nutrient diffusion rates : (a) : u=1 (b) : u=2 (c) : u=2.5 (d) : u=2.7 (e) : u=3 (f) : u=4 (g) : u=5 (h) : u=6. (i) : u=8 (j) : u=9.</figcaption></figure>
6.10. SI J. Distribution and fitness of the floating macrophytes population in case of evolution when $u = 7$

<figure><img src="image_18.png" /><figcaption>FIG. S10 – The example of nutrient diffusion rate $u = 7$ studied in Fig. 3.4 when **S** and **F** evolve where were observed the emergence of phenotypes deeper than 0 in the floating macrophytes population which prevented the occurrence of cyclic behaviour. **Left**: Distribution of the floating macrophytes density at different time steps; initially we observe the presence of seeding uniformly distributed along the trait set which disappear at time $t = 600$ except at depth $z_f \sim 1.5$. At a longer time step we can see that only the population at the surface remained. **Right**: fitness of the population **F** at the same time steps for every trait $z_f \in [0, 10]$ which represents their growth rate.</figcaption></figure>