Structural stability estimated through critical perturbation determines evolutionary persistence in mutualistic model ecosystems

TL;DR

This study shows how evolutionary adaptation in mutualistic ecosystems affects structural stability and biodiversity under environmental change.

Contribution

The paper introduces evolutionary adaptation into mutualistic models to link ecological stability with evolutionary persistence.

Findings

Evolution increases critical competition, enhancing structural stability in mutualistic networks.

Natural selection steers ecosystems toward reduced effective competition, fostering stability.

Structural stability predicts biodiversity loss under evolutionary and environmental changes.

Abstract

Understanding the factors that influence the persistence and stability of complex ecological networks is a central focus of ecological research. Recent research into these factors has predominantly attempted to unveil the ecological processes and structural constraints that influence network stability. Comparatively little attention has been given to the consequences of evolutionary events, despite the fact that the interplay between ecology and evolution has been recognized as fundamental to understand the formation of ecological communities and predict their reaction to change. We extend existing mutualistic population dynamical models by incorporating evolutionary adaptation events to address this critical gap. We relate ecological aspects of mutualistic community stability to the stability of persistent evolutionary pathways. Our findings highlight the significance of the structural stability of ecological systems in predicting biodiversity loss under both evolutionary and environmental changes, particularly in relation to species-level selection. Notably, our simulations reveal that the evolution of mutualistic networks tends to increase a network-dependent parameter termed critical competition, which places systems in a regime in which mutualistic interactions enhance structural stability and, consequently, biodiversity. This research emphasizes the pivotal role of natural selection in shaping ecological networks, steering them towards reduced effective competition below a critical threshold where mutualistic interactions foster stability in the face of environmental change.

Figures

Eco-evolutionary model and critical competition ρc. (a) Schematic representation of the model. At each step of the simulation, a randomly selected species (green square) is subjected to an evolutionary event (box ‘evolutionary dynamics’). Evolutionary events are of two types: link swap (E1) or creation/loss (E2). The probabilities of events (p, q respectively) are indicated. The effect of the event on the network is evaluated by running the ecological dynamics, where environmental perturbations may be added. Acceptance of the event depends on whether selection is considered. If selection is considered, the change is accepted only if it enhances the fitness of the modified species (i.e. there is an increase in its biomass). If there is no selection all changes are accepted, regardless of their effect on species or the community. (b) Simulated time series of species abundances illustrating the trade-off between the interspecific competition parameter ρ and the critical competition ρc. (c) Example of eco-evolutionary dynamics showing six trajectories for no selection (blue) and selection (orange) scenarios. The evolutionary equilibrium is achieved when no extinctions are observed throughout 350 substitutions (i.e. accepted evolutionary events). (d) log–log representation of the critical competition against the total biodiversity for a set of real plant–pollinator networks extracted from the Web of Life. The power–law relationships for both plants and animals have an exponent approximately −0.5. By means of the trade-off illustrated in (b), observing positive mutualistic effects in large communities demands decreasing the intraguild interspecific competition to favour coexistence.

Structural stability explains diversity at the evolutionary equilibrium. (a) Biodiversity of the eco-evolutionary process at equilibrium for different values of the interspecific competition ρ and environmental perturbation Δ. Each panel represents the outcome of the evolutionary process for different values of Δ. In boxplots, box limits represent interquartile range (IQR) and the horizontal line the median value of the final species richness (α-diversity) across 50 simulations. Vertical lines (i.e. whiskers) show the ±1.5 IQR values. Points represent the outliers. (b) Predicted structural stability (critical environmental perturbation Δc) for the starting network for each value of interspecific competition ρ. (c) Relationship between structural stability and final biodiversity, normalized by the starting diversity before the evolutionary process. (d) Slope of the relationship shown in (c) for different values of environmental perturbation. The closer the slope is to one the more accurate the prediction. Ordinates are not significantly different from zero for Δ=0.1.

Disruptive evolution prompts biodiversity loss and biomass increase. Mean α-diversity (i.e. species richness) (left), effective diversity (middle) and biomass (right) throughout eco-evolutionary trajectories (halted after stability in number of species is achieved over 350 evolutionary steps) for different values of environmental perturbation Δ (rows) and interspecific competition ρ (columns). Mean values (points) and standard deviation (vertical lines) are shown every 25 evolutionary steps (i.e. fixed mutation or substitutions) across 50 independent simulations for selection (orange) and no-selection (blue) scenarios.

Evolution results in more connected networks that are less nested. Mean connectance (left), normalized nestedness (quantified as NODFc, middle) and standard deviation of species degrees (right) throughout eco-evolutionary trajectories (halted after stability in number of species is achieved over 350 evolutionary steps) for different values of environmental perturbation Δ (rows) and interspecific competition ρ (columns). Mean values (points) and standard deviation (vertical lines) are shown every 25 evolutionary steps (i.e. fixed mutation or substitutions) across 50 independent simulations for selection (orange) or no-selection (blue) scenarios.

Evolution increases critical competition above effective competition. Mean effective competition ρeff, critical competition ρc (middle) and their relative distance D (right) throughout eco-evolutionary trajectories (halted after stability in number of species is achieved over 350 evolutionary steps) for different values of environmental perturbation Δ (rows) and interspecific competition ρ (columns). For reference, ρ is also displayed with dotted lines in left and middle columns, and D=0 in the right column. Mean values (points) and standard deviation (vertical lines) are shown every 25 evolutionary steps (i.e. fixed mutation or substitutions) across 50 independent simulations for selection (orange) or no-selection (blue) scenarios. Evolution under selection always leads to systems with ρeff<ρ and ρc>ρ, implying D>0 (ρeff<ρc).

Selection increases structural stability across evolutionary time. (Left) Predicted structural stability, Δc, for two interspecific competition values and without environmental fluctuations (upper row) and with fluctuations of amplitude Δ=0.1 (lower row). (Right) Test of structural stability with simulations. Systems at steady state throughout the evolutionary simulations were tested numerically for structural stability every 25 evolutionary steps (x-axis) in the evolutionary scenarios of ρ and Δ combinations indicated in the right hand side. Each system was subjected to perturbations of amplitude Δ (shown on the top) and their dynamics investigated. Fifty replicates of ecological dynamics, corresponding to the 50 evolutionary replicates at the corresponding evolutionary time-points (i.e. number of substitutions), were simulated for each system and amplitude, and the fraction of simulations with at least one extinction is represented in the y-axis.

Equations

Full text

Introduction

Unveiling the drivers of persistence and stability of complex ecological networks remains a primary focus of ecological research. Seminal work by MacArthur in the 1960s suggested a positive relation between the complexity of ecological interactions and stability [1]. A decade later, May challenged this notion with mathematical proof of an opposing result, showing the existence of a limit in the size and connectivity of an ecological system for it to be stable [2]. May’s work also suggested that mutualistic interactions are less compatible with the stability of large ecosystems compared with other interaction types: ‘It is tempting to speculate that such community stability considerations may play a role in explaining why mutualism is a fascinating biological topic but its importance in populations in general is small’ (May in [3] quoting Williamson). This claim seems paradoxical considering the significant biomass hosted by flowering plants across the world’s biomes and the extensive biodiversity of their pollinators [4]. These systems are fundamentally supported by mutualistic interactions. Furthermore, mutualistic interactions in macroscopic [5], as well as in microscopic organisms [6] are widespread across many ecological contexts.

Lotka–Volterra (LV) models of interacting populations dynamics have played a pivotal role in our theoretical understanding of community persistence. A central focus of research grounded on these models has been on comparing the stability of purely competitive communities with those incorporating mutualistic interactions. Reported negative effects of mutualisms on community persistence have been linked to the absence of saturating responses for mutualistic effects. Ironically, as early as in 1979 Goh showed that nonlinear models of mutualism were more tractable than competitive or prey–predator ones [7]. Nowadays, it seems clear that nonlinear responses lead to stable systems sustaining biodiversity [8,9], a result that is robust even in the face of invasions by alien species [10].

