The consequences of gene flow for local adaptation and differentiation: a two-locus two-deme model

We consider a population subdivided into two demes connected by migration in which selection acts in opposite direction. We explore the effects of recombination and migration on the maintenance of multilocus polymorphism, on local adaptation, and on differentiation by employing a deterministic model with genic selection on two linked diallelic loci (i.e., no dominance or epistasis). For the following cases, we characterize explicitly the possible equilibrium configurations: weak, strong, highly asymmetric, and super-symmetric migration, no or weak recombination, and independent or strongly recombining loci. For independent loci (linkage equilibrium) and for completely linked loci, we derive the possible bifurcation patterns as functions of the total migration rate, assuming all other parameters are fixed but arbitrary. For these and other cases, we determine analytically the maximum migration rate below which a stable fully polymorphic equilibrium exists. In this case, differentiation and local adaptation are maintained. Their degree is quantified by a new multilocus version of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_\mathrm{ST}$$\end{document} and by the migration load, respectively. In addition, we investigate the invasion conditions of locally beneficial mutants and show that linkage to a locus that is already in migration-selection balance facilitates invasion. Hence, loci of much smaller effect can invade than predicted by one-locus theory if linkage is sufficiently tight. We study how this minimum amount of linkage admitting invasion depends on the migration pattern. This suggests the emergence of clusters of locally beneficial mutations, which may form ‘genomic islands of divergence’. Finally, the influence of linkage and two-way migration on the effective migration rate at a linked neutral locus is explored. Numerical work complements our analytical results.


Introduction
Migration in a geographically structured population may have opposing effects on the genetic composition of that population and, hence, on its evolutionary potential. On the one hand, gene flow caused by migration may be so strong that it not only limits but hinders local adaptation by swamping the whole population with a genotype that has high fitness in only one or a few demes. On the other hand, if migration is sufficiently weak, gene flow may replenish local populations with genetic variation and contribute to future adaptation. In this case, locally adapted genotypes may coexist in the population and maintain high levels of genetic variation as well as differentiation between subpopulations. For reviews of the corresponding, well developed one-locus theory, see Karlin (1982), Lenormand (2002), and Nagylaki and Lou (2008).
If selection acts on more than one locus, additional questions arise immediately. For instance, what are the consequences of the genetic architecture, such as linkage between loci, relative magnitude of locus effects, or epistasis, on the degree of local adaptation and of differentiation achieved for a given amount of gene flow? What are the consequences for genetic variation at linked neutral sites? What genetic architectures can be expected to evolve under various forms of spatially heterogeneous selection?
For selection acting on multiple loci, the available theory is much less well developed than for a single locus. One of the main reasons is that the interaction of migration and selection, even if the latter is nonepistatic, leads to linkage disequilibrium (LD) between loci (Li and Nei 1974;Christiansen and Feldman 1975;Slatkin 1975;Barton 1983). LD causes substantial, often insurmountable, complications in the analysis of multilocus models. Therefore, many multilocus studies are primarily numerical and focus on quite specific situations or problems. For instance, Spichtig and Kawecki (2004) investigated numerically the influence of the number of loci and of epistasis on the degree of polymorphism if selection acts antagonistically in two demes. Yeaman and Whitlock (2011) showed that concentrated genetic architecture, i.e., clusters of linked, locally beneficial alleles, evolve if stabilizing selection acts on a trait such that the fitness optima in two demes differ.
Linkage disequilibrium is also essential for the evolution of recombination. The evolution of recombination in heterogeneous environments has been studied by a number of authors (e.g., Charlesworth and Charlesworth 1979;Pylkov et al. 1998;Lenormand and Otto 2000), and the results depend strongly on the kind of variability of selection across environments, the magnitude of migration, and the sign and strength of epistasis.
Recent years have seen some advances in developing general theory for multilocus migration-selection models. The focus of this work was on the properties of the evolutionary dynamics and the conditions for the maintenance of multilocus polymorphism in limiting or special cases, such as weak or strong migration (Bürger 2009a,b), or in the Levene model (Nagylaki 2009;Bürger 2009cBürger , 2010Barton 2010;Chasnov 2012). This progress was facilitated by the fact that in each case, LD is weak or absent.
Using a continent-island-model framework, Bürger and Akerman (2011) and Bank et al. (2012) analyzed the effects of gene flow on local adaptation, differentiation, the emergence of Dobzhansky-Muller incompatibilities, and the maintenance of polymorphism at two linked diallelic loci. They obtained analytical characterizations of the possible equilibrium configurations and bifurcation patterns for wide ranges of parameter combinations. In these models, typically high LD is maintained. In particular, explicit formulas were derived for the maximum migration rate below which a fully polymorphic equilibrium can be maintained, as well as for the minimum migration rate above which the island is swamped by the continental haplotype.
Here, we explore the robustness of some of these results by admitting arbitrary (forward and backward) migration between two demes. This generalization leads to substantial mathematical complications, but also to new biological insight. Because our focus is on the consequences of gene flow for local adaptation and differentiation, we assume divergent selection among the demes, i.e., alleles A 1 and B 1 are favored in deme 1, and A 2 and B 2 are favored in deme 2. The loci may recombine at an arbitrary rate. By ignoring epistasis and dominance, we assume genic selection. Mutation and random drift are neglected. Because we assume evolution in continuous time, our model also describes selection on haploids.
The model is set up in Sect. 2. In Sect. 3, we derive the equilibrium and stability structure for several important special cases. These include weak, strong, highly asymmetric, and super-symmetric migration, no or weak recombination, independent or strongly recombining loci, and absence of genotype-environment interaction. In Sect. 4, we study the dependence of the equilibrium and stability patterns on the total migration rate while keeping the ratio of migration rates, the recombination rate, and the selection coefficients constant (but arbitrary). In particular, we derive the possible bifurcation patterns for the cases of independent loci (linkage equilibrium) and for completely linked loci. With the help of perturbation theory, we obtain the equilibrium and stability configurations for weak or strong migration, highly asymmetric migration, and weak or strong recombination. For these cases, we determine the maximum migration rate below which a stable, fully polymorphic equilibrium is maintained, and the minimum migration rate above which the population is monomorphic. Numerical work complements our analytical results.
The next four sections are devoted to applications of the theory developed in Sects. 3 and 4. In Sects. 5 and 6, we use the migration load and a new, genuine multilocus, fixation index (F ST ), respectively, to quantify the dependence of local adaptation and of differentiation on various parameters, especially, on the migration and the recombination rate. In Sect. 7, we investigate the invasion conditions for a mutant of small effect ( A 1 ) that is beneficial in one deme but disadvantageous in the other deme. We assume that the mutant is linked to a polymorphic locus which is in selection-migration balance. We show that linkage between the loci facilitates invasion. Therefore, in such a scenario, clusters of locally adapted alleles are expected to emerge (cf. Yeaman and Whitlock 2011;Bürger and Akerman 2011). In Sect. 8, we study the strength of barriers to gene flow at neutral sites linked to the selected loci by deriving an explicit approximation for the effective migration rate at a linked neutral site. Our results are summarized and discussed in Sect. 9. Several purely technical proofs are relegated to the Appendix.

The model
We consider a sexually reproducing population of monoecious, diploid individuals that is subdivided into two demes connected by genotype-independent migration. Within each deme, there is random mating. We assume that two diallelic loci are under genic selection, i.e., there is no dominance or epistasis, and different alleles are favored in different demes. We assume soft selection, i.e., population regulation occurs within each deme. We ignore random genetic drift and mutation and employ a deterministic continuous-time model to describe evolution. A continuous-time model is obtained from the corresponding discrete-time model in the limit of weak evolutionary forces (here, selection, recombination, and migration).
We denote the rate at which individuals in deme 1 (deme 2) are replaced by immigrants from the other deme by m 1 ≥ 0 (m 2 ≥ 0). Then m = m 1 + m 2 is the total migration rate. The recombination rate between the two loci is designated by ρ ≥ 0.
Alleles at locus A are denoted by A 1 and A 2 , at locus B by B 1 and B 2 . We posit that A 1 and B 1 are favored in deme 1, whereas A 2 and B 2 are favored in deme 2. In deme k (k = 1, 2), we assign the Malthusian parameters 1 2 α k and − 1 2 α k to A 1 and A 2 , and 1 2 β k and − 1 2 β k to B 1 and B 2 . Because we assume absence of dominance and of epistasis, the resulting fitness matrix for the genotypes reads (2.1) By relabeling alleles, we can assume without loss of generality α 1 > 0 > α 2 and β 1 > 0 > β 2 . Hence, A 1 B 1 and A 2 B 2 may be called the locally adapted haplotypes in deme 1 and deme 2, respectively. By relabeling loci, we can assume β 1 ≥ α 1 . We define θ = α 1 β 2 − α 2 β 1 .
(2.2) By exchanging demes, i.e., by the transformationα k = −α k * andβ k = −β k * (where k * denotes the deme = k), or by exchanging loci, i.e., by the transformationα k = β k andβ k = α k , we can further assume θ ≥ 0 without loss of generality, cf. Appendix A.1. The fitness matrix (2.1) is also obtained if the two loci contribute additively to a quantitative trait that is under linear directional selection in each deme (Bürger 2009c). Then θ = 0 if the genotypic values are deme independent, i.e., if there is no genotype-environment interaction on the trait level.
(2.4) Therefore, locus A is under weaker selection than locus B in both demes, i.e., |α k | ≤ |β k | for k = 1, 2, if and only if β 2 < α 2 holds. The population can be described by the gamete frequencies in each of the demes. We denote the frequencies of the four possible gametes A 1 B 1 , A 1 B 2 , A 2 B 1 , and A 2 B 2 in deme k by x 1,k , x 2,k , x 3,k , and x 4,k . Then the state space is S 4 × S 4 , where S 4 = {(x 1 , x 2 , x 3 , x 4 ) : x i ≥ 0 and 4 i=1 x i = 1} is the simplex. The following differential equations for the evolution of gamete frequencies in deme k can be derived straightforwardly: (2.5) Here the marginal fitness w i,k of gamete i and the mean fitnessw k in deme k are calculated from (2.1), η 1 = η 4 = −η 2 = −η 3 = 1, and D k = x 1,k x 4,k − x 2,k x 3,k is the linkage-disequilibrium (LD) measure. We note that D k > 0 corresponds to an excess of the locally adapted haplotypes in deme k. The equations (2.5) also describe the dynamics of a haploid population if in deme k we assign the fitnesses α k , −α k , β k , −β k to the alleles A 1 , A 2 , B 1 , B 2 , respectively. Instead of gamete frequencies it is often more convenient to work with allele frequencies and the LD measures D k . We write p k = x 1,k + x 2,k and q k = x 1,k + x 3,k for the frequencies of A 1 and B 1 in deme k. Then the gamete frequencies x i,k are calculated from the p k , q k , and D k by (2.6b) The constraints x i,k ≥ 0 and 4 i=1 x i,k = 1 for i = 1, 2, 3, 4, and k = 1, 2 It follows that p k , q k , and D k evolve according to (2.7c) We emphasize that, because we are treating a continuous-time model, the parameters ρ, m k , α k , and β k are rates (of recombination, migration, growth), whence they can be arbitrarily large. Their magnitude is determined by the time scale. By rescaling time, for instance to units of ρ or m, the number of independent parameters could be reduced by one without changing the equilibrium properties.

