Clines in quantitative traits: the role of migration patterns and selection scenarios.

The existence, uniqueness, and shape of clines in a quantitative trait under selection toward a spatially varying optimum is studied. The focus is on deterministic diploid two-locus n-deme models subject to various migration patterns and selection scenarios. Migration patterns may exhibit isolation by distance, as in the stepping-stone model, or random dispersal, as in the island model. The phenotypic optimum may change abruptly in a single environmental step, more gradually, or not at all. Symmetry assumptions are imposed on phenotypic optima and migration rates. We study clines in the mean, variance, and linkage disequilibrium (LD). Clines result from polymorphic equilibria. The possible equilibrium configurations are determined as functions of the migration rate. Whereas for weak migration, many polymorphic equilibria may be simultaneously stable, their number decreases with increasing migration rate. Also for intermediate migration rates polymorphic equilibria are in general not unique, however, for loci of equal effects the corresponding clines in the mean, variance, and LD are unique. For sufficiently strong migration, no polymorphism is maintained. Both migration pattern and selection scenario exert strong influence on the existence and shape of clines. The results for discrete demes are compared with those from models in which space varies continuously and dispersal is modeled by diffusion. Comparisons with previous studies, which investigated clines under neutrality or under linkage equilibrium, are performed. If there is no long-distance migration, the environment does not change abruptly, and linkage is not very tight, populations are almost everywhere close to linkage equilibrium.

Introduction

Strength and patterns of migration in a spatially structured population are important determinants of the degree of local adaptation and the genetic variation that can be maintained in a heterogeneous environment. If there is an environmental gradient, clines in the gene frequencies or in the mean or other characteristics of a quantitative trait may be established. Such clines have been an important topic of both empirical and theoretical research since the pioneering work of  Haldane (1948). The existence and shape of clines depends on the strength and patterns of migration, the properties of spatially varying selection, and the underlying genetics. In this work we assume that genetic variation is maintained by selection and migration, and ignore mutation and random genetic drift.

For populations subdivided into discrete demes, migration is frequently modeled by an island or a stepping-stone model. The former assumes that outbreeding individuals disperse uniformly to all other demes, whereas the latter assumes that the probability of migration decreases with distance, i.e., there is isolation by distance. For populations occupying a continuous habitat, migration is usually approximated by diffusion. Diffusion models, as well as certain generalizations, are derived by assuming that large migration steps are unlikely in short time intervals and selection is weak (Nagylaki, 1975, Nagylaki, 1989). Naturally, such models exhibit isolation by distance.

The large majority of previous theoretical investigations assumes that selection acts on a single diallelic locus. For discrete demes and numerous types of migration patterns and selection schemes,  Karlin (1982) performed a comprehensive investigation on the maintenance of protected polymorphisms (corresponding to the existence of clines). Although his results show that more mixing (e.g., by a higher migration rate or by migration to more distant demes) tends to restrict the conditions for a protected polymorphism, he also gave examples where less mixing inhibits a polymorphism. We shall compare the propensity of frequently employed migration patterns in maintaining clines at two recombining loci. Reviews of the extensive literature on one-locus migration–selection models may be found in  Lenormand (2002),Nagylaki and Lou (2008), and  Bürger (2014).

Also for the diffusion model, as well as more general forms of dispersal in continuous space, a wealth of results about existence, uniqueness, and properties of polymorphic equilibria and clines in gene frequencies at a single locus has accumulated. This literature is reviewed by  Nagylaki and Lou (2008) and  Lou et al. (2013). The maintenance of a cline is facilitated by reducing the ratio of diffusion rate to selection intensity, and it is impeded if long-distance dispersal is incorporated into the diffusion model (Nagylaki, 2012,  Su and Nagylaki, in press).

Due to its complexity, multilocus theory is much less developed. This holds both for models with discrete or continuous space. Most investigations make rather restrictive assumptions, such as absence of epistasis or of linkage disequilibrium (LD), or assume two demes. We relax these assumptions and, additionally, provide a comparison of multi-deme models with diffusion models. The available theory for discrete demes is reviewed in  Bürger (2014). Multilocus or quantitative–genetic models with diffusion in a spatially varying environment have been studied by  Slatkin, 1975, Slatkin, 1978,  Felsenstein (1977),  Barton, 1983, Barton, 1999,  Kruuk et al. (1999), and  Hu (2005).

In the present work we consider a quantitative trait that is subject to selection toward a phenotypic optimum in each location. The trait is determined by two diallelic, recombining loci. The diploid sexual population may be subdivided into a finite number of demes or occupy a continuous domain. Mating is random in each location. The phenotypic optimum varies in space. If it is close to the middle of available phenotypes, the trait is under stabilizing selection; if it is close to or at an extreme phenotype, the trait is under directional selection. We impose selection scenarios that differ in the way the optimum changes across space. This change may be gradual, occur in several steps of moderate size, or abruptly in one big step such that there are only two different environments. Such selection scenarios have been discussed in the literature on hybrid zones (e.g.,  Barton, 1999,  Kruuk et al., 1999,  Kawakami and Butlin, 2012). Spatially uniform stabilizing selection is also investigated. We study the following migration patterns: (i) the island model in which migrating individuals reach every deme (island) with the same probability, (ii) stepping-stone models in which individuals migrate either only to next neighbors or to demes in the vicinity such that the probability decreases with distance, and (iii) a diffusion model for a population that inhabits a continuous bounded one-dimensional habitat.

Our main goal is to investigate how the conditions for the existence of polymorphic stationary solutions, or clines, and their properties (e.g., spatial shape) depend on the number of demes, the rate and pattern of migration, the selection scenario, and recombination. Most analysis is dedicated to models with a finite number of demes. However, an essential component will be the comparison of 12-deme models with diffusion models. To make useful comparisons between different patterns or scenarios, the equilibrium configurations and bifurcation patterns are described as functions of the migration rate. The important limiting cases of weak and of strong migration are treated in Section  3. For two demes and loci of equal effect an almost complete mathematical analysis is obtained (Section  4). It complements previous analyses assuming absence of epistasis (Akerman and Bürger, 2014a, Akerman and Bürger, 2014b) or selection on haploids (Geroldinger and Bürger, 2014). The analysis of the two-deme case is not only an important guide to the, mainly numerical, analysis of models with a higher number of demes (Section  5), but also helps to establish analytical results for the island model.

We describe the spatial dependence of the distribution of the trait by its mean phenotype, its genetic variance, and the LD between loci. In Section  6, the properties and shapes of the corresponding clines are compared for the different migration patterns and selection scenarios. Section  7 is dedicated to the comparison of our results with those from previous multilocus analyses of neutral clines (Feldman and Christiansen, 1975, Christiansen, 1986) and analyses of multilocus or quantitative–genetic diffusion models, in particular those of  Slatkin, 1975, Slatkin, 1978,  Felsenstein (1977), and  Barton, 1983, Barton, 1999. In Section  8, our main results are summarized and discussed.

Model

We study a deterministic migration–selection model in which a sexually reproducing, diploid population is subdivided into demes connected by genotype-independent migration. It is assumed that the genotypic fitnesses are uniquely determined by the genotypic value of a quantitative trait. We posit that in each deme , fitness is given by the quadratic function where the phenotypic optimum depends on , and measures the strength of selection. It is assumed that is sufficiently small such that on the range of genotypic values (also called phenotypes). If the optimum is close to the middle of the phenotypic range, the trait is under stabilizing selection in deme ; if it is close to the boundary, it is under directional selection.

The trait is determined additively by two diallelic loci, and , which recombine at rate . We assign the genotypic contributions , , , and to the four alleles , , , and , respectively. The genotypic values of all 16 genotypes are obtained by adding all allelic contributions. Without loss of generality, we use a scale such that . Then the phenotypic range is , the two double homozygotes and have the (extreme) phenotypes −1 and 1, respectively, and all four double heterozygotes have phenotype 0. We restrict the phenotypic optima to this range, i.e., we assume . Finally, we introduce the ratio of locus effects Unless mentioned otherwise we assume , i.e.,  .

The frequencies of the four gametes, , , , , in deme are designated , , , , respectively. The fitness of zygotes consisting of gametes and in deme is , where is the genotypic value. The mean fitness in deme is given by , where denotes the marginal fitness of haplotype in deme .

We assume equivalent sexes, random mating within demes, and that population regulation occurs within each deme (soft selection). We denote linkage disequilibrium in deme by . Then the change of gamete frequencies in deme due to selection and recombination is where , , and denotes LD after selection.

Let denote the backward-migration matrix, i.e.,  denotes the probability that an individual in deme immigrated from deme . After migration random mating and reproduction occur within demes. Therefore, the frequency of the th gamete in deme in the next generation is: Eqs. (2.3) define a discrete dynamical system on the -fold Euclidean product of the simplex

For convenience we introduce the allele frequencies and of alleles and in deme . Then the gamete frequencies are given by the relations

We shall use the notation . See Table 1 for a glossary of symbols.

Table 1

Glossary of symbols. We define the symbols in the main text that occur in more than one paragraph. Roman and Greek letters are listed separately. Uppercase letters precede lower case ones and listing is in order of appearance in the text. The references are to the position of first appearance in the text. Reference (2.1)−, refers to the text above Eq. (2.1), whereas (2.1)+refers to the text below Eq. (2.1).