Another critical factor influencing the relationship between mutualisms and community stability is the consideration that, although models with random parameters are extremely valuable as null models (see e.g. [11]), explanatory models often require non-random parameters. Indeed, interest in the study of mutualistic interaction networks surged partly from field experiments reporting non-random nested patterns in the organization of their links [12]. These observations prompted the development of theoretical approaches investigating the importance of network topology for the stability of simulated mutualistic communities.

While some researchers have reported a positive relationship between nestedness and the stability of mutualistic communities [9,13], others have suggested a more prominent role of other topological properties such as the connectance or the degree distribution within these networks [14]. These discrepancies have led to the suggestion that fair comparisons of the stabilizing effects of different topological properties should only be made after ensuring the feasibility of the systems being compared [15–17]. This is particularly important for comparisons between systems showing substantial differences in the number of parameters, e.g. purely competitive communities versus those considering mutualistic and competitive interactions together. Random parameter value selection confounds feasibility and stability, a distinction not always acknowledged [18], making systems with fewer parameters appear more stable.

Since feasibility is a necessary condition for species coexistence, there is renewed interest in the study of the effect of environmental perturbations on ecosystems’ feasibility, under the umbrella of structural stability analysis. A system is said to be structurally stable if small perturbations in the system’s parameters do not change the topology of the dynamical phase space. Introduced by Pontrjagin & Andronov in 1937 [19] and generalized to $[eqn]$ -dimensions by Peixoto [20], this concept slowly gained momentum in ecology through Levins’ [21] and Gilpin's [22] feasibility analysis alongside Vandermeer [23], and Svirezhev & Logofet mathematical developments [24]. Later Chesson’s coexistence conditions [25,26] inspired the application of structural analysis to connect competition to maximum biodiversity limits [9,27], and the introduction of the feasibility domain volume—estimating coexistence probability for a given interaction matrix—provided a more explicit quantitive link to structural stability [15,17].

However, the volume of the feasibility domain does not assess the structural stability of a specific system, namely one in which the abundances and growth rates at equilibrium are known. To tackle this question, other structural stability measures have been proposed. Among proposals in this direction we find the critical perturbation, $[eqn]$ , a measure that we will exploit in this paper, which quantifies the maximum perturbation of the environment that affects the growth rates that a specific system can withstand before at least one extinction is statistically observed [16]. Another example is a structural analog of the fitness difference, quantifying how far a given system is from the most structurally stable state [28]. Similar metrics and refinements have been later proposed, combining the volume and shape of the feasibility domain alongside the position of the system’s growth rates [29–31].

Still, most structural stability metrics focus on the effect of perturbations to growth rates rather than on interaction parameters, which may have favoured an interest in the study of the effect of random modifications in the interactions or in species assembly processes. Recent work investigating the evolutionary assembly of mutualistic networks argues for the existence of mutualistic regimes favouring biodiversity. These efforts have revealed the existence of structural patterns such as nestedness [32] or modularity [33] which, as was the case for dynamical stability, may also depend on the assumptions made. For example, the absence of explicit ecological dynamics [34], a focus on indirect mutualism [35], the explicit consideration of demographic stochasticity [36] or species traits [37] have led to different conclusions on the emergence of structural properties in mutualistic networks. These potential sources of discrepancy across studies highlight the need for a unified eco-evolutionary framework considering the minimal set of ingredients capable of generating expected structural outcomes of mutualistic networks depending on different choices of mechanisms and parameters.

From the ecological side, analysis of the critical perturbation in mutualistic systems made it possible to find that mutualistic interactions favour biodiversity only under certain conditions. The strength of interspecific competitive interactions among species within the same trophic level (e.g. plants or pollinators) must be smaller than a threshold termed critical competition, $[eqn]$ , for a mutualistic system to be more structurally stable than its purely competitive system counterpart [16]. The value of metrics describing what topological properties make these systems more structurally stable also changes depending on whether the strength of interspecific competition is above or below this critical value. Therefore, the critical competition sets a reference to explain how different parameterizations may lead to different outcomes, explaining previous discrepancies regarding the effect of mutualisms on complex ecosystems and the topological properties relevant for biodiversity maintenance [38,39].

From the evolutionary side, a key difficulty comes from combining frequency-dependent evolution within populations with realistic ecological dynamics requiring an appropriate control of species densities [40]. A potential simplifying assumption comes from considering monomorphic populations, with low mutation rates [41,42], facilitating a separation in time-scales between the ecological and evolutionary dynamics, assumptions that led to recent advances in eco-evolutionary LV models [43].

In this work, we depart from these ecological and evolutionary results to explore whether it is possible to predict the evolutionary outcome (expected biodiversity) from the properties of a starting unperturbed system. In particular, how different considerations regarding the connectivity and endurance to change of these systems may influence such predictions. There is increasing interest in relating environmental conditions and network properties (see e.g. [44,45]), given the rapidly changing environmental conditions we are witnessing.

We explored the effects that the combined action of evolutionary events (in the form of changes in species interactions and subsequent species-level selection) and environmental fluctuations have on the biodiversity of mutualistic systems (figure 1a). More specifically, we asked: (i) whether it is possible to predict biodiversity loss without species reintroduction; (ii) to what extent the behaviour observed for the structural stability of the mutualistic network, i.e. mutualistic benefits depending on the value of intraguild competition relative to a critical threshold, still holds when evolutionary events are considered (illustrated in figure 1b); and, as a corollary, (iii) what mutualistic regime (beneficial or detrimental for stability) systems tend to evolve towards (interspecific competition below or above critical competition).

Eco-evolutionary model and critical competition ρc. (a) Schematic representation of the model. At each step of the simulation, a randomly selected species (green square) is subjected to an evolutionary event (box ‘evolutionary dynamics’). Evolutionary events are of two types: link swap (E1) or creation/loss (E2). The probabilities of events (p, q respectively) are indicated. The effect of the event on the network is evaluated by running the ecological dynamics, where environmental perturbations may be added. Acceptance of the event depends on whether selection is considered. If selection is considered, the change is accepted only if it enhances the fitness of the modified species (i.e. there is an increase in its biomass). If there is no selection all changes are accepted, regardless of their effect on species or the community. (b) Simulated time series of species abundances illustrating the trade-off between the interspecific competition parameter ρ and the critical competition ρc. (c) Example of eco-evolutionary dynamics showing six trajectories for no selection (blue) and selection (orange) scenarios. The evolutionary equilibrium is achieved when no extinctions are observed throughout 350 substitutions (i.e. accepted evolutionary events). (d) log–log representation of the critical competition against the total biodiversity for a set of real plant–pollinator networks extracted from the Web of Life. The power–law relationships for both plants and animals have an exponent approximately −0.5. By means of the trade-off illustrated in (b), observing positive mutualistic effects in large communities demands decreasing the intraguild interspecific competition to favour coexistence.

Our simulations demonstrate that less competitive systems are not only more structurally stable but also more robust against changes brought about by evolutionary events. Crucially, estimation of the structural stability of the system before the evolutionary process is started is a good predictor of the evolutionary outcome in the presence of selection, especially when nuanced environmental perturbations are considered. We show that the main mechanism sustaining biodiversity throughout evolutionary time is an increase in the system’s critical competition. Consequently, the system increases its chances of stability in the face of environmental perturbations while maintaining mutualistic interactions in the regime where mutualism positively contributes to structural stability, promoting biodiversity maintenance. We would thus expect ecological networks shaped by natural selection to exhibit a tendency to keep effective competition below the critical value of competition to maintain biodiversity and withstand perturbations.

Results

The interplay between effective and critical competition determines stability and diversity