Equilibria and their stability
We distinguish three types of equilibria: (i) monomorphic equilibria (ME), (ii) singlelocus polymorphisms (SLPs), and (iii) full (two-locus) polymorphisms (FPs). The first two types are boundary equilibria, whereas FPs are internal equilibria (except when ρ = 0). The stability properties of the ME and the coordinates and conditions for admissibility of the SLPs can be derived explicitly. However, the stability conditions for the SLPs and the conditions for existence or stability of FPs could be derived only for a number of limiting cases. These include strong recombination, weak or no recombination, weak, strong, or highly asymmetric migration.

Existence of boundary equilibria
The four ME, corresponding to fixation of one of the gametes, exist always. Their coordinates are as follows: where aˆsignifies an equilibrium. There are up to four SLPs, one in each marginal one-locus system. We denote the SLPs where B 1 or B 2 is fixed by P A,1 or P A,2 , respectively, and the SLPs where A 1 or A 2 is fixed by P B,1 or P B,2 . Their coordinates and the conditions for their admissibility can be calculated explicitly (Eyland 1971). We define (3.1) By (2.3a), we have In addition, it is easy to show that the assumptions (2.3) imply: If locus B is fixed (for B 1 or B 2 ), the equilibrium allele frequencies at locus A arê If locus A is fixed, the equilibrium allele frequencies at locus B are given bŷ Thus, the four SLPs have the following coordinates: The equilibria P A,1 and P A,2 are admissible if and only if |σ 1 + σ 2 | < 1, (3.7) and the equilibria P B,1 and P B,2 are admissible if and only if The SLPs leave the state space through one of their 'neighboring' ME if |σ 1 + σ 2 | or |τ 1 + τ 2 | increases above 1. In particular, we find (3.9c) τ 1 + τ 2 ↑ 1 ⇐⇒ P B,1 → M 2 and P B,2 → M 4 . (3.9d) Throughout, we use ↓ to indicated convergence from above and ↑ to indicate convergence from below. Figure 1 illustrates the location of the possible equilibria. The SLPs are asymptotically stable within their marginal one-locus system if and only if they are admissible. Then they are also globally asymptotically stable within their marginal system (Eyland 1971). (We use globally stable in the sense that at least all trajectories from the interior of the designated set converge to the equilibrium.) The Fig. 1 Location of equilibria. In terms of gamete frequencies, the state space is S 4 × S 4 , where each S 4 corresponds to one deme. This figure shows (schematically) the location in S 4 of all boundary equilibria and of the stable internal equilibrium F. F converges to F ∞ if ρ → ∞ and to F 0 if ρ → 0. The LE manifold is indicated by hatching reader may notice that (3.7) and (3.8) are precisely the conditions for maintaining a protected polymorphism at locus A and B, respectively.

Stability of monomorphic equilibria
At each monomorphic equilibrium, the characteristic polynomial factors into three quadratic polynomials. Two of them determine stability with respect to the marginal one-locus systems, whereas the third determines stability with respect to the interior of the state space. The stability properties of the monomorphic equilibria are as follows. The proof is given in Appendix A.2. Proposition 3.1 M 1 is asymptotically stable if σ 1 + σ 2 < −1 and τ 1 + τ 2 < −1 (3.10) and one of the following conditions hold: M 2 is always unstable. M 3 is asymptotically stable if σ 1 + σ 2 > 1 and τ 1 + τ 2 < −1. (3.12) M 4 is asymptotically stable if σ 1 + σ 2 > 1 and τ 1 + τ 2 > 1 (3.13) and one of the following conditions hold: If one of the inequalities in (3.10), (3.12), or (3.13), or one of the inequalities for m 2 in (3.11b) or (3.14b) is reversed, the respective equilibrium is unstable.
The above result shows that each of M 1 , M 3 , or M 4 can be stable, but never simultaneously. For sufficiently loose linkage, the stability of a ME is determined solely by its stability within the two marginal one-locus systems in which it occurs. Stability of M 3 is independent of the recombination rate. For given migration rates, the equilibria M 1 and M 4 may be stable for high recombination rates but unstable for low ones. For a low total migration rate (m 1 + m 2 ), no ME is stable. For a sufficiently high total migration rate, there is a globally asymptotically stable ME (Sect. 4.5).

Stability of single-locus polymorphisms
As already mentioned, a single-locus polymorphism is globally attracting within its marginal one-locus system whenever it is admissible. Although the coordinates of the SLPs are given explicitly, the conditions for stability within the full, six-dimensional system on S 4 × S 4 are uninformative because the four eigenvalues that determine stability transversal to the marginal one-locus system are solutions of a complicated quartic equation.
In the following we treat several limiting cases in which the conditions for stability of the SLPs and for existence and stability of FPs can be obtained explicitly.

