# Adaptive diversification in heterogeneous environments

Olivier David, Christian Lannou, Herve Monod, Julien Papaix, Djidi Traore

► To cite this version:

Olivier David, Christian Lannou, Herve Monod, Julien Papaix, Djidi Traore. Adaptive diversification in heterogeneous environments. Theoretical Population Biology, 2017, 114, pp.1-9.
10.1016/j.tpb.2016.11.003. hal-01508877

HAL Id: hal-01508877

https://hal.science/hal-01508877

Submitted on 17 Nov 2021

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

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

<figure><img src="image_1.png" /></figure>
Adaptive diversification in heterogeneous environments

Olivier David¹*, Christian Lannou², Hervé Monod¹, Julien Papaïx³, Djidi Traore¹

1 MaIAGE, INRA, Université Paris-Saclay, 78350 Jouy-en-Josas, France

2 UMR 1290 BIOGER, INRA, Thiverval Grignon, France

3 UR 546 BioSP, INRA, Avignon, France

* Correspondence: O. David, Unité MaIAGE, Bât. 210, INRA, Domaine de Vilvert, 78352 Jouy-en-Josas Cedex, France

E-mail: Olivier.David@inrae.fr

Abbreviations: ES (evolutionary singularity).

This article is an author version of the article: Olivier David, Christian Lannou, Hervé Monod, Julien Papaïx, Djidi Traore, Adaptive diversification in heterogeneous environments, Theoretical Population Biology, volume 114, 2017, pages 1-9.

https://doi.org/10.1016/j.tpb.2016.11.003

©2017. This manuscript version is made available under the CC-BY-NC-ND 4.0 license

http://creativecommons.org/licenses/by-nc-nd/4.0/
Abstract

The role of environmental heterogeneity in the evolution of biological diversity
has been studied only for simple types of heterogeneities and dispersals. This arti-
cle broadens previous results by considering heterogeneities and dispersals that are
structured by several environmental factors. It studies the evolution of a metapop-
ulation, living in a network of patches connected by dispersal, under the effects
of mutation, selection and migration. First it is assumed that patches are equally
connected and that they carry habitats characterized by several factors exerting
selection pressures on several individual traits. Habitat factors may vary in the
environment independently or they may be correlated. It is shown that correlations
between habitat factors promote adaptive diversification and that this effect may
be modified by trait interactions on survival. Then it is assumed that patches are
structured by two crossed factors, called the row and column factors, such that
patches are more connected when they occur in the same row or in the same col-
umn. Environmental patterns in which each habitat appears in each row the same
number of times and appears in each column the same number of times are found
to hinder adaptive diversification.

*Keywords:* Adaptation; Adaptive dynamics; Migration; Multivariate evolution; Spatial heterogeneity; Structured populations.
# 1 Introduction

Adaptive diversification is the evolution of different types of organisms under the effects of evolutionary forces such as selection and mutation (Doebeli, 2011). In natural populations, it may play a role in the evolution of biological diversity (Kocher, 2004; Johnson, 2007). In agrosystems, it may allow some genotypes to have higher performances in some environments (Rhoné et al., 2008; Gautier et al., 2009) but also some pathogens to better develop on some hosts and possibly to be more damaging (Pariaud et al., 2009; REX Consortium, 2013; Fabre et al., 2015). In medicine, it may be involved in the emergence of drug resistance (REX Consortium, 2013; Roemhild et al., 2015). It is therefore important to identify the circumstances that promote diversification.

Biotic and abiotic conditions are often variable within environments. Various features of this environmental heterogeneity may influence adaptive diversification. The overall level of heterogeneity generally promotes diversification, although intermediate levels sometimes maximize diversification propensity (Meszéna et al., 1997; Geritz et al., 1998; Doebeli and Dieckmann, 2003; Débarre et al., 2013; Haller et al., 2013). When geographical sites are equally connected by dispersal, balancing the frequencies of the habitats resulting from environmental heterogeneity promotes diversification (Geritz et al., 1998). When dispersal is more important between neighbouring sites, gradients and landscapes isolating populations into large spatially continuous habitats promote diversification (Débarre and Gandon, 2010; Birand et al., 2012; Haller et al., 2013; Papaïx et al., 2013). When dispersal is more important within groups clustering patches than between groups, diversification is hindered by making within-group habitat distributions equal (Papaïx et al., 2013).

Most theoretical works on the role of environmental heterogeneity in adaptive diversification have studied the influence of a single environmental factor. However the complexity of environments makes it likely that many factors are involved in adaptive diversification (Poisot et al., 2011; Laughlin and Messier, 2015). For example such factors could be related to temperature, rainfall, soil type, host genotype etc. Gavrilets and Vose (2005) studied how populations adapted to their environment when habitats were characterized by several environmental factors and were assigned to patches with equal probabilities. However habitat factors may be associated in a landscape in various ways: they may vary independently between patches or they may be correlated. An exploration of the consequences of the relationships between habitat factors is needed.

Most theoretical studies of the diversifying effects of environmental heterogeneity have considered simple models of dispersal, in particular models in which dispersal rates are constant, decay with spatial separation or depend on an environmental factor. However dispersal may be more complex in reality (Karlin, 1982). For example, plant pathogens may be dispersed by air, rain, water, soil, or by vectors such as animals, pollen, microbes, people, machinery (West, 2014). The spread of human diseases may depend on the clustering of hosts into species or age classes (Sloan et al., 2011; Johnson et al., 2015).
Animal diseases may spread between farms at a regional scale through animal movements due to animal trade (Beaunée et al., 2015). An exploration of the interactions between environmental heterogeneity and various forms of dispersal would be welcome.

In this article, we study how environmental heterogeneity promotes the gradual evolution of polymorphism when this heterogeneity and dispersal depend on several factors. We consider a population of individuals living on a network of patches holding different habitats (Section 2). We restrict ourselves to models that are sufficiently simple to be analyzed analytically. We study how the allocation of habitats to patches influences diversification when habitats are characterized by several factors and when patches are equally connected by dispersal (Section 3). Then we study how it influences diversification when dispersal depends on two crossed factors (Section 4).

# 2 Model

## 2.1 Individuals and environment

Individuals are characterized by a vector $x = (x_1, \dots, x_d)^T$ of $d \ge 1$ continuous adaptive traits, where $T$ denotes transposition. They are assumed to reproduce asexually and to have the same trait vector $x$ as their parent, unless mutations occur.

The environment comprises $P$ patches, each holding $K$ individuals, with $K$ sufficiently large that demographic stochasticity can be neglected. Each patch belongs to one habitat among $H$ possible ones. Each habitat is characterized by a vector $\beta = (\beta_1, \dots, \beta_d)^T$, that corresponds to the optimal phenotype for this habitat. As will be detailed below, the growth of a juvenile in this habitat is optimal if $x = \beta$. The allocation of habitats to patches is called an environmental pattern.

## 2.2 Life cycle

