Eco-evolutionary model on spatial graphs reveals how habitat structure affects phenotypic differentiation

Eco-evolutionary model on spatial graphs

We establish an individual-based model (IBM) where individuals are structured over a trait space and a graph representing a landscape. For the sake of simplicity, we consider the case of asexual reproduction and haploid genetics²⁹. Individuals die, reproduce, mutate and migrate in a stochastic fashion, which together results in macroscopic properties. The formulation of the stochastic IBM allows an analytical description of the dynamics at the population level, which links emergent properties to the elementary processes that generate them.

The trait space ({{{{{{{mathcal{X}}}}}}}}subseteq {{mathbb{R}}}^{d}) is continuous and can be split into a neutral trait space ({{{{{{{mathcal{U}}}}}}}}) and an adaptive trait space ({{{{{{{mathcal{S}}}}}}}}). We refer to neutral traits (uin {{{{{{{mathcal{U}}}}}}}}) as traits that are not under selection, in contrast to adaptive traits (sin {{{{{{{mathcal{S}}}}}}}}), which experience selection. The graph denoted by G is composed of a set of vertices {v₁,v₂,…,v_M} that correspond to habitat patches (suitable geographical areas), and a set of edges that constrain the movement of individuals between the habitat patches. We use the original measure of genetic differentiation for quantitative traits Q_ST (standing for Q-statistics) in the case of haploid populations^45,46. We denote the neutral trait value of the kth individual on v_i as ({u}_{k}^{(i)}), the number of individuals on v_i as N⁽ⁱ⁾, the mean neutral trait on v_i as ({overline{u}}^{(i)}), and the mean neutral trait in the metapopulation as (overline{u}). It follows that we quantify neutral differentiation Q_ST,u as

where ({sigma }_{B,u}^{2}={mathbb{E}}[frac{1}{M}{sum }_{i}{left({overline{u}}^{(i)}-overline{u}right)}^{2}]) denotes the expected neutral trait variance between the vertices and ({sigma }_{W,u}^{2}=frac{1}{M}mathop{sum }nolimits_{i}^{M}{mathbb{E}}left[frac{1}{{N}^{(i)}}{sum }_{k}{left({u}_{k}^{(i)}-{overline{u}}^{(i)}right)}^{2}right]) denotes the average expected neutral trait variance within vertices. We similarly quantify adaptive differentiation Q_ST,s.

Following the Gillespie update rule⁴⁷, individuals with trait ({x}_{k}in {{{{{{{mathcal{X}}}}}}}}) on vertex v_i are randomly selected to give birth at rate b⁽ⁱ⁾(x_k) and die at rate d(N⁽ⁱ⁾) = N⁽ⁱ⁾/K, where K is the local carrying capacity. The definition of d therefore captures competition, which is proportional to the number of individuals on a vertex and does not depend on the individuals’ traits (we relax this assumption later on). The offspring resulting from a birth event inherits the parental traits, which can independently be affected by mutations with probability μ. A mutated trait differs from the parental trait by a random change that follows a normal distribution with variance ({sigma }_{mu }^{2}) (corresponding to the continuum of alleles model⁴⁸). The offspring can further migrate to neighbouring vertices by executing a simple random walk on G with probability m. A schematic overview of the two different settings considered is provided in Fig. 1. Under the setting with no selection, individuals are only characterised by neutral traits so that ({{{{{{{mathcal{X}}}}}}}}={{{{{{{mathcal{U}}}}}}}}). For individuals on a vertex with trait x_k ≡ u_k we define b⁽ⁱ⁾(x_k) ≡ b, so that the birth rate is constant. This ensures that neutral traits do not provide any selective advantage. Under the setting with heterogeneous selection, each vertex of the graph v_i is labelled by a habitat type with environmental condition Θ_i that specifies the optimal adaptive trait value on v_i. It follows that, for individuals with traits ({x}_{k}=({u}_{k},{s}_{k})in {{{{{{{mathcal{U}}}}}}}}times {{{{{{{mathcal{S}}}}}}}}) on v_i, we define

where p is the selection strength⁴¹. This ensures that the maximum birth rate on v_i is attained for s_k = Θ_i, which results in a differential advantage that acts as an evolutionary stabilising force. In the following we consider two habitat types denoted by I and II with symmetric environmental conditions θ_I and θ_II, so that Θ_i ∈ {θ_I, θ_II} and θ_II = − θ_I = θ, where θ can be viewed as the habitat heterogeneity⁴¹.

Deterministic approximation of the population dynamics under no selection