Weak migration
The equilibrium and stability structure for weak migration can be deduced from the model with no migration by perturbation theory. In the absence of migration (m 1 = m 2 = 0), the two subpopulations evolve independently. Because selection is nonepistatic and there is no dominance, in each deme the fittest haplotype becomes eventually fixed. In fact, mean fitness is nondecreasing (Ewens 1969). Our assumptions about fitness, i.e., (2.1) and (2.3a), imply that in deme 1 the equilibrium witĥ p 1 =q 1 = 1 andD 1 = 0 is globally attracting, and in deme 2 the equilibrium with p 2 =q 2 = 0 andD 2 = 0 is globally attracting. Therefore, in the combined system, i.e., on S 4 × S 4 , but still with m 1 = m 2 = 0, the (unique) globally attracting equilibrium is given bŷ All other equilibria are on the boundary and unstable. Because, generically, all equilibria in the system without migration are hyperbolic and it is a gradient system (Shahshahani 1979;Bürger 2000, p. 42), Theorem 5.4 in Bürger (2009a applies and shows that the perturbation F of the equilibrium (3.15) is globally asymptotically stable for sufficiently small migration rates m 1 and m 2 . Boundary equilibria remain unstable for sufficiently small migration rates. It is straightforward to calculate the coordinates of the perturbed equilibrium to leading order in m 1 and m 2 . They are given bŷ Therefore, we conclude Proposition 3.2 For sufficiently weak migration, there is a unique, globally attracting, fully polymorphic equilibrium F. To leading order in m 1 and m 2 , its coordinates are given by (3.16).
Proposition 3.2 remains valid if the assumptions (2.3b) and (2.3c) are dropped. Apart from the obvious fact that migration reduces differences between subpopulations, the above approximations show that the lower the recombination rate, the smaller is this reduction. Thus, for given (small) migration rates, differentiation between subpopulations is always enhanced by reduced recombination. Linkage disequilibria within subpopulations are always positive.
In (3.17), the differential equations for the two loci are decoupled, i.e., (3.17a) and (3.17b) as well as (3.17c) and (3.17d) form closed systems. Thus, the dynamics of the full system is a Cartesian product of the two one-locus dynamics. Therefore, in addition to the ME and to the SLPs determined above, the following internal equilibrium, denoted by F ∞ , may exist Because in the one-locus model the FP is globally asymptotically stable (hence, it attracts all trajectories from the interior of the state space) whenever it is admissible (Eyland 1971;Hadeler and Glas 1983, Theorem 2;Nagylaki and Lou 2008, Section 4.3.2), and because the full dynamics is the Cartesian product of the one-locus dynamics, the fully polymorphic equilibrium F ∞ is globally asymptotically stable whenever it is admissible. Similarly, we conclude that a boundary equilibrium is globally asymptotically stable whenever it is asymptotically stable in the full system. These results in combination with those in Sects. 3.1 and 3.2 yield the following proposition. Proposition 3.3 Assume (3.17). Then a globally asymptotically stable equilibrium exists always. This equilibrium is internal, hence equals F ∞ (3.18), if and only if (3.7) and (3.8) hold. It is a SLP if one of (3.7) or (3.8) is violated, and a ME if both (3.7) and (3.8) are violated.
We note that F ∞ → P A,2 does not occur because it requires τ 1 + τ 2 ↑ 1 and |σ 1 + σ 2 | < 1, which is impossible by (3.3). We leave the simple determination of the conditions for bifurcations of F ∞ with one of the ME to the interested reader. Proposition 3.3 can be extended straightforwardly to an arbitrary number of loci because the dynamics at each locus is independent of that at the other loci. This decoupling of loci occurs because there is no epistasis.

Strong recombination: quasi-linkage equilibrium
If recombination is strong, a regular perturbation analysis of the internal equilibrium F ∞ of (3.17) can be performed. The allele frequencies and linkage disequilibria can be calculated to order 1/ρ. Formally, we set (3.20) keep σ k and τ k constant, and let ρ → ∞. Then, we obtain and analogous formulas hold for the second deme. Because LD is of order 1/ρ, this approximation may be called the quasi-linkage equilibrium approximation of the fully polymorphic equilibrium (Kimura 1965;Turelli and Barton 1990;Nagylaki et al. 1999). We note that LD is positive in both demes and increases with increasing differentiation between the demes, increasing migration, or decreasing recombination. Proposition 5.1 in Bürger (2009a) shows that in every small neighborhood of an equilibrium of the model with LE (3.17), there is one equilibrium of the perturbed system, and it has the same stability properties as the unperturbed equilibrium. Because of the simple structure of (3.17), a stronger result can be obtained. In an isolated onelocus system on [0, 1] 2 (e.g., (3.17a) and (3.17b)), every trajectory from the interior converges to the unique asymptotically stable equilibrium (Sect. 3.5), and the chainrecurrent points (Conley 1978) are the equilibria. Therefore, the same holds for the LE dynamics (3.17), and the regular global perturbation result of Nagylaki et al. (1999) (the proof of their Theorem 2.3) applies for large ρ. Hence the dynamical behavior with strong recombination is qualitatively the same as that under LE. We conclude that for sufficiently strong recombination every asymptotically stable equilibrium is globally asymptotically stable.

No recombination
Let recombination be absent, i.e., ρ = 0. Then, effectively, we have a one-locus model in which the four alleles correspond to the four gametes In deme k, they have the selection coefficients 1 respectively. According to Theorem 2.4 of Nagylaki and Lou (2001), generically, no more than two gametes can be present at an equilibrium. We will prove a stronger result and characterize all possible equilibria and their local stability.
Because ρ = 0, there may be a polymorphic equilibrium at which only the gametes A 1 B 1 and A 2 B 2 are present. We call it F 0 and set Then one-locus theory (Sect. 3.1) informs us that F 0 is admissible if and only if |κ 1 + κ 2 | < 1. (3.23) Its coordinates are given bŷ . Within the subsystem in which only the gametes A 1 B 1 and A 2 B 2 are present, F 0 is asymptotically stable whenever it is admissible. One-locus theory implies that (3.25b) A simple application of Corollary 3.9 of Nagylaki and Lou (2007) shows that the gamete A 1 B 2 will always be lost (to apply their result, recall assumptions (2.3) and use γ 22 = γ 23 = 0, α 1 α 1 +β 1 < γ 21 < α 2 α 2 +β 2 , γ 24 = 1 − γ 21 ). This strengthens the result in Sect. 3.2 that M 2 is always unstable. Thus, we are left with the analysis of the tri-gametic system consisting of A 1 B 1 , A 2 B 1 , and A 2 B 2 . (If θ < 0, then gamete A 2 B 1 is lost.) In Appendix A.3 it is proved that F 0 is the only equilibrium at which both loci are polymorphic except when If (3.26) holds, then there is a line of internal equilibria connecting F 0 with P A,1 or P B,2 (or M 3 ); see Appendix A.3. We find that F 0 is asymptotically stable if m 1 m 2 <m, (3.28) and unstable if the inequality is reversed (Appendix A.4). For sufficiently small migration rates, Proposition 3.2 implies that F 0 = F and F 0 is globally asymptotically stable. If the inequality in (3.28) is reversed, F 0 may or may not be admissible. Of course, if F 0 is asymptotically stable, then the equilibria M 1 and M 4 are unstable; cf. (3.25). The following argument shows that M 3 cannot be simultaneously stable with F 0 . We rewrite (3.28) as (3.29) Because Since M 3 is asymptotically stable if (3.12) holds and because, as is easy to show, (3.12) and (3.31) are incompatible, the assertion follows. It can also be shown from (3.12) and (3.31) that M 3 cannot become stable when F 0 loses its stability except in the degenerate case when σ 1 + σ 2 = 1 and τ 1 + τ 2 = −1.
In our tri-gametic system, P A,1 and P B,2 are the only possible SLPs. They may exist simultaneously with F 0 if (3.28) holds, i.e., if F 0 is stable, but not otherwise (Appendix A.5). If (3.28) holds, both are unstable (if admissible). P A,1 or P B,2 have an eigenvalue 0 if and only if (3.26) holds or if they leave or enter the state space through a ME. In Appendix A.5 it is shown that P A,1 is asymptotically stable if and only if τ 1 + τ 2 < −1 and |σ 1 + σ 2 | < 1 and m 1 m 2 >m, (3.32) and P B,2 is asymptotically stable if and only if 1 < σ 1 + σ 2 and |τ 1 + τ 2 | < 1 and m 1 m 2 >m. (3.33) Hence, if m 1 m 2 increases abovem, the SLP that is admissible becomes asymptotically stable. Upon collision of the stable SLP with one of the adjacent ME, the corresponding ME becomes stable and remains so for all higher migration rates. We summarize these findings as follows: The proposition shows that, except for the nongeneric case when (3.26) holds, there is always precisely one stable equilibrium point. Numerical results support the conjecture that the stable equilibrium is globally asymptotically stable. Bifurcation patterns as functions of m are derived in Sect. 4.8.
In addition to F 0 , there exists a second FP on the edge connecting M 2 and M 3 . Although its coordinates can be calculated easily, it is not of interest here as it is unstable for every choice of selection and migration parameters. This unstable equilibrium leaves the state space under small perturbations, i.e., if ρ > 0.

Highly asymmetric migration
All special cases treated above suggest that there always exists a globally asymptotically stable equilibrium. This, however, is generally not true as was demonstrated by the analysis of the two-locus continent-island (CI) model in Bürger and Akerman (2011). There, all possible bifurcation patterns were derived and it was shown that the fully polymorphic equilibrium can be simultaneously stable with a boundary equilibrium. For highly asymmetric migration rates, the equilibrium and stability structure can be obtained by a perturbation analysis of this CI model. Therefore, we first summarize the most relevant features of the analysis in Bürger and Akerman (2011). Because in that analysis the haplotype A 2 B 2 is fixed on the continent (here, deme 2) and there is no back migration (m 2 = 0), it is sufficient to treat the dynamics on the island (here, deme 1) where immigration of A 2 B 2 occurs at rate m 1 . Thus, the state space is S 4 .
It was shown that up to two internal (fully polymorphic) equilibria, denoted by E + and E − , may exist. Only one (E + ) can be stable. Two SLPs, E A and E B , may exist. At E A , locus A is polymorphic and allele B 2 is fixed; at E B , locus B is polymorphic and allele A 2 is fixed. E A (E B ) is admissible if and only if m 1 < α 1 (m 1 < β 1 ). E A is always unstable. Finally, there always exists the monomorphic equilibrium E C at which the haplotype A 2 B 2 is fixed on the island. The equilibrium coordinates of all equilibria were obtained explicitly. In addition, it was proved (see also Bank et al. 2012, Supporting Information, Theorem S.4) that precisely the following two types of bifurcation patterns can occur: Type 1 There exists a critical migration rate m • > 0 such that: • If 0 < m 1 < m • , a unique internal equilibrium, E + , exists. It is asymptotically stable and, presumably, globally asymptotically stable.
Type 2 There exist critical migration rates m • and m • satisfying m • > m • > 0 such that: . It is asymptotically stable and, presumably, globally stable. • At m 1 = m • , an unstable equilibrium (E − ) enters the state space by an exchangeof-stability bifurcation with a boundary equilibrium (E B or E C ). • If m • < m 1 < m • , there are two internal equilibria, one asymptotically stable (E + ), the other unstable (E − ), and one of the boundary equilibria (E B or E C ) is asymptotically stable.
• At m 1 = m • , the two internal equilibria merge and annihilate each other by a saddle-node bifurcation. • If m 1 > m • , a boundary equilibrium (E B or E C ) is asymptotically stable and, presumably, globally stable.
For sufficiently large migration rates (m > m •• ≥ m • ), E C is globally asymptotically stable in both cases. Bifurcation patterns of Type 2 occur only if the recombination rate is intermediate, i.e., if ρ is about as large as α 1 .
By imbedding the CI model into the two-deme dynamics, (2.5) or (2.7), perturbation theory can be applied to obtain analogous results for highly asymmetric migration, i.e., for sufficiently small m 2 /m 1 (Karlin and McGregor 1972). This is so because all equilibria in the CI model are hyperbolic except when collisions between equilibria occur (Bürger and Akerman 2011). Since the coordinates of the internal equilibria E + and E − were derived, the perturbed equilibrium frequencies can be obtained. Because they are too complicated to be informative, we do not present them. The perturbation of E + , denoted by F, is asymptotically stable. As E − is internal, it cannot be lost by a small perturbation. Also the boundary equilibria and their stability properties are preserved under small perturbations. In particular, E C gives rise to M 4 , and the SLPs E A and E B give rise to P A,2 and P B,2 , respectively, If recombination is intermediate, (at least) under highly asymmetric two-way migration, one stable and one unstable FP can coexist. In this case the stable FP, F, is simultaneously stable with either M 4 or P B,2 . Although there is precisely one (perturbed) equilibrium in a small neighborhood of every equilibrium of the CI model, we can not exclude that other internal equilibria or limit sets are generated by perturbation.
3.9 The case θ = 0 The analyses in the previous sections are based on the assumptions (2.3), in particular, on θ > 0. However, many of the results obtained above remain valid if θ = 0. Here, we point out the necessary adjustments.
Without loss of generality, we can assume in addition to θ = 0 and (2.3a). Then we observe that Therefore, either applies, where equality in (3.36a) and (3.36b) holds if α k = β k (k = 1, 2). In addition, With these preliminaries, we can treat the changes required in the above propositions if θ = 0.
As already noted, Proposition 3.2 remains valid independently of the value of θ . In Proposition 3.3, the only SLPs through which the internal equilibrium F ∞ can leave the state space are P B,1 and P B,2 ; see (3.19c) and (3.19d). The reason is that, except when σ 1 +σ 2 = 0 (and (3.37) applies), P A,1 and P A,2 are only admissible if P B,1 and P B,2 are. Thus, the locus under weaker selection always becomes monomorphic at lower rates of gene flow than the locus under stronger selection.
In the highly symmetric case of (3.36c), SLPs cannot be lost. Thus, F ∞ is always admissible and globally stable, cf. Proposition 3.3. If ρ = 0, (3.37) implies that F 0 exists always (and is stable). In the next section we show that in this highly symmetric case the FP is always admissible for arbitrary recombination rates.

The super-symmetric case
In many, especially ecological, applications highly symmetric migration-selection models are studied. Frequently made assumptions are that the migration rates between the demes are identical (m 1 = m 2 ), selection in deme 2 mirrors that in deme 1 (α k = −α k * ), and the loci are equivalent (α k = β k ). Thus, θ = 0 and (3.36c) holds, which we assume now.
Conditions (3.7) and (3.8) imply that all four SLPs are admissible. Hence, all monomorphisms are unstable. In addition, it can be proved that all SLPs are unstable (Appendix A.6). If migration is weak, a globally asymptotically stable, fully polymorphic equilibrium (F) exists (Proposition 3.2).
Because every boundary equilibrium is hyperbolic for every parameter choice, the index theorem of Hofbauer (1990) can be applied. Since none of the boundary equilibria is saturated, it follows that an internal equilibrium with index 1 exists. For small migration rates, this is F because it is unique. Since the boundary equilibria are always hyperbolic, no internal equilibrium can leave the state space through the boundary.
However, we cannot exclude that the internal equilibrium undergoes a pitchfork or a Hopf bifurcation. Numerical results support the conjecture that the internal equilibrium is unique and globally attracting, independently of the strength of migration. This is a very special feature of this super-symmetric case; cf. Proposition 4.3.

General case
Because a satisfactory analysis for general parameter choices seems out of reach, we performed extensive numerical work to determine the possible equilibrium structures. In no case did we find more complicated equilibrium structures than indicated above, i.e., apparently there are never more than two internal equilibria. If there is one internal equilibrium, it appears to be globally asymptotically stable. If there are two internal equilibria, then one is unstable and the other is simultaneously stable with one boundary equilibrium (as in the CI model). Apparently, two internal equilibria occur only for sufficiently asymmetric migration rates and only if the recombination rate is of similar magnitude as the selection coefficients.
For low migration rates as well as for high recombination rates, there is a unique, fully polymorphic equilibrium which is globally asymptotically stable and exhibits positive LD (Sects. 3.4 or 3.6). We denote the (presumably unique) asymptotically stable, fully polymorphic equilibrium by F. If migration is weak, or recombination is weak, or recombination is strong, we have proved that F is unique. Useful approximations are available for weak migration or strong recombination; see (3.16) or (3.21). Finally, for sufficiently high migration rates one of the monomorphic equilibria is globally asymptotically stable.

Bifurcation patterns and maintenance of polymorphism
Here we study how genetic variation and polymorphism depend on the strength and pattern of migration. In particular, we are interested in determining how the maximum migration rate that permits genetic polymorphism depends on the other parameters. For this end, we explore properties of our model, such as the possible bifurcation patterns, as functions of the total migration rate m. We do this by assuming that α 1 , α 2 , β 1 , β 2 , ρ, and the migration ratio where m > 0 and 0 ≤ φ ≤ 1, are constant. The values φ = 0 and φ = 1 correspond to one-way migration, as in the CI model. If φ = 1 2 , migration between the demes is symmetric, an assumption made in many studies of migration-selection models. Fixing φ and treating m as the only migration parameter corresponds to the migration scheme introduced by Deakin (1966).

Important quantities
We define several important quantities that will be needed to describe our results and summarize the relevant relations between them. Let and The quantities m A , m B , and m F 0 yield the bounds for the intervals of total migration rates m in which the SLPs at A, B, and the polymorphic equilibrium F 0 , respectively, are admissible: Here, the left and the right inequalities correspond, and we have (4.7) The quantities m M 1 and m M 4 occur in the stability conditions of the monomorphic equilibria M 1 and M 4 (Proposition 4.1), and m * determines the range of stability of F 0 ; see (4.56). They satisfy (4.8) We note that m A , m B , m F 0 , and m M 4 assume their minima if φ = 0 and their maxima if φ = 1, whereas m M 1 assumes its minimum or maximum at φ = 1 or φ = 0, respectively. m * is a convex function of φ, and symmetric around its minimum φ = 1/2. The definitions of (several of) the quantities φ X are motivated by the following relations: where we have (4.10) The following relations apply to m * : where we derived (4.11b) and (4.11c) from (4.9c) and (4.9f) using (4.10).
In the following, we summarize the most important inequalities between the quantities φ X : (4.14) They can be derived straightforwardly from their definitions and our general assumption (2.3). Finally, if ρ = 0 andθ = α 1 α 2 − β 1 β 2 , the following relations hold: Additional relations that are needed only in the proofs may be found in Appendix A.7.

Admissibility of SLPs
We begin by expressing the conditions for admissibility of the SLPs in terms of the total migration rate m and the migration ratio φ. Since, by (3.7), (3.8), and (4.4), every SLP is admissible if m is sufficiently small and leaves the state space at a uniquely defined critical migration rate, it is sufficient to determine this critical rate and the monomorphism through which it leaves the state space. Using (4.3a), (4.3b), (4.5), and (4.4), we infer from (3.9) that We observe that locus A is polymorphic and locus B is monomorphic if and only if If β 2 < α 2 , we infer from (10.18c) and (10.18d) that (4.18) holds if and only if Therefore, (2.4) implies that if locus A is under weaker selection than locus B in both demes (|β k | > |α k |), then there is a range of values φ and m such that locus A is polymorphic whereas B is monomorphic. This is in contrast to the CI model or highly asymmetric migration rates or θ = 0, where it is always the locus under weaker selection that first loses its polymorphism while m increases. This is a pure one-locus result and a consequence of the classical condition for a protected polymorphism, e.g., (3.7). With two-way migration, a locus with alleles of small and similar (absolute) effects in the demes (α 1 ≈ −α 2 ) may be maintained polymorphic for higher migration rates than a locus with alleles of large and very different (absolute) effects.

Stability of monomorphic equilibria
Here, we reformulate the stability conditions of the ME derived in Sect. 3.2 in terms of m and φ.
If in these conditions one inequality is reversed, the corresponding equilibrium is unstable.
Proof We prove only that the statement about M 1 is equivalent to that in Proposition 3.1. The others follow analogously or are immediate.
). Invoking (10.31), we can combine conditions (a), (b), and (c) to obtain (4.25a). Statement (ii) follows directly from (10.18f) and (10.35a). Statement (iii) follows by observing that only internal equilibria in LD will depend on ρ, the factor t 3 (10.4) in the characteristic polynomial at M 1 is the only one that depends on ρ, and t 3 gives rise to an eigenvalue zero if and only if m = m M 1 . An analogous argument holds for M 4 .
The asymmetry between (4.25) and (4.26) results from the fact that α 1 < β 1 is assumed, whereas β 2 < α 2 or β 2 ≥ α 2 is possible. The reader may recall the comments made below Proposition 3.1. In addition, we note that if the fitness parameters and ρ and φ are fixed, a stable ME remains stable if m is increased. This is not necessarily so if m and φ are varied simultaneously. For related phenomena in the one-locus case, see Karlin (1982) and Nagylaki (2012). In Sect. 4.5, we will prove global convergence to one of the asymptotically stable ME if m is sufficiently large.