Symbol	Reference	Definition
A	(2.1)+	First locus
A	(2.1)+	First allele at locus A
a	(2.1)+	Second allele at locus A
B	(2.1)+	Second locus
B	(2.1)+	First allele at locus B
b	(2.1)+	Second allele at locus B
c_1	(2.1)+	Substitution effect at A
c_2	(2.1)+	Substitution effect at B
D_k	(2.3a)−	Linkage disequilibrium in deme k
EkA,∗, EkB,∗	(3.1)	SLPs in deme k for m=0, where ∗∈{0,1}
F_k	(3.1)+	Internal equilibrium in deme k for m=0
G	(2.1)−	Genotypic value of the trait
G¯_k	(5.6)+	Genotypic mean in deme k
I_m(G)	(3.3)+	Weak-migration perturbation of the equilibrium G
I^j	(4.2)−	Internal equilibria (0≤j≤5)
I	(5.1)−	Migration matrix of the island model
M	(2.3c)−	Backward-migration matrix
Mki	(3.1)−	Equilibrium in deme k corresponding to fixation of gamete i
m_kl	(2.3c)−	Probability that an individual in deme k immigrated from deme l
m	(3.2)	Migration rate
m_max	(4.9)	Maximum migration rate below which a stable polymorphic equilibrium can occur
m_st(G)	(4.1)−	Migration rate at which the equilibrium G gets stable for n=2
m_un(G)	(4.1)−	Migration rate at which the equilibrium G gets unstable for n=2
m_ad(G)	(4.1)−	Migration rate at which the equilibrium G gets admissible for n=2
m_na(G)	(4.1)−	Migration rate at which the equilibrium G loses admissibility for n=2
m¯	(5.2)	Rescaled migration rate
m∗X,M(G)	(5.8)+	Migration rate at which the state of G changes, where ∗∈{st,un,ad,na}. X indicates the selection scenario and M the migration matrix
N_i	(3.4)	Partitions of the set of demes {1,…,n}
n	(2.1)−	Number of demes
P_k	(2.1)	Phenotypic optimum in deme k
P	(4.1)−	Phenotypic optimum for two environments (P=−P_1=P_2)
P_c	(4.5)−	Critical value for P for n=2
p_k	(2.4)+	Frequency of allele A in deme k
q_k	(2.4)+	Frequency of allele B in deme k
r	(2.1)+	Recombination rate
S_4	(2.4)	Simplex
S	(5.3)−	Migration matrix of the stepping-stone model
S_2	(5.3)+	Migration matrix of the generalized stepping-stone model
s	(2.1)	Selection intensity
s~	(5.16)−	Rescaled selection intensity
t	(7.4a)	Time
w_k(G)	(2.1)	Fitness of genotypic value G in deme k
w_ij,k	(2.2)+	Fitness of genotype ij in deme k
w_i,k	(2.2)+	Fitness of gamete i in deme k
w¯_k	(2.2)+	Mean fitness in deme k
V_k	(5.6)+	Phenotypic variance in deme k
V_T	(6.2)+	Phenotypic variance in the entire population
x_i,k	(2.2)+	Frequency of gamete i in deme k
x_k	(2.5)+	Vector of gamete frequencies in deme k
y	(7.4a)	Spatial variable in a continuous domain
κ	(2.2)	Ratio of locus effects
η_i	(2.3a)	Constants
σ^2	(7.4a)	Diffusion rate in a continuous domain
^(s)	(2.3a)	Indicates haplotype (or gene) frequencies after selection and recombination
^^′	(2.3c)	Indicates haplotype (or gene) frequencies in the next generation
ˆ	(3.1)+	Indicates an equilibrium value

Limiting cases

We determine equilibria and their stability properties analytically for the limiting cases of no migration, weak migration, and strong migration.

No migration

For panmictic populations the model has been analyzed previously (reviewed in  Bürger, 2000, Chap. VI.2). We recapitulate the relevant results. Because in the absence of migration the dynamics of the demes are decoupled, we describe the equilibrium configuration for a single deme .

Three types of equilibria may exist: (i) monomorphic equilibria, (ii) single-locus polymorphisms (SLPs), and (iii) fully polymorphic equilibria. The monomorphic equilibrium corresponding to fixation of gamete in deme is denoted by . Four SLPs, corresponding to the fixation of one allele at one locus, exist. Their coordinates are where the superscript or of indicates the polymorphic locus, and the superscript 0 or 1 which allele is fixed at the other locus. The hat, , signifies an equilibrium.

The equilibria and are admissible if and only if ; and are admissible if and only if . Stability conditions of all boundary equilibria are available for arbitrary locus effects (Bürger, 2000). Because , the SLPs are asymptotically stable when they are admissible. If or , the trait is under directional selection and or , respectively, is globally asymptotically stable (Appendix A.1). If , then and are simultaneously asymptotically stable. If (and ), there exists a unique internal equilibrium which is always unstable (Appendix A.1).

These results show that, depending on the phenotypic optimum , one of the following five qualitatively different equilibrium configurations occurs.

(i) If , is globally asymptotically stable.

(ii) If , and are simultaneously asymptotically stable and is unstable.

(iii) If , and are simultaneously asymptotically stable and is unstable.

(iv) If , and are simultaneously asymptotically stable and is unstable.

(v) If , is globally asymptotically stable.

Weak migration

We apply the perturbation theory developed by  Karlin and McGregor, 1972a, Karlin and McGregor, 1972b to infer existence and local stability of equilibria for weak migration from the model with no migration. The migration matrix is supposed to satisfy where the are constants that satisfy .

If , the dynamics (2.3) on is given by the Euclidean product of the single-deme dynamics on . Every equilibrium is of the form where is an equilibrium in deme . If all components of the equilibrium (3.3) are identical, i.e.,  for , we denote the equilibrium (3.3) by .

If in the absence of migration every equilibrium is hyperbolic, perturbation theory shows that the following holds for sufficiently small (Karlin and McGregor, 1972b): (i) in the neighborhood of each asymptotically stable equilibrium for , there exists exactly one equilibrium for and it is asymptotically stable; (ii) in the neighborhood of each unstable internal equilibrium for , there exists exactly one equilibrium for and it is unstable; (iii) in the neighborhood of each unstable boundary equilibrium for , there exists at most one equilibrium for , and if it exists, it is unstable. If we denote the perturbation of by , then as . The proof of Theorem 4.1 in  Karlin and McGregor (1972b) shows that an equilibrium may leave the state space after perturbation only if it is transversally unstable.

The following proposition combines the results for panmictic populations summarized in Section  3.1 with the perturbation theory outlined above.

Assume   (2.3), let be sufficiently small, and define the setsThe following asymptotically stable equilibria exist:whereThe following are unstable equilibria:whereand

The proof is given in Appendix A.2.

We note that not all unstable equilibria are given by (3.6): e.g., if for some , then and  are unstable. In general, it has to be checked separately whether the perturbation of an unstable boundary equilibrium leaves the state space. Notably, Proposition 3.1 holds independently of the migration matrix .

The values and are excluded in the above Proposition, because then not every equilibrium is hyperbolic and separate treatment is needed. Numerical work suggests that Proposition 3.1 remains valid if is added to , and are added to , and is added to .

We conjecture that almost all trajectories converge to one of the equilibria in (3.5) if is small. If for every either or holds, this conjecture follows from global perturbation theory (Bürger, 2009, Section 5; or  Bürger, 2014, Theorem 7.7 and Remark 7.8). The reason is that if and , there is global convergence to an asymptotically stable equilibrium (Section  3.1); if and , the Lyapunov function establishes (exponential) convergence to for every trajectory with , and to for every trajectory with ; trajectories satisfying converge to . In both cases the convergence patterns persist for small .

Strong migration

With the special migration schemes of Section  5 in mind, we assume an even number of demes and posit that selection and migration satisfy the following symmetry conditions for every : These conditions describe a mirror symmetry between demes and , such that in deme selection acts on the haplotypes , , , in the same way as selection in deme on , , , , respectively. In particular, selection on the trait occurs in opposite direction in these two sets of demes.

Assume   (3.7).

1.  If migration is sufficiently strong, i.e., and are sufficiently large, then and are asymptotically stable and no other equilibrium is stable. The equilibrium exists, is the unique internal equilibrium, and is unstable.

2.  The critical migration rate at which and become asymptotically stable, denoted (Section   5.2), is independent of . and are stable if and only if they are stable with respect to their marginal one-locus systems.

The proof is given in Appendix A.3. There it is also shown that if , the migration rates at which and become stable depend on .

The analysis below will show that the equilibrium configuration of the strong-migration limit, as defined by Proposition 3.4.1, may apply only for much larger migration rates than .

Two demes

In this section we assume and (3.7), i.e.,  and . Then the strength of divergent selection between the demes increases with increasing . Our goal here is to describe the equilibrium configurations and bifurcation patterns as the migration rate increases.

We find the equilibria of (2.3) by using the algorithm NSolve of Mathematica (Wolfram Research, Inc., 2010) and determine their local stability properties by calculating the eigenvalues numerically. Global stability results are inferred from forward iterations of (2.3). They were performed with Mathematica and the following adjustments: In each deme, 1000 initial values from the interior were chosen as , where the are independent and uniformly distributed in . Iterations were stopped if the Euclidean distance between successive values declined below . Two equilibrium values were considered as equal if their Euclidean distance (in ) was less than .