The model can be formulated as a measure-valued point process (³⁰ and Supplementary Note). Under this formalism, we demonstrate in the Supplementary Note how the population size and the trait dynamics show a deterministic behaviour when a stabilising force dampens the stochastic fluctuations. This makes it possible to express the dynamics of the macroscopic properties with deterministic differential equations, connecting emergent patterns to the processes that generate them. In particular, in the setting of no selection, competition stabilises the population size fluctuations, and the dynamics can be considered deterministic and expressed as

where (A={({a}_{i,j})}_{1le i,jle M}) is the adjacency matrix of the graph G and D = (d₁,d₂,…,d_M) is a vector containing the degree of each vertex (number of edges incident to the vertex). The first term on the right-hand side corresponds to logistic growth, which accounts for birth and death events of non-migrating individuals. The second term captures the gains due to migrations, which depend on the graph topology. Assuming that all vertices with the same degree have an equivalent position on the graph, corresponding to a “mean field” approach (see Methods), one can obtain a closed-form solution from Eq. (3) (see Eq. (12)), which shows that the average population size (overline{N}) scales with ({langle sqrt{k}rangle }^{2}/langle krangle), where 〈k〉 is the average vertex degree and (langle sqrt{k}rangle) is the average square-rooted vertex degree. The quantity ({langle sqrt{k}rangle }^{2}/langle krangle), denoted as h_d, relates to the homogeneity in vertex degree of the graph and can therefore be viewed as a measure negatively associated with heterogeneity in connectivity. Simulations of the IBM illustrate that h_d can explain differences in population size for complex graph topologies with varying migration regimes (Fig. 2a for graphs with M = 7 vertices and Supplementary Fig. 1a for M = 9). This analytical result is connected to theoretical work on reaction-diffusion processes⁴⁹ and highlights that irregular graphs (graphs whose vertices do not have the same degree) result in unbalanced migration fluxes that affect the ecological balance between births and deaths. Highly connected vertices present an oversaturated carrying capacity (N⁽ⁱ⁾ > bK, see Methods), increasing local competition and lowering total population size compared with regular graphs (Fig. 2a). Because populations with small sizes experience more drift (³¹ and Supplementary Fig. 2), this result indicates that graph topology affects neutral differentiation not only through population isolation, but also by affecting population dynamics.

<div class="c-article-section__figure js-c-reading-companion-figures-item" data-test="figure" data-container-section="figure" id="figure-2" data-title="Effect of and hd on average population size

$$overline{N}$$

N

¯

and neutral differentiation QST,u in the setting with no selection.”>

Nonetheless, the stochasticity of the processes at the individual level can propagate to the population level and substantially affect the macroscopic properties. In particular, neutral differentiation emerges from the stochastic fluctuations of the populations’ neutral trait distribution. These fluctuations complicate an analytical underpinning of the dynamics, and in this case simulations of IBM offer a straightforward approach to evaluate the level of neutral differentiation.

Effect of graph topology on neutral differentiation under no selection

We study a setting with no selection and investigate the effect of the graph topology on neutral differentiation. When migration is limited, individuals’ traits are coherent on each vertex but stochastic drift at the population level generates neutral differentiation between the vertices. Migration attenuates neutral differentiation because it has a correlative effect on local trait distributions. Following^21,22,26, we expect that the intensity of the correlative effect depends on the average path length of the graph 〈l〉, defined as the average shortest path between all pairs of vertices⁵⁰. For a constant number of vertices, 〈l〉 is strictly related to the mean betweenness centrality and quantifies the graph connectivity⁵⁰. High 〈l〉 implies low connectivity and greater isolation of populations, and hence we expect that graphs with high 〈l〉 are associated with high differentiation levels. We consider various graphs with an identical number of vertices and run simulations of the IBM to obtain the neutral differentiation level Q_ST,u attained after a time long enough to discard transient dynamics (see Methods). We then interpret the discrepancies in Q_ST,u across the simulations by relating them to the underlying graph topologies.

We observe strong differences in Q_ST,u across graphs for varying m, and find that 〈l〉 explains at least 55% of the variation in Q_ST,u across all graphs with M = 7 vertices for (Fig. 2b). Nonetheless, some specific graphs, such as the star graph, present higher levels of Q_ST,u than expected by their average path length. To explain this discrepancy, we explore the effect of homogeneity in vertex degree h_d, as we showed in Eq. (12) that it decreases population size, which should in turn increase Q_ST,u by intensifying stochastic drift. We find that h_d explains 57% of the variation for low m (Fig. 2c). However, the fit remains similar after correcting for differences in population size (see Supplementary Table 1), indicating that irregular graphs structurally amplify the isolation of populations. Unbalanced migration fluxes lead central vertices to host more individuals than allowed by their carrying capacity. This causes increased competition that results in a higher death rate, so that migrants have a lower probability of further spreading their trait. Highly connected vertices therefore behave as bottlenecks, increasing the isolation of peripheral vertices and consequently amplifying Q_ST,u.