2.1.

We developed an eco-evolutionary model for mutualistic networks by coupling an ecological community dynamics model with an evolutionary model of species interactions, illustrated in figure 1a. Our approach makes two important assumptions. First, we assumed that evolutionary changes occurring at the level of species interactions are representative of the underlying changes in the corresponding species traits that enable these ecological interactions. Second, we also assumed a separation of time-scales between the ecological and evolutionary rates (see §4), with evolutionary changes happening only once the ecological equilibrium is achieved—hence omitting changes in shorter time-scales that could be interpreted as ‘adaptive rewiring’ (e.g. [46]).

The model proceeds as follows: we begin by setting the system at a feasible and dynamically stable state. Next, we introduce an evolutionary event by randomly selecting a species that will either: (i) swap one of its interactions with another randomly selected species with probability $[eqn]$ ; (ii) establish a new interaction with probability $[eqn]$ ; or (iii) remove it with probability $[eqn]$ . Unless otherwise stated, we present simulations for $[eqn]$ (see §4 and electronic supplementary material for other scenarios). Once the change is introduced, small perturbations emulating environmental noise might be introduced as uniform changes on the growth rates $[eqn]$ of all species in the system. The parameter $[eqn]$ controls the perturbations strength, by determining the maximum relative amplitude that growth rates of species will experience upon perturbation. More precisely, the perturbed growth rate $[eqn]$ is related to the unperturbed rate $[eqn]$ by the formula $[eqn]$ , where $[eqn]$ is a random number from a uniform distribution. At the maximum perturbation ( $[eqn]$ ), a species’ growth rate can be reduced to zero.

The resulting model is then integrated until equilibrium. Once the new equilibrium is achieved, the process is iterated by introducing a new evolutionary event, environmental perturbation and integrating. In addition, we may consider the presence of selection if the species selected to introduce an evolutionary change must increase in biomass for the change to be accepted (the mutation becomes a substitution). Our choice may be viewed as a numerical approximation to the ‘invasion implies substitution’ criteria found by Cressman et al. for LV systems [43], detailed in §4. If there is no selection, all mutations become substitutions independently of the fate of the species targeted. The simulations continue until an evolutionary equilibrium is achieved. Evolutionary stability is usually defined as a community-level state in which no mutants can invade the community [47]. However, this criteria is not valid for simulations under no selection, because all mutations are accepted, hindering the comparison between both regimes. Since it was shown that structural stability determines the biodiversity the system can host [9], we focused on evolutionary persistence as an alternative criteria that stops the simulations when no extinctions are observed over a large number of substitutions (350)—a criteria valid in both selection and non-selection scenarios (see §4).

The ecological model considered that mutualistic interactions between two pools of species (e.g. plants indexed $[eqn]$ and their animal pollinator species indexed $[eqn]$ ) were encoded in a network of interactions $[eqn]$ with a mean interaction strength $[eqn]$ (see §4 for a full description and figure 1 for a schematic representation). Competitive interactions within each pool were encoded in mean-field matrices $[eqn]$ defined with a parameter $[eqn]$ for intra-specific and $[eqn]$ for inter-specific competition ( $[eqn]$ ) (see figure 1a and §4). It is possible to rewrite the system containing both competitive and mutualistic interactions (and/or predatory, see [27]) as a purely competitive system, in which non-competitive interactions are effectively incorporated as variations in the bare interspecific competition $[eqn]$ , becoming an effective competition parameter $[eqn]$ . The most important advantage of the simplification provided by the formalism derives from its ecological interpretability, partly thanks to the relationship between the stability of the effective system and the whole system [48]. For instance, for mutualistic systems, it has been verified that $[eqn]$ , meaning that mutualistic interactions effectively reduce the interspecific competition [9]. The effect that specific network topologies have on the stability of the system can also be explored similarly [16]. Crucially, mutualistic benefits on structural stability were observed only if the effective interspecific competition parameter $[eqn]$ was lower than a threshold, $[eqn]$ , termed critical competition (figure 1b) [16].

In this work, and consistent with previous results, we observed that the relationship between $[eqn]$ and $[eqn]$ determined different scenarios in which evolutionary changes showed a distinct behaviour depending on whether selection was considered or not. In figure 1c, we illustrate this for different runs of a simulation in which the same starting network was evolved under no selection and selection scenarios. Evolution under selection consistently converged earlier to systems with a higher diversity. Since in real networks a negative relationship between diversity and critical competition has been observed (figure 1d), an open challenge is to derive the expected relationship between their effective and critical competitions to theoretically understand the mechanisms behind their stability. We aim to tackle this challenge by investigating the effect that evolutionary events and environmental perturbations have on $[eqn]$ and $[eqn]$ , which influence the structural stability and, in turn, the diversity of the system.

Structural stability explains diversity at evolutionary persistent states

2.2.

We started by selecting a real network from the dataset available from the Web of Life (www.web-of-life.es/) to be evolved. We chose a network studied in detail in previous work [16] with approximately the same number of plant (47) and animal species (46) and connectance close to the median of the entire set $[eqn]$ [49]. This choice allowed us to start with a critical competition value approximately in the middle of the range observed for the set of real networks ( $[eqn]$ in figure 1d). In this way, changes in these quantities throughout the simulation will more likely fall within ranges observed for other real networks. In addition, using the same initial mutualistic network across model simulations we could fairly compare the effect that the different ecological and evolutionary parameters had on the evolutionary outcome. We present in electronic supplementary material, figure S1 results for the evolution of randomly generated networks, finding similar results.

We ran evolutionary simulations varying the value of intraguild interspecific competition $[eqn]$ and the level of environmental perturbations from $[eqn]$ (no perturbations: any species loss should be attributed to evolutionary events alone) to $[eqn]$ (10% of the maximum perturbation, see §4 for details). The maximum value of $[eqn]$ (0.1) was chosen because it determines the structural stability of the most vulnerable system (the one with highest competition, $[eqn]$ ). Structural stability was defined as the minimum amplitude $[eqn]$ of the perturbation on species’ growth rates leading to extinctions. We termed this value of the amplitude the critical perturbation $[eqn]$ which can be estimated analytically, and verified numerically with simulations (see §4 and [16] for details). Unless otherwise stated, for computational convenience we will use the analytical calculation. Therefore, $[eqn]$ was close to the $[eqn]$ of the less structurally stable system.

We found that the lower the interspecific competition, the higher the evolutionary persistence, measured as species $[eqn]$ −diversity (figure 2a). Interestingly, the specific relationship between diversity and competition changes depending on the presence and magnitude of environmental perturbations and on the existence of selection. In the absence of perturbations, the dependence of the final diversity of the system on the interspecific competition $[eqn]$ is weaker than if perturbations are present (figure 2a). This dependence is, however, more nuanced under selection, with interspecific competition having little effect on biodiversity in no-selection simulations with strong environmental perturbations. The relative values between selection and no-selection simulations also change depending on the specific value of the environmental perturbations, $[eqn]$ . Selection is capable of preserving more diversity than no-selection random change except for scenarios of relatively high interspecific competition in the absence of perturbations ( $[eqn]$ ). For high levels of perturbation, the loss of species observed for no-selection simulations was up to approximately fivefold higher than in the presence of selection (for $[eqn]$ .1 and $[eqn]$ ).

Structural stability explains diversity at the evolutionary equilibrium. (a) Biodiversity of the eco-evolutionary process at equilibrium for different values of the interspecific competition ρ and environmental perturbation Δ. Each panel represents the outcome of the evolutionary process for different values of Δ. In boxplots, box limits represent interquartile range (IQR) and the horizontal line the median value of the final species richness (α-diversity) across 50 simulations. Vertical lines (i.e. whiskers) show the ±1.5 IQR values. Points represent the outliers. (b) Predicted structural stability (critical environmental perturbation Δc) for the starting network for each value of interspecific competition ρ. (c) Relationship between structural stability and final biodiversity, normalized by the starting diversity before the evolutionary process. (d) Slope of the relationship shown in (c) for different values of environmental perturbation. The closer the slope is to one the more accurate the prediction. Ordinates are not significantly different from zero for Δ=0.1.