In combination with our analytical results for weak and for strong migration, we obtain a presumably complete classification of bifurcations in which the stable equilibria are involved.

For any equilibrium , we designate by or the critical migration rate at which becomes stable or unstable, respectively, as increases above this value. Analogously, we write or for the critical migration rate at which gains or loses admissibility, respectively. These critical migration rates turn out to be unique.

Proposition 3.4 shows that in the limit of strong migration, and are simultaneously stable. A linear stability analysis of and reveals that these equilibria are stable if and only if , where We note that is independent of (cf. Section  3.3) and if and only if .

From one-locus theory (Karlin and Campbell, 1980, Bürger, 2014) and our symmetry assumptions ( and ), we infer that four SLPs exist if and only if all monomorphic equilibria are unstable. Otherwise, no SLP exists. The allele frequency at an SLP is a zero of a cubic polynomial which does not have simple form. Numerical investigations suggest that the SLPs are always unstable. (They are stable within their marginal one-locus system but unstable with respect to the interior of the state space). They play no role in the further analysis.

At several instances we define internal equilibria by weak-migration perturbations, e.g.  . Then we use the notation for the whole range of parameters where this equilibrium exists. The following equilibrium plays a central role in the subsequent analysis Its coordinates are continuous in since and as (Appendix A.1).

Because Proposition 3.1 shows that the equilibrium configuration for weak migration depends on , we distinguish three cases according to increasing strength of divergent selection.

Case I

Let . Then there is stabilizing selection in each deme, and divergent selection between demes is weak. According to (4.1), and are asymptotically stable for every . For sufficiently weak migration, and are the only internal stable equilibria (Proposition 3.1). The equilibria , , , , and are admissible and unstable if . As increases, the following three bifurcations1 occur which reduce the number of equilibria and, eventually, yield the equilibrium configuration of the strong-migration limit.

The equilibrium collides with the two unstable equilibria and in a subcritical pitchfork bifurcation in which loses its stability but persists, and the unstable equilibria are annihilated. Analogously, collides with the two unstable equilibria and in a subcritical pitchfork bifurcation. The value at which and loose their stability is denoted by . At the value , a third subcritical pitchfork bifurcation occurs in which the three unstable internal equilibria , , and collide, and are annihilated, and remains admissible and unstable. In this case the sequence of bifurcation points is (Fig. 1a). If , the equilibrium configuration of the strong-migration limit applies.

Fig. 1

Bifurcation patterns for two demes. The functions , , , and provide two-dimensional projections of the six-dimensional coordinates and are given in Appendix A.5. Solid and dotted lines represent stable and unstable equilibria, respectively. The equilibrium is displayed in green, the equilibria and are displayed in orange, and and in red. Gray dotted lines in panel a show the equilibria , , , and , whereas in the other panels gray lines show , , , . Panels d and f are zoomed-in versions of panels c and e, respectively. In Case II.a, the bifurcations can occur in different orders; see (4.6). The SLPs are not shown because they are always unstable and bifurcate only with the monomorphic equilibria when they leave the state space. The asymmetries in panels d and f result from the nonlinear projections and , respectively. Parameters are and . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

For each of the equilibria , , and , the equilibrium allele frequencies and LD in demes 1 and 2 satisfy the symmetry relation Eqs. (A.10), (A.10) and (A.11) in Appendix A.4 provide approximations for , , and by assuming weak evolutionary forces and linkage equilibrium (see Fig. B.13).

Case II

Let . Then there is (asymmetric) stabilizing selection in each deme, and divergent selection between demes is moderately strong. We recall from Section  3.1 that if , the equilibria and ( and ) are simultaneously asymptotically stable in deme 1 (deme 2). Additionally, there is the unstable internal equilibrium in each deme. If migration is weak, there are nine internal equilibria. Among them, , , , and are asymptotically stable. The definitions of and extend those in Case I, because and converge to , and and converge to as (Section  3.1).

As the migration rate increases, the stable equilibrium collides with the two unstable equilibria and in a subcritical pitchfork bifurcation, i.e.,  becomes unstable and and are annihilated. Analogously, collides with the two unstable equilibria and in a subcritical pitchfork bifurcation. Both bifurcations occur at the same migration rate . As the migration rate increases further, a third subcritical pitchfork bifurcation occurs at which the three unstable equilibria , and collide, and are annihilated and remains admissible and unstable (Figs. 1d,f).

For larger we distinguish two subcases, depending on whether the two stable equilibria and do or do not collide with . The equilibria and collide if and only if , where is an increasing function of which cannot be calculated explicitly. In the continuous-time approximation (Appendix A.4, Fig. B.13), however, the coordinates of and can be calculated (A.13), and is given by

Case II.a

If , the two stable equilibria and do not collide and leave the state space at the migration rate by transcritical bifurcations with the two boundary equilibria and , respectively. The equilibrium is unstable for all migration rates. The bifurcations can occur in three different orders: The first order is the most common (Fig. 2) and is displayed in Figs. 1c,d. Above the highest indicated bifurcation point, the equilibrium configuration of the strong-migration limit applies.

Fig. 2

Panel a shows the regions of stability of the equilibria , , , , and as a function of . The red line shows and ; the black line (4.1). The orange line shows and is obtained by numerical calculation of the bifurcation point. Panel b shows the fraction of trajectories converging to one of the four simultaneously stable equilibria if (Case I). Initial values were chosen as described at the beginning of Section  4. In both panels, and . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Case II.b

If , a supercritical pitchfork bifurcation occurs at when the three equilibria , , and collide. The equilibria and annihilate each other and becomes stable. At the critical migration rate , a second supercritical pitchfork bifurcation occurs, in which becomes unstable and and are re-established. As the migration rate increases further, the stable equilibria and leave the state space at by transcritical bifurcations with and , respectively. The sequence of bifurcation events is given by (Figs. 1e,f).

Case III

Let . Then there is directional selection in each deme, and divergent selection between demes is strong. If migration is weak, is the globally attracting internal equilibrium. At , the equilibrium becomes unstable and the two stable equilibria and are established in a supercritical pitchfork bifurcation. As the migration rate increases, and leave the state space at by transcritical bifurcations with and , respectively. The sequence of bifurcation events is given by (Fig. 1b).

If , Case I applies, and if , Case III applies. Because they are degenerate (Remark 3.2), they require separate treatment.

The above results are related to our previous work (Geroldinger and Bürger, 2014), where we mainly studied a haploid model and explored the influence of unequal locus effects () and of the recombination rate on the maximum migration rates admitting polymorphism. Case I exhibits the same bifurcation pattern as Pattern I.sr.0 in  Geroldinger and Bürger (2014). Case II does not have an analogue in the haploid model. Since the SLPs are not admissible in the haploid model for , at most two internal equilibria can be stable, whereas in the present diploid model the four internal equilibria , , , and may be simultaneously stable. Case III is identical to Pattern D.sr.1 with in  Geroldinger and Bürger (2014). The equilibria and correspond to and in  Geroldinger and Bürger (2014).

If , then for every stable equilibrium and every . For weak migration, approximations of the internal equilibria show that and (if ), (if ), and (if ). For the continuous-time model it can be shown that for all migration rates; see (A.13). Therefore, if migration is weak, at every stable equilibrium we have if (stabilizing selection) and if (directional selection). Numerical work suggests that this also holds for intermediate migration rates. For the haploid model, it could be proved that in the case analogous to Case I (i.e., Pattern I.sr.0 in  Geroldinger and Bürger, 2014), LD at and is negative for all migration rates.

The above analysis shows that for sufficiently strong divergent selection () there is an interval of migration rates for which a unique asymptotically stable internal equilibrium () exists which, presumably, is globally attracting. This interval increases with (Fig. 2a) and includes 0 if . If , there are always multiple simultaneously stable equilibria. For Case I, Fig. 2b shows the fraction of trajectories converging to one of the stable equilibria , , and as a function of the migration rate.

The maximum migration rate up to which a stable polymorphic equilibrium can occur is given by In Fig. 2a these migration rates are displayed as functions of . Whereas is increasing with the recombination rate (results not shown), is independent of (4.1). Therefore, the critical ratio above which the equilibrium configuration of the strong-migration limit applies is independent of if is sufficiently large.

Island and stepping-stone models

In this section we investigate the influence of the migration pattern, of the number of demes, and of different selection scenarios on the equilibrium configurations. In particular, we shall compare migration patterns exhibiting different degrees of mixing and different degrees of isolation by distance. Our selection scenarios include models in which there is one major step-like change in the environment, models in which the environment changes (more) gradually, and a model with uniform stabilizing selection.

Migration patterns and selection scenarios

We investigate the island model and two stepping-stone models. Whereas the former has no geographic structure, the latter exhibit isolation by distance. Two versions of the stepping-stone model will be considered. In the first, individuals migrate only to neighboring demes, whereas in the second migration to more distant demes, or islands, is admitted but occurs with decreasing probability.

The island model

The (forward and backward) migration matrix of the island model is given by Proposition A.1 demonstrates that in each deme the coordinates of the equilibria depend on only through the position of the optimum . This holds for every choice , , and , and confirms that the island model exhibits no spatial structure. Also the following relation between an island model with an even number of islands to the two-deme model is notable.