Weak migration
We recall from Proposition 3.2 that for sufficiently weak migration, there is a fully polymorphic equilibrium, it is globally asymptotically stable, and exhibits positive LD in both demes.

Proposition 4.3 For sufficiently large m, one of the monomorphic equilibria
Proof The proof is based on the perturbation results about the strong-migration limit in Section 4.2 of Bürger (2009a). The strong-migration limit is obtained if max k=1,2 {|α k |, |β k |, ρ}/m → 0. In this limit, the demes become homogeneous and the system of differential equations (2.7) converges to a system, where in each demė are the spatially averaged selection coefficients and averaging is performed with respect to the Perron-Frobenius eigenvector (1 − φ, φ) of the migration matrix (see Section 4.2 in Bürger 2009a for a much more general treatment starting with a multilocus model in discrete time). Therefore, Proposition 4.10 in Bürger (2009a) applies and, provided m is sufficiently large, all trajectories of (2.7) converge to a manifold on which the allele frequencies and the linkage disequilibria in both demes are nearly identical. In addition, in the neighborhood of each hyperbolic equilibrium of (4.27) there is exactly one equilibrium of (2.7), and it has the same stability.
In the present case, the conclusion of Proposition 4.10 in Bürger (2009a) can be considerably strengthened. Because the system (4.27) describes evolution in an ordinary two-locus model under genic selection, the ME representing the gamete of highest fitness is globally asymptotically stable. In fact, (4.27) is also a generalized gradient system for which Lemma 2.2 of Nagylaki et al. (1999) holds. Therefore, the analog of statement (c) in Theorem 4.3 of Bürger (2009a) applies and yields global convergence to the unique stable equilibrium.
Finally, it is an easy exercise to show that, in the strong-migration limit, i.e., with fitnesses averaged according to (4.28), gamete Since there is no dominance, the corresponding ME is the unique stable equilibrium.