In electronic supplementary materials, we extended the values of $[eqn]$ explored by setting perturbations relative to the $[eqn]$ of each system (i.e. $[eqn]$ with $[eqn]$ ). We found that more competitive systems still lose more species for the same relative perturbation (i.e. same value of $[eqn]$ ), even though the absolute perturbation $[eqn]$ is much higher for less competitive systems (electronic supplementary material, figure S2a). Moreover, we found that for each interspecific competition $[eqn]$ value there is a characteristic environmental perturbation $[eqn]$ in which the relative difference in final $[eqn]$ diversity between simulations under selection and no-selection scenarios peaked, with both quantities related as $[eqn]$ (electronic supplementary material, figure S2b–d). Therefore, given $[eqn]$ it is possible to estimate the region of environmental perturbations amplitudes in which selection effects are stronger.

To quantify the agreement between evolutionary persistence and structural stability, we compared the $[eqn]$ diversity persistence (our proxy of evolutionary stability) against the predicted structural stability, $[eqn]$ , shown in figure 2c (see §4). Notably, evolutionary persistence goes hand in hand with structural stability showing a similar trend, with less competitive networks (i.e. those with higher $[eqn]$ ; figure 2b) losing fewer species than those with high competition.

Interestingly, representing the final diversity against the structural stability (figure 2c) we observed that structural stability had a good qualitative agreement in all scenarios, following significant linear trends in all cases. Moreover, as soon as environmental perturbations were present ( $[eqn]$ ), structural stability accurately explained the final diversity, with slopes approaching one (note that final diversity was normalized by the starting diversity to get a more direct interpretation of the slopes, shown in figure 2c). In figure 2d, we summarize variations in the slopes of these relationships across the gradient of environmental perturbation $[eqn]$ . In the absence of environmental perturbation, structural stability underestimates the final expected diversity independently of the presence of selection, especially for highly competitive systems (lower structural stability; figure 2b). This is illustrated by the slope of the relationship in figure 2c being lower than one (figure 2d). For intermediate levels of perturbation, structural stability tends to overestimate the final diversity in the scenario under selection, and to underestimate it in the no-selection case. For high levels of perturbation, on the other hand, it tends to accurately explain the final diversity (slope tends to one) in the presence of selection, and strongly underestimates it for no-selection simulations. Importantly, in electronic supplementary material, figure S3 we show that our results are robust for networks with different ratios in the number of species of plants and pollinators, or with a similar ratio and different number of total species, with slopes $[eqn]$ close to one in all cases : $[eqn]$ with $[eqn]$ (104 species); $[eqn]$ with $[eqn]$ (97 species); and $[eqn]$ with $[eqn]$ (78 species).

In summary, structural stability qualitatively describes the evolutionary outcome, and it can be used as a predictor of biodiversity in simulations under selection, in particular in the presence of higher levels of environmental perturbations.

It is important to note that structural stability predictions consider how environmental perturbations acting on the growth rates of the starting system affect species persistence. On the other hand, evolutionary events affect the interactions between species, either swapping a randomly selected interaction between two species, introducing a new one, or removing a randomly selected interaction. Moreover, structural stability is defined on the starting system but, during the evolutionary process, the system will lose species, thus changing the structural stability of the system as evolution progresses. Therefore, the starting structural stability of the system sets the boundaries within which evolutionary biodiversity loss will lie.

Biodiversity decline follows increasing biomass and connectance

2.3.

To gather deeper insights into the impacts of evolutionary change on community and network structure, we selected specific simulations for further analysis. We focus on the two extreme scenarios of environmental perturbations ( $[eqn]$ ) and of interspecific competition ( $[eqn]$ ), and assess their evolutionary trajectories in detail.

In the absence of environmental perturbation ( $[eqn]$ ), the qualitative behaviour for the loss of biodiversity was different for low and high competition (figure 3). For low competition under selection the system stabilized very quickly, conserving most species, while no-selection simulations took longer times to stabilize, losing more species. On the contrary, for high competition the system under selection lost species faster but stabilized after approximately 1500 substitutions, while no-selection simulations had a slower decay, and they were not stabilized after the limit of 2000 substitutions. For high environmental perturbations ( $[eqn]$ ), the loss of species was exacerbated for both low and high competition, with both no-selection and selection simulations having a very similar behaviour for high competition (figure 3). This is due to the fact that, for this parameterization, the environmental perturbation corresponds to the critical threshold of environmental perturbation, $[eqn]$ , hence being highly disruptive.

Disruptive evolution prompts biodiversity loss and biomass increase. Mean α-diversity (i.e. species richness) (left), effective diversity (middle) and biomass (right) throughout eco-evolutionary trajectories (halted after stability in number of species is achieved over 350 evolutionary steps) for different values of environmental perturbation Δ (rows) and interspecific competition ρ (columns). Mean values (points) and standard deviation (vertical lines) are shown every 25 evolutionary steps (i.e. fixed mutation or substitutions) across 50 independent simulations for selection (orange) and no-selection (blue) scenarios.

In the presence of selection, the loss of biodiversity went hand in hand with a strong increase in total community biomass, except for scenarios of high environmental perturbations and high competition. An increase was expected, as the condition for selection is that the species subjected to evolutionary change must increase in biomass for the change to be accepted. Therefore, without perturbations, a faster loss of biodiversity under selection may be due to a destabilizing effect of mutualistic interactions when biomass increases, as shown below. Notably, biomass was distributed unevenly across species in all systems, with some species concentrating more biomass, as indicated by a lower effective diversity compared with species $[eqn]$ diversity (figure 3).

Biodiversity loss and biomass increase were accompanied by a significant increase in the connectance of the mutualistic network (see figure 3). In the absence of environmental perturbations the increase was clearly a selected feature, as evidenced by the well-separated trend with respect to no-selection simulations. When high environmental perturbations were present, however, the increase in connectance seemed to be a collateral effect arising from the loss of species under high interspecific competition, with weak evidence for selection when competition was low. To confirm this interpretation, we ran simulations for two additional values of the parameter $[eqn]$ , which controls the probability of performing an interaction insertion or deletion versus a swap in the evolutionary events. When the probability of inserting new links was low and interactions rewiring (i.e. swaps) was favoured ( $[eqn]$ ), there was a more dramatic loss of biodiversity under selection (electronic supplementary material, figure S4a). Therefore, under this scenario, the increase in connectance was mainly driven by species loss (electronic supplementary material, figure S4). On the contrary, for high probability of interaction insertion/deletion $[eqn]$ ) similar levels of connectance were achieved without losing species (electronic supplementary material, figure S4b), pointing to the role of the establishment of new interactions in driving increases in connectivity. Interestingly, trajectories stabilized for similar levels of connectance across different values of the parameter $[eqn]$ , but the number of species varied. These results suggest that a central mechanism to sustain biodiversity possibly relies on building a core of well-connected species. This is also confirmed by the fact that the number of connected components in the network (the number of groups of nodes that did not share links with nodes outside of it) consistently converged to a single group (electronic supplementary material, figure S5).