If demes have optimum and demes have optimum , then for every equilibrium in the two-deme model with migration rate there is an equilibrium in the island model with migration rate (Proposition A.2). This rescaling of the migration rate is a consequence of the following argument. In the two-deme model, denotes the probability that an individual breeds in the other deme. This coincides with the probability that an individual migrates to a deme with a different environment. In the island model, the second interpretation of does not hold if . Instead, the probability of switching the selective environment is . Therefore, (5.2) transforms critical migration rates at which the equilibrium structure changes for the two-deme model to analogous critical migration rates for the island model. It is useful even if because spatially heterogeneous equilibria may exist in the two-deme model (e.g.,  , ).

Stepping-stone models

The backward-migration matrix of the (single-step) stepping-stone model is given by In this migration pattern individuals can migrate only to neighboring demes. Alternatively, we consider a generalized stepping-stone model, where migration to more distant demes is possible but its rate decreases with distance. The matrix of this generalized stepping-stone model is given in Appendix A.7 for and . In all our migration patterns, may be interpreted as the probability of outbreeding.

Obviously, the statements of Proposition A.1, Proposition A.2 do not apply to the stepping-stone models. In the stepping-stone models an increasing number of demes increases isolation by distance. Therefore, equilibrium frequencies change gradually in space even if the environment changes sharply.

In Section  7, we will compare our results on the stepping-stone models to previous investigations using diffusion approximations in continuous time and space in an unbounded domain.

Selection scenarios

For each of the migration patterns, we consider the following selection scenarios (Fig. 3): Scenario A models a sharp change, or single step, in the phenotypic optimum from −1 to 1. In one half of the demes () genotype is the best adapted, whereas is the best adapted in the other half (). Scenario B assumes two steps in the phenotypic optimum, from −1 to 0 and from 0 to 1. Therefore, heterozygotes and the (repulsion) genotypes and are selectively favored in the center of the domain. Scenario C assumes that the phenotypic optimum changes linearly in space, which ensures that each genotype is well adapted in some deme if the number of demes is large enough. In Scenario D there is uniform stabilizing selection toward in all demes.

Fig. 3

The selection scenarios (5.4) for and .

Simple calculations or a glance at Fig. 3 reveal that for fixed selection intensity , the maximum fitness difference between genotypes in each deme, , varies among the selection scenarios.

Recalling (3.1), we note that except for Scenario C the relative sizes of the sets () are independent of . In Scenario A, we have and ; in Scenario B, and hold; and in Scenario D, . However, in Scenario C we have and if , but , , and if . This fact is responsible for some peculiar dependencies of critical migration rates on in Scenario C. To avoid this phenomenon in Scenario B, and also to keep the number of demes even, we assumed .

Equilibrium configurations

We start by noting that the migration patterns (5.1), (5.3), (A.22), and the selection scenarios (A.38) satisfy (3.7). Therefore, the following proposition follows immediately from (2.3).

1.  Equilibria that do not satisfyoccur in pairs, and . Each pair satisfies the relationsThe equilibria of each pair have the same stability properties and satisfy and , where and denote the mean genotypic value and the genetic variance in deme .

2.  Equilibria that do not satisfyoccur in pairs, and . Each pair satisfies the relationsThe equilibria of each pair have the same stability properties, the same mean genotypic value, and the same variance.

Whereas the first statement also holds if , the second statement requires . Eq. (5.5) generalizes (4.4). In (A.23), the equilibria , , , , and are defined for weak migration and demes. As in the two-deme model, the equilibria , , and satisfy (5.5). In addition, fulfills (5.7), and and fulfill (5.8). The equilibria and satisfy (5.6) and (5.8).

The coordinates of the stable equilibria were calculated from forward iterations of (2.3) (see Section  4). Because several equilibria lie on the manifold given by the symmetry relation (5.5), their coordinates could be computed efficiently by iteration of (2.3) on this manifold. Local stability was determined by numerical evaluation of the eigenvalues of the Jacobian of (2.3).

For increasing migration rate, the number of stable equilibria decreases from its usually high value for weak migration (Proposition 3.1). The numerical computations suggest that, in close analogy to the two-deme model, the reduction of internal equilibria is always due to pitchfork bifurcations.

In this section we investigate the number of stable internal equilibria and the migration rates at which the bifurcations occur. These migration rates depend on the migration pattern, the selection scenario, the number of demes , as well as on and . They are denoted by , where indicates whether the equilibrium changes admissibility or stability (as in Section  4), indicates the selection scenario, and the migration matrix.

Numerical work suggests the following: Internal equilibria never enter the state space through the boundary, and SLPs are never stable if . There is at least one internal equilibrium () satisfying (5.5). Proposition 3.1 implies that for every migration pattern , we have if and . In Scenarios A, B, and C, the equilibrium configuration of the strong-migration limit applies if .

Scenario A

Proposition 3.1 implies that for weak migration there is a unique stable internal equilibrium which we denote by (A.23a). In the absence of migration every trajectory converges to (if ) or (if ) (Section  3.1). Therefore, is globally asymptotically stable for weak migration (Section  3.2). As the migration rate increases, the equilibrium becomes unstable and two stable equilibria and are established in a supercritical pitchfork bifurcation at . The equilibria and leave the state space through and , respectively, at . Therefore, the bifurcation pattern is analogous to that of Case III in the two-deme model (Fig. 1b). The critical migration rates and depend on the number of demes , the migration pattern , and the selection intensity ; the former depends also on .

For the island model, the migration rates and are obtained from the two-deme model by rescaling according to (5.2): where and are the critical migration rates from the two-deme model; see (4.1), (4.1) and (A.14), respectively.

The migration rates and in the stepping-stone models cannot be determined analytically and are evaluated numerically in Table A.1 and Fig. 4. It is important to note that for large they may exceed (our maximum migration rate). Indeed, in the stepping-stone models rather small is required such that . Then always seems to hold for . Therefore, in the island model the equilibrium configuration of the strong-migration limit, hence a homogeneous population, is reached at lower migration rates than in the stepping-stone models. The reason is that short-range migration has a weaker homogenizing effect than distance-independent migration. For the same reason, , , , and increase with the number of demes in the stepping-stone models. Thus, both sets of inequalities support the notion that increasing isolation by distance facilitates the maintenance of genetic variation.

Table A.1

Critical migration rates for different selection scenarios and migration patterns. The symbol ‘-’ indicates that the critical migration rate does not exist. The symbol ‘*’ indicates that the migration rate would exist for smaller selection intensities . Equilibrium configurations were calculated in steps of . The values indicate the lowest migration rate for which a different equilibrium configuration was observed. The data for and are visualized in Fig. 4.

	Island				Generalized stepping-stone				Stepping-stone
	A	B	C	D	A	B	C	D	A	B	C	D
s=0.1,r=0.5
n=6
munX,M(I^2,3)	0	0.009	0.008	0.008	0	0.008	0.007	0.008	0	0.007	0.007	0.035
mstX,M(I^1)	0	0.021	0.029	–	0	0.029	0.058	–	0	0.047	–	–
munX,M(I^1)	0.220	0.138	0.063	–	*	0.396	0.159	–	*	*	–	–
mstX,M(M^2,3)	0.245	0.187	0.140	0	*	*	0.378	0	*	*	*	0
n=12
munX,M(I^2,3)	0	0.010	0.005	0.008	0	0.010	0.010	0.042	0	0.014	0.015	0.203
mstX,M(I^1)	0	0.023	–	–	0	0.044	0.117	–	0	0.163	–	–
munX,M(I^1)	0.242	0.152	–	–	*	*	0.303	–	*	*	–	–
mstX,M(M^2,3)	0.269	0.206	0.133	0	*	*	*	0	*	*	*	0
s=0.2,r=0.5
n=6
munX,M(I^2,3)	0	0.018	0.015	0.014	0	0.016	0.014	0.035	0	0.013	0.014	0.067
mstX,M(I^1)	0	0.041	0.056	–	0	0.058	0.111	–	0	0.095	0.269	–
munX,M(I^1)	0.372	0.262	0.131	–	*	*	0.337	–	*	*	0.382	–
mstX,M(M^2,3)	0.403	0.335	0.265	0	*	*	*	0	*	*	*	0
n=12
munX,M(I^2,3)	0	0.020	0.013	0.015	0	0.019	0.018	0.08	0	0.027	0.028	0.393
mstX,M(I^1)	0	0.045	–	–	0	0.087	0.225	–	0	0.332	–	–
munX,M(I^1)	0.409	0.289	–	–	*	*	*	–	*	*	–	–
mstX,M(M^2,3)	0.444	0.369	0.254	0	*	*	*	0	*	*	*	0
s=0.2,r=0.01
n=6
munX,M(I^2,3)	0	0.007	0.005	0.003	0	0.005	0.004	0.006	0	0.003	0.003	0.009
mstX,M(I^1)	0	0.039	0.050	–	0	0.058	0.111	–	0	0.095	0.288	–
munX,M(I^1)	0.372	0.264	0.138	–	*	*	0.332	–	*	*	0.352	–
mstX,M(M^2,3)	0.403	0.335	0.265	0	*	*	*	0	*	*	*	0
n=12
munX,M(I^2,3)	0	0.007	0.003	0.003	0	0.005	0.003	0.008	0	0.006	0.006	0.04
mstX,M(I^1)	0	0.042	–	–	0	0.087	0.226	–	0	0.333	–	–
munX,M(I^1)	0.409	0.290	–	–	*	*	*	–	*	*	–	–
mstX,M(M^2,3)	0.444	0.369	0.254	0	*	*	*	0	*	*	*	0