We then evaluate the concurrent effect of 〈l〉 and h_d on Q_ST,u with a multivariate regression model that we fit independently for low and high migration regimes (Fig. 2d). The multivariate regression model explains at least 70% of the variation in Q_ST,u for the migration regimes considered and for graphs with M = 7 vertices (see Supplementary Table 2 for details). Moreover, we find that 〈l〉 and h_d have akin contributions to neutral differentiation for low m, but the effect of 〈l〉 increases for higher migration regimes while the effect of h_d decreases. To ensure that these conclusions can be generalised to larger graphs, we conduct the same analysis on a subset of graphs with M = 9 vertices and find congruent results (Supplementary Fig. 1). In the absence of selection and with competitive interactions, graphs with a high average path length 〈l〉 and low homogeneity in vertex degree h_d, or similarly graphs with low connectivity and high heterogeneity in connectivity, show high levels of neutral differentiation.

Deterministic approximation of the population dynamics and adaptation under heterogeneous selection

We next consider heterogeneous selection and investigate the response of adaptive differentiation to the spatial distribution of habitat types, denoted as the Θ-spatial distribution. Adaptive differentiation emerges from local adaptation, but migration destabilises adaptation as a result of the influx of maladaptive migrants. We expect that higher connectivity between vertices of similar habitat type increases the level of adaptive differentiation, because it increases the proportion of well-adapted migrants. Local adaptation can be investigated by approximating the stochastic dynamics of the trait distribution with a deterministic partial differential equation (PDE). We demonstrate under mean-field assumption how the deterministic approximation can be reduced to an equivalent two-habitat model. We analyse the reduced model with the theory of adaptive dynamics^36,41 and find a critical migration threshold m^⋆ that determines local adaptation. m^⋆ depends on a quantity coined the habitat assortativity r_Θ, and we demonstrate with numerical simulations that r_Θ determines the overall adaptive differentiation level Q_ST,s reached at steady state in the deterministic approximation.

Heterogeneous selection, captured by the dependence of the birth rate on Θ_i, generates a stabilising force that dampens the stochastic fluctuations of the adaptive trait distribution. The dynamics of the adaptive trait distribution consequently shows a deterministic behavior and we demonstrate in the Supplementary Note and Supplementary Figs. 3 and 4 that the number of individuals on v_i with traits (sin {{Omega }}subset {{{{{{{mathcal{S}}}}}}}}) can be approximated by the quantity ∫_Ωn⁽ⁱ⁾(s)ds, where n⁽ⁱ⁾ is a continuous function solution of the PDE

Equation (4) is similar to Eq. (3), except that it incorporates an additional term corresponding to mutation processes and that the birth rate is trait-dependent. We show how Eq. (4) can be reduced to an equivalent two-habitat model under mean-field assumption. The mean-field approach differs slightly from the setting with no selection because vertices are labelled with Θ_i. Here we assume that vertices with similar habitat types have an equivalent position on the graph (see Supplementary Fig. 5 for a graphical representation), so that all vertices with habitat type I are characterised by the identical adaptive trait distribution that we denote by ({overline{n}}^{{{{{{{{bf{I}}}}}}}}}), and are associated with the birth rate ({b}^{{{{{{{{bf{I}}}}}}}}}(s)=b(1-p{(s-{theta }_{{{{{{{{bf{I}}}}}}}}})}^{2})). Let P(I, II) denote the proportion of edges connecting a vertex v_i of type II to a vertex v_j of type I, and let P(I) denote the proportion of vertices v_i of type I. By further assuming that habitats are homogeneously distributed on the graph so that (P({{{{{{{bf{I}}}}}}}})=P({{{{{{{bf{II}}}}}}}})=frac{1}{2}), Eq. (4) transforms into

(see Methods), where we define

as the habitat assortativity of the graph, which ranges from −1 to 1. When r_Θ = − 1, all edges connect dissimilar habitat types (disassortative graph), while as r_Θ tends towards 1 the graph is composed of two clusters of vertices with identical habitat types (assortative graph). Eq. (5) can be analysed with the theory of adaptive dynamics^36,38,41, a mathematical framework that provides analytical insights by assuming a “trait substitution process”. Following this assumption, the mutation term in Eq. (5) is omitted and the phenotypic distribution results in a collection of discrete individual types that are gradually replaced by others until evolutionary stability is reached (see Methods and^36,38,41 for details). By applying the theory of adaptive dynamics, we find a critical migration rate m^⋆

so that when m > m^⋆, a single type of individual exists with adaptive trait ({s}^{* }=left({theta }_{{{{{{{{bf{II}}}}}}}}}+{theta }_{{{{{{{{bf{I}}}}}}}}}right)/2=0) in the steady-state (see Methods for the derivation of Eq. (7)). In this case, adaptive differentiation Q_ST,s is nil and the average population size is given by (overline{N}=bK{(1-ptheta )}^{2}). In contrast, when m = 0 and/or r_Θ = 1, all individuals are locally well-adapted with trait Θ_i on v_i, and it follows that the average population size is higher and equal to (overline{N}=bK), while adaptive differentiation is maximal and equal to ({Q}_{ST,s}={{{{{{{rm{Var}}}}}}}}({{Theta }})/left({{{{{{{rm{Var}}}}}}}}({{Theta }})+0right)=1). When 0 < m < m^⋆, the coexistence of two types of individuals on each vertex v_i is predicted but the calculation of the trait values is more subtle. To understand the effect of m and r_Θ on the local trait distributions and on Q_ST,s, we therefore leave behind the adaptive dynamics framework and numerically solve Eq. (5) by including the mutation term. When 0 < m < m^⋆, the local trait distributions are bimodal with peaks corresponding to the two types of individuals predicted by the adaptive dynamics. The highest peak corresponds to the well-adapted individuals, whose adaptation is destabilised by the influx of maladaptive migrants (Fig. 3a). This phenomenon is dampened as r_Θ increases, since the proportion of maladaptive migrants is reduced in assortative graphs (Fig. 3b). As a consequence, the habitat assortativity r_Θ increases the differentiation Q_ST,s when 0 < m < m^⋆ (Fig. 3c). The simulations further confirm that the adaptive dynamics prediction given by Eq. (7) is still valid when the continuous accumulation of mutations is considered, so that for m > m^⋆ the local trait distributions obtained from Eq. (5) are unimodal and Q_ST,s vanishes (Fig. 3a,c). Our analysis of the mean-field deterministic approximation Eq. (5) therefore demonstrates that assortative graphs present high levels of adaptive differentiation Q_ST,s. On the other hand, the analysis shows that Q_ST,s rapidly declines with increasing m on disassortative graphs, until Q_ST,s vanishes when m > m^⋆.

Effect of graph topology on adaptive differentiation under heterogeneous selection

To generalise the conclusions drawn from the mean-field deterministic approximation Eq. (5), we generate different Θ-spatial distributions for varying graph topology, and compare outputs of the IBM simulations with those of Eq. (5) (see Methods for the details of the simulations). For each combination of Θ-spatial distribution and graph, we compute the habitat assortativity r_Θ, since r_Θ can be generalised from Eq. (6) to any graph topology following the original definition of⁵¹ as

where Θ_× and Θ_∧ denote the sets of habitats found at the toe and tip of each directed vertex of graph V, and 〈Θ_×〉, 〈Θ_∧〉 and ({sigma }_{{{{Theta }}}_{times }},{sigma }_{{{{Theta }}}_{wedge }}) denote their respective means and standard deviations (see Supplementary Note). The mean-field deterministic approximation Eq. (5) is in very good agreement with the IBM simulations for general graph ensembles at low migration regimes, and captures the response of (overline{N}) and Q_ST,s to r_Θ (Fig. 4). Nonetheless, under high migration regimes, higher levels of Q_ST,s are observed in the stochastic simulations compared with the mean field deterministic approximation (Supplementary Fig. 6). We hypothesize that this reinforcement is generated by stochastic drift, which must become the main driver of differentiation when local adaptation is lost for m > m^⋆, and perform a multivariate regression analysis to investigate the additional effect of 〈l〉 and h_d on Q_ST,s. As expected, the analysis highlights that the effect of 〈l〉 and h_d are substantial and complement the effect of r_θ for high m (Fig. 5c for graphs with M = 7 vertices and Supplementary Fig. 7a for M = 9), further explaining the discrepancies observed (see Supplementary Table 3).