Differences between trajectories under selection and no-selection seem to be common across most topological properties analysed: in the absence of environmental perturbations there are significant differences while, if there are environmental perturbations, these differences are maintained only when interspecific competition is low. To provide further evidence, all networks under selection, except the one at high competition and high environmental perturbations, evolved their connectances towards a more skewed distribution of the degrees ( $[eqn]$ figure 4, right column), being more disassortative for selection when environmental perturbations were absent and having similar values when perturbations were present (electronic supplementary material, figure S5). Nestedness is expected to increase for more disassortative networks and, indeed, when it was estimated with the measure of ecological overlap (namely, the mean fraction of shared mutualistic partners across species, see §4 [9]) it strongly increased along evolutionary time (electronic supplementary material, figure S5). However, when it was controlled by the expected nestedness of random networks with the same size and connectance using the NODFc measure (see §4; [50]), it decreased for all simulations unless competition was high and environmental perturbations present (figure 4, middle column). Therefore, the disruptive nature of the evolutionary process we proposed limits our capacity to elucidate the role of nestedness on the final outcome of network assembly, a question discussed in more detail below.

Evolution results in more connected networks that are less nested. Mean connectance (left), normalized nestedness (quantified as NODFc, middle) and standard deviation of species degrees (right) throughout eco-evolutionary trajectories (halted after stability in number of species is achieved over 350 evolutionary steps) for different values of environmental perturbation Δ (rows) and interspecific competition ρ (columns). Mean values (points) and standard deviation (vertical lines) are shown every 25 evolutionary steps (i.e. fixed mutation or substitutions) across 50 independent simulations for selection (orange) or no-selection (blue) scenarios.

In summary, the disruptive evolutionary processes we built led to networks with fewer species, higher and unevenly distributed biomass, and higher, unevenly distributed connectances in the presence of selection, with nestedness having a secondary role in these networks.

Evolution increases critical competition, favouring positive mutualistic effects

2.4.

As shown above, the sign of mutualistic effects depends on the relative value of $[eqn]$ versus $[eqn]$ . Mutualistic interactions are detrimental for stability if $[eqn]$ and beneficial if $[eqn]$ . Therefore, we considered the normalized distance of $[eqn]$ to $[eqn]$ , $[eqn]$ , defined between −1 and 1, to monitor the mutualistic regime, with $[eqn]$ favouring and $[eqn]$ being detrimental to stability (see §4). Since evolutionary changes modify both $[eqn]$ and $[eqn]$ , we are interested in the mutualistic regimes that are favoured throughout the eco-evolutionary process we considered.

We observed that, independently of the environmental perturbations, the critical competition of the system tended to increase (figure 5, middle column). However, while for no environmental perturbations this increase was a consequence of the selected topological properties (with no overall increase for no-selection simulations), in the presence of environmental perturbations it seemed to be a byproduct of the loss of species. For the effective competition, selection generated an overall reduction which, for no-selection simulations, depended on the environmental perturbations. For both critical and effective competition we observed a shift within the first approximately 500 substitutions aligned to changes in the topology of the mutualistic network, in particular to the variance of the degrees (compare figure 5 with figure 4).

Evolution increases critical competition above effective competition. Mean effective competition ρeff, critical competition ρc (middle) and their relative distance D (right) throughout eco-evolutionary trajectories (halted after stability in number of species is achieved over 350 evolutionary steps) for different values of environmental perturbation Δ (rows) and interspecific competition ρ (columns). For reference, ρ is also displayed with dotted lines in left and middle columns, and D=0 in the right column. Mean values (points) and standard deviation (vertical lines) are shown every 25 evolutionary steps (i.e. fixed mutation or substitutions) across 50 independent simulations for selection (orange) or no-selection (blue) scenarios. Evolution under selection always leads to systems with ρeff<ρ and ρc>ρ, implying D>0 (ρeff<ρc).

The complex interplay between both critical and effective competition explain differences in species loss between different parameter and selection regimes. In particular, for high competition we observed that the relative distance between effective and critical competition ( $[eqn]$ ) run from negative to positive values (figure 5). When it was negative, mutualism had destabilizing effects and the strong increase in biomass observed for simulations under selection led to a dramatic loss of species. The system became stabilized, with no more species lost, when $[eqn]$ became positive (crossing the dotted line in figure 5).

Evolution increases structural stability, favouring biodiversity maintenance

2.5.

We showed that the structural stability of the system at the beginning of the simulations sets the limits in which species loss will lie in an evolutionary process in the presence of environmental fluctuations. The structural stability, however, depends on the number of species and topological properties of the system, which vary throughout evolutionary trajectories. We were therefore interested in understanding how the structural stability itself evolved.

Notably, the structural stability of the system, estimated by the predicted critical perturbation, $[eqn]$ , experienced a long-term increase in the presence of selection (figure 6, left). While in some cases (e.g. no-selection, high competition and high environmental perturbations) this was a trivial consequence of having few highly connected species, this was not the case in other simulations (e.g. under selection, low competition, high environmental perturbations) where the number of species was above 50 (figure 3) and the connectance lower than 0.5, with a threefold increase in biomass (figure 6, left). Importantly, in some cases, the long-term increase occurred only after a nuanced decrease, consistent with a similar shift observed for the effective competition, and suggesting that a topological rearrangement of the system that weakens its structural stability is a prerequisite to later enabling its increase.

Selection increases structural stability across evolutionary time. (Left) Predicted structural stability, Δc, for two interspecific competition values and without environmental fluctuations (upper row) and with fluctuations of amplitude Δ=0.1 (lower row). (Right) Test of structural stability with simulations. Systems at steady state throughout the evolutionary simulations were tested numerically for structural stability every 25 evolutionary steps (x-axis) in the evolutionary scenarios of ρ and Δ combinations indicated in the right hand side. Each system was subjected to perturbations of amplitude Δ (shown on the top) and their dynamics investigated. Fifty replicates of ecological dynamics, corresponding to the 50 evolutionary replicates at the corresponding evolutionary time-points (i.e. number of substitutions), were simulated for each system and amplitude, and the fraction of simulations with at least one extinction is represented in the y-axis.

To validate structural stability predictions, we considered networks evenly distributed throughout the evolutionary trajectories and tested their structural stability by perturbing the steady state found at the corresponding evolutionary time-point (see §4). We considered different amplitudes of the perturbations and calculated the fraction of simulations in which at least one species went extinct.

Simulations were in agreement with predictions across all combinations of selection regime, competition and environmental perturbations (figure 6) with increasing $[eqn]$ associated with a decrease in the fraction of simulations reporting extinctions. Importantly, we again observed a shift in which the system was less stable (increasing fractions of simulations with observed extinctions) to then become more structurally stable, especially for scenarios of high environmental perturbation and low competition (figure 6b), where the predicted structural stability of the system takes the larger relative fall (figure 6a).

Our results showed that, in the face of disruptive eco-evolutionary changes under selection, mutualistic systems will tend to lose species and to rearrange their topology to increase the critical competition above the effective one. This ensures that mutualism operates in the competition regime in which it favours stability.

Discussion

Ecological and evolutionary perspectives on the importance of mutualisms for the emergence and persistence of biodiversity have revealed different facets of mutualistic networks assembly that influence their stability (reviewed in [38,39]). Studies focusing on the evolution of mutualistic interactions have shown that interaction rewiring and speciation with inheritance are key ingredients to the emergence of structure in mutualistic networks ([32,34]). On the other hand, those focusing on the stability effects of mutualisms on system dynamics at ecological time-scales, excluding evolutionary or assembly events (sometimes called internal stability [10]), have highlighted the crucial role of mutualistic interactions for community stability [8,9,15,16,51,52]. These two angles have been, to a great extent, so far disconnected, possibly due to a lack of theoretical framework unifying ecological dynamics, evolutionary changes and invasions.

Indeed, most theoretical studies either explore the formation of patterns without taking into account the dynamics (e.g. [34]), or excluding the effects of environmental perturbations (e.g. [32]). Here, we contributed to fill this gap by studying how community coexistence becomes affected when evolutionary events are introduced and the ecological dynamics are explicitly simulated in the presence of environmental perturbations.