Fig. 4

Intervals of the migration rate in which the equilibrium configurations of the various selection scenarios and migration patterns occur. Colors indicate the equilibrium configurations. Orange: more than two equilibria are stable (). Red: the two internal equilibria and are stable. Green: is globally stable (). Gray: and are stable. The numbers give the critical migration rate at which the corresponding configuration emerges, provided it is non-zero. The recombination rate is . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Scenario B

In Scenario B, in which there are three environments (5.4b), we assume . The stable internal equilibria for weak migration are obtained from Proposition 3.1 and are given in (A.24). There are such equilibria. If , the number of stable equilibria quickly reduces to four (the number of stable equilibria if ) as increases from zero. These four equilibria are , , , and (A.23).

Numerical work suggests that holds always. Therefore, the number of stable internal equilibria is greater than or equal to four if . If is slightly larger than , and are the only stable equilibria. Except for and all stable internal equilibria get annihilated by bifurcations with unstable internal equilibria.

The equilibria and may either leave the state space through and , respectively, or collide with the internal unstable equilibrium (A.23a). The first case is analogous to Case II.a in the two-deme model (Fig. 1c) and the second case is analogous to Case II.b in the two-deme model (Fig. 1e). If and collide with , the two stable equilibria and are annihilated and becomes stable. As the migration rate increases, gets unstable, and the two stable equilibria and are re-established. Finally, and leave the state space by transcritical bifurcations with and , respectively, at .

In the island model, the migration rate can be calculated using an argument analogous to that for Scenario A by invoking Proposition A.2.2: where is derived from a linear stability analysis in the three-island model with , , and . The scaling factor of arises because in Scenario B the probability of switching the selective environment is .

Comparing the numerically evaluated critical migration rates in Table A.1 for the different migration patterns, we observe that, in addition to (5.10), (5.11) with , holds for ; see also Fig. 4.

For the critical migration rate , both (Fig. 4c, Table A.1, ) and (Fig. 4d, Table A.1, ) may hold. Therefore, in contrast to , , (see (5.10), (5.11), (5.13)), is not necessarily increasing with isolation by distance. The source of this ambiguous dependence is the following. On the one hand, strong migration homogenizes the spatial genetic differences and depletes genetic variation (this effect is determining all other critical migration rates, which are higher). On the other hand, immigrants from demes with different selective environments aid within-deme variation. The second effect becomes very weak with increasing isolation by distance because neighboring demes tend to have the same environment. It is weak if , but it is dominating if .

Scenario C

In Scenario C, the environment changes steadily (5.4c). For sufficiently weak migration the stable equilibria are given by (A.25) and their number by (A.26), which gives for , 12, respectively. The qualitative dependence of the equilibrium configurations on is similar to Scenario B, except that for very small there are more equilibria. However, the bifurcation pattern corresponding to Case II.a of the two-deme model occurs much more often than that of Case II.b. Fig. 4 and Table A.1 also show that in several cases, never becomes stable (eg., the green region is missing in Fig. 4f). Finally, the inequalities (5.10), (5.11), and (5.13) hold for if the corresponding migration rates are between 0 and , which is not always the case.

Fig. 4f, shows that the migration pattern may affect the establishment of a globally attracting equilibrium in a non-intuitive way. Whereas the equilibrium becomes stable for , it does not for or (see also Fig. B.1 in Appendix B, Online Supplement).

Scenario D

If there is uniform stabilizing selection toward 0, there are stable equilibria for weak migration of which are internal; and are stable for every ; see Proposition 5.2 and (A.27). The equilibria and do not exist and is never stable. In a series of pitchfork bifurcations, these stable internal equilibria are reduced to the stable internal equilibria and , which are obtained from (A.23) with . Similar to Case I of the two-deme model (Fig. 1a), the four equilibria , , , and are stable up to . As in that case, the equilibrium configuration of the strong-migration limit applies if . Thus, in contrast to Scenarios A, B, and C, the strong-migration limit does not apply for every . For the island model we infer the critical migration rate above which no (stable) polymorphism is possible from the two-deme model (Remark 5.1): where is the corresponding migration rate in the two-deme model. In contrast to Scenarios B and C, the influence of the different migration patterns on is simple, i.e.,  holds (Figs. 4g,h).

Comparison and summary

If migration is sufficiently weak, the equilibrium configuration depends on the number of demes and the selection scenario, but is independent of the migration pattern (Proposition 3.1). For sufficiently strong migration, the equilibrium configuration of the strong-migration limit applies (Proposition 3.4). It is independent of the migration pattern, the number of demes, and the selection scenario. The equilibrium configurations in the parameter range where migration and selection are intermediate can be described with the help of the critical migration rates , , , and . They partition the interval in up to five parts:

(i) If , more than two equilibria (internal or monomorphic) are stable;

(ii) if , two internal equilibria (, ) are stable;

(iii) if , one internal equilibrium () is stable;

(iv) if , two internal equilibria (, ) are stable;

(v) if , the monomorphisms and are stable.

Fig. 4 displays these intervals for every selection scenario and migration model. For Scenarios B and C, all five types may occur; for Scenario A only (iii), (iv), and (v) occur; for Scenario D only (i) and (v) occur.

Comparison of the first bar with the second and third in panels a–f of Fig. 4 shows that in the stepping-stone models genetic variation is lost (gray regions) for higher migration rates than in the island model. Often these migration rates exceed 0.5; then there is no gray region. Clearly, this reflects the fact that gene flow has a stronger homogenizing effect in the absence of isolation by distance than in its presence. To demonstrate the ubiquity of this finding and to compare patterns and in more detail, we proceed as follows. Start with for given , , , and X. Denote by the selection intensity in the corresponding stepping-stone model such that , i.e., such that with and the transition to the strong-migration limit occurs at the same as with and . In particular, and get stable at the same migration rate in both migration patterns.

The fourth bar in panels a–f of Fig. 4 shows the intervals in which the different equilibrium configurations occur for the stepping-stone model with selection intensity . By definition of , the gray regions occur above the same migration rate as for the island model (first bar). Comparison of the first and the fourth bar in panels a–f shows that in the island model with selection intensity , the regions where is stable are larger than in the stepping-stone model with selection intensity . This appears to reflect the greater importance of initial conditions in migration patterns involving isolation by distance. The number of demes has two effects on the equilibrium configuration. First, for weak migration, the number of stable internal equilibria increases with the number of demes if . Second, in the stepping-stone models the degree of isolation by distance increases with . Therefore, critical migration rates in the stepping-stone models increase with (Fig. 5). For , the role of is well understood in Scenarios A and B; see (5.9) and (5.12), (5.9) and (5.12). However, in Scenario C, decreases from for to for . This is due to the variation of the relative sizes of (3.1) as explained below Eq. (5.4).

Fig. 5

The migration rate shown as a function of for Scenario A (red), B (green), and C (blue). The selection intensity is . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

A comparison of Scenarios A, B, and C shows that the parameter range where is stable decreases from A to C, i.e.,  where (in Fig. 4, compare the green stacks among panels a, c, and e, as well as among b, d, and f). Also the maximum migration rate below which polymorphism is possible, decreases from Scenario A to Scenario D, i.e.,  (Fig. 4, Table A.1). In this sense, a single abrupt change in the environment is more favorable for the maintenance of genetic variation than a more gradual change.

For weak migration the number of coexisting stable internal equilibria increases from Scenario A to D. Hence, in a more gradually changing environment, initial conditions affect evolution much more than in an environment that changes abruptly.

Similar to the two-deme model (Remark 4.5), is increasing in if ; see Table A.1. The reason is that LD at the equilibria and is negative in the demes under stabilizing selection (Eq. (A.34), Fig. 8c). Therefore, more recombination increases genetic variance because it reduces the negative LD (Bürger, 2000, p. 74). The critical migration rates and may increase or decrease with but depend only very weakly on (Table A.1). However, is independent of (Proposition 3.4).

Fig. 8

Clines in the mean phenotypic value (a), in the genetic variance (b), and in LD (c) at simultaneously stable equilibria in the stepping-stone model with Scenario B and . The parameters are , , and which is smaller than . Ten equilibria (out of the 16 for weak migration) are stable. Five pairs exhibit different mean, variance and LD. Each color corresponds to a pair of simultaneously stable equilibria. The red dashed lines show and and correspond to the red lines in Figs. 6d,e,f. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Clines in the mean phenotype, genetic variance, and LD

Here we investigate how the migration patterns and selection scenarios determine the spatial distribution of the population across demes. In particular, we are interested in how mean phenotype, genetic variance, and LD vary in space. We focus on the range and briefly treat the case of very small migration rates, when four or more equilibria may be stable simultaneously, further below. For every migration pattern and selection scenario, as well as for representative values of and , we calculated the mean phenotype, the genetic variance, and the measure of LD at equilibrium in every deme and displayed them as functions of the deme number . This was done for a fine grid of admissible migration rates (Figs. B.2–B.7 in Appendix B, Online Supplement). Fig. 6 displays representative results for one migration rate.