*Juvenile dispersal.* Juveniles are assumed to disperse after birth. A proportion $m_{ii'}$ of juveniles move from patch $i'$ to patch $i$ during a life cycle, with $\sum_i m_{ii'} = 1$ and $m_{ii'} = m_{i'i}$. All patches are assumed to be connected.

*Growth and selection.* After dispersal, a juvenile either develops into an adult or dies. The proportion $g(x, \beta)$ of juveniles with trait vector $x$ in a habitat selecting for phenotype $\beta$ that become adults is assumed to have a $d$-dimensional normal shape (Geritz et al., 1998; Doebeli and Ispolatov, 2010):

$$g(x, \beta) = \alpha \exp \left( -(x - \beta)^T B (x - \beta) / 2 \right),$$

with $d \times d$ inverse variance-covariance matrix $B$, that is symmetric and positive definite. An individual is more adapted to a habitat when its phenotype is closer to the optimal phenotype of this habitat. The normal shape of $g(x, \beta)$ imposes a trade-off between the
growth in the habitats: adaptation to a habitat causes maladaptation to the other habitats. The non-diagonal elements of B quantify the interactions between traits on juvenile development. When B is diagonal, ln(g(x, β)) is simply a sum of trait contributions:

$$ \ln(g(x, \beta)) = \ln(\alpha) - \sum_{j=1}^{d} B_{jj}(x_j - \beta_j)^2/2, $$

where $B_{jj'}$ is the element of B in row j and column j'. When B has some non-zero non-diagonal elements, ln(g(x, β)) also involves contributions from pairs of traits, i.e. the sum $-\sum_{j<j'} B_{jj'}(x_j - \beta_j)(x_{j'} - \beta_{j'})$, so that traits interact.

*Density regulation.* After growth, density dependence makes the size of each patch equal to K.

*Reproduction.* Finally, each adult gives birth to f juveniles before dying. Thus, generations are non-overlapping.

No other assumptions are done at this stage but later when we address some specific issues.

## 2.3 Adaptive dynamics

We use the adaptive dynamics framework (Geritz et al., 1998; Diekmann, 2004). Let us consider a monomorphic resident population, i.e. a population in which all the individuals have the same trait vector $x$. The population evolves thanks to the recurrent fixation of mutants. Mutations are assumed to be rare so that the population has time to reach its demographic equilibrium between mutation events. They are also assumed to have small effects so that evolution is gradual.

Adaptive dynamics is based on the definition of an invasion fitness that indicates if a mutant can invade a resident population. To define fitness, the initial demography of a mutant population is approximated by a matrix model (Appendix A). Invasion fitness is then defined as the dominant eigenvalue of the projection matrix. The mutant is assumed to replace the resident when its fitness is larger than one and to go extinct when its fitness is smaller than one.

Evolutionary singularities (ESs) are evolutionary equilibria. It is shown in Appendix D that our model has a unique ES that is equal to:

$$ \hat{x} = \sum_{h} p_{1...d}(h) \beta(h), $$

where $p_{1...d}(h)$ is the frequency of habitat $h$, i.e. the proportion of patches with this habitat, and $\beta(h)$ is the optimal phenotype of habitat $h$. This ES depends on the environmental pattern through habitat frequencies only.

An ES $\hat{x}$ is attracting (or convergence stable) if the resident trait converges to $\hat{x}$ by gradual evolutionary changes. It is shown in Appendix E that $\hat{x}$ is an attractor. In the
following sections, we study how environmental heterogeneity destabilizes $\hat{x}$ when habitats and dispersal are patterned by factors.

## 3 Diversification in multi-dimensional environments

*Multi-dimensional habitats.* Environments are likely to exert selection pressures through several variables in nature. Thus, we assume that habitats are multi-dimensional in the sense that they are characterized by *d* factors. Here each factor has two levels, so that there are $2^d$ possible habitats, that are denoted by the row vectors *h* = (*h*₁, ..., *h*₉), where *h*ᵢ = 1, 2 is the level of the *j*th factor for *j* = 1, ..., *d*. The optimal phenotype of habitat *h* is assumed to be equal to $\beta(h) = ((-1)^{h_1}, ..., (-1)^{h_d})^T \theta$, where $\theta$ is a positive scalar. The *j*th factor influences the *j*th component of $\beta(h)$ and thus exerts a selection pressure on the *j*th trait, but it may also affect other traits through trait interactions. The discrepancy between habitats *h* and *h'* can be quantified by the distance:

$$ \sqrt{(\beta(h) - \beta(h'))^T B (\beta(h) - \beta(h'))}. \quad (1) $$

The more distant habitats *h* and *h'*, the stronger the selection pressure exerted on an individual adapted to habitat *h* when it is in habitat *h'*. For example when *d* = 2, the distance between habitats (1, 1) and (2, 2) (respectively (1, 2) and (2, 1)) is equal to $2\theta\sqrt{B_{11} + 2B_{12} + B_{22}}$ (respectively $2\theta\sqrt{B_{11} - 2B_{12} + B_{22}}$). A positive interaction ($B_{12} > 0$) increases (respectively decreases) the distance between habitats (1, 1) and (2, 2) (respectively (1, 2) and (2, 1)). A negative interaction has the opposite effect.

*Juvenile dispersal.* Dispersal parameters are assumed to be equal to $m_{ii'} = \mu_1 + \mu_2$ if *i* = *i'* (philopatry) and $m_{ii'} = \mu_1$ otherwise (dispersal), where $0 < \mu_1 \le 1$ to ensure that all patches are connected and $0 \le \mu_2 \le 1$. Thus, patches are assumed to be equally connected and a juvenile has more chance to remain in its birth patch than to move to another patch when $\mu_2 > 0$. This dispersal is called a Deakin's dispersal in this article (Karlin, 1982). As patches are equally connected, the spatial distribution of habitats is not relevant and an environmental pattern can be defined by its habitat frequencies only.

*Factor correlation.* Habitat factors may vary in the environment independently or they may be correlated. Factors *j* and *j'* are independent when they satisfy:

$$ p_{jj'}(l, l') = p_j(l) p_{j'}(l'), \quad (2) $$

for *l*, *l'* = 1, 2, where $p_j(l)$ (respectively $p_{jj'}(l, l')$) denotes the frequency of level *l* of factor *j* (respectively pair of levels (*l*, *l'*)) of factors *j* and *j'*), i.e. the proportion patches with this level (respectively pair of levels). For example, factors are independent for patterns $D_{(0)}$ in Tables 1 and 2. The dependence between factors *j* and *j'* is quantified through the correlation coefficient:

$$ \rho_{jj'} = \frac{p_{jj'}(1, 1) - p_j p_{j'}}{\sqrt{p_j (1 - p_j) p_{j'} (1 - p_{j'})}}, $$
<table><caption>Table 1: Frequencies of two-dimensional habitats. There are four habitats that are denoted by row vectors ($h_1, h_2$), $1 \le h_1, h_2 \le 2$. Patterns $D_{(1)}$ and $D_{(-1)}$ involve a fraction of these habitats whereas $D_{(0)}$ involves these habitats in equal frequency. For each habitat factor, the frequency of level 1 is 1/2. The subscript of a pattern denotes the value of factor correlation.</caption><thead><tr><th>Habitat</th><th>$D_{(-1)}$</th><th>$D_{(-1/3)}$</th><th>$D_{(0)}$</th><th>$D_{(1)}$</th></tr></thead><tbody><tr><td>(1, 1)</td><td>0</td><td>1/6</td><td>1/4</td><td>1/2</td></tr><tr><td>(1, 2)</td><td>1/2</td><td>2/6</td><td>1/4</td><td>0</td></tr><tr><td>(2, 1)</td><td>1/2</td><td>2/6</td><td>1/4</td><td>0</td></tr><tr><td>(2, 2)</td><td>0</td><td>1/6</td><td>1/4</td><td>1/2</td></tr></tbody></table>

for $0 < p_j, p_{j'} < 1$, where $p_j$ denotes $p_j(1)$. The frequency of a pair of levels can be written as follows:

$$p_{jj'}(l, l') = p_j(l) p_{j'}(l') + (-1)^{l+l'} \rho_{jj'} \sqrt{p_j(1-p_j) p_{j'}(1-p_{j'})}, \quad (3)$$

for $1 \le j \ne j' \le d$ and $l, l' = 1, 2$. According to (3), factors $j$ and $j'$ are independent iff $\rho_{jj'} = 0$. A positive value of $\rho_{jj'}$ increases (respectively decreases) the frequencies of the pairs of levels (1, 1) and (2, 2) for factors $j$ and $j'$ (respectively (1, 2) and (2, 1)). A negative value has the opposite effect. When factor correlations $\rho_{jj'}$ are all equal, their common value is denoted by $\rho$: $\rho_{jj'} = \rho$ for $1 \le j \ne j' \le d$. An environmental pattern such that factor correlation $\rho_{jj'}$ is equal to $\rho$ whatever the pair of factors is denoted by $D_{(\rho)}$. For example, factor correlations are all equal to 0 for pattern $D_{(0)}$ in Table 2. Several patterns may have the same value of $\rho$. Pattern $D_{(\rho)}$ may not exist for given values of $\rho$, $d$ and $p_j$.

*Example.* In plant epidemiology, a patch may be an agricultural plot and a habitat may be the variety grown in a plot. A plant leaf may be divided into healthy sites that can be colonized by a pathogen (Papaïx et al., 2014). It may be assumed that $K$ healthy sites become infected in each plot during a pathogen life cycle, while the sites infected during the previous life cycle become removed, thereby making patches of size $K$ and generations non-overlapping. A habitat factor may be a plant locus, whose levels are different alleles, and in this case factor correlation $\rho_{jj'}$ corresponds to linkage disequilibrium (Hartl and Clark, 1997). Our framework could give insight into how the deployment of some varieties in an agricultural landscape influences pathogen adaptation to these varieties (Papaïx et al., 2011; REX Consortium, 2013; Fabre et al., 2015).

## 3.1 Influence of factor correlations

The ES is uninvadable if no nearby mutant can invade a monomorphic resident population with trait $\hat{x}$. In this case, selection is stabilizing at $\hat{x}$ and if $\hat{x}$ is also attracting, the
<table><caption>Table 2: Frequencies of four-dimensional habitats. There are eight habitats that are denoted by row vectors ($h_1, \dots, h_4$), $1 \le h_1, \dots, h_4 \le 2$. Patterns $D_{(1)}$ and $D_{(-1/3)}$ involve a fraction of these habitats whereas $D_{(0)}$ involves these habitats in equal frequency. For each factor, the frequency of level 1 is 1/2. The subscript of a pattern denotes the value of factor correlation.</caption><thead><tr><th>Habitat</th><th>$D_{(-1/3)}$</th><th>$D_{(0)}$</th><th>$D_{(1)}$</th></tr></thead><tbody><tr><td>(1, 1, 1, 1)</td><td>0</td><td>1/8</td><td>1/2</td></tr><tr><td>(1, 1, 2, 2)</td><td>1/6</td><td>1/8</td><td>0</td></tr><tr><td>(1, 2, 1, 2)</td><td>1/6</td><td>1/8</td><td>0</td></tr><tr><td>(1, 2, 2, 1)</td><td>1/6</td><td>1/8</td><td>0</td></tr><tr><td>(2, 1, 1, 2)</td><td>1/6</td><td>1/8</td><td>0</td></tr><tr><td>(2, 1, 2, 1)</td><td>1/6</td><td>1/8</td><td>0</td></tr><tr><td>(2, 2, 1, 1)</td><td>1/6</td><td>1/8</td><td>0</td></tr><tr><td>(2, 2, 2, 2)</td><td>0</td><td>1/8</td><td>1/2</td></tr></tbody></table>

population evolves into a monomorphic population with trait $\hat{x}$. We first study how factor correlations influence ES invadability when there are no trait interactions on juvenile development. It is shown in Appendix F that $\hat{x}$ is then uninvadable if the matrix with jth diagonal element:

$$4 \left( 1 + 2 \frac{\mu_2}{1 - \mu_2} \right) p_j (1 - p_j) B_{jj}^2 \theta^2 - B_{jj} \quad (4)$$

and with element in row j and column $j'$ ($j \neq j'$):

$$4 \left( 1 + 2 \frac{\mu_2}{1 - \mu_2} \right) B_{jj} B_{j'j'} \theta^2 \rho_{jj'} \sqrt{p_j (1 - p_j) p_{j'} (1 - p_{j'})} \quad (5)$$

is negative definite. Correlations between factors introduce non-diagonal elements in the matrix defined by (4) and (5), which tends to increase the eigenvalues of this matrix and to destabilize evolutionary equilibria (Débarre et al., 2014; Svardal et al., 2014).

To study the influence of factor correlations on diversification further, we focus on the case when the frequency of level 1 ($p_j$) is equal to $p$ whatever the factor, factor correlations ($\rho_{jj'}$) are all equal to $\rho$ and the diagonal elements of $B$ are all equal to $\tau$. It is shown in Appendix F that $\hat{x}$ is then uninvadable if:

$$4 \left( 1 + 2 \frac{\mu_2}{1 - \mu_2} \right) p(1-p) \tau \theta^2 \max(1-\rho, 1+(d-1)\rho) - 1 < 0. \quad (6)$$

This condition is the same as that when a single trait evolves ($d = 1, \rho = 0$) iff $\rho = 0$. Thus, although traits do not interact, they do not evolve in isolation from one another
unless $\rho = 0$. When $d = 2$, factor correlation increases the left-hand side of (6) whatever its sign ($\gamma = 0$, Fig. 1). Increasing (respectively decreasing) $\rho$ increases the frequency of habitats (1, 1) and (2, 2) (respectively (1, 2) and (2, 1)), which strengthens the selection exerted by these habitats, destabilizing $\hat{x}$. When $d = 4$, factor correlation has the same effect on invadability when $\rho = -1/3$ but has a stronger effect when $\rho = 1$ ($\gamma = 0$, Fig. 2).

## 3.2 Joint influence of factor correlations and trait interactions

The traits of an organism may show interactions in nature (Hartl and Clark, 1997; Barton et al., 2007). Thus, we now let traits interact, but we focus on the case when $B$ is completely symmetric to simplify, i.e. its diagonal elements are equal to $\tau > 0$ and its non-diagonal elements are equal to $\tau\gamma$. The parameter $\gamma$ quantifies trait interactions and satisfies $-1/(d-1) < \gamma < 1$ to ensure that $B$ is positive definite, which introduces a dissymmetry in its domain when $d > 2$. We also assume that the frequencies of level 1 are all equal to $p$ and that factor correlations are all equal to $\rho$ to simplify. It is shown in Appendix F that $\hat{x}$ is then uninvadable if:

$$4 \left( 1 + 2 \frac{\mu_2}{1 - \mu_2} \right) p(1-p)\tau\theta^2 \max\left( (1-\gamma)(1-\rho), (1+(d-1)\gamma)(1+(d-1)\rho) \right) - 1 < 0. \quad (7)$$

Thus, trait interactions and factor correlations jointly affect invadability. When $d = 2$, increasing $\gamma$ increases the left-hand side of (7) when $\rho = 1$ but have the opposite effect when $\rho = -1$ (Fig. 1). This is because increasing (respectively decreasing) both $\rho$ and $\gamma$ increases the frequency of and the distance between habitats (1, 1) and (2, 2) (respectively (1, 2) and (2, 1)), which strengthens selection, destabilizing $\hat{x}$. Conversely increasing $\rho$ while decreasing $\gamma$ increases the frequency of habitats (1, 1) and (2, 2) but decreases the distance between these habitats, which weakens selection, stabilizing $\hat{x}$. Thus, trait interactions and factor correlations may jointly stabilize the ES (negative invadability criterion). When $d = 4$, the effect of trait interactions is stronger for $\rho = 1$ and is monotonous for $\rho = -1/3$ (Fig. 2).

It can be verified that pattern $D_{(\rho)}$ with $\rho = -\frac{\gamma}{1+(d-2)\gamma}$, when it exists, minimizes the left-hand side of (7) for a given value of $p$ and thus hinders diversification. For example patterns $D_{(0)}$ in Tables 1 and 2, that balance habitat frequencies, limit diversification when trait interactions are low (Fig. 1 and 2). Alternatively when traits show a moderate positive interaction and $d = 2$, pattern $D_{(-1/3)}$ in Table 1 limits diversification. Pattern $D_{(1)}$ maximizes the left-hand side of (7) for a given value of $p$ when $\gamma \ge \frac{2-d}{2-d+d^2}$ and thus promotes diversification. This result is in agreement with intuition since $D_{(1)}$ maximizes the frequencies of habitats (1, ..., 1) and (2, ..., 2), and a large $\gamma$ value increases the distance between these habitats. When $d = 2$, $D_{(-1)}$ maximizes the left-hand side of (7) for a given value of $p$ when $\gamma \le 0$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Influence of trait interaction ($\gamma$) and of factor correlation ($\rho$) on the invadability of the ES when two traits are under selection ($d = 2$). The figure shows the invadability criterion (left-hand side of (7)) for the patterns of Table 1. The horizontal line locates the threshold value of 0: the ES is invadable when the invadability criterion is positive and uninvadable when the invadability criterion is negative. Parameter values were chosen so that the ES is on the brink of invadability when a single trait is under selection (the left-hand side of (7) vanishes when $\gamma = \rho = 0$). Parameter values: $p = 1/2$, $P\mu_1 = \mu_2 = 1/2$, $\tau\theta^2 = 1/3$.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Influence of trait interaction ($\gamma$) and of factor correlation ($\rho$) on the invadability of the ES when four traits are under selection ($d = 4$). The figure shows the invadability criterion (left-hand side of (7)) for the patterns of Table 2. The horizontal line locates the threshold value of 0: the ES is invadable when the invadability criterion is positive and uninvadable when the invadability criterion is negative. Parameter values were chosen so that the ES is on the brink of invadability when a single trait is under selection (the left-hand side of (7) vanishes when $\gamma = \rho = 0$). The axes have the same scale as in Fig. 1 to facilitate comparison. The parameter $\gamma$ is larger than -1/3 to ensure that B is positive definite. Parameter values: $p = 1/2$, $P\mu_1 = \mu_2 = 1/2$, $\tau\theta^2 = 1/3$.</figcaption></figure>
# 4 Diversification with a structured dispersal

*Juvenile dispersal.* Patches are not always equally connected in the real world (Karlin, 1982). To gain insight into the effects of a nonuniform dispersal, we now assume that patches are structured by two crossed factors affecting dispersal, that are called the row and column factors. It is assumed that there is a single patch at the intersection of a row and a column, thus patches can be arranged in a lattice with *R* rows and *C* columns (*P* = *RC*). Dispersal parameters are assumed to be equal to:

$m_{ii'} = \mu_1 + \mu_3$ if *i* and *i'* are different patches in the same row,

$m_{ii'} = \mu_1 + \mu_4$ if *i* and *i'* are different patches in the same column,

$m_{ii'} = \mu_1 + \mu_2 + \mu_3 + \mu_4$ if *i* = *i'* (philopatry),

$m_{ii'} = \mu_1$ otherwise,

where $0 \le \mu_1, \mu_2, \mu_3, \mu_4 \le 1$ and $\mu_1 > 0$ or $\mu_3 \mu_4 > 0$ to ensure that all patches are connected. When $\mu_3 > 0$ and $\mu_4 > 0$, a juvenile that leaves its birth patch has more chance to move to a patch in the same row or in the same column. This dispersal is called a row-column dispersal in this article.

For example, individuals could disperse according to independent between-row and between-column Deakin's dispersals with parameters $\mu_1^{(r)}$ and $\mu_2^{(r)}$ (respectively $\mu_1^{(c)}$ and $\mu_2^{(c)}$) for the between-row (respectively between-column) dispersal. It can be verified that the resulting dispersal is a row-column dispersal with parameters $\mu_1 = \mu_1^{(r)} \mu_1^{(c)}$, $\mu_2 = \mu_2^{(r)} \mu_2^{(c)}$, $\mu_3 = \mu_2^{(r)} \mu_1^{(c)}$ and $\mu_4 = \mu_1^{(r)} \mu_2^{(c)}$. It also corresponds to a Kronecker product dispersal, as the matrix of dispersal parameters can be written as the Kronecker product of a between-row dispersal matrix and a between-column dispersal matrix (Karlin, 1982).

*Habitats.* We assume that habitats are one-dimensional and that a single trait evolves (*d* = 1) to simplify. Two habitats are present in the environment, their optimal phenotypes are equal to $\pm\theta$ and *B* is a scalar. We denote by *p* the frequency of habitat 1.

*Example.* In human epidemiology, some pathogens may spread more easily between persons in the same geographical region and in the same age class, as individuals may have more contacts within such classes (Karlin, 1982; Watts et al., 2005; Wallinga et al., 2006; Sloan et al., 2011; Liccardo and Fierro, 2013). Thus, a patch may be the combination of a geographical region and an age class. As before, it may be assumed that *K* individuals become infected in each patch during a pathogen life cycle while the individuals infected during the previous life cycle become removed. The habitats may be different medicines. Our framework could give insight into how treatment strategies influence pathogen adaptation to some medicines (REX Consortium, 2013; Roemhild et al., 2015).

*Evolutionary branching.* An attracting and invadable ES is a branching point when some phenotypes on either side of the ES can mutually invade each other (mutual invadability) and when the distance between these phenotypes increases over time (evolutionary divergence) (Geritz et al., 1998). We compare pattern $D_1$, in which patches in the same row have the same habitat, and pattern $D_2$, in which each habitat appears in each row
<table><caption>Table 3: Environmental patterns for a population living on a lattice with four rows and four columns. The habitats are denoted by 1 and 2. Pattern $D_1$ is rowwise whereas pattern $D_2$ is balanced in the sense that habitat 1 appears in each row and in each column once and habitat 2 appears in each row and in each column thrice. The frequency of habitat 1 is 1/4.</caption><thead><tr><th colspan="4">$D_1$</th><th colspan="4">$D_2$</th></tr></thead><tbody><tr><td>2</td><td>2</td><td>2</td><td>2</td><td>1</td><td>2</td><td>2</td><td>2</td></tr><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>2</td><td>1</td><td>2</td><td>2</td></tr><tr><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>1</td><td>2</td></tr><tr><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>2</td><td>1</td></tr></tbody></table>

the same number of times and appears in each column the same number of times, for a given value of $p$ (Table 3). In Appendix G, it is shown that the ES is a branching point for the rowwise pattern $D_1$ when:

$$4 \left( 1 + 2 \frac{\mu_2 + C \mu_3}{1 - (\mu_2 + C \mu_3)} \right) p(1-p) B\theta^2 - 1 > 0, \quad (8)$$

whereas the population evolves into a monomorphic population with phenotype $\hat{x}$ when the left-hand side of this equation is negative. The branching condition becomes:

$$4 \left( 1 + 2 \frac{\mu_2}{1 - \mu_2} \right) p(1 - p) B\theta^2 - 1 > 0, \quad (9)$$

for the balanced pattern $D_2$.

It follows that $D_2$ is less favorable to branching than $D_1$ (i.e. left-hand side smaller for (9) than for (8)) (Fig. 3). It can be verified that (i) $D_2$ is the least favorable pattern to branching among all the patterns with the same value of $p$ and that (ii) $D_1$ is the most favorable pattern to branching among all the patterns with the same value of $p$ when $C \mu_3 \ge R \mu_4$. Pattern $D_2$ reaches the same branching criterion value as that when patches are equally connected with philopatry parameter $\mu_2$; $D_1$ reaches the same branching criterion value as that when rows form large patches that are equally connected with philopatry parameter $\mu_2 + C \mu_3$. By symmetry, columnwise patterns have similar properties to those of rowwise patterns.

# 5 Discussion

This article studies how environmental heterogeneity influences adaptive diversification in a multivariate framework. Previous studies on the diversification of several traits found that interactions between traits and the number of resources promoted diversification
<figure><img src="image_4.png" /><figcaption>Figure 3: Influence of frequency of habitat 1 ($p$) and of environmental pattern on the diversification of a lattice population. The figure shows the branching criterion (left-hand sides of (8) and (9)) for a rowwise pattern $D_1$ and a balanced pattern $D_2$ (Table 3). The horizontal line locates the threshold value of 0: the population is polymorphic when its branching criterion is positive and monomorphic when its branching criterion is negative. Parameter values were chosen so that the population is on the brink of diversification for pattern $D_2$ and $p = 1/2$ (the left-hand side of (9) vanishes). Parameter values: $RC \mu_1 = C \mu_3 = R \mu_4 = \mu_2 = 1/4$, $B\theta^2 = 3/5$.</figcaption></figure>
(Doebeli and Ispolatov, 2010; Débarre et al., 2014; Svardal et al., 2014). This article complements these studies by considering multi-dimensional habitats. It shows that associating the levels of some habitat factors across patches (factor correlation) promotes diversification in the absence of trait interactions. Conversely mixing the levels of such factors across patches makes these factors independent, hindering diversification when trait interactions are low. In addition, factor correlations interact with trait interactions, so that they may jointly either increase or decrease the strength of selection. As a result, the patterns that are optimal for hindering or promoting diversification depend on trait interactions.

Our study on multi-dimensional habitats provides explicit invadability criteria, making the influence of environmental heterogeneity easier to study than using a sensitivity analysis. We restricted ourselves to study invadability rather than evolutionary branching because an analytic theory of multi-dimensional branching seems to be currently lacking (Vukics et al., 2003; Durinx et al., 2008). Considering invadability only left many issues unsolved, in particular if evolutionary branching does occur when the ES is invadable, how many types coexist at the end of the diversification process, how many traits become polymorphic and how multi-dimensional habitats influence evolutionary speed.

It has recently been shown that evolutionary trajectories are less likely to converge to an ES and are more likely to be chaotic in multi-dimensional trait spaces (Doebeli and Ispolatov, 2014; Ispolatov et al., 2016). Our model has a unique ES that is always attracting, so that its trajectories do not seem to be chaotic. It has also recently been shown that diversification may occur from other points than ESs in multi-dimensional trait spaces (Ito and Dieckmann, 2014; Ispolatov et al., 2016). For example traits may evolve at different speeds so that the fast traits may branch while the slow traits may still be under directional evolution. Our study considered diversification at the ES only, but diversification is perhaps possible from other points of evolutionary trajectories.

The influence of environmental heterogeneity on diversification depends on the manner individuals disperse. Previous studies found that when dispersal rates decreased with distance, landscapes showing large spatially continuous habitat clusters promoted diversification (Débarre and Gandon, 2010; Birand et al., 2012; Haller et al., 2013; Papaïx et al., 2013). When patches clustered into groups such that patches in the same group were more connected, diversification was found to be hindered by making the within-group habitat distributions equal (Papaïx et al., 2013). Our study extends the latter result to a dispersal involving two crossed factors, called the row and column factors. Diversification is found to be hindered by making habitat distributions equal both within rows and within columns. These environmental patterns are reminiscent of Latin squares in the domain of experimental designs (Bailey, 2008).

In our study, habitat differentiation promotes diversification (invadability and branching criteria are increasing functions of θ), in agreement with some earlier results (see e.g. Geritz et al., 1998; Débarre and Gandon, 2010; Papaïx et al., 2013). This result presumably partly comes from the assumption that patch sizes are fixed and equal. When
population size is variable, intermediate levels of heterogeneity may maximize branching
propensity (Doebeli and Dieckmann, 2003; Débarre et al., 2013; Haller et al., 2013). A
possible explanation is that when habitat differenciation is too pronounced, dispersers
cannot persit long enough in hostile habitats to adapt. However, polymorphism may
evolve even if population size is variable and heterogeneity is high, if mutations have
large effects, the initial population is polymorphic or demography is stochastic (Débarre
et al., 2013).

Our results are valid when patches are large. When patch number and size are both
small, demographic stochasticity, that was not taken into account by our deterministic
model, may delay or even prevent branching and reduce diversity within diverging pop-
ulations when branching occurs (Claessen et al., 2007; Débarre et al., 2013; Wakano and
Iwasa, 2013). When patches are numerous, demographic stochasticity at the level of the
whole metapopulation may be safely ignored. However if in addition patches are small
and dispersal is low, interactions between mutant individuals within patches during the
initial growth phase of the mutant population, that were neglected in our quantification
of fitness (Appendix A), may become influential. In our model, they could lead to kin
competition and could hinder branching (Wakano and Lehmann, 2014; Mullon et al.,
2016).

We focused on a discrete-time model in which dispersal occurs between density reg-
ulation and selection. It corresponds to a soft selection model with philopatry (Débarre
and Gandon, 2011). Our results are presumably highly dependent on this choice. The
evolution of polymorphism in discrete-time models is influenced by the life cycle chosen
and in particular the order of demographic events (Débarre and Gandon, 2011). For ex-
ample placing dispersal between selection and density regulation leads to a hard selection
model with philopatry in which polymorphism is unlikely to evolve.

Acknowledgements

This study was funded by Agence Nationale de la Recherche (project GESTER ANR-11-
AGRO-0003). We thank two anonymous reviewers for helpful comments.

Appendix

A Invasion fitness of a mutant in a resident popula-
tion

Let a mutant with trait vector *y* appear in the resident population with trait *x*. To model the dynamics of the mutant's population, we take account of the patch structure and consider that time is discrete with one time step corresponding to a life cycle. The
number $n_i(t + 1)$ of mutant individuals in patch $i$ at time $t + 1$ satisfies (Papaïx et al., 2013):

$$n_i(t + 1) = \frac{\sum_{i'} g(y, \beta(h(i))) m_{ii'} n_{i'}(t)}{\sum_{i'} g(x, \beta(h(i))) m_{ii'} (K - n_{i'}(t)) + \sum_{i'} g(y, \beta(h(i))) m_{ii'} n_{i'}(t)} K,$$

where $h(i)$ is the habitat of patch $i$. As $n_{i'}(t) \ll K$ when the mutant is rare and $\sum_{i'} m_{ii'} = 1$, the initial demography of the mutant population in patch $i$ is approximated by:

$$n_i(t + 1) \approx \frac{\sum_{i'} g(y, \beta(h(i))) m_{ii'} n_{i'}(t)}{g(x, \beta(h(i)))},$$

for $i = 1, \dots, P$. This system of equations is written as a matrix model (Caswell, 2001; Ellner and Guckenheimer, 2006):

$$n(t + 1) \approx A(x, y) n(t),$$

where $n(t) = (n_1(t), \dots, n_P(t))^T$ is the vector of population sizes and $A(x, y)$ is a $P \times P$ projection matrix with non-negative and constant coefficients (given $x$ and $y$). The entry in row $i$ and column $i'$ of the projection matrix is equal to:

$$[A(x, y)]_{ii'} = \frac{g(y, \beta(h(i)))}{g(x, \beta(h(i)))} m_{ii'}. \qquad (A.1)$$

The matrix $A(x, y)$ is assumed to be diagonalizable with eigenvalues, right and left eigenvectors denoted by $\lambda^{(i)}(x, y)$, $r^{(i)}(x, y)$ and $l^{(i)}(x, y)$, $i = 1, \dots, P$, respectively. The eigenvectors of $A(x, y)$ are assumed to be scaled so that $l^{(i)}(x, y)^* r^{(i)}(x, y) = 1$ and $l^{(i)}(x, y)^* r^{(i')}(x, y) = 0$ for $i \neq i'$, where $l^{(i)}(x, y)^*$ denotes the complex conjugate trans- pose of $l^{(i)}(x, y)$. As every patch can contribute to every other patch (Section 2), $A(x, y)$ is irreducible and has a dominant eigenvalue $\lambda^{(1)}(x, y)$ (Caswell, 2001). Since some of the eigenvalues of $A(x, y)$ may be complex, $\lambda^{(1)}(x, y)$ is defined as the eigenvalue with the largest real part. It is real, simple and its associated left and right eigenvectors are posi- tive. As it quantifies the growth of the mutant's population (Caswell, 2001), the invasion fitness $s(x, y)$ of a mutant with trait vector $y$ in a resident population with trait vector $x$ is defined as follows:

$$s(x, y) = \lambda^{(1)}(x, y).$$

The mutant is assumed to replace the resident when its fitness is larger than one and to go extinct when its fitness is smaller than one.

## B Diagonalization of $A(x, x)$

From (A.1), the entry in row $i$ and column $i'$ of $A(x, x)$ is equal to:

$$[A(x, x)]_{ii'} = m_{ii'}. $$
For the row-column dispersal, the diagonalization of $A(x, x)$ gives:

$$A(x, x) = (RC\mu_1 + C\mu_3 + R\mu_4 + \mu_2)S_1 + (C\mu_3 + \mu_2)S_2 + (R\mu_4 + \mu_2)S_3 + \mu_2S_4,$$

where $S_1 = (J_R/R) \otimes (J_C/C)$, $S_2 = (I_R - J_R/R) \otimes (J_C/C)$, $S_3 = (J_R/R) \otimes (I_C - J_C/C)$, $S_4 = (I_R - J_R/R) \otimes (I_C - J_C/C)$, $I_R$ is the $R \times R$ identity matrix, $J_R$ is the $R \times R$ matrix of ones and $\otimes$ is the Kronecker product. Thus, the eigenvalues of $A(x, x)$ are $\lambda^{(1)}(x, x) = 1$ with multiplicity 1, $\lambda^{(2)}(x, x) = C\mu_3 + \mu_2$ with multiplicity $R-1$, $\lambda^{(3)}(x, x) = R\mu_4 + \mu_2$ with multiplicity $C-1$ and $\lambda^{(4)}(x, x) = \mu_2$ with multiplicity $(R-1)(C-1)$. The dominant right and left eigenvectors satisfy:

$$[r^{(1)}(x, x)]_i = 1/P, \qquad (B.1)$$

$$[l^{(1)}(x, x)]_i = 1, \qquad (B.2)$$

for $i = 1, \dots, P$. Fixing $\mu_3 = \mu_4 = 0$ yields the diagonalization of $A(x, x)$ for the case when patches are equally connected.

## C Fitness derivatives

Let $\phi_u(x, y)$ denote the first-order derivative of the function $\phi(x, y)$ with respect to (w.r.t.) the scalar (i.e. one-dimensional) trait $u$; for example when $\phi = s$ and $u = y_1$, we obtain $s_{y_1}(x, y) = \frac{\partial s}{\partial y_1}(x, y)$. The first-order derivative $s_u(x, y)$ of fitness w.r.t. the scalar trait $u$, or fitness gradient, is equal to (Caswell, 2001):

$$s_u(x, y) = l^{(1)}(x, y)^T A_u(x, y) r^{(1)}(x, y), \qquad (C.1)$$

for $u = x_1, \dots, x_d, y_1, \dots, y_d$. In (C.1), $A_u(x, y)$ is the $P \times P$ matrix whose element in row $i$ and column $i'$ is equal to $\frac{\partial[A]_{ii'}}{\partial u}(x, y)$.

Let $\phi_{uv}(x, y)$ denote the second-order derivative of $\phi(x, y)$ w.r.t. the scalar traits $u$ and $v$; for example when $\phi = s$, $u = y_1$ and $v = y_2$, we obtain $s_{y_1 y_2}(x, y) = \frac{\partial^2 s}{\partial y_1 \partial y_2}(x, y)$. Using (C.1), we find:

$$s_{uv} = l_v^{(1)T} A_u r^{(1)} + l^{(1)T} A_{uv} r^{(1)} + l^{(1)T} A_u r_v^{(1)},$$

where we have suppressed the dependence of quantities such as $s_{uv}(x, y)$ on $(x, y)$ for notational convenience and $A_{uv}$ is the $P \times P$ matrix whose element in row $i$ and column $i'$ is equal to $\frac{\partial^2[A]_{ii'}}{\partial u \partial v}$. According to equations 9.131 and 9.132 in Caswell (2001), we have:

$$r_v^{(1)} = \sum_{i=2}^{P} \frac{l^{(i)\star} A_v r^{(1)}}{\lambda^{(1)} - \lambda^{(i)}} r^{(i)},$$

$$l_v^{(1)} = \sum_{i=2}^{P} \frac{r^{(i)\star} A_v^* l^{(1)}}{\bar{\lambda}^{(1)} - \bar{\lambda}^{(i)}} l^{(i)},$$
where $\bar{\lambda}^{(i)}$ is the complex conjugate of $\lambda^{(i)}$ and $r^{(i)\star}$ is the complex conjugate transpose of $r^{(i)}$. As $l^{(1)}$ is real, we have:

$$l_v^{(1)T} = l_v^{(1)\star} = \sum_{i=2}^{P} \frac{l^{(1)T} A_v r^{(i)}}{\lambda^{(1)} - \lambda^{(i)}} l^{(i)\star},$$

so that the derivative $s_{uv}(x, y)$ of fitness w.r.t. the scalar traits $u$ and $v$ is equal to:

$$s_{uv} = l^{(1)T} A_{uv} r^{(1)} + \sum_{i=2}^{P} \frac{(l^{(1)T} A_u r^{(i)})(l^{(i)\star} A_v r^{(1)}) + (l^{(1)T} A_v r^{(i)})(l^{(i)\star} A_u r^{(1)})}{\lambda^{(1)} - \lambda^{(i)}}, \quad (C.2)$$

for $u, v = x_1, \dots, x_d, y_1, \dots, y_d$.

## D Evolutionary singularities

ESs $\hat{x}$ are the solutions of the equation (Meszéna et al., 2002):

$$s_y(x, x) = 0,$$

where $s_y$ is the vector of dimension $d$ equal to the derivative of $s$ w.r.t. the vector $y$. We calculate $s_y(x, x)$ using (C.1). As the elements of $A(x, y)$ satisfy (A.1) and the proportion $g(x, \beta)$ of juveniles that become adults can be written as follows:

$$g(x, \beta) = \alpha \exp \left( - \sum_j B_{jj}(x_j - \beta_j)^2 / 2 - 2 \sum_{j<j'} B_{jj'}(x_j - \beta_j)(x_{j'} - \beta_{j'}) / 2 \right),$$

the derivatives of $A$ w.r.t. the mutant's traits are equal to:

$$A_{y_j}(x, y) = - \left( \sum_{j'} B_{jj'} \Delta^{(j')}(y) \right) A(x, y), \quad j = 1, \dots, d, \quad (D.1)$$

where $\Delta^{(j)}(y)$ is the $P \times P$ diagonal matrix whose $i$th diagonal element is equal to $y_j - \beta(h(i))_j$, where $h(i)$ is the habitat of patch $i$. Using (C.1), we obtain that ESs are the solutions of the system of equations:

$$\sum_{j'} B_{jj'} l^{(1)}(x, x)^T \Delta^{(j')}(x) r^{(1)}(x, x) = 0, \quad j = 1, \dots, d.$$

As $B$ is invertible, we find that ESs satisfy:

$$l^{(1)}(\hat{x}, \hat{x})^T \Delta^{(j)}(\hat{x}) r^{(1)}(\hat{x}, \hat{x}) = 0, \quad j = 1, \dots, d.$$
Using (B.1) and (B.2), we obtain:

$$ \hat{x} = \frac{1}{P} \sum_{i} \beta(h(i)) = \sum_{h} p_{1...d}(h) \beta(h). $$

As $\beta(h) = ((-1)^{h_1}, \dots, (-1)^{h_d})^T \theta$, the ES is equal to:

$$ \hat{x} = (1 - 2p_1, \dots, 1 - 2p_d)^T \theta, \qquad (D.2) $$

where $p_j$ denotes $p_j(1)$.

## E Attractivity

For a multi-dimensional trait, the attractivity of an ES depends on the variance-covariance matrix of mutational steps. For clonally reproducing organisms, $\hat{x}$ is attracting whatever the mutational variance-covariance matrix (strong convergence stability) if:

$$ s_{xx}(\hat{x}, \hat{x}) - s_{yy}(\hat{x}, \hat{x}) > 0, $$

where $s_{xx}$ (respectively $s_{yy}$) denotes the $d \times d$ matrix equal to the second-order derivative of $s$ w.r.t. the vector $x$ (respectively $y$) and the symbol $>$ means positive definite (Meszéna et al., 2002; Leimar, 2005).

We calculate $s_{xx}(\hat{x}, \hat{x}) - s_{yy}(\hat{x}, \hat{x})$ using (C.2). The first-order derivatives of $A$ w.r.t. the mutant's traits are given by (D.1). Those w.r.t. the resident's traits are equal to:

$$ A_{x_j} = -A_{y_j}, \quad j = 1, \dots, d. $$

Thus, the second term of the right-hand side of (C.2) is identical for $s_{x_j x_{j'}}(\hat{x}, \hat{x})$ and $s_{y_j y_{j'}}(\hat{x}, \hat{x})$, $1 \le j, j' \le d$, so that:

$$ s_{x_j x_{j'}}(\hat{x}, \hat{x}) - s_{y_j y_{j'}}(\hat{x}, \hat{x}) = l^{(1)}(\hat{x}, \hat{x})^T A_{x_j x_{j'}}(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}) - l^{(1)}(\hat{x}, \hat{x})^T A_{y_j y_{j'}}(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}). $$

Using (D.1), we find that the second-order derivatives of $A$ are equal to:

$$ A_{y_j y_{j'}}(x, y) = \left( \left( \sum_k B_{jk} \Delta^{(k)}(y) \right) \left( \sum_k B_{j'k} \Delta^{(k)}(y) \right) - B_{jj'} I_P \right) \times A(x, y), \qquad (E.1) $$

$$ A_{x_j x_{j'}}(x, y) = \left( \left( \sum_k B_{jk} \Delta^{(k)}(x) \right) \left( \sum_k B_{j'k} \Delta^{(k)}(x) \right) + B_{jj'} I_P \right) A(x, y), $$
whatever $1 \le j, j' \le d$. Thus, we obtain:

$$s_{x_j x_{j'}}(\hat{x}, \hat{x}) - s_{y_j y_{j'}}(\hat{x}, \hat{x}) = 2B_{jj'} l^{(1)}(\hat{x}, \hat{x})^T A(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}) = 2B_{jj'}.$$

It follows that:

$$s_{xx}(\hat{x}, \hat{x}) - s_{yy}(\hat{x}, \hat{x}) = 2B > 0,$$

and that $\hat{x}$ is attracting whatever the mutational variance-covariance matrix, the matrix $B$, the dispersal parameters and the environmental pattern.

# F Diversification in multi-dimensional environments

If the fitness Hessian $s_{yy}(\hat{x}, \hat{x})$ is negative definite:

$$s_{yy}(\hat{x}, \hat{x}) < 0,$$

then the function $s(\hat{x}, y)$ attains a local maximum at $y = \hat{x}$ and $\hat{x}$ is uninvadable (Meszéna et al., 2002; Leimar, 2005).

Using (C.2), we obtain:

$$\begin{aligned} s_{y_j y_{j'}}(\hat{x}, \hat{x}) = & l^{(1)}(\hat{x}, \hat{x})^T A_{y_j y_{j'}}(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}) & \quad (\text{F.1}) \\ & + l^{(1)}(\hat{x}, \hat{x})^T A_{y_j}(\hat{x}, \hat{x}) \left( \sum_{i=2}^{P} \frac{r^{(i)}(\hat{x}, \hat{x}) l^{(i)}(\hat{x}, \hat{x})^*}{\lambda^{(1)}(\hat{x}, \hat{x}) - \lambda^{(i)}(\hat{x}, \hat{x})} \right) A_{y_{j'}}(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}) \\ & + l^{(1)}(\hat{x}, \hat{x})^T A_{y_{j'}}(\hat{x}, \hat{x}) \left( \sum_{i=2}^{P} \frac{r^{(i)}(\hat{x}, \hat{x}) l^{(i)}(\hat{x}, \hat{x})^*}{\lambda^{(1)}(\hat{x}, \hat{x}) - \lambda^{(i)}(\hat{x}, \hat{x})} \right) A_{y_j}(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}). \end{aligned}$$

Using (E.1), we find that the first term of the right-hand side of (F.1) is equal to:

$$l^{(1)}(\hat{x}, \hat{x})^T A_{y_j y_{j'}}(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}) = \sum_h p_{1 \dots d}(h) \left( \sum_k B_{jk}(\hat{x}_k - \beta(h)_k) \right) \left( \sum_k B_{j'k}(\hat{x}_k - \beta(h)_k) \right) - B_{jj'}.$$

Using Appendix B, we obtain:

$$\sum_{i=2}^{P} \frac{r^{(i)}(\hat{x}, \hat{x}) l^{(i)}(\hat{x}, \hat{x})^*}{\lambda^{(1)}(\hat{x}, \hat{x}) - \lambda^{(i)}(\hat{x}, \hat{x})} = \frac{1}{1 - \mu_2} (I_P - J_P / P).$$

Using (D.1), we obtain that the second and third terms of the right-hand side of (F.1) are equal to:

$$\left(\frac{\mu_2}{1 - \mu_2}\right) \sum_h p_{1 \dots d}(h) \left(\sum_k B_{jk}(\hat{x}_k - \beta(h)_k)\right) \left(\sum_k B_{j'k}(\hat{x}_k - \beta(h)_k)\right).$$
It follows that:

$$ s_{y_j y_j'}(\hat{x}, \hat{x}) = \left(1 + 2\frac{\mu_2}{1 - \mu_2}\right) \times \sum_h p_{1\dots d}(h) \left(\sum_k B_{jk}(\hat{x}_k - \beta(h)_k)\right) \left(\sum_k B_{j'k}(\hat{x}_k - \beta(h)_k)\right) - B_{jj'}. \quad (F.2) $$

## F.1 Influence of factor correlations

Equations (F.2), (D.2) and (3) lead to (4) and (5) when traits do not interact (B is diagonal).

## F.2 Joint influence of factor correlations and trait interactions

We focus on the case when the matrix B is completely symmetric to simplify: its diagonal elements are equal to $\tau$ and its non-diagonal elements are equal to $\tau\gamma$. We also focus on environmental patterns such that the frequencies of level 1 are all equal to $p$ and such that factor correlations are all equal to $\rho$. Using (F.2), it can be verified that the elements of $s_{yy}(\hat{x}, \hat{x})$ are equal to:

$$ s_{y_j y_j'}(\hat{x}, \hat{x}) = 4 \left(1 + 2\frac{\mu_2}{1 - \mu_2}\right) \tau^2 \theta^2 p(1 - p) \times \left((2\gamma + (d-2)\gamma^2)(1 - \rho) + (1 + (d-1)\gamma)^2\rho\right) - \tau\gamma, \\ s_{y_j y_j'}(\hat{x}, \hat{x}) = 4 \left(1 + 2\frac{\mu_2}{1 - \mu_2}\right) \tau^2 \theta^2 p(1 - p) \times \left((1 + (d-1)\gamma^2)(1 - \rho) + (1 + (d-1)\gamma)^2\rho\right) - \tau, $$

for $1 \le j \ne j' \le d$. Thus, $s_{yy}(\hat{x}, \hat{x})$ is completely symmetric and its eigenvalues are equal to:

$$ \omega^{(1)} = 4 \left(1 + 2\frac{\mu_2}{1 - \mu_2}\right) \tau^2 \theta^2 p(1 - p)(1 - \gamma)^2 (1 - \rho) - \tau(1 - \gamma), \\ \omega^{(2)} = 4 \left(1 + 2\frac{\mu_2}{1 - \mu_2}\right) \tau^2 \theta^2 p(1 - p)(1 + (d-1)\gamma)^2 (1 + (d-1)\rho) - \tau(1 + (d-1)\gamma). $$

The matrix $s_{yy}(\hat{x}, \hat{x})$ is negative definite iff both eigenvalues are negative, i.e. iff (7) is satisfied. When $\gamma = 0$, we obtain (6).

## G Diversification with a structured dispersal

When a single trait evolves and when $s_{yy}(\hat{x}, \hat{x}) > 0$, an attracting ES is a branching point (Geritz et al., 1998). Thus, we study when this condition is satisfied.
Using (C.2), we obtain:

$$ s_{yy}(\hat{x}, \hat{x}) = l^{(1)}(\hat{x}, \hat{x})^T A_{yy}(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}) \qquad (G.1) \\ + 2 l^{(1)}(\hat{x}, \hat{x})^T A_y(\hat{x}, \hat{x}) \left( \sum_{i=2}^{P} \frac{r^{(i)}(\hat{x}, \hat{x}) l^{(i)}(\hat{x}, \hat{x})^*}{\lambda^{(1)}(\hat{x}, \hat{x}) - \lambda^{(i)}(\hat{x}, \hat{x})} \right) A_y(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}). $$

Using similar calculations as in Appendix F, we find that:

$$ A_y(x, y) = -B \Delta(y) A(x, y), \qquad (G.2) $$

$$ A_{yy}(x, y) = B(B \Delta^2(y) - I_P) A(x, y), \qquad (G.3) $$

where $\Delta(y)$ is the $P \times P$ diagonal matrix whose $i$th diagonal element is equal to $y - \beta(h(i))$ and $h(i)$ is the habitat of patch $i$. Using (D.2), (B.1), (B.2) and (G.3), we find that the first term of the right-hand side of (G.1) is equal to:

$$ l^{(1)}(\hat{x}, \hat{x})^T A_{yy}(\hat{x}, \hat{x}) r^{(1)}(\hat{x}, \hat{x}) = B(4p(1-p)B\theta^2 - 1). \qquad (G.4) $$

Using (G.2), we find that the second term of the right-hand side of (G.1) is equal to:

$$ 2B^2 l^{(1)}(\hat{x}, \hat{x})^T \Delta(\hat{x}) \left( \sum_{i=2}^{P} \frac{\lambda^{(i)}(\hat{x}, \hat{x}) r^{(i)}(\hat{x}, \hat{x}) l^{(i)}(\hat{x}, \hat{x})^*}{\lambda^{(1)}(\hat{x}, \hat{x}) - \lambda^{(i)}(\hat{x}, \hat{x})} \right) \Delta(\hat{x}) r^{(1)}(\hat{x}, \hat{x}). \qquad (G.5) $$

Using Appendix B, we obtain:

$$ \begin{aligned} \sum_{i=2}^{P} \frac{\lambda^{(i)}(\hat{x}, \hat{x}) r^{(i)}(\hat{x}, \hat{x}) l^{(i)}(\hat{x}, \hat{x})^*}{\lambda^{(1)}(\hat{x}, \hat{x}) - \lambda^{(i)}(\hat{x}, \hat{x})} &= \eta_2 S_2 + \eta_3 S_3 + \eta_4 S_4, \\ &= (\eta_2 - \eta_4)S_2 + (\eta_3 - \eta_4)S_3 \\ &\quad + \eta_4(S_2 + S_3 + S_4), \end{aligned} \qquad (G.6) $$

where $\eta_i = \frac{\lambda^{(i)}(\hat{x}, \hat{x})}{\lambda^{(1)}(\hat{x}, \hat{x}) - \lambda^{(i)}(\hat{x}, \hat{x})}$, $i = 2, 3, 4$, $\eta_4 = \frac{\mu_2}{1-\mu_2}$, $\eta_2 - \eta_4 = \frac{C \mu_3}{(1-\mu_2)(RC \mu_1 + R \mu_4)}$ and $\eta_3 - \eta_4 = \frac{R \mu_4}{(1-\mu_2)(RC \mu_1 + C \mu_3)}$. Using (D.2), (B.1) and (B.2), we obtain:

$$ l^{(1)}(\hat{x}, \hat{x})^T \Delta(\hat{x}) (S_2 + S_3 + S_4) \Delta(\hat{x}) r^{(1)}(\hat{x}, \hat{x}) = 4p(1-p)\theta^2, \qquad (G.7) $$

since $S_2 + S_3 + S_4 = I_P - S_1$. Using (G.4)-(G.7), we obtain:

$$ s_{yy}(\hat{x}, \hat{x}) = B \left( 4 \left( 1 + 2 \frac{\mu_2}{1 - \mu_2} \right) p(1 - p)B\theta^2 - 1 \right) \\ + 2B^2 l^{(1)}(\hat{x}, \hat{x})^T \Delta(\hat{x}) ((\eta_2 - \eta_4)S_2 + (\eta_3 - \eta_4)S_3) \Delta(\hat{x}) r^{(1)}(\hat{x}, \hat{x}). $$
Using (D.2), (B.1) and (B.2), we obtain for the rowwise pattern $D_1$:

$$l^{(1)}(\hat{x}, \hat{x})^T \Delta(\hat{x}) S_2 \Delta(\hat{x}) r^{(1)}(\hat{x}, \hat{x}) = 4p(1-p)\theta^2,$$

$$l^{(1)}(\hat{x}, \hat{x})^T \Delta(\hat{x}) S_3 \Delta(\hat{x}) r^{(1)}(\hat{x}, \hat{x}) = 0,$$

which leads to (8). For the balanced pattern $D_2$, we have:

$$l^{(1)}(\hat{x}, \hat{x})^T \Delta(\hat{x}) S_2 \Delta(\hat{x}) r^{(1)}(\hat{x}, \hat{x}) = 0,$$

$$l^{(1)}(\hat{x}, \hat{x})^T \Delta(\hat{x}) S_3 \Delta(\hat{x}) r^{(1)}(\hat{x}, \hat{x}) = 0,$$

which leads to (9).

## References

Bailey, R. A., 2008. Design of comparative experiments. Vol. 25. Cambridge University Press.

Barton, N. H., Briggs, D. E. G., Eisen, J. A., Goldstein, D. B., Patel, N. H., 2007. Evolution. Cold Spring Harbor Laboratory Press, NY.

Beaunée, G., Vergu, E., Ezanno, P., 2015. Modelling of paratuberculosis spread between dairy cattle farms at a regional scale. Veterinary Research 46 (1).

Birand, A., Vose, A., Gavrilets, S., 2012. Patterns of species ranges, speciation, and extinction. Am Nat 179 (1), 1–21.

Caswell, H., 2001. Matrix Population Models. Sinauer Associates Inc. Publishers, Sunderland (USA).

Claessen, D., Andersson, J., Persson, L., de Roos, A. M., 2007. Delayed evolutionary branching in small populations. Evolutionary Ecology Research 9 (1), 51–69.

Débarre, F., Gandon, S., 2010. Evolution of specialization in a spatially continuous environment. Journal of Evolutionary Biology 23 (5), 1090–1099.

Débarre, F., Gandon, S., 2011. Evolution in heterogeneous environments: between soft and hard selection. The American Naturalist 177 (3), E84–E97.

Débarre, F., Nuismer, S. L., Doebeli, M., 2014. Multidimensional (co)evolutionary stability. Am Nat 184 (2), 158–171.

Débarre, F., Ronce, O., Gandon, S., 2013. Quantifying the effects of migration and mutation on adaptation and demography in spatially heterogeneous environments. J Evol Biol 26 (6), 1185–1202.
Diekmann, O., 2004. A beginner's guide to adaptive dynamics. Banach Center Publications 63, 47-86.

Doebeli, M., 2011. Adaptive Diversification. Princeton University Press.

Doebeli, M., Dieckmann, U., 2003. Speciation along environmental gradients. Nature 421, 259-264.

Doebeli, M., Ispolatov, I., 2010. Complexity and diversity. Science 328 (5977), 494-497.

Doebeli, M., Ispolatov, I., 2014. Chaos and unpredictability in evolution. Evolution 68 (5), 1365-1373.

Durinx, M., Metz, J. A. J. H., Meszéna, G., 2008. Adaptive dynamics for physiologically structured population models. J Math Biol 56 (5), 673-742.

Ellner, S. P., Guckenheimer, J., 2006. Dynamic Models in Biology, 1st Edition. Princeton University Press.

Fabre, F., Rousseau, E., Mailleret, L., Moury, B., 2015. Epidemiological and evolutionary management of plant resistance: optimizing the deployment of cultivar mixtures in time and space in agricultural landscapes. Evolutionary Applications 8 (10), 919-932.

Gautier, M., Flori, L., Riebler, A., Jaffrézic, F., Laloé, D., Gut, I., Moazami-Goudarzi, K., Foulley, J.-L., 2009. A whole genome Bayesian scan for adaptive genetic divergence in West African cattle. BMC genomics 10 (1), 550.

Gavrilets, S., Vose, A., 2005. Dynamic patterns of adaptive radiation. Proceedings of the National Academy of Sciences of the United States of America 102 (50), 18040-18045.

Geritz, S. A. H., Kisdi, E., Meszéna, G., Metz, J. A. J. H., 1998. Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree. Evolutionary Ecology 12, 35-57.

Haller, B. C., Mazzucco, R., Dieckmann, U., 2013. Evolutionary branching in complex landscapes. Am Nat 182 (4), E127-E141.

Hartl, D. L., Clark, A. G., 1997. Principles of population genetics. Vol. 116. Sinauer associates Sunderland.

Ispolatov, I., Madhok, V., Doebeli, M., 2016. Individual-based models for adaptive diversification in high-dimensional phenotype spaces. Journal of Theoretical Biology 390, 97-105.

Ito, H. C., Dieckmann, U., 2014. Evolutionary branching under slow directional evolution. Journal of Theoretical Biology 360, 290-314.
Johnson, M. T. J., 2007. Genotype-by-environment interactions leads to variable selection on life-history strategy in common evening primrose (Oenothera biennis). Journal of evolutionary biology 20 (1), 190–200.

Johnson, P. T. J., de Roode, J. C., Fenton, A., 2015. Why infectious disease research needs community ecology. Science 349 (6252).

Karlin, S., 1982. Classifications of selection-migration structures and conditions for a protected polymorphism. Evol. Biol 14, 61–204.

Kocher, T. D., 2004. Adaptive evolution and explosive speciation: the cichlid fish model. Nature Reviews Genetics 5 (4), 288–298.

Laughlin, D. C., Messier, J., 2015. Fitness of multidimensional phenotypes in dynamic adaptive landscapes. Trends Ecol Evol 30 (8), 487–496.

Leimar, O., 2005. The evolution of phenotypic polymorphism: randomized strategies versus evolutionary branching. American Naturalist 165 (6), 669–681.

Liccardo, A., Fierro, A., 2013. A lattice model for influenza spreading. PloS one 8 (5), e63935.

Meszéna, G., Czibula, I., Geritz, S., 1997. Adaptative dynamics in a 2-patch environment: a toy model for allopatric and parapatric speciation. Journal of Biological Systems 5, 265–284.

Meszéna, G., Kisdi, E., Dieckmann, U., Geritz, S. A. H., Metz, J. A. J., 2002. Evolutionary optimisation models and matrix games in the unified perspective of adaptive dynamics. Selection 2 (1), 193–220.

Mullon, C., Keller, L., Lehmann, L., Nuismer, S. L., Bronstein, J. L., 2016. Evolutionary stability of jointly evolving traits in subdivided populations. The American Naturalist 188 (2), 175–195.

Papaïx, J., David, O., Lannou, C., Monod, H., 2013. Dynamics of adaptation in spatially heterogeneous metapopulations. PLoS One 8 (2), e54697.

Papaïx, J., Monod, H., Goyeau, H., du Cheyron, P., Lannou, C., 2011. Influence of cultivated landscape composition on variety resistance: an assessment based on wheat leaf rust epidemics. New Phytologist 191, 1095–1107.

Papaïx, J., Touzeau, S., Monod, H., Lannou, C., 2014. Can epidemic control be achieved by altering landscape connectivity in agricultural systems? Ecological Modelling 284, 35–47.
Pariaud, B., Ravigné, V., Halkett, F., Goyeau, H., Carlier, J., Lannou, C., 2009. Ag-gressiveness and its role in the adaptation of plant pathogens. Plant Pathology 58 (3), 409-424.

Poisot, T., Bever, J. D., Nemri, A., Thrall, P. H., Hochberg, M. E., 2011. A conceptual framework for the evolution of ecological specialisation. Ecology Letters 14 (9), 841-851.

REX Consortium, 2013. Heterogeneity of selection and the evolution of resistance. Trends in Ecology & Evolution 28 (2), 110-118.

Rhoné, B., Remoue, C., Galic, N., Goldringer, I., Bonnin, I., 2008. Insight into the genetic bases of climatic adaptation in experimentally evolving wheat populations. Molecular Ecology 17 (3), 930-943.

Roemhild, R., Barbosa, C., Beardmore, R. E., Jansen, G., Schulenburg, H., 2015. Temporal variation in antibiotic environments slows down resistance evolution in pathogenic Pseudomonas aeruginosa. Evolutionary Applications 8 (10), 945-955.

Sloan, C., Moore, M. L., Hartert, T., 2011. Impact of pollution, climate, and sociode-mographic factors on spatiotemporal dynamics of seasonal respiratory viruses. Clinical and Translational Science 4 (1), 48-54.

Svardal, H., Rueffler, C., Doebeli, M., 2014. Organismal complexity and the potential for evolutionary diversification. Evolution 68 (11), 3248-3259.

Vukics, A., Asbóth, J., Meszéna, G., 2003. Speciation in multidimensional evolutionary space. Phys Rev E Stat Nonlin Soft Matter Phys 68, 041903.

Wakano, J. Y., Iwasa, Y., 2013. Evolutionary branching in a finite population: deterministic branching vs. stochastic branching. Genetics 193 (1), 229-241.

Wakano, J. Y., Lehmann, L., 2014. Evolutionary branching in deme-structured populations. Journal of theoretical biology 351, 83-95.

Wallinga, J., Teunis, P., Kretzschmar, M., 2006. Using data on social contacts to estimate age-specific transmission parameters for respiratory-spread infectious agents. American journal of epidemiology 164 (10), 936-944.

Watts, D. J., Muhamad, R., Medina, D. C., Dodds, P. S., 2005. Multiscale, resurgent epidemics in a hierarchical metapopulation model. Proceedings of the National Academy of Sciences of the United States of America 102 (32), 11157-11162.

West, J. S., 2014. Plant pathogen dispersal. In: eLS. John Wiley & Sons Ltd, Chichester.
URL http://www.els.net