Our results showed that the equilibrium diversity of the eco-evolutionary process can be explained from the critical perturbation that the starting system can withstand, to the point it could be used in some cases as a predictor. We observed predictability increased under selection, a result reported in evolutionary research [53]. To understand this result we should keep in mind that the critical perturbation is the maximum amplitude of environmental perturbations the system can cope with until at least one extinction is observed. On the other hand, the equilibrium diversity in the eco-evolutionary process was obtained after introducing repeated perturbations well below the critical value. Therefore, our results suggest a connection between repeated small perturbations and single large ones. A minimal amplitude of environmental perturbations seems to be a prerequisite to our ability to accurately explain outcomes of biodiversity, possibly due to its role in overriding the effects of evolutionary changes.

We built on a well-studied population dynamic model for which the relationship between mutualism and biodiversity depends on the effective competition that species experience after accounting for all indirect pathways. By simulating evolutionary events in different environmentally changing scenarios, we aimed to unveil the mutualistic regime (positive versus negative for sustaining biodiversity) in which we could expect mutualistic systems to operate. We believe our approach provides an understanding of the ecological baseline on which other processes such as invasions or adaptation would play out (see below).

Our results suggest that, especially in the presence of selection, disruptive evolutionary change can drive mutualistic systems towards a regime in which mutualism supports biodiversity. We observed an increase in the total biomass of the system unevenly distributed across species, an increase in network connectance and, for selection, a decrease in the evenness of mutualistic degrees. Some of these patterns were a direct consequence of the loss of species (e.g. connectance in the presence of high environmental perturbations) but, in some simulation regimes, they were selected for.

The systematic increase of a well-connected core of species observed in our simulations seems to be aligned with results suggesting that connectance is underestimated in field experiments due to the typically short observations time-scales [54], and with the importance of core species in nested patterns [55]. It seems, however, to contradict findings pointing to the importance of an increase in modularity [33]. In our view, the discrepancy stands from the type of perturbations considered. While Sheykhali et al. [33] considered targeted attacks on specific species, we considered perturbations on model parameters, reflecting stochastic variation of environmental conditions. Directly targeting species abundances may make sense under some circumstances (e.g. human activity such as in fisheries or the use of certain types of antibiotics). In other scenarios, however, such as environmental conditions influencing species performance and, in turn, their abundances, the perturbations used in our study represent a more natural approach to model these conditions. Another explanation may come from differences in the selection strength since an increase in network connectance was an observed evolutionary outcome when selection was strong, while modularity was observed when selection was weak [56].

The main consequence of the topological and biomass changes we observed was an increase in the critical ( $[eqn]$ ), and an overall decrease in the effective ( $[eqn]$ ), competition. A reduction of the effective competition can be interpreted as a decrease in the species niche overlap, a result reported in an evolutionary model for purely competitive systems [57], which was also associated to an increase in the structural stability as measured by the structural resistance [29], conceptually similar to the one we used. In addition, the relationship between $[eqn]$ and $[eqn]$ determine the mutualistic regime, stabilizing when $[eqn]$ (figure 1b) and destabilizing when the opposite is true. Importantly, in previous work we observed that in real networks a negative relationship between the critical competition and the number of species exists (figure 1d). Considering our results and assuming that real networks also operate in a regime in which mutualism positively influences stability, would imply that the larger the network is, the smaller the effective competition must be. This is an important finding, given the lack of experimental data to determine specific parameter values in theoretical modelling, a question that so far remains unsolved, particularly for competitive interactions [58]. This highlights the need to design sound experimental settings to measure and quantify specific demographic, environmental and ecological (e.g. interaction strengths) dimensions of ecological systems. This research avenue should be prioritized in future studies aimed at understanding the drivers of the emergence and maintenance of biodiversity.

Our results also suggest that determining the critical competition of mutualistic networks would allow theoreticians to determine regions of the parameter space describing their different behaviours, which can then be confronted with experimental data. For example, Lampo et al. recently reported that seasonal changes may have an influence in the relative dominance of different structural properties such as the nestedness and the modularity (a property we did not explore here) [59]. Interestingly, the authors found a non-linear relationship between the connectance and the feasibility of synthetic systems, which could be explained by the difference between the effective and critical competition, both changing with the connectance. Therefore, linking system properties with environmental disturbances can have important consequences for the monitoring of ecosystems’ health and the design of management interventions [60,61]. As an illustration, we noted that, given an estimation of the interspecific competition, we could predict the region of environmental disturbance values where selection would be expected to have the strongest effects (see electronic supplementary material, figure S2).

Despite their simplicity, our results suggest connections with more complicated mechanisms. For example, the reduction in the effective competition seems to be aligned with the idea of mutualism having stabilizing effects within evolutionary time-scales [62]. This also allows us to shed light on how additional mechanisms might explain coexistence when they operate by modulating competition, such as mutualistic-mediated adaptation to new niches [63], or competition–colonization trade-offs [64]. Further, we also observed that, although systems tend to become more structurally stable than the original network in the long-term, they achieved this state after a nuanced decrease in the structural stability. Therefore, environmental perturbations and evolutionary changes drive the system towards more vulnerable configurations, where some species were purged. This result provides a quantitative example of a qualitative metaphor described in the panarchy cycle [65], where it is suggested there is a need for a stage of ‘destruction’ to withstand perturbations. We observe, however, differences in the expected role of the connectance (destabilizing versus stabilizing), the discussion of which remains beyond the scope of this work.

Our model has some limitations. Firstly, increasing the number of mutualistic interactions (i.e. beneficial effects) has no cost to a species, which may disproportionally favour more connected networks. A simple way to explore the effect of this was to reduce the probability of introducing/removing a link. In this setting, we found that increasing the probability of swapping interactions (at the expense of creating new ones) also induced an increase in network connectivity but through a different route: a decline in biodiversity. In our simulations, the sometimes abrupt decrease in the number of species and the increase in the connectance rapidly reduces the sparsity of the interactions matrix, and hence the size of the ensemble of possible networks with similar levels of connectivity, making it unlikely that significantly nested patterns will be found.

Secondly, we did not consider the introduction of new species [36,46,66]. The intended aim for this decision was to have a one-to-one comparison between evolutionary and structural stability. Since structural stability is defined as the maximum environmental perturbation a system can sustain before a species goes extinct, a direct comparison is to explore the evolutionary loss of species. We acknowledge that introducing new species may lead to other network topologies. For example, Cai et al. found similar results to ours with increased nestedness or modularity depending on model parameters [46]. Moreover, two mutualistic regimes were found in the comparison of stability between evolved and random networks when competition was varied, that could likely be explained by the critical competition threshold. Nevertheless, the pattern of variation between nestedness and modularity they observed when varying strength of mutualistic interactions may also be explained by the relative importance of nestedness and connectance in different parameter regimes [16]. Direct comparison with our results is not straightforward because they did not vary the connectance of the system. Combining their results with ours, we conjecture that each assembly event may reduce the distance between the effective and the critical competition, positioning the system at the edge of critical competition. Therefore, in our view, how a nested pattern emerges in an eco-evolutionary process with interspecific competition and environmental perturbations present is far from being understood. It may be that processes such as adaptive foraging [67] are needed to favour specialist strategies, thus leading to more sparse and nested networks [66].

A final caveat of our work comes from the evolutionary assumptions underlying our simulations, partly motivated by its computational and analytical tractability. Perhaps the most important one is that phenotypic frequencies were neglected to have control on species densities, central in the evolution of multispecies systems, and which is not satisfactorily handled in frameworks such as evolutionary games theory [40]. To continue, although we considered the set of ecological interactions at the species level as the traits subject to evolutionary change, we note that our results can be interpreted from other non-evolutionary perspectives of interaction modification. For example, as soon as the separation of time-scales assumed in our model is relaxed, our implementation is compatible with other interaction changes paradigms such as behavioural modifications or other forms of ‘adaptive rewiring’ that permit species-level adaptation in non-evolutionary contexts (e.g. [46]). We believe that the separation of time-scales assumed in our model is justified because even though rapid evolution is potentially overlooked in empirical systems, it is still rare in macroscopic organisms [68]. It is nevertheless important to remind ourselves that a suite of environmental drivers have the potential to alter ecological and evolutionary rates and, in turn, may modify this picture [69].