Fig. 6

Clines in the mean phenotype (left), the genetic variance (middle), and LD (right) for different selection scenarios and migration patterns. Blue lines indicate the island model, red lines the stepping-stone model, and green lines the generalized stepping-stone model. Magenta lines show the stepping-stone model with . The corresponding equilibrium configuration for each migration patterns can be inferred from Fig. 4. Solid lines indicate that is the unique stable equilibrium, whereas dashed lines indicate that and are simultaneously stable (they exhibit the same mean, variance, and LD). In the left column, dots mark the positions of the optimum. The parameters are , , , and . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Clines in the mean phenotype and local adaptation

The panels in the left column of Fig. 6 display the clines in the mean phenotype. The degree of local adaptation in deme is measured by . In most cases, the stepping-stone models favor local adaptation compared to the island model (Figs. B.2–B.4). However, the relation is valid in every deme only in Scenario A. In Scenario B, it is the island model that maximizes local adaptation in the demes under stabilizing selection because if ; cf.  (A.32), (A.35). However, it leads to poor adaptation in demes under directional selection (Fig. 6d). For Scenario C, counter examples to (6.1) occur in demes with stabilizing selection, e.g., if , , , , and (Fig. B.4).

Genetic variance

In a step environment (Scenario A) and with stepping-stone migration (, ), the within-deme variance is always maximized in the center of the cline (Fig. 6b), and it decreases toward the boundaries. For Scenarios B and C, this does not hold: the variance may be maximized in the center or elsewhere (Figs. 6e,h). The bimodal patterns occur mainly for weak single-step migration (Figs. B.3, B.4). In the absence of migration, different haplotypes are fixed in demes with directional or stabilizing selection. Therefore, in Scenario B weak migration induces substantial variance in the demes adjacent to an environmental change. In Scenario C, the following arguments show that for weak migration the variance is bimodal. If , then for , and otherwise (Section  3.1). Therefore, if migration is weak, the variance in the demes with is higher than in the demes in the center of the range or close to the boundary.

If in Scenarios B or C, migration rates are such that is the unique stable equilibrium, i.e.,  (whence migration is no longer weak), the genetic variance decreases from the center of the cline to its boundaries (Figs. B.3, B.4).

For the island model, the genetic variance is either spatially uniform (Scenario A) or weakly dependent on space (Scenarios B and C). In the latter case, it may be maximized or minimized in the center, or it may be bimodal (Figs. B.3, B.4).

Although Fig. 6 suggests the simple relation for the variances (in deme ) maintained by the three migration patterns, it does not hold in general. Obviously, (6.2) is violated if , but it may also be violated if (Figs. B.2–B.4).

Finally we consider the genetic variance in the entire population. It is calculated from the spatially averaged gamete frequencies and displayed in Fig. 7 as a function of . Whereas in Scenario A, is monotone decreasing in , weak migration may increase in Scenario B and Scenario C. The variance decreases rapidly when approaches , i.e., when the cline starts to collapse. For given , decreases from Scenario A to B to C. Further, the effect of linkage on decreases from Scenario A to B to C because the absolute magnitude of LD decreases from Scenario A to B to C (Section  6.3).

Fig. 7

The genetic variance in the entire population as a function of the migration rate. The island model (blue) is shown for (dark) and (light). Red and orange lines show the stepping-stone model for and , respectively. Green and magenta lines display the generalized stepping-stone model and the stepping-stone model with , respectively (). At dashed lines, equilibria are simultaneously stable. For reasons of visibility only is shown for , whereas at the other stable equilibria is not displayed. Parameters are and . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Linkage disequilibrium

Linkage disequilibrium depends strongly on the selection scenario, the migration pattern, and the spatial location. In Scenario A, the situation is simple. For the stepping-stone models, assumes its maximum in the center of the cline and decays monotonically to a very small positive or negative value at the boundaries (e.g.,  Fig. 6c). A similar pattern was reported by  Slatkin (1975), who modeled dispersal in continuous space by diffusion and assumed nonepistatic directional selection at every location. At the boundaries of the cline, LD may be negative (). This peculiar phenomenon is likely due to the fact that in the demes at the boundary, migration is unidirectional. In an infinite domain, LD will approach zero in increasingly distant demes. For the island model with weak migration, LD is positive and the same in all islands (A.28).

In Scenarios B and C, LD may be a complicated function of the distance from the center (Figs. 6f,i). It tends to be positive in some demes and negative in others. In Scenario B with the stepping-stone models and weak to moderate migration (Fig. 6f), LD is maximized in demes and , which are the demes under directional selection next to the environmental step. For higher migration rates, LD is usually maximized in the center of the cline (Figs. B.3, B.4). In Scenario B with , LD assumes the same positive value in all demes under directional selection and the same positive or negative value elsewhere.

In Scenario C with stepping-stone migration, each deme is close to linkage equilibrium for a wide range of migration rates (Fig. B.4). However, the island model exhibits deviations from linkage equilibrium. They are not negligible if linkage is tight (Fig. B.7).

There are two general conclusions that can be drawn. (i) For every investigated migration pattern, Scenario A is the one in which the highest LD occurs (in the demes next to the environmental step), and Scenario C is the one in which the maximum (absolute) LD is the lowest. This does not mean that in Scenario C, LD is everywhere lower than in Scenario A. (ii) For weak and intermediate migration and each of the selection scenarios A, B, or C, the average absolute amount of LD is highest with and lowest with .

In order to explain the patterns of LD in Scenario B and C, we recall from Remark 4.3 that in the two-deme model migration induces negative LD if (stabilizing selection) and positive LD if (directional selection). Proposition A.3 and Remark A.4 partially generalize this result: In Scenario B, weak migration induces non-positive LD in the demes under stabilizing selection and non-negative LD in the demes under directional selection. If migration connects environments under stabilizing selection with environments under directional selection (as in Scenarios B and C), negative and positive LD may offset each other. It is apparent from Fig. 6 that LD in Scenario C is much lower than in Scenario B which, in turn, is lower than in Scenario A. Its magnitude depends on the migration pattern, , and the deme (Figs. 6f,i and Figs. B.3, B.4, B.8). The degree of isolation by distance can have an ambiguous effect on the sign of LD. In Fig. 6f (), LD is negative for and , but positive for .

The parameter range

In this usually very small range of migration rates (Fig. 4), at least four equilibria are simultaneously stable in Scenarios , and . These equilibria may exhibit different means and variances (Fig. 8). The maximum variance (among stable equilibria) in the center of the habitat is of the same magnitude as for intermediate migration rates; compare Fig. 6e () with Fig. 8b (). The maximum LD (among stable equilibria) in the center of the habitat may be much higher than for intermediate migration rates; compare Fig. 6f () with Fig. 8c ().

Comparison with other multilocus models

Here, we compare our results of Section  6 to previously investigated clinal multilocus models.

A neutral model

Feldman and Christiansen (1975) studied a model without selection in which two ‘continents’ are fixed for different genetic backgrounds and are connected by demes with (single-step) stepping-stone migration into which they feed their genotypes. For two neutral loci, with fixed in deme 1 and fixed in deme , there is a unique cline which is linear in the allele frequencies and globally asymptotically stable. Linkage disequilibrium is unimodal with a maximum value of in the center of the cline. For a generalization to multiple loci, see  Christiansen (1986).

Comparing Fig. 6c with the approximation shows that LD in the neutral cline is usually much lower than in Scenario A with strong selection. However, it tends to be higher than in Scenario C (Fig. 6i), in which the cline in the mean is nearly linear. The variance in the center of the neutral cline is of the same order of magnitude as in the cline under migration–selection balance. At the boundaries, though, it is (fixed at) zero (compare Fig. B.9 with Figs. B.2–B.4).

The approximation (7.1) uses that in the model of  Feldman and Christiansen (1975) allele frequencies in adjacent demes differ by .  Kruuk et al. (1999, eq. (A.4)) generalized (7.1) to which assumes Fig. B.10 shows the accuracy of approximation (7.2) for each of the selection scenarios A, B, or C. If (7.3) is approximately satisfied, as in Scenario C, (7.2) approximates LD well for small migration rates. In the center of the cline with Scenario A (7.3) is obviously violated and (7.2) performs poorly. For the performance of an extension of (7.2) to weak selection by  Barton and Shpak (2000) see Fig. B.10.

Continuous space

Several models have been set up to describe clines at multiple loci or in polygenic traits in a continuous domain. Slatkin (1975) studied the effects of linkage on the clines in allele frequency at two loci and the associated LD. He assumed a step environment on the real line, analogous to our Scenario A, and used partial differential equations of the form where is the diffusion rate, and are the frequency and marginal fitness, respectively, of gamete at position at time , is the mean fitness, and denotes LD (Appendix A.11). These equations can be deduced as an approximation to the discrete-time model (2.3) with stepping-stone migration among a large number of demes in the same way as in  Nagylaki, 1975, Nagylaki, 1989, and need to be complemented by appropriate boundary conditions (usually, zero-flux conditions). Then is the (scaled) variance in dispersal distance. By assuming absence of epistasis, dominance, and LD, the cline in allele frequencies, i.e., the non-trivial equilibrium solution, can be calculated explicitly (Slatkin, 1975). However, he also showed numerically that, even in the absence of epistasis, LD is positive, maximized in the center of the cline, and decaying to zero with increasing distance from the center. If the loci are tightly linked, may approach its maximum value of at the center. In addition, a decreasing recombination rate steepens the cline in allele frequencies because stronger LD strengthens selection (Barton and Shpak, 2000).