Linkage equilibrium
We shall establish all possible equilibrium configurations and their dependence on the parameters under LE. In Fig. 2, the equilibrium configurations are displayed as schematic bifurcation diagrams with the total migration rate m as the bifurcation parameter. In Theorem 4.4, we assign to each diagram its pertinent parameter combinations.
In order to have only one bifurcation diagram covering cases that can be obtained from each other by simple symmetry considerations but are structurally equivalent otherwise, we use the sub-and superscripts X and Y in the labels of Fig. 2. For an efficient presentation of the results, we define  C. Figure 3 shows the order in which the bifurcation diagrams of Fig. 2 arise if φ is increased from 0 to 1.
Proof We prove parts A and B simultaneously, essentially by rewriting the conditions in Proposition 3.3 on admissibility and stability of the equilibria in terms of m, m A , and m B (4.3).  19) and (4.4), we infer easily: (4.39g) Invoking the relations (10.18), we can rewrite conditions (4.39a), (4.39c)-(4.39f) in the form (4.40a) F ∞ → P B,1 ⇐⇒ m ↑ −m A and β 2 < α 2 and φ <φ AB , (4.40b) (4.40e) We conclude immediately that (4.40b) applies in case ( The bifurcations of equilibria that cannot be stable can be derived easily from Sects. 4.2 and 4.3 and the above theorem by noting that these are boundary equilibria and corresponding pairs of SLPs are admissible for the same parameters; see (3.7) and (3.8). Inclusion of these bifurcations would require the introduction of subcases.

Corollary 4.5 Under the assumption of LE, the maximum migration rate, below which a stable two-locus polymorphism exists, is given by
The corollary is a simple consequence of Proposition 3.3 and (4.40).

Strong recombination: quasi-linkage equilibrium
We recall from Sect. 3.6 that for sufficiently strong recombination, global convergence to the unique stable equilibrium occurs. From the coordinates (3.21) of the perturbed internal equilibrium, which is in quasi-linkage equilibrium, approximations could be derived for the critical migration rates at which the internal equilibrium collides with a boundary equilibrium and leaves the state space. It is not difficult to check with Mathematica that for large ρ, F collides with We note that m ρ max (P B,2 ) > 0 if and only if φ > φ AB , as is expected from (4.40c). Closer examination of (4.44) reveals that both m ρ max (P B,2 ) > m A and m ρ max (P B,2 ) < m A may hold.
Thus, the fully polymorphic equilibrium may be maintained for higher or lower migration rates than in the case of LE. This does not conform with the intuitive expectation that for reduced recombination m ρ max (P B,2 ) > m ∞ max should hold because the locally adapted haplotypes (A k B k in deme k) are less frequently broken apart. However, numerical evaluation of (4.44) shows that m ρ max (P B,2 ) < m ∞ max occurs only for about 3% of the admissible parameter combinations and if it holds, m ρ max (P B,2 ) is only very slightly less than m ∞ max (results not shown). If ρ is about as large as the largest selection coefficient or smaller, m max increases with decreasing ρ. Expressions analogous to (4.44) can be obtained for collisions of F with the other equilibria.

No recombination
Our aim is to establish all possible equilibrium configurations and their dependence on the parameters if recombination is absent. In Fig. 4, the equilibrium configurations are displayed as schematic bifurcation diagrams with the total migration rate m as the bifurcation parameter. In Theorem 4.6, we assign to each diagram its pertinent parameter combinations.
In order to have only one bifurcation diagram covering cases that can be obtained from each other by simple symmetry considerations but are structurally equivalent otherwise, we use the sub-and superscripts X and Y in the labels of Fig. 4. For an efficient presentation of the results, we define Theorem 4.6 Let ρ = 0. Figure 4 shows all possible bifurcation diagrams that involve bifurcations with equilibria that can be stable for some m given the other parameters. 1. Diagram (a) in Fig. 4 applies if one of the following two cases holds: Fig. 4 applies if one of the following two cases holds: (4.46b) Fig. 4 applies if one of the following two cases holds:

Diagram (e) in
(4.47b) Fig. 4 applies if one of the following four cases holds:  Fig. 4 applies if one of the following two cases holds:

Diagram (b) in
(4.49b) Fig. 4 applies if one of the following two cases holds:

Diagram (d) in
(4.50b) Fig. 4 applies if one of the following two cases holds:

Diagram (h) in
(4.52) Fig. 4 applies if one of the following two cases holds: Fig. 4 applies if one of the following two cases holds:
Next, we treat the bifurcations of the SLPs. The SLPs are admissible in intervals of the form 0 < m < |m A | or 0 < m < |m B | and leave the state space upon collision with a ME (Sect. 4.2). From (3.32) we conclude by simple calculations that P A,1 is asymptotically stable if and only if as is the case in (4.46a), (4.47a), (4.48a) (if min{φ F 0 , φ AB } = φ AB ), and (4.48b), as well as in (4.50a), (4.51a), and (4.53b). From (3.33), we conclude that P B,2 is asymptotically stable if and only if as is the case in (4.46b), (4.47b), (4.48c), and (4.48d) (if max{φ F 0 , φ AB } = φ AB ), as well as in (4.50b), (4.51b), and (4.53a). It remains to study the stability of the ME. For ρ = 0, we infer from Sect. 4.1 and Proposition 4.1: In conjunction with the above results on F 0 and the SLPs, this shows that, except in the degenerate cases (4.49), (4.50), (4.52), and (4.54), a ME becomes stable through a transcritical bifurcation with either F 0 , P A,1 , or P B,2 . In particular, M 1 becomes asymptotically stable for large m if (4.45a), (4.46a), or (4.49a) applies, M 3 becomes asymptotically stable if one of (4.47), (4.48), (4.51), (4.52), (4.53), or (4.54) applies, and M 4 becomes asymptotically stable if (4.45b), (4.46b), or (4.49b) applies. If φ = φ A or φ = φ B (4.50), then P A,1 or P B,2 , respectively, is admissible and asymptotically stable for every m, and every ME is unstable. This finishes the proof of parts A and B. Part C of Theorem 4.6 follows immediately from parts A and B by applying the relations in (4.15).
This theorem demonstrates that, for given selection parameters, the equilibrium structure, hence also the evolutionary dynamics, depends strongly on the degree φ of asymmetry of the migration rates. However, it is also important to note (and maybe counter intuitive) that for symmetric migration (φ = 1/2) any of the ten possible bifurcation diagrams may apply, simply by choosing the selection parameters accordingly.
The bifurcations of equilibria that cannot be stable can be derived easily from Sects. 4.2, 4.3, 3.7, and the above theorem by noting that these are boundary equilibria and corresponding pairs of SLPs are admissible for the same parameters. Inclusion of these bifurcations would require the introduction of subcases. In particular, P A,2 , P B,1 , and M 2 are always unstable because gamete A 1 B 2 is eventually lost. We observe from (10.20) and (10.38) that at most one pair of SLPs can be admissible if F 0 is either unstable or not admissible. If this is the case, then one of these SLPs is asymptotically stable (Fig. 4).
Corollary 4.7 If ρ = 0, the maximum migration rate, below which a stable two-locus polymorphism exists, is given by The corollary follows from the arguments surrounding (4.56) and (4.57).

Weak recombination
If m = m * , a regular perturbation analysis of F 0 yields the coordinates of a fully polymorphic (internal) equilibrium to leading order in ρ. This equilibrium, F, is asymptotically stable (Karlin and McGregor 1972). We denote the first-order approximation of F by F ρ . Therefore, we have Because the coordinates of F ρ are much too complicated to be informative, we refrain from presenting them. For sufficiently small ρ, the following properties of F ρ (hence, of F) can be inferred from Proposition 4.1, Remark 4.2, and Theorem 4.6, Part A.1: (4.62b) The above perturbation analysis can not be used to investigate the properties of the internal equilibrium F for given small positive ρ when m is varied in the proximity of m * . Therefore, we performed numerical calculations to study the fate of F when ρ is small and fixed, and m increases. It suggests the following: where m * A and m * B are close to m * . Thus, if ρ is small, F stays close to F 0 as m increases from 0 until a value close to m * is reached. Then, within a very short interval of m, F moves 'quickly' along the manifold given by (10.8) and (10.11) to one of the boundary equilibria (P A,1 , P B,2 , or M 3 ) on the 'opposite' side of the state space, where it exchanges stability upon collision with the respective equilibrium (at m * A , m * B , or m * ). F appears to be asymptotically stable whenever it is admissible.
If one of the cases in (4.62) applies, then F 0 can be maintained for higher migration rates than F because m M 1 and m M 4 are decreasing functions in ρ. Numerical investigations support the conjecture that F 0 can be maintained for higher migration rates than F whenever recombination is weak but positive. Thus, when recombination is weak, decreasing ρ increases the maximum migration rate below which a stable, fully polymorphic equilibrium can be maintained.

Highly asymmetric migration
As already discussed in Sect. 3.8, by introducing weak back migration (i.e., φ close to 0 or 1) to the CI model, every equilibrium in the CI model gives rise to a unique equilibrium in a small neighborhood. This (perturbed) equilibrium has the same stability as the unperturbed. For weak or strong recombination, we can strengthen this conclusion. Because the CI model with ρ = 0 is a generalized gradient system (Bürger and Akerman 2011, Section 3.4.4) and the LE dynamics (3.17) has a globally asymptotically stable equilibrium (Theorem 4.4), the proof of Theorem 2.3 of Nagylaki et al. (1999) applies and shows that in both cases the global dynamics remains qualitatively unchanged under small perturbations. In particular, no new equilibria or limit sets are generated by a small perturbation.
Therefore, if ρ is sufficiently small and φ is sufficiently close to 0 or 1, we infer from Section 3.8 and Theorem 2 in Bürger and Akerman (2011) that the following bifurcation pattern applies (where i = 1 or 4): • If 0 < m < m M i , a unique internal equilibrium, F, exists. It is globally asymptotically stable. • At m = m M i , F leaves the state space through the ME M i by an exchange-ofstability bifurcation.
This pattern is displayed in diagram (a) of Fig. 4, where F 0 needs to be substituted by F. We conjecture that it applies whenever ρ is sufficiently small and either φ < φ M 1 or φ > φ M 4 holds. The bounds φ M 1 and φ M 4 follow from Remark 4.2 becauseφ M 1 is not needed if ρ is sufficiently small; see (10.32b). However, the upper bounds for ρ given in Remark 4.2 are, in general, too large to guarantee the above bifurcation pattern. This is known from the CI model in which the monomorphic equilibrium (M i ) may be simultaneously stable with the internal equilibrium F because an unstable internal equilibrium enters the state space at m = m M i through M i . If φ = 1, this may occur if 1 3 (α 1 + β 1 ) < ρ < 3α 1 − β 1 , cf. (4.65b). For φ = 0 or φ = 1, we have not been able to determine the upper bound for ρ below which F indeed leaves the state space through M i . Now we treat large ρ. Proposition 4.1 and Remark 4.2 show that if ρ > max{−α 2 , −β 2 }, then M 1 is asymptotically stable if and only if φ < φ A and m > max{−m A , −m B }, and if ρ > α 1 , then M 4 is asymptotically stable if and only if φ > φ B and m > max{m A , m B }.
If, in addition to ρ being sufficiently large, φ is small or large, then Theorem 4.4 implies that the internal equilibrium (F) leaves the state space through P A,1 , P B,1 , or P B,2 . The respective conditions are small perturbations of those given in (4.40a), (4.40b), of (4.40c), respectively. Combining theses conditions with those for the stability of the ME and observing (4.12) and (4.13), we conclude that the following bifurcation pattern applies if one of the conditions (a) α 2 ≤ β 2 and φ < φ A , or (b) β 2 < α 2 and φ <φ AB , or (c) φ > φ B holds approximately: • If 0 < m < m • , a unique internal equilibrium, F, exists. It is asymptotically stable.
• At m = m • , F leaves the state space through a SLP by an exchange-of-stability bifurcation.

Maintenance of polymorphism
As already noted in Sect. 3.11, for general parameters the equilibrium configurations could not be determined analytically. To explore the potential of spatially heterogeneous selection in maintaining genetic variation in the presence of gene flow, we investigate the maximum total migration rate, m max , that admits a stable, fully polymorphic equilibrium. We have already shown that m max = m ∞ max holds in the LE approximation (Corollary 4.5), and m max = m 0 max holds if ρ = 0 (Corollary 4.7). From (10.20) and (10.38) we conclude that where, as is not difficult to show, equality holds if and only if φ = φ AB . For the CI model with φ = 1, Proposition 1 in Bürger and Akerman (2011) yields In this case, the fully polymorphic equilibrium is globally asymptotically stable if (4.65a) or (4.65c) apply, but only locally stable if (4.65b) holds and m is close to m max . A formula analogous to (4.65), but with −α 2 and −β 2 instead of α 1 and β 1 , holds if φ = 0.
In general, we have no explicit formula for m max . However, extensive numerical work, as well as (4.65)  Here, m max is obtained by determining numerically the critical migration rate when the stable internal equilibrium hits the boundary. This is done by computing when the leading eigenvalue at the boundary equilibrium is zero and by calculating the coordinates of the fully polymorphic equilibrium in a small neighborhood. In a and b, we haveφ AB = 1 4 (indicated by the kink in the dashed line in a), φ M 1 = 5 17 , φ A = 1 3 , φ AB = 3 8 , φ B = 1 2 , φ M 4 = 5 8 . In c and d, we have holds always. This is illustrated by Fig. 6, which displays the dependence of m max on the migration ratio φ (Fig. 6a, c) and on the recombination rate ρ (Fig. 6b, d) for two selection regimes. In Fig. 6a, b, locus B is under stronger selection in both demes. In Fig. 6c, d, each locus is under stronger selection in one deme. In Fig. 6a, c, m ∞ max and m 0 max are shown as functions of φ. The inequality (4.64) is a conspicuous feature in both cases. Also the shapes of m ∞ max and m 0 max are conspicuous. The following properties are easy to prove: m ∞ max is not differentiable at φ = φ AB and φ =φ AB , and m 0 max is not differentiable at The latter condition is equivalent to (3.36c). Therefore, the analysis in Sect. 3.10 applies and shows that an internal equilibrium, which presumably is globally asymptotically stable, exists always, i.e., m max = m ∞ max = m 0 max = ∞. Figure 6b, d illustrates the effect of recombination on m max for different values of φ. In all cases investigated, m max decreased monotonically with increasing ρ. These findings support the conjecture that (4.66) is always valid. Therefore, m 0 max − m ∞ max serves as a useful estimate for the sensitivity of m max to variation in ρ. We can prove that m 0 max − m ∞ max is maximized at φ = φ M 1 , or at φ = φ M 4 , or at φ = 0 if β 2 < α 2 . Although we proved that m max < m ∞ max can occur (Sect. 4.7), all numerical examples showed that m max is only very slightly smaller than m ∞ max in this case (results not shown). Therefore, our results suggest that the cases of LE (infinitely strong recombination) and of no recombination 'essentially' bracket the range of parameters for which both loci can be maintained polymorphic.
As Fig. 6a, c shows, the range of values φ for which the equilibrium structure can be expected to be similar to the CI model, i.e., φ < min{φ M 1 ,φ AB } or φ > φ M 4 (Sect. 4.10), can vary considerably.
Finally, we infer from Proposition 4.1 that none of the ME is stable if m < max{|m A |, |m B |}. Hence, in this case at least one locus is maintained polymorphic. By contrast, we have shown in Sect. 4.2 that no SLP is admissible if m > max{|m A |, |m B |}. However, as demonstrated by our results for ρ = 0, an internal equilibrium may be asymptotically stable if max{|m A |, |m B |} < m < m 0 max . These results suggest that no genetic variability can be maintained if This bound is best possible if ρ = 0. For sufficiently large ρ, the corresponding bound is max{|m A |, |m B |}.

Migration load and local adaptation
Here, we briefly investigate some properties of the migration load of the subpopulations and of the total population. We use these migration loads as simple measures for local adaptation (but see Blanquart et al. 2012). Mean fitness in deme k is given bȳ w k = α k (2 p k − 1) + β k (2q k − 1), with its maximum at α k + β k . Therefore, the migration loads in demes 1 and 2, defined as the deviation ofw k from its maximum, are given by Assuming that the subpopulations are of equal size, we define the load of the total population by L = 1 2 (L 1 + L 2 ). If migration is weak, we can calculate the migration load in each deme at the fully polymorphic equilibrium F (Proposition 3.2) to leading order in m 1 and m 2 . For deme 1, we obtain and an analogous formula holds for deme 2. Obviously, the migration load increases with increasing migration rates m 1 or m 2 , hence with m, in each of the demes and in the total population. Simple calculations show that each of the loads also increases with increasing recombination rate ρ if migration is weak. In general, however, the load in each deme does not always increase with increasing m. The reason is that for sufficiently strong migration, generically, first one locus, then one of the haplotypes becomes fixed (Proposition 4.3). If this is either A 1 B 1 or A 2 B 2 , then the load in the corresponding deme will vanish for high migration rates, whereas that in the other deme will be very high. In such a case, the load of the total population may also decrease with increasing m. This occurs for large migration rates (not far below m max ) and it can occur for completely linked loci as well as for loci in LE. In the CI model, the load always increases with the migration rate (Bürger and Akerman 2011) Finally, although L is increasing in ρ if migration is weak, this is not necessarily so if migration is strong. By using a grid of parameter combinations, we showed numerically that in about 0.34% of more than 10 6 combinations of α 1 , α 2 , β 1 , β 2 , m, and φ, the total load L at the equilibrium F ∞ is lower than that at F 0 (results not shown). Again, this occurs for high migration rates, not far below the value m ∞ max at which F ∞ leaves the state space. Then a population maintained fully polymorphic by tight linkage may have a higher total load than a population in which fixation of a locus or a haplotype is facilitated by high recombination. In all such cases, selection in one deme was (considerably) stronger than in the other, and in more than 70% of the cases, a specialist haplotype became fixed at very high migration rates. In summary, under a wide range of conditions in this model, reduced recombination is favored, but there are instances where increased recombination is favored (cf. Pylkov et al. 1998;Lenormand and Otto 2000).

F ST and differentiation
The most commonly used measure for quantifying differentiation in spatially structured populations is F ST . For diallelic loci, F ST can be defined as F ST = Var( p) p(1−p) , where Var( p) is the variance of the allele frequencies in the total population andp is the allele frequency averaged over the demes. Estimators of multilocus F ST are usually defined as weighted averages of one-locus F ST estimators (e.g., Weir and Cockerham 1984;Leviyang and Hamilton 2011). Here, we extend Nagylaki's (1998) approach and define a genuine multilocus version of F ST that measures the covariance of the frequencies of (multilocus) haplotypes. We restrict attention to the diallelic two-locus case, but the extension to multiple multiallelic loci is evident. A general multilocus theory of fixation indices will be developed elsewhere.
Let c k denote the proportion of the population in deme k, so that k c k = 1. Then the frequency of haplotype i in the entire population isx i = k c k x i,k . Because our subpopulations are randomly mating, the frequency of genotype i j in the entire population is given by x i x j = k c k x i,k x j,k . Following eqs. (6a) and (6b) in Nagylaki (1998), we define F ST,i j as a standardized measure of the covariance of the frequencies of haplotypes i and j: The multilocus, or haplotype, heterozygosity in the entire population can be defined ash where i runs over all haplotypes. If the entire population were panmictic, its multilocus heterozygosity would be Thus, 1 − h T is the probability that two gametes chosen at random from the entire population are the same haplotype. Following eq. (32) in Nagylaki (1998), we define F ST by Then F ST can be written as (6.5) in direct generalization of the classical formula given above. We focus on the dependence of the equilibrium value of F ST on the migration parameters m and φ and on the recombination rate ρ. Because we obtained the coordinates of the stable, fully polymorphic equilibrium equilibrium F explicitly only in special or limiting cases, explicit formulas for F ST can be derived only in these cases. For instance, if migration is weak, we obtain from (3.16) that, to leading order in m, Here, F ST increases with decreasing ρ, and decreases with increasing m. Thus, stronger linkage leads to increased differentiation if migration is weak. Figure 7 illustrates for two selection scenarios how F ST , evaluated at the stable, fully polymorphic equilibrium F, depends on the total migration rate m and the recombination rate ρ. In diagrams (a) and (c) of Fig. 7, it is assumed that locus B is under stronger selection than locus A in both demes. It shows that F ST usually declines with increasing migration rate. However, there are a few instances, where F ST increases if m is slightly below the migration rate at which the fully polymorphic equilibrium In panels a and c, locus B is under stronger selection in both demes (α 1 = 1 2 , α 2 = −1, β 1 = −β 2 = 2, θ = 1). In panels b and d, locus A is under stronger selection than B in deme 2, and locus B is under stronger selection than A in deme 1 (α 1 = −β 2 = 0.4, β 1 = −α 2 = 2, θ = 3.84). Note that in all cases, F ST is also monotone decreasing in ρ. For ρ = 0 and ρ = ∞ (LE), the lines are from numerical evaluation of (6.5) by substitution of the coordinates of F 0 (3.24) or F ∞ (3.18). For other values of ρ, the numerically determined coordinates of the internal equilibrium are used loses admissibility. In diagrams (a) and (c) of Fig. 7, differentiation between the populations experiences the fastest decline for weak migration (relative to the selection parameters), whereas this is not necessarily so in diagrams (b) and (d). There, F ST may experience its strongest decrease if migration is strong. Figure 7 also shows that at large migration rates, F ST may increase if the recombination rate increases, i.e., F ST is not minimized under linkage equilibrium. However, this occurs only for large recombination rates, i.e., larger than the largest selection coefficient. This is compatible with the finding in Sect. 4.11 that at high recombination rates, m max may (slightly) increase in ρ, and the finding in Sect. 5 that the load L may decrease with increasing m. We note that this 'aberrant' behavior of m max , L, and F ST does not necessarily occur for the same parameter combinations. Among more than 10 6 parameter combinations of α 1 , α 2 , β 1 , β 2 , m, and φ, we found no instance where F ST evaluated at the equilibrium F ∞ was higher than that at F 0 (results not shown). Importantly, if recombination is weak or migration is weak then F ST apparently always increases with tighter linkage.
Comparison of our multilocus F ST with averages of single-locus F ST values showed that the multilocus F ST declines somwehat faster at small migration rates than the averaged single-locus F ST . For large parameter regions, the qualitative behavior of these measures of differentiation is the same. Differences occur only for a subset of selection coefficients at high migration rates and high recombination rates. Finally, we mention that our multilocus F ST is a sensitive measure of differentiation only if the effective number of haplotypes is low. This parallels the well known fact that the classical F ST is a sensitive measure of differentiation only if the effective number of alleles is low (e.g., Nagylaki 1998Nagylaki , 2011. Thus, our multilocus F ST may be most useful if applied to short sequences of DNA. A thorough and more general study is in preparation.

Invasion of a locally beneficial mutant
Differentiation between subpopulations can be increased by the invasion of mutants that establish a stable polymorphism at their locus. Therefore, we consider a locus (A) at which a new mutant A 1 arises that is advantageous relative to the wild type A 2 in deme 1, but disadvantageous in deme 2. In terms of our model, we assume α 1 > 0 > α 2 . If locus A is isolated, this mutant can invade and become established in a stable polymorphism if and only if |σ 1 + σ 2 | < 1; cf. (3.7) and (3.9). Using m and φ, this condition can be rewritten as m < |m A |, (7.1) see (4.3a) and (4.4a), or We restrict attention to the case φ > φ A (4.2a) when the influx of the deleterious allele A 2 into deme 1 is sufficiently strong such that A 2 is protected. (The case φ < φ A is symmetric and more suitable to study invasion of A 2 under influx into deme 2 of A 1 which is deleterious there.) Then the mutant A 1 can invade if any of the following equivalent conditions hold: where φ inv > 1 if and only if m < α 1 . Thus, A 1 can always invade if m < α 1 . For the CI model (φ = 1), each of the conditions in (7.3) simplifies to the well known invasion condition m < α 1 (Haldane 1930). The conditions (7.3) show that invasion is facilitated whenever back migration is increased, either by keeping m 1 constant and increasing m 2 , or by fixing m and decreasing φ. For the CI model it was proved that invasion of a locally beneficial mutant is always facilitated by increased linkage to a locus in migration-selection balance (Bürger and Akerman 2011). In fact, mutants of arbitrarily small effect can invade provided they are sufficiently tightly linked to this polymorphic locus which may be considered as the background in which the new mutant appears.
Here, we investigate whether this is also the case with two-way migration. Thus, we assume that locus B is in migration-selection balance (which requires that analogs of (7.3) are satisfied for β 1 and β 2 ) and a locally beneficial mutant A 1 arises at the linked locus A. Hence, the model in Sect. 2 applies and we assume (2.3).
Because we are mainly interested in the invasion properties of mutants of small effect, we assume that locus B is under stronger selection than A, i.e., |α k | < |β k | in deme k = 1, 2. Before the mutant A 1 arises, the population is at the equilibrium P B,2 (where |τ 1 + τ 2 | < 1 must hold for admissibility; see Sect. 3). A 1 can invade if P B,2 is unstable. Since the eigenvalues determining external stability are zeros of a complicated quartic equation, the stability of P B,2 cannot be determined analytically. We expect that the new stable equilibrium that will be reached is the fully polymorphic equilibrium F. For the CI model, this was be proved in (Bürger and Akerman 2011). For the case of LE, it follows from Theorem 4.4. Figure 8 displays typical results about the invasion of the mutant A 1 . In Fig. 8a, the maximum recombination rate admitting invasion, denoted by ρ max , is shown as a function of φ. In the shaded region, A 1 can invade. If φ ≤ φ inv = 0.55, (7.3c) implies that A 1 can always invade. If φ > φ inv , there exists ρ max < ∞, such that A 1 can invade only if ρ < ρ max , i.e., if A 1 is sufficiently tightly linked to locus B. In Fig. 8b, the minimum selection coefficient α 1 necessary for invasion of A 1 is shown as a function of ρ/m for various values of φ. These values are obtained by computing when the leading eigenvalue that determines external stability of P B,2 equals zero.
We conclude that, as in the CI model, mutants of arbitrarily small effect can invade provided they are sufficiently tightly linked to a locus that is already maintained in migration-selection balance. In addition, as shown by both panels in Fig. 8, increasingly symmetric migration facilitates the invasion and establishment of locally beneficial alleles.

The effective migration rate at a linked neutral site
Linkage to loci under selection may impede or enhance gene flow at a neutral marker locus. In the first case, linkage may act as a barrier to gene flow. This was shown by the work of Petry (1983), Bengtsson (1985), Barton and Bengtsson (1986), and Charlesworth et al. (1997), who developed and studied the concept of the effective migration rate as a measure of the 'effective' gene flow at a neutral site. More recently, the effective migration rate was studied for CI models with selection on a single locus in a class-structured population (Kobayashi et al. 2008) or with selection on two linked loci (Bürger and Akerman 2011). Fusco and Uyenoyama (2011) investigated the consequences of a selectively maintained polymorphism on the rate of introgression at a linked neutral site under symmetric migration between two demes.
Here, we derive an explicit expression for the effective migration rate at a neutral locus (N) that is located between the two selected loci, A and B. Recombination between locus A (B) and the neutral locus occurs with rate ρ AN (ρ NB ) such that ρ = (a) (b) Fig. 8 Invasion properties of locally beneficial alleles. In a, the maximum recombination rate between loci A and B, below which invasion of A 1 can occur, is displayed as a function of φ. The parameters α 1 = −α 2 = 0.1, β 1 = −2β 2 = 2, and m = 1 are fixed. Therefore, φ A = 1 2 and φ inv = 0.55. In b, the minimum selective advantage α 1 required for invasion of A 1 is shown as a function of ρ for different values of φ. The parameters α 2 = −0.1, β 1 = −2β 2 = 2, and m = 1 are fixed ρ AN + ρ NB . Thus, only one crossover event occurs in a sufficiently small time interval. We assume that ρ AN and ρ NB are positive, i.e., the neutral locus is not completely linked to a selected site. We consider two variants at the neutral locus, N 1 and N 2 , each with arbitrary, positive initial frequency in at least one deme. The frequency of N 1 in deme k(= 1, 2) is denoted by n k . We model evolution at the three loci by a system of 7 × 2 ordinary differential equations for the allele frequencies and linkage disequilibria ( p 1 , p 2 , q 1 , q 2 , D AB 1 , D AB 2 , n 1 , n 2 , ). We refrain from presenting the equations for the allele frequencies at the neutral locus and the associated linkage disequilibria because they are a straightforward extension of those in Section 4.6 of Bürger and Akerman (2011).
Obviously, the equilibrium allele frequencies at the neutral locus are the same in each deme and given by the initial allele frequencies averaged over the two demes: The equilibrium frequencies at the two selected loci are independent of the neutral locus and, thus, the same as in the two-locus model treated above. The linkage disequilibria involving the neutral locus (D AN k , D NB k , and D ANB k ) are zero at equilibrium. By (8.1), there is a one-dimensional manifold of equilibria resulting from the absence of selection at the neutral locus.
We assume that parameters are such that the fully polymorphic equilibrium F is admissible and globally asymptotically stable. Using the above order for the allele frequencies and linkage equilibria, the Jacobian at the equilibrium F has block structure, where J S is the Jacobian describing convergence of ( p 1 , p 2 , q 1 , q 2 , D AB 1 , D AB 2 ) to F, and J N is the Jacobian describing convergence of (n 1 , n 2 , ) to (n,n, 0, 0, 0, 0, 0, 0). Because zero is the leading eigenvalue of J N , the rate of convergence to equilibrium at the neutral locus is determined by the second largest eigenvalue of J N , which we denote by λ N . We define the effective (total) migration rate by m eff = −λ N (Bengtsson 1985;Kobayashi et al. 2008;Bürger and Akerman 2011). It can be checked that under weak migration, i.e., to leading order in m 1 and m 2 , one obtains (a Mathematica notebook is available on request). If the neutral site is linked only to one selected locus (e.g., because β 1 = β 2 = 0), then is obtained. Thus, two linked selected loci act as a much stronger barrier to gene flow than a single selected locus, especially if the recombination rate between the two loci is not much larger than the selective coefficients. In Fig. 9, the approximation (8.3) of the effective migration rate m eff is displayed as a function of m for various parameter combinations and compared with the exact value obtained by numerical evaluation of λ N . We note that m eff is (approximately) the sum of the two effective one-way migration rates (Bürger and Akerman 2011) and closely related to Kobayashi and Telschow's (2011) effective recombination rate. Our result complements their explicit example on two-locus incompatibilities. We refer to their paper for the discussion of the relation of this concept of an effective migration rate to that of Bengtsson (1985) and for applications in the context of speciation theory.

Discussion
The purpose of this investigation was to improve our understanding of how genetic architecture, in particular recombination and locus effects, as well as the pattern and amount of migration determine polymorphism, local adaptation, and differentiation in a subdivided population inhabiting a heterogeneous environment. For simplicity, we restricted attention to two linked, diallelic loci and to migration between two demes. The study of diversifying selection in just two demes may also shape our intuition about clinal variation if the two subpopulations are from different ends of the cline. If alleles are beneficial in only one environment and detrimental in the other, local adaptation of subpopulations and differentiation between them can be obtained only if a (multilocus) polymorphism is maintained. Therefore, most of our mathematical results focus on existence and stability of polymorphic equilibria and on the dependence of the equilibrium configurations on the model parameters (migration rates, selection coefficients, recombination rate).
The model is introduced in Sect. 2. Sections 3 and 4 are devoted to the derivation of the possible equilibrium configurations and bifurcation patterns. They contain our main mathematical results. Explicit analytical results about existence and stability of equilibria were obtained for several limiting or special cases and are complemented by numerical work.
The conditions for admissibility of all single-locus polymorphisms (SLPs) are given in Sect. 3.1, those for asymptotic stability of the monomorphic equilibria (ME) in Proposition 3.1 in Sect. 3.2. The stability of SLPs could not generally be determined (Sect. 3.3). Weak migration is treated by perturbation methods in Sect. 3.4. For sufficiently weak migration, there exists a globally attracting fully polymorphic equilibrium, F (Proposition 3.2). Its approximate coordinates are given by (3.16).
The complete equilibrium and stability structure could be derived under the assumption of linkage equilibrium (Sect. 3.5). The unique, fully polymorphic equilibrium F = F ∞ is admissible and globally attracting if and only if all four SLPs are admissible. Otherwise, one boundary equilibrium (SLP or ME) is globally asymptotically stable (Proposition 3.3). These results extend straightforwardly to an arbitrary number of diallelic loci. Based on these results, nonlinear perturbation theory establishes the existence of a globally stable, fully polymorphic equilibrium in a perturbed parameter range if recombination is sufficiently strong (Sect. 3.6). This equilibrium is in quasi-linkage equilibrium and given by (3.21).
Also for completely linked loci all equilibria and their local stability properties could be derived (Sect. 3.7). In this case, the fully polymorphic equilibrium F 0 (3.24) may lose stability while it is admissible (3.28). At this threshold a boundary equilibrium becomes stable by a 'jump bifurcation' (Proposition 3.4). In general, however, more complicated equilibrium patterns than determined by Propositions 3.3 and 3.4 can occur, in particular, multiple stable equilibria.
In Sect. 3.8, we apply perturbation theory to infer the equilibrium properties under highly asymmetric migration from those derived for the continent-island model in Bürger and Akerman (2011) and Bank et al. (2012). There, a stable (F) and an unstable fully polymorphic equilibrium may exist if recombination is intermediate, and F is simultaneously stable with a boundary equilibrium. In general (Sect. 3.11), we cannot exclude the existence of more than two internal equilibria or complicated dynamical behavior. Numerical searches produced no such instances. What can be shown easily is that, if ρ < ∞, any fully polymorphic equilibrium exhibits LD. In all cases, where an internal equilibrium was calculated (numerically or analytically), it exhibited positive LD.
In the super-symmetric case, in which selection in deme 2 mirrors that in deme 1 and migration is symmetric, an assumption made in several applications, a fully polymorphic equilibrium exists always and, presumably, is stable (Sect. 3.10). This is a highly degenerate situation because if θ = 0, only a monomorphic equilibrium can be stable for sufficiently large migration rates (Proposition 4.3). If θ = 0 (Sect. 3.9), then a fully polymorphic equilibrium can exist for arbitrarily large migration rates if φ = φ AB (see also Sect. 4.11).
Whereas in Sect. 3 the focus was on the efficient presentation of the existence and stability results of equilibria, in Sect. 4 these results are used to derive the possible bifurcation patterns with the total migration rate m as the bifurcation parameter. All possible bifurcation patterns could be derived under the assumption of LE (Theorem 4.4,Figs. 2,3), and under the assumption of complete linkage (Theorem 4.6,Figs. 4,5). The latter case is considerably more complex. Interestingly, in each case, every bifurcation pattern can occur for every ratio φ = m 1 /m of migration rates by choosing the selection coefficients appropriately. Hence, the assumption of symmetric migration does not yield simpler equilibrium configurations than general migration if arbitrary selection coefficients are admitted.
In each of these cases (LE or ρ = 0), we determined the maximum migration rate m max admitting an asymptotically stable, fully polymorphic equilibrium (Corollaries 4.5 and 4.7). The maximum migration rate m 0 max for ρ = 0 always exceeds or equals that (m ∞ max ) for LE, i.e., m ∞ max ≤ m 0 max . Although for strong recombination, m max can be very slightly smaller than m ∞ max (Sect. 4.7), in the vast majority of investigated cases, m max is bracketed by m ∞ max and m 0 max (Fig. 6, Sect. 4.11).
Proposition 4.3 demonstrates that a ME is globally attracting if migration is sufficiently strong (except in the degenerate case noted above). If we interpret the equilibria M 2 and M 3 as fixation of a generalist ( A 1 B 2 and A 2 B 1 are haplotypes of intermediate fitness), and M 1 and M 4 as fixation of a specialist (A 1 B 1 and A 2 B 2 are the locally adapted haplotypes), then depending on the sign of θ one of the generalists becomes fixed for high but note that, depending on the selection coefficients, both φ A and φ B can be arbitrarily close to 0 or 1.). The critical value m as well as φ A and φ B are independent of ρ. Otherwise, one of the specialists becomes fixed for large m.
The fact that a generalist becomes fixed for strong migration is a distinct feature of (balanced) two-way migration: in the CI model or if migration is sufficiently asymmetric (φ < φ A or φ > φ B if θ > 0), one of the specialist haplotypes swamps the populations and becomes fixed. Another difference between highly asymmetric and more symmetric migration patterns is that in the first case, it is always the locus under weaker selection that first loses its polymorphism while m increases, whereas this not necessarily so in the latter case (see Sect. 4.2 and Theorem 4.6,cases A3 and A4).
In summary, we determined quantitatively when the following three evolutionarily stable states discussed by Kawecki and Ebert (2004) occur: (i) existence of a single specialist optimally adapted to one deme and poorly to the other, (ii) existence of a single generalist type which has higher average fitness in the whole population than than any of the specialists, and (iii) existence of a set of specialists each adapted to its deme, i.e., coexistence in a polymorphism. Local adaptation and differentiation occur only in case (iii).
In Sect. 5, we used the migration load in each deme to quantify the degree of local adaptation. In Sect. 6 we introduced a new multilocus version of F ST to measure differentiation. If migration is weak, then local adaptation and differentiation decrease with increasing migration rate and increase with increasing linkage between the loci (Fig. 7). In particular, for given (small) migration rate, local adaptation and differentiation are maximized if the fitness effects are concentrated on a single locus (corresponding to ρ = 0 in our model). However, as discussed in Sect. 5, for high migration rates, the migration load of the total population can decrease with increasing recombination or migration rate. Similarly, at high recombination and migration rates, F ST can increase with increasing migration or recombination rate. Thus, for given, relatively high migration rate, F ST may be minimized at intermediate recombination rates. Apparently, it is always maximized in the absence of recombination.
In Sect. 7, we investigated the conditions for invasion of locally beneficial mutants. At an isolated locus, such a mutant can invade and become established in a migrationselection equilibrium if and only if its advantage exceeds a threshold that increases with the immigration rate of the wild type; see (7.3b). If, however, this mutant occurs at a locus that is linked to a locus that is already in migration-selection balance, then its invasion is facilitated, i.e., its local selective advantage can be smaller (Fig. 8b). Equivalently, for given selection coefficients and total migration rate, the minimum recombination rate needed for invasion increases if φ, or the influx of the (deleterious) wild type relative to the efflux of the new mutant, increases (Fig. 8a). For the extreme case of one-way migration from a 'continental' population to an 'island' population that is adapting to a new environment, Bürger and Akerman (2011) proved that invasion of a locally beneficial mutant is always facilitated by increased linkage to a locus in migration-selection balance.
Thus, our results complement the numerical finding by Yeaman and Whitlock (2011) for a multilocus quantitative-genetic model that clusters of locally adaptive mutations, or concentrated genetic architectures, build up in spatially structured populations with opposing selection pressures in two demes. Because tighter linkage is required for invasion under increasingly asymmetric migration rates, more concentrated architectures and a greater advantage for recombination-reducing mechanisms (such as chromosome inversions) should be expected for highly asymmetric migration. In finite populations, invasion of new mutants occurs only with a certain probability, and genetic drift may erase polymorphism. Numerical work, supported by analytical methods, has already shed some light on the dependence of the probability of establishment of new, locally adaptive mutations on the recombination rate and other factors (Yeaman and Otto 2011;Feder et al. 2012). Analytical work on the role of genetic drift and finite population size on these issues is in progress.
Our results also show that, in the absence of epistasis and under the present form of balancing selection, reduced recombination between selected loci is favored, except when migration rates are sufficiently symmetric and high (Sect. 5). Selection inducing certain forms of epistasis may favor high recombination in structured populations more easily (Pylkov et al. 1998;Lenormand and Otto 2000;Bank et al. 2012). Therefore, general predictions about the emergence of clusters of locally adaptive mutations in regions of reduced recombination, or of genomic islands of speciation (Wu and Ting 2004) or of differentiation (Feder et al. 2012), can not be made in the absence of detailed information about epistasis and the spatial pattern of selection and migration. At least in the absence of epistasis, the most favorable situation for the emergence of such clusters should occur in populations that are adapting to a new environment, still receiving maladaptive gene flow but sending out only very few or no migrants (corresponding to a continent-island model).
In Sect. 8, we derived the approximation (8.3) for the effective migration rate at a linked neutral locus that is located between the selected loci. This approximation is simply the sum of the two effective migration rates under one-way migration (Bürger and Akerman 2011). Because in the present model, polymorphism at the selected loci is maintained by balancing selection, the effective migration rate may be greatly reduced compared with the actual migration rate (see Fig. 9). Thus, strong barriers against gene flow may build up at such neutral sites and enhance (neutral) differentiation (see Charlesworth and Charlesworth 2010, Chap. 8.3). Future work will have to study the actual amount and pattern of neutral diversity at such sites in finite populations.
In the following, we derive the stability conditions (3.10) and (3.11) for M 1 . Those for M 4 can be deduced analogously or by symmetry considerations by taking into account that (2.3b) implies min{α 1 , β 1 } = α 1 . The stability analysis of M 2 and M 3 is much simpler and left to the reader.
A.3 Calculation of equilibria with two polymorphic loci if ρ = 0 As shown in the main text, by Corollary 3.9 of Nagylaki and Lou (2007) it is sufficient to assume that A 1 B 2 is absent, which implies D k = p k (1 − q k ) and p k ≤ q k . Setting ρ = 0, we find from the equationsṗ 1 = 0 andq 1 = 0 (2.7) that holds at equilibrium. Substituting (10.8) intoṗ 2 andq 2 , we obtain at equilibrium, where g 1 and g 2 are quadratic polynomials in ( p 1 , q 1 ). The obvious substitution results in the equilibrium condition (10.10a) It is easy to check that F 0 always fulfills this condition and it is the only solution satisfying 0 ≤ p 1 = q 1 ≤ 1. Hence, unless there is curve ( p 1 , q 1 ) of solutions of (10.10) that passes through F 0 and through either a point on p 1 = 0 with 0 < q 1 ≤ 1 or on q 1 = 1 with 0 ≤ p 1 < 1, F 0 is the unique admissible solution of (10.10).
Because F 0 has an eigenvalue 0 only if either (3.26) is satisfied or if |κ 1 + κ 2 | = 1 (which occurs if and only if F 0 collides with either M 1 or M 4 ), F 0 is the only equilibrium with both loci polymorphic, except when (3.26) is satisfied. In the latter case, a line of equilibria exists, as we show now.
We calculate m 2 from (3.26) and substitute into (10.10). The right-hand side factorizes into two linear terms. Only one of them gives rise to admissible equilibria and, in fact, yields the manifold: where 0 ≤ q 1 ≤ 1. The allele frequencies in the other deme are obtained from (10.8).
It is straightforward to check that not only F 0 , but also the equilibria P A,1 and P B,2 lie on this manifold. In terms of the gamete frequencies, this manifold is a straight line.
A.4 Stability of F 0 In this section we derive the stability of F 0 .
As A 1 B 2 is lost if ρ = 0 and (2.3) hold, it is sufficient to consider the dynamics (2.5) in S 3 × S 3 . In this case, the characteristic polynomial at F 0 factors into two quadratic polynomials, P(λ) = t 1 (λ)t 2 (λ). These are given by The polynomial t 1 determines the stability with respect to the (effectively one-locus) system where only 'alleles' A 1 B 1 and A 2 B 2 are present. It is convex with t 1 (0) ≥ 0 if and only if |κ 1 + κ 2 | ≤ 1 (where the equalities correspond), i.e., whenever F 0 is admissible, cf. (3.23). If |κ 1 + κ 2 | < 1, t 1 (0) > 0 and t 1 attains a negative value at its minimum (as can be shown easily). Therefore, all eigenvalues emanating from t 1 are real and negative whenever F 0 is admissible. The polynomial t 2 determines stability with respect to the interior of S 3 × S 3 . It is convex and attains its minimum at where λ min < 0 by (2.3a) and (3.22). As t 2 (λ min ) < 0, the eigenvalues emanating from t 2 are real. As t 2 (0) ≥ 0 ⇐⇒ m 1 m 2 ≤m, (10.14) where the equalities correspond andm is defined in (3.27), and because t 2 (0) > 0, we conclude that the two eigenvalues emanating from t 2 are negative if and only if (3.28) holds.
A.5 Stability of SLPs if ρ = 0 For ρ = 0 it is sufficient to study the dynamics (2.5) in S 3 × S 3 . SLPs wherex k,2 > 0 (k = 1, 2), i.e., P A,2 and P B,1 , are unstable. It remains to study the stability of P A,1 and P B,2 . We present the analysis for P A,1 in detail, as results for P B,2 follow analogously. At P A,1 the characteristic polynomial factors into two quadratic polynomials, P(λ) = t 1 (λ)t 2 (λ), given by t 1 determines stability with respect to the one-locus system where B 1 is fixed. t 1 (0) = 0 if and only if |σ 1 + σ 2 | = 1, i.e., whenever P A,1 collides with a ME according to (3.9a) and (3.9b). Whenever |σ 1 + σ 2 | < 1, i.e., P A,1 is admissible (3.7), t 1 (0) > 0 and t 1 (0) > 0. As t 1 (λ) > 0 for every λ, t 1 attains a minimum, where it is straightforward to show that t 1 takes a negative value at its minimum. Thus, all eigenvalues emanating from t 1 are real and negative whenever P A,1 is admissible. t 2 determines stability with respect to the interior of S 3 × S 3 . t 2 (0) ≥ 0, if and only if m 1 m 2 ≥m, cf. (3.27), where the equalities correspond. Whenever m 1 m 2 >m, t 2 (0) > 0. As t 2 (λ) > 0 for every λ, t 2 attains a minimum, where it is straightforward to show that t 2 takes a negative value at its minimum. Thus, all eigenvalues emanating from t 1 are real and negative whenever m 1 m 2 >m holds. Otherwise, at least one eigenvalue is positive.
A.6 The super-symmetric case We prove that in the super-symmetric case of Sect. 3.10, all SLPs are unstable. We assume symmetric migration rates (m 1 = m 2 = m), equivalent loci (α k = β k = a), and selection in deme 2 mirrors that in deme 1 (α k = −α k * ). Thus, θ = 0. Equilibria may collide (thus leave or enter the state space) if and only if at least one of their eigenvalues is zero. Eigenvalues are zeros of the characteristic polynomial, which has the form P(λ) = c 6 λ 6 + · · · + c 1 λ + c 0 . If zero is an eigenvalue at an equilibrium, i.e., P(0) = 0, the constant term c 0 must vanish. In the super-symmetric case every characteristic polynomials at an SLP has the same constant term c 0 =−a 2 ρ 2a 2 a 2 +m 2 +m(3m −ρ) a 2 +m 2 −(a 2 +m 2 )(3m −ρ) . (10.17) One can show that c 0 = 0 is impossible if m > 0.

A.7 Important quantities and relations
The following section complements Sect. 4.1. Here, we derive all relations of φ X (4.2) and m X (4.3) needed inSects. 4.2 to 4.8 and in the proofs of the theorems there.
Our approach to derive the possible relations between m F 0 and m A or m B is as follows: First, we derive all relevant relations of φ X (4.2) for arbitrary recombination ρ. We use these relations to determine the required relations between m A , m B , m M 1 and m M 4 for arbitrary ρ. By setting ρ = 0 in the results obtained and by the equivalence given in (4.10), the possible relations between m F 0 and m A or m B follow immediately.