Recent findings have highlighted the potential impact that phenotypic variability may have on ecosystem stability [70], especially when rapid evolution is considered [71]. However, models simulating eco-evolutionary feedbacks have so far considered small systems comprising no more than a handful of species (see e.g. [72–74]). Further efforts are hence needed to investigate eco-evolutionary dynamics in polymorphic populations within large ecosystems. Recent frameworks proposed to investigate the (co)evolution of interactions in mutualistic communities driven by changes in species morphological [75] and phenological [76] traits provide a promising avenue for future research aimed at relaxing the assumptions made by our model. This can extend our findings to the explicit consideration of phenotypic traits in the evolutionary process driving ecological interactions.

In summary, our work contributes towards an eco-evolutionary theory of ecological networks combining environmental fluctuations and evolutionary events that considers mutualistic and competitive interactions together. Our results can have implications for conservation initiatives where the targets for conservation are not only species but their interactions. This is especially since species interactions, and the corresponding effects on effective competition, can be strongly affected by different aspects of global change such as habitat loss and climate warming [77–79].

Methods

Ecological community dynamics

4.1.

Following previous work [9,16,48], we model two groups of species: $[eqn]$ plants and $[eqn]$ animals. Species within each group compete with all others through fully connected competition matrices $[eqn]$ of Lotka–Volterra (LV) type. Species between groups interact via mutualistic interactions $[eqn]$ , that are non-zero only if there is a link between $[eqn]$ and $[eqn]$ in the mutualistic network. Mutualistic benefits saturate at large abundances of the mutualistic partners providing the benefits [8]. Dynamics are hence governed by the following equations:

$$[eqn]$$

We only present the equations for plants, since those for animals can be exactly obtained by interchanging the superscripts P and A. In equation (4.1), $[eqn]$ denotes the abundance of species $[eqn]$ , $[eqn]$ its intrinsic growth rate, and $[eqn]$ the maximum mutualistic growth rate of species $[eqn]$ , whose inverse, $[eqn]$ , can be interpreted as the handling time.

Parameterization of the system

4.2.

Competition and mutualistic matrices

4.2.1.

In our model, we consider a competition matrix with different parameters for intra- ( $[eqn]$ ) and inter-specific competitions ( $[eqn]$ ). The competitive matrix for plant species can be written as:

$$[eqn]$$

where $[eqn]$ is Kronecker’s delta and $[eqn]$ are dimensionless parameter values uniformly distributed in $[eqn]$ . We parameterize mutualistic interactions as

$$[eqn]$$

where $[eqn]$ measures the strength of mutualism with respect to competition, $[eqn]$ is the adjacency matrix of the mutualistic network, and $[eqn]$ are dimensionless parameter values uniformly distributed between $[eqn]$ and $[eqn]$ if $[eqn]$ and are zero if $[eqn]$ .

Empirical network and metaparameters values

4.2.2.