Felsenstein (1977),  Slatkin (1978), and  Barton, 1983, Barton, 1999 investigated models of stabilizing selection on a quantitative trait by approximating gene flow by diffusion. Since the models of  Felsenstein (1977) and  Slatkin (1978) occur as limiting or special cases of Barton’s models (and are discussed there), we focus on two of Barton’s models but ignore mutation. Under the assumptions of a Gaussian distribution of allelic effects at each of loci and of linkage equilibrium,  Barton (1999, eqs. (4) and (5)) obtained Eqs. (A.38) for the evolution of the mean and the variance of the trait. As discussed there, the assumption of a Gaussian allelic distribution is rather restrictive; however, if it holds and the optimum changes gradually, the assumption of linkage equilibrium is supported by Felsenstein’s (1977) analysis.  Barton (1999, eq. (10)) also investigated a model in which diallelic loci of equal effect and in linkage equilibrium contribute to the trait. For two loci, his ‘rare-allele model’ is specified in (A.39). It is equivalent to (7.4) if is assumed in (7.4).

In Figs. 9, B.11, our results for with 12 demes are compared with the diffusion approximation (7.4) and with Barton’s models (A.38) and (A.39). To compare the diffusion approximations with our discrete model, we assumed that the habitat is the interval . Therefore, the boundary conditions are imposed. The diffusion rate , calculated as the variance in dispersal distance, depends on the position of the demes; see (A.41), (A.42), (A.43). Because holds to a close approximation (A.44), we use this as the uniform value. Equations (A.38), (A.39), and (7.4) were solved by using the Mathematica routine NDSolve and assuming spatially uniform initial conditions.

Fig. 9

Clines in the mean phenotype (left), the genetic variance (middle), and LD (right) for the Gaussian PDE model (A.38) (green lines), Barton’s PDE model (A.39) which assumes linkage equilibrium (blue dashed lines), and our model (7.4). The model (7.4) is shown for (yellow lines) and (red lines). Because (A.38) and (A.39) assume linkage equilibrium, LD is shown only for (7.4). Yellow and red dots show clines from the 12-deme stepping-stone model for and , respectively. In panels d, e, g, and h, blue dashed, yellow, and red lines overlap. Other parameters are and . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Figs. 9, B.11 show that the 12-deme and the PDE models yield similar clines in the mean phenotype. In Scenario A, the clines for linked loci are slightly steeper in the center of the cline than for unlinked ones. This was already predicted by  Slatkin (1975) for a model with linear directional selection. In Scenario C, the Gaussian model (A.38) leads to less adaptation than all other models near the boundaries of the habitat. The Gaussian model also exhibits large deviations from the variance maintained in all other models. This may not be too surprising because a Gaussian distribution of allelic effects is not suitable to approximate the distribution in a diallelic model. This could be different for models with many or a continuum of alleles.

In general, mean, variance, and LD in the 12-deme stepping-stone model are very accurately approximated by the PDE model (7.4), although with stronger dispersal the approximation for LD may become slightly less accurate (Fig. B.11 f, j). If LD is low, whether recombination is strong or weak, Barton’s linkage equilibrium approximation (A.39) for the mean and the variance is essentially indistinguishable from that based on (7.4). If LD is large, which is the case only in Scenario A near the environmental step, it affects the genetic variance to a notable extent. Then it leads to an elevated variance near the environmental step which, in turn, entails a slightly steeper gradient of the cline in the mean.

The clines in the mean, variance, and LD displayed in Figs. 9, B.11 are unique although in Scenarios B and C the underlying genotype-frequency equilibria are not. In fact, there are pairs of stable equilibria that have the same mean, variance, and LD (cf.  Fig. 6 and Proposition 5.2). If migration is much weaker than selection, the clines in the variance and in LD are no longer uniquely determined as already shown for the stepping-stone model (Fig. 8).

Loci with unequal effects

If unequal locus effects () and arbitrary linkage are admitted, there is a much larger number of possible equilibrium configurations. Already for , the analysis is much more intricate than that in Section  4 (cf.  Geroldinger and Bürger, 2014). There are several reasons for these complications.

(i) Single-locus polymorphisms may be stable for and their stability depends in a complicated way on and .

(ii) Four different equilibrium configurations can occur in the limit of strong migration, one with a globally asymptotically stable internal equilibrium, one with two locally stable internal equilibria, one with two locally stable SLPs, and the one with and locally stable (see  Bürger, 2000, p. 207). The latter applies if and only if and , where . In this case a numerical linear stability analysis shows that decreases in .

(iii) The definition (A.23) of the equilibria , , , , and can be extended to . However, the movement of the equilibria in the state space with increasing migration rate is much more complicated. Nevertheless, for a large set in the parameter space either two equilibria ( and ) are simultaneously stable or a unique stable equilibrium () exists. In contrast to , where pairs of equilibria (such as , or , ) have the same mean, variance, and LD, this not so if . A frequently occurring analogue of the patterns in Fig. 6 is shown in Fig. B.12.

Discussion

Here we recapitulate our model and results in a non-technical way and discuss the relation to previous literature. The purpose of this work was to investigate the effects of migration patterns and selection scenarios on the maintenance and the properties of clines in a quantitative trait. We assumed that the trait is determined additively by two diallelic, recombining loci. Fitness decays quadratically from a phenotypic optimum. Because the position of the optimum, which depends on space, may be anywhere within the range of possible phenotypes, the trait may be under stabilizing or directional selection.

The population is subdivided into discrete demes, representing different locations in space (Section  2), or inhabits continuous space, i.e., a bounded one-dimensional interval (Section  7.2). One of the advantages of using discrete demes is that different migration patterns can be modeled easily, whereas diffusion models are based on the assumption that there is mainly short-distance migration and evolutionary forces are weak. In particular, we used the island model (denoted by ), in which there is no distance because every island (deme) is reached with the same probability, the single-step stepping-stone model (), and a multi-step stepping-stone model () in which more distant demes can be reached with reduced probability (Section  5.1). Demes are ordered from 1 to , although this (spatial) order is irrelevant in the island model. If , all three migration patterns coincide. Therefore, the analysis of the two-deme model in Section  4 is of central importance.

In addition, we employed four different selection scenarios. They are illustrated in Fig. 3. In Scenario A, the phenotypic optimum is at the left boundary of the phenotypic range in demes ( even) and at the right boundary in the others. Thus, there is a single, abrupt environmental change in the middle of the spatial domain. In Scenario B, there is directional selection toward the extreme phenotypes (as in Scenario A) in the left and the right third of demes, whereas the trait is under stabilizing selection in the middle third. Thus, there are two environmental steps and hybrids are favored in the middle of the spatial range. In Scenario C, the phenotypic optimum changes linearly from the left to the right boundary of the phenotypic range, thus reflecting a gradual change of the environment. In Scenario D, there is spatially uniform stabilizing selection toward the middle of the phenotypic range.

In the hybrid zone literature, selection schemes which disfavor hybrids everywhere (Scenario A) are sometimes referred to as ‘ecotone zone’, whereas selection schemes which favor intermediate genotypes (Scenarios B and C) are also known as ‘hybrid superiority zones’ (Kawakami and Butlin, 2012). Although we assume exogenous selection, selection is not purely exogenous in the sense of  Kruuk et al. (1999) because our fitnesses display dominance and epistasis.

Our study complements or extends previous investigations that were based on diffusion models in several ways. For instance,  Felsenstein (1977) and  Slatkin (1978) assumed a multivariate normal distribution of allelic effects. Felsenstein assumed a linearly changing optimum, as in our Scenario C, whereas Slatkin considered Scenarios A and C.  Slatkin (1975) studied a diallelic two-locus diffusion model with a step environment (Scenario A) that can be interpreted as a model of a quantitative trait under linear directional selection. Except for different assumptions about selection, Slatkin’s (1975) model is identical to our diffusion model (7.4).  Barton (1983) and  Barton and Shpak (2000) assumed multiple loci and spatially independent (endogenous) selection against hybrids (without or with epistasis, respectively).  Kruuk et al. (1999) compared aspects of models of endogenous selection with a model similar to our Scenario C.  Barton (1999) investigated a diffusion-based model that is equivalent to our diffusion model (7.4), except that his ignores LD. He compared his so-called rare-alleles model to several of the above mentioned models (which are stated in Appendix A.11).

Because clines result from polymorphic equilibria, we first summarize the results about equilibrium configurations and bifurcation patterns. The important limiting cases of weak and of strong migration are analyzed in Section  3. They apply to every migration pattern and selection scenario.