We extend our analyses to realistic landscapes with a continuum of habitat types by running simulations on graphs obtained from real spatial habitat datasets and by considering mean annual temperature as a proxy for habitat type (see Supplementary Fig. 8 and Supplementary Table 4). We also consider simulations accounting for trait-dependent competition to test whether our results hold under more complex ecological processes (see Supplementary Note for the implementation details and Supplementary Table 5 for the results). The simulations are congruent and show that the effects of r_Θ, h_d, and 〈l〉 are similar under these alternative settings, underlining the robustness of these metrics and the generality of our conclusions. Taken together, these results indicate that under sufficiently strong selection and sufficiently high habitat heterogeneity, adaptive differentiation Q_ST,s is mainly driven by habitat assortativity r_Θ. Nonetheless, local adaptation is lost in disassortative graphs when m > m⋆, such that 〈l〉 and h_d become complementary determinants of Q_ST,s for high migration regimes.

Effect of habitat assortativity on neutral differentiation under heterogeneous selection

We finally consider a setting with heterogeneous selection where individuals carry both neutral and adaptive traits. With distinct habitat types, selection promotes neutral differentiation by reducing the birth rate of maladaptive migrants, reinforcing the isolation of local populations. We have shown above that adaptive differentiation Q_ST,s is driven by habitat assortativity r_Θ, so we expect r_Θ, together with the topological metrics found in the setting with no selection, to influence the level of neutral differentiation Q_ST,u. We first investigate how the response of Q_ST,u to migration compares between the setting with no selection and the setting with heterogeneous selection for graphs with an identical topology. We then examine how the response compares between graphs with an identical topology but different r_Θ. We finally consider simulations on different graphs with varying r_Θ to assess the concurrent effect of 〈l〉, h_d, and r_Θ on Q_ST,u.

Migration has a fitness cost because maladaptive migrants present lower fitness. Under an equivalent migration regime, migrants therefore have a lower probability of reproduction, increasing the populations’ isolation compared with a setting without selection. Simulations with varying m on the complete graph confirm that selection in heterogeneous habitats reinforces Q_ST,u compared with a setting without selection (Fig. 5a). Nonetheless, previous results show that adaptive differentiation Q_ST,s vanishes on a disassortative graph when m > m^⋆, implying that individuals become equally fit in all habitats. In this case, the isolation effect of heterogeneous selection is lost and Q_ST,u reaches a similar level as in the setting with no selection for m > m^⋆ (Fig. 5a), although Q_ST,u is slightly higher in the setting with heterogeneous selection due to lower population size ((overline{N}=bK(1-ptheta )) vs. (overline{N}=bK), see section above and Methods). This suggests that r_Θ reinforces Q_ST,u, as assortative graphs sustain higher levels of adaptive differentiation (Figs. 3 and 4). Simulations on the path graph with varying Θ-spatial distribution support this conclusion for high migration regimes, but show the opposite relationship under low migration regimes, where the habitat assortativity r_Θ decreases Q_ST,u (Fig. 5b). Assortative graphs are composed of large clusters of vertices with similar habitats, within which migrants can circulate without fitness losses. Local neutral trait distributions become more correlated within these clusters, resulting in a decline in Q_ST,u for assortative graphs compared with disassortative graphs. Figure 5b therefore highlights the ambivalent effect of r_Θ on Q_ST,u. r_Θ reinforces Q_ST,u by favouring adaptive differentiation, but also decreases Q_ST,u by decreasing population isolation within clusters of vertices with the same habitat type.

We compare the effect of r_Θ on Q_ST,u to the effect of the topology metrics 〈l〉 and h_d found in the setting with no selection using multivariate regression analysis on simulation results obtained for different graphs with varying Θ-spatial distribution (Fig. 5d for graphs with M = 7 vertices and Supplementary Fig. 7b for M = 9). The multivariate model explains the discrepancies in Q_ST,u across the simulations for low and high migration regimes (see Supplementary Table 3 for details), and we find that r_Θ, 〈l〉, and h_d contribute similarly to neutral differentiation. Hence, the effects of r_Θ and the topology metrics 〈l〉 and h_d add up under heterogeneous selection. A change in sign of the standardized effect of r_Θ on Q_ST,s for low and high migration regimes verifies that the ambivalent effect of r_Θ on Q_ST,u found on the path graph holds for general graph ensembles. Simulations with trait-dependent competition and simulations on realistic graphs with a continuum of habitat types equally confirm the ambivalent effect of r_Θ and further support the complementary effect of 〈l〉 and h_d on Q_ST,u (see Supplementary Fig. 8). 〈l〉 and h_d therefore drive neutral differentiation with and without heterogeneous selection. r_Θ becomes an additional determinant of neutral differentiation under heterogeneous selection. In contrast to the non-ambivalent, positive effect of habitat assortativity on adaptive differentiation, r_Θ can amplify or depress neutral differentiation depending on the migration regime considered.

Source: Ecology - nature.com