In our analyses, we consider as starting network a plant–pollinator mutualistic network with $[eqn]$ animal and $[eqn]$ plant species describing a rainforest mutualistic community comprising plants and their animal pollinators (insects and birds) from Cainama, Venezuela [1]. This network features approximately the same number of animal and plant species and connectance $[eqn]$ close to the median of a large set ( $[eqn]$ ) of empirical networks available from the Web of Life database (https://www.web-of-life.es/).

The dynamical equations for all species in the system depend on meta-parameters: the direct interspecific competition parameters between plant (animal) $[eqn]$ ( $[eqn]$ ) species ( $[eqn]$ ), the mutualistic handling times for plants $[eqn]$ (animals, $[eqn]$ ), the mutualistic strength $[eqn]$ , and the equilibrium abundances $[eqn]$ . We used parameter values within the regime that allow us to simulate a weak-mutualistic facultative system, the regime in which topological trade-offs were clearly found in previous work [16]. In particular, we set $[eqn]$ , $[eqn]$ and $[eqn]$ for all species in the system. For each simulation, we parameterize the initial system randomly drawing a value for each parameter following a uniform distribution around its metaparameter value $[eqn]$ : $[eqn]$ with $[eqn]$ , except for $[eqn]$ which remains constant, since we are interested in understanding the behaviour of the system across a wide range of interspecific competition values (see below). For simplicity, we set $[eqn]$ . Finally, the growth rates ( $[eqn]$ ) are determined once all the other parameters are drawn by solving the system at equilibrium, ensuring that they are positive and hence that the fixed point is feasible.

Eco-evolutionary simulations

4.3.

Evolutionary events and selection

4.3.1.

To simulate eco-evolutionary dynamics, we performed a two-step process in which ecological dynamics play out according to the dynamical equations (equation (4.1)) until the ecosystem reached equilibrium. Evolutionary events were performed over the equilibrium state of the system, and subsequently interrogated running a new ecological dynamics. Our approach implies considering a N-species monomorphic population in which mutation rates are rare [41,42], which justifies treating each mutation as an independent event without overlapping dynamics.

To support the biological plausibility of this separation of time-scales, let us consider, as an example, the per-site mutation rate in bumblebees (Bombus terrestris), estimated as $[eqn]$ mutations/site/generation [80]. Assuming a locus of 1 kb (1000 base pairs), the per-individual probability of a mutation at that locus is approximately $[eqn]$ . For a population of $[eqn]$ individuals, the per-generation probability that at least one individual mutates at the locus is on the order of $[eqn]$ .

This scenario defines a geometric distribution of waiting times for success in a Bernoulli trial, with an expected waiting time of approximately 100 generations (around 30 generations in the bumblebee case, given a slightly higher locus-wide mutation rate). In our parameterization, the starting intrinsic growth rates of species range from 1 to 5 depending on interspecific competition, and equilibrium is typically reached in fewer than 25 time units. Assuming one generation per time unit, this corresponds to 25−100 generations, aligning well with the expected evolutionary waiting times and thus justifying our modelling assumptions.

Evolutionary events at an equilibrium were introduced randomly selecting a species $[eqn]$ (including plant and animal species), on which random changes will be performed. Three possible random mutation events (interaction created, lost, or swapped, illustrated in figure 1) were considered. In the main text, we chose $[eqn]$ (see electronic supplementary material, figure S4 for other scenarios). The strength of new interactions $[eqn]$ was drawn in the same way as the initial mutualistic interaction coefficients 3. Interspecific competitive interactions were not affected by the evolutionary process. In the context of LV models, the fitness of monomorphic species at equilibrium are given by $[eqn]$ for all resident species [40]. When a mutant $[eqn]$ arises from a resident species $[eqn]$ , its abundance must be low ( $[eqn]$ ) and, since mutations only occur on traits involving mutualistic interactions, we also have that $[eqn]$ , $[eqn]$ and $[eqn]$ . Under these conditions, the mutant can invade the population if $[eqn]$ [40].

To simulate the effect of selection, we considered a computationally convenient approach to estimate whether a mutation fixates in a population (termed ‘invasion implies substitution’ theorem [81,82]). To illustrate it in the context of LV-systems, let us consider two N-species equilibria: the equilibrium before the mutated species appears, $[eqn]$ , and another one, $[eqn]$ , in which the mutant replaces the resident species with all remaining species having positive biomasses. It has been shown that, if $[eqn]$ exists and it is globally stable, $[eqn]$ implies replacement of the resident species by the mutant [43]. This theorem may, however, not hold when environmental perturbations are introduced, as we do here. As an approximation, in scenarios under selection mutations were accepted if the outcome of the ecological dynamics at the new equilibrium $[eqn]$ yielded a higher biomass for the mutated species than that of the resident species at $[eqn]$ . In the context of our model, since the complete model can be approximated at equilibrium by an effective-competition LV model ( $[eqn]$ ), a higher mutant abundance implies that the effective competition that the mutant species feels is lower than the one from the resident species, implying higher fitness. Our criteria therefore acts as a proxy for Cressman et al.’s rigorous LV substitution criteria, where invasion and replacement (stable mutant-dominated equilibrium) are implicitly ensured through numerical integration and biomass comparisons, although we should note that our criteria does not necessarily imply that the final equilibrium has the same number of species than the original one. To have a fair comparison between evolutionary trajectories with and without selection, we considered trajectories with the same number of substitutions (i.e. rejected mutations are excluded).

Eco-evolutionary process

4.3.2.

To explore the effect of the strength of competition on the evolutionary trajectories of our networks, we performed simulations across a range of values of interspecific competition: $[eqn]$ . Intra-specific competition $[eqn]$ was kept the same and equal to 1 throughout the duration of each single simulation. Further, to investigate the interplay between network evolution and the level of environmental perturbation, we explored the effect of varying the environmental perturbation on the system (i.e. $[eqn]$ ) by setting it to different values of $[eqn]$ , where 0.1 corresponds to the lowest $[eqn]$ value predicted for the network used in our simulations across all values of $[eqn]$ used (i.e. $[eqn]$ , see figure 2b).

We ran 50 simulations for each combination of environmental perturbation ( $[eqn]$ ) and interspecific competition ( $[eqn]$ ) values for both no-selection and selection scenarios. During each simulation, at each evolutionary timestep (i.e. after ecological dynamics were integrated), species whose abundance fell below $[eqn]$ were considered extinct and removed from the system. Since we are considering monomorphic populations all individuals within each species have the same phenotype. Therefore, evolutionary stability would be achieved when no species can be successfully invaded by a mutant phenotype [47]. However, this criteria cannot be fulfilled if there is no selection, because all mutations are accepted. We thus considered as an alternative criteria evolutionary persistence, which more naturally links to structural stability and its known relationship with the number of species it can host [9]. Simulations halted when the number of species was maintained across at least 350 evolutionary timesteps. This procedure yielded a total of 2800 evolutionary trajectories (7 $[eqn]$ ’s × 4 $[eqn]$ ’s × 2 selection regimes × 50 replicates) overall comprising approximately $[eqn]$ population dynamics simulations (plus those performed for additional experiments considering different values of p, different initial networks, and structural stability analyses).

Numerical simulations of ecological dynamics were performed in Fortran, by integrating the equations (equation (4.1)) using the Bulirsch–Stoer algorithm with adaptive step until convergence. This code was wrapped in R [83] to perform the evolutionary dynamics. Analyses of the simulation outputs were performed in R and plots generated with ggplot2 [84].

Network properties

4.3.3.

Full analytical expressions of the network properties used are found in electronic supplementary materials. To assess changes in the state of the system along evolutionary trajectories we quantified total community biomass as the sum of the biomass of all species present in the community alongside a suite of network properties on the mutualistic interaction network. The $[eqn]$ -diversity was estimated calculating the richness (total number of species) and the effective diversity (sometimes termed true diversity), which is the exponential of the Shannon diversity. The connectance was calculated as $[eqn]$ (with $[eqn]$ the total number of interactions in the network and $[eqn]$ and $[eqn]$ the number of species of plants and animals, respectively). Nestedness, the extent to which the set of partners of specialist species are proper subsets of the partners of more generalist ones, was quantified using two related metrics, the ecological overlap $[eqn]$ , which is the mean fraction of mutualistic partners that are shared with other species [9], and the $[eqn]$ metric [85]:

$$[eqn]$$

To enable comparison among networks of different sizes and number of interactions, we used a normalized version of $[eqn]$ proposed by Song et al. [50] (termed $[eqn]$ , see electronic supplementary material, Methods). Heterogeneity in the interactions across species was quantified calculating the standard deviation of the degree of the species in the network ( $[eqn]$ ) and the assortativity ( $[eqn]$ )—a network property that assesses the tendency for nodes (i.e. species) to be connected to others with similar degrees [86]. The degree of a species in the network is simply the number of mutualistic links it possesses.

Quantification of structural stability

4.4.

Effective and critical competition

4.4.1.

The linearization of the system governed by the system of equation (4.1) at equilibrium leads to the definition of the effective competition matrix, $[eqn]$ (see electronic supplementary materials and [9,27] for a more general motivation of the approach). This matrix ( $[eqn]$ ) encapsulates both the competitive interactions between species in the same group, either P or A, and their mutualistic interaction with species in the other group (see electronic supplementary materials). The diagonalization of this matrix allows us to estimate the effective interspecific competition parameter for plant species (symmetric expressions are derived for animal species) as:

$$[eqn]$$

This parameter represents the mean-field competition experienced by a plant species in the system. Comparing $[eqn]$ with the ‘bare’ interspecific competition parameter, $[eqn]$ , allows us to assess the effect of mutualistic interactions on species. For example, if $[eqn]$ mutualism effectively reduces interspecific competition. In subsequent work [16], it was shown that mutualistic interactions reduce the interspecific competition of the system only if $[eqn]$ is smaller than a threshold, termed critical competition, $[eqn]$ :

$$[eqn]$$

where $[eqn]$ , i.e. it is proportional to the mutualistic overlap between species in the same group. A full analytical derivation can be found in [16]. In figure 1, we recovered the $[eqn]$ values obtained in [16] for a set of empirical networks retrieved from the web of life (https://www.web-of-life.es). The relationship between both quantities can be compactly expressed through the relative distance between the effective competition and the critical competition:

$$[eqn]$$

which is positive in the regime in which mutualism favours biodiversity, and negative otherwise.

Analytical prediction of structural stability

4.4.2.

We follow the definition of structural stability proposed in [16], which considers the structural stability of a specific fixed point of a system. This differs from other definitions (e.g. [17]) that estimate the structural stability as the volume of the parameter space leading to feasible systems compatible with a given matrix of interactions, i.e. an estimation of the volume of all possible fixed points. Therefore, the definition we follow requires considering not only the interaction matrix but also the growth rates of all species and their realized abundances.

We model the effect of environmental perturbations of amplitude $[eqn]$ on species growth rates as $[eqn]$ , with $[eqn]$ being a random number drawn from a uniform distribution. We are interested in predicting the critical perturbation $[eqn]$ in which at least one extinction was observed. This is our estimate of the structural stability of the system. In previous work [16], we showed that a good approximation of $[eqn]$ can be obtained with the expression:

$$[eqn]$$

with symmetric equations for animal species. Hence, the system’s critical perturbation is the minimum found across both pools of species $[eqn]$ The term $[eqn]$ represents a natural scale of biodiversity set by the effective competition $[eqn]$ . The quantities $[eqn]$ and $[eqn]$ , quantify the propagation of perturbations and the vulnerability of the system at the unperturbed state, respectively (see electronic supplementary material). Full details of the derivation of this expression are provided in [16].

Numerical quantification of structural stability

4.4.3.

We verified that the analytical prediction of structural stability was consistent with the persistence observed in simulations. To this end, we considered systems at steady states every 25 substitutions throughout the eco-evolutionary trajectories. For each system at a steady state, we performed a random perturbation on the growth rates of amplitude $[eqn]$ , and integrated the population dynamics of the system with the perturbed growth rates until a new steady state is reached, recording the number of observed extinctions. We repeated the same procedure for 50 random realizations of noise, corresponding to the 50 evolutionary replicates, and reported the fraction of simulations in which at least one extinction was observed. We investigated how this fraction changed for different values of $[eqn]$ .