If migration is sufficiently weak relative to selection, then in selection scenarios B, C, and D there are multiple, simultaneously stable polymorphic equilibria for every migration pattern. Their number increases (approximately) exponentially in and, for given , from B to C to D (Section  5.2). For Scenarios B and C, the critical migration rate below which more than two polymorphic equilibria are simultaneously stable is usually one or two orders of magnitude smaller than the selection parameter . This range is indicated by the orange bars in Fig. 4. Its upper bound is the critical migration rate (see also Table A.1). For Scenario A and weak (or moderate) migration, there is always a unique fully polymorphic equilibrium (), i.e.,  .

If migration is sufficiently strong relative to selection, then no polymorphism is maintained for any selection scenario or migration pattern because one of the haplotypes with intermediate phenotype ( or ) swamps the whole population. In this case, the monomorphic equilibria and are the only stable equilibria. The critical migration rate above which no stable polymorphism can be maintained is given by for Scenarios and by for Scenario ; see (5.17). The gray regions in Fig. 4 show . Notably, for every given migration pattern, decreases from Scenario A to B to C to D (5.18). Hence, in a step environment stable polymorphic equilibria can be maintained for much higher gene flow than in a gradually changing environment.

As the number of demes increases, increases slowly (in proportion to ) for the island model and Scenarios A, B, and D; Eqs.  (5.9), (5.12), (5.14). In Scenario C, this can be violated for reasons explained in Section  5.3. For the stepping-stone model, increases faster than linear in for Scenarios A, B, and C (Fig. 5). This much faster increase is not surprising because isolation by distance increases with the number of demes. Finally, for any given selection scenario, increases from to to , again supporting the intuition that increasing isolation by distance facilitates the maintenance of polymorphism.

The range of migration rates between the critical values and can be partitioned into up to three different intervals in which there is either a unique, globally asymptotically stable internal equilibrium, , or a pair of asymptotically stable internal equilibria, and (Section  5.3). For Scenario D, such an intermediate range does not exist because . The range of migration rates for which there is a unique stable equilibrium decreases from Scenario A to B to C to D, for which it vanishes (Fig. 4, Table A.1). Interestingly, every bifurcation pattern that was found for in any of the selection scenarios occurs in essentially the same form for a certain range of positions of the optima ( and ) in the two-deme model of Section  4. They are displayed in Fig. 1. The only qualitative difference is that for two demes and weak migration, at most two internal equilibria are stable instead of many (up to in Scenario D).

Stable polymorphic equilibria give rise to (stable) clines. Because our interest is in quantitative traits and how local adaptation and genetic variation depend on migration patterns and selection scenarios, we studied clines in the mean phenotype and in the (total) genetic variance. In addition, we investigated LD and its spatial dependence. Except for very weak migration (), when there are many simultaneously stable equilibria, the clines in the mean, variance, and LD are unique even if the underlying polymorphic equilibria differ. This is due to the symmetry assumptions of the model (cf. Fig. B.11).

In the language of hybrid zones, the results discussed above show that in a hybrid superiority zone clines exist only for lower migration rates than in an ecotone zone and, for very low migration rates, initial conditions play a more important role because of the existence of multiple clines. The reason for the former finding is that in a hybrid superiority zone, the haplotypes and swamp the entire population easier than in an ecotone zone.

The shape of the clines is strongly influenced by both the migration pattern and by the selection scenario. This is exemplified by Fig. 6 and documented extensively by Figs. B.2–B.4. In most cases, the degree of local adaptation (as measured by the deviation of the mean from the optimum) increases from migration pattern to to , i.e., with increasing isolation by distance. This seems to be universally true for Scenario A, but is violated for Scenario B in the demes with stabilizing selection (e.g.,  Fig. 6d). In these demes, the island model provides maximum adaptation, whereas in the demes with directional selection the stepping-stone model () maximizes local adaptation. There are also rare exceptions in Scenario C.

For the island model, the genetic variance is spatially uniform in Scenario A and weakly dependent on space in Scenarios B and C. In the latter, it may be maximized or minimized in the center, or it may be bimodal (Figs. B.3, B.4). For the stepping-stone models and a step environment (Scenario A), the genetic variance is always maximized in the center of the cline and decreases toward its boundaries. For Scenarios B and C, the variance may be maximized in the center or elsewhere (Figs. 6e,h). Distinctive bimodal patterns occur mainly for weak single-step migration (Figs. B.3, B.4). The modes occur in demes at the boundary between regions of stabilizing and of directional selection.

An increase in has a simple effect on local adaptation: it is progressively reduced until the cline collapses at . Its effects on the genetic variance are more complex, as is documented by Figs. B.2–B.4. However, the variance of the total population is rather insensitive to changes of over a wide range (Fig. 7). It may be slowly decreasing in or be maximized at intermediate values. In Scenario A, tight linkage may substantially increase , whereas it is almost independent of in Scenario B and C.

Next, we discuss LD and the role of recombination. The examples presented in Fig. 6 are representative for a large range of parameters. A much more complete picture is obtained from Figs. B.2–B.4 for and Figs. B.5–B.7 for . Although the details are complex, some general conclusions emerge.

(i) The highest linkage disequilibria occur in Scenario A in the demes adjacent to the environmental step. LD is always positive, as is expected under a balance between directional selection and migration (Li and Nei, 1974,  Christiansen and Feldman, 1975,  Slatkin, 1975,  Bürger and Akerman, 2011,  Akerman and Bürger, 2014b), although this is not universally true in the presence of epistasis (Geroldinger and Bürger, 2014). Of course, LD increases with tighter linkage.

(ii) In Scenario C, LD is very weak under stepping-stone migration. This is in line with Felsenstein’s (1977) result that for normally distributed allelic effects and a linearly changing optimum, LD is absent at equilibrium. For the island model, small positive LD is maintained in the demes under directional selection and small negative LD otherwise.

(iii) The most complex patterns occur for Scenario B because in the central demes there is stabilizing selection which induces negative LD. In general, the absolute magnitude of LD is between those of Scenarios A and C, and stronger recombination obviously reduces LD. The typical spatial patterns are displayed in Fig. 6f. Interestingly, for the stepping-stone models, LD is nearly absent for weak or moderately strong migration, but becomes appreciable for strong migration. For the island model, essentially the opposite is true; LD is relatively high for low migration and vanishes for large (see Fig. B.8 for details). The reason for this finding is the different degree of mixing exhibited by the migration patterns.

Because LD is low in Scenario B and almost absent in Scenario C, the clines in the mean and the variance are hardly affected by recombination or LD. In Scenario A, recombination and LD affect the clines as follows. Because lower induces higher LD, the variance is somewhat inflated if the increase in LD is sufficiently high. For the stepping-stone model, this occurs near the center of the cline, and for the island model it is a spatially universal effect (Fig. B.5). As in Slatkin’s (1975) model, reduced leads to a slightly steeper cline and to a slight increase in local adaptation. The reason is that stronger positive LD strengthens selection at each locus (cf.  Barton, 1983).

The approximations for LD based on the assumptions of neutrality (7.3) or quasi-linkage equilibrium (Barton and Shpak, 2000, Eq. (14)) perform well in Scenarios B and C if migration is sufficiently weak, so that LD is very small (Fig. B.10). In Scenario A, but not otherwise, a variant of (7.3) performs very well over a wide range of migration rates. If migration is not weak, the neutral approximation (7.3) tends to overestimate LD, whereas the quasi-linkage-equilibrium approximation tends to underestimate it.

The majority of our numerical results is based on the assumption of strong selection. The choice in many of the figures implies that in the demes under directional selection, the fitness of the least fit phenotype is only 20% of that of the optimum phenotype. In Scenario A, this applies to every deme. Nevertheless, comparison of the 12-deme model to the diffusion model (7.4), whose derivation is based on the assumption of weak evolutionary forces (Nagylaki, 1975, Nagylaki, 1989), shows excellent concordance (Figs. 9, B.11). Therefore, most of our discussion above carries over to the corresponding diffusion models. These figures (as well as the discussion above) also show that Barton’s (1999) ‘rare-alleles’ diffusion model (Eq. (A.39) in Appendix A), which ignores LD, provides accurate approximations to the clines in the mean and the variance unless loci are tightly linked. The Gaussian model (Eq. (A.38)) yields almost accurate clines in the mean, but distinctively deviant ones in the variance.

Finally, most of our analysis is based on symmetry assumptions. Throughout, we assumed a one-to-one correspondence of demes in which the phenotypic optimum is or , and we assumed symmetric migration (3.7). Most of the analysis is also based on the assumption of loci of equal effects. Deviation from any of these assumptions will have multiple consequences. First, most pitchfork bifurcations will be replaced by (pairs of) saddle–node bifurcations. Second, different polymorphic equilibria will give rise to different clines, hence clines in the mean, variance, and LD will no longer be unique, unless there is a unique polymorphic equilibrium (corresponding to ). Third, stability of single-locus polymorphisms will be facilitated. Fourth, even in the limit of strong migration, a globally stable fully polymorphic equilibrium (hence a cline) can be maintained if locus effects are sufficiently different and linkage is tight (Section  7.3). Therefore, can be infinite. As demonstrated by Geroldinger and Bürger (2014) for a haploid model, even if is finite, a reduction of the ratio of locus effects can lead to an increase or a decrease of , depending on whether recombination is low or high. Fifth, deviation from the symmetry assumptions about selection or migration will, in general, lead to a reduction of by facilitating fixation of the haplotype with the highest mean fitness, i.e., averaged across demes and weighted by the principal eigenvector of the migration matrix.