# Eco-evolutionary model on spatial graphs reveals how habitat

14 min read### 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^{29}. 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 \(\mathcalX\subseteq \mathbbR^d\) is continuous and can be split into a neutral trait space \(\mathcalU\) and an adaptive trait space \(\mathcalS\). We refer to neutral traits \(u\in \mathcalU\) as traits that are not under selection, in contrast to adaptive traits \(s\in \mathcalS\), which experience selection. The graph denoted by *G* is composed of a set of vertices *v*_{1},*v*_{2},…,*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 *k*th individual on *v*_{i} as \(u_k^(i)\), the number of individuals on *v*_{i} as *N*^{(i)}, the mean neutral trait on *v*_{i} as \(\overlineu^(i)\), and the mean neutral trait in the metapopulation as \(\overlineu\). It follows that we quantify neutral differentiation *Q*_{ST,u} as

$$Q_ST,u=\sigma _B,u^2/(\sigma _B,u^2+\sigma _W,u^2)$$

(1)

where \(\sigma _B,u^2=\mathbbE[\frac1M\sum _i\left(\overlineu^(i)-\overlineu\right)^2]\) denotes the expected neutral trait variance between the vertices and \(\sigma _W,u^2=\frac1M\mathop\sum \nolimits_i^M\mathbbE\left[\frac1N^(i)\sum _k\left(u_k^(i)-\overlineu^(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^{47}, individuals with trait \(x_k\in {{{{{\mathcalX}}}}}\) on vertex *v*_{i} are randomly selected to give birth at rate *b*^{(i)}(*x*_{k}) and die at rate *d*(*N*^{(i)}) = *N*^{(i)}/*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^{48}). 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 \({{{{{{{\mathcalX}}}}}}}={{{{{{{\mathcalU}}}}}}}\). For individuals on a vertex with trait *x*_{k} ≡ *u*_{k} we define *b*^{(i)}(*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 {{{{{{{\mathcalU}}}}}}}\times {{\mathcalS}}\) on *v*_{i}, we define

$$b^(i)(x_k)\equiv b^(i)(s_k)=b(1-p(s_k-\Theta _i)^2)$$

(2)

where *p* is the selection strength^{41}. 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^{41}.

### Deterministic approximation of the population dynamics under no selection

The model can be formulated as a measure-valued point process (^{30} 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

$$\partial _tN_t^(i)=N_t^(i)\left[b(1-m)-\fracN_t^(i)K\right]+mb\mathop\sum\limits_j\ne i\fraca_i,jd_jN_t^(j)$$

(3)

where \(A=(a_i,j)_1\le i,j\le M\) is the adjacency matrix of the graph *G* and *D* = (*d*_{1},*d*_{2},…,*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 \(\overlineN\) scales with \(\langle \sqrtk\rangle ^2/\langle k\rangle\), where 〈*k*〉 is the average vertex degree and \(\langle \sqrtk\rangle\) is the average square-rooted vertex degree. The quantity \(\langle \sqrtk\rangle ^2/\langle k\rangle\), 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^{49} 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*^{(i)} > *b**K*, see Methods), increasing local competition and lowering total population size compared with regular graphs (Fig. 2a). Because populations with small sizes experience more drift (^{31} and Supplementary Fig. 2), this result indicates that graph topology affects neutral differentiation not only through population isolation, but also by affecting population dynamics.

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^{50}. For a constant number of vertices, 〈*l*〉 is strictly related to the mean betweenness centrality and quantifies the graph connectivity^{50}. 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 \(s\in \Omega \subset {{{{{{{\mathcalS}}}}}}}\) can be approximated by the quantity ∫_{Ω}*n*^{(i)}(*s*)*d**s*, where *n*^{(i)} is a continuous function solution of the PDE

$$\partial _tn_t^(i)(s)= \, n_t^(i)(s)\left[b^(i)(s)(1-m)-\frac1K\int_{{{{{{{{\mathcalS}}}}}}}}n_t^(i)(\bfs)d{{{\bfs}}}\right]\\ +m\mathop\sum\limits_j\ne ib_j(s)\fraca_i,jd_jn_t^(j)(s)+\frac12\mu \sigma _\mu ^2\Delta _s\left[b^(i)(s)n_t^(i)(s)\right]$$

(4)

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 \(\overlinen^\bfI\), and are associated with the birth rate \(b^\bfI(s)=b(1-p(s-\theta _\bfI)^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(\bfI)=P(\bfII)=\frac12\), Eq. (4) transforms into

$$\partial _t\overlinen_t^\bfI(s)= \,\overlinen_t^\bfI(s)\left[b^{{{{\bfI}}}}(s)(1-m)-\frac1K\int_{{{{{{{{\mathcalS}}}}}}}}\overlinen_t^{{{{{{{{\bfI}}}}}}}}({{{{{{{\bfs}}}}}}})d{{{{{{{\bfs}}}}}}}\right]+\frac12\mu \sigma _\mu ^2(\Delta _sb^{{{{{{{{\bfI}}}}}}}}\overlinen_t^{{{{{{{{\bfI}}}}}}}})(s)\\ +\fracm2\,[(1-r_\Theta )b^{\bfII}(s)\overlinen_t^{{{{{{{{\bfII}}}}}}}}(s)+(1+r_\Theta )b^{{{{{{{{\bfI}}}}}}}}(s)\overlinen_t^{{{{{{{{\bfI}}}}}}}}(t)]$$

(5)

(see Methods), where we define

$$r_\Theta =2\left(P({{{{{{{\bfI}}}}}}},{{{{{{{\bfI}}}}}}})-P({{{{{{{\bfI}}}}}}},{{{{{{{\bfII}}}}}}})\right)$$

(6)

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*^{⋆}

$$m^\star =\frac1(1-r_\Theta )\frac4p\theta ^2(1+3p\theta ^2)$$

(7)

so that when *m* > *m*^{⋆}, a single type of individual exists with adaptive trait \(s^* =\left(\theta _{{{{{{{{\bfII}}}}}}}}+\theta _{{{{{{{{\bfI}}}}}}}}\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 \(\overlineN=bK(1-p\theta )^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 \(\overlineN=bK\), while adaptive differentiation is maximal and equal to \(Q_ST,s={\rmVar}(\Theta )/\left({{{{{{{\rmVar}}}}}}}(\Theta )+0\right)=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^{51} as

$$r_\Theta =\frac{{{{{\rmCov}}}}(\Theta _\times ,\Theta _\wedge )}{\sigma _{\Theta _\times }\sigma _{\Theta _\wedge }}$$

(8)

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 \(\overlineN\) 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 (\(\overlineN=bK(1-p\theta )\) vs. \(\overlineN=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.