A dynamical theory of speciation on holey adaptive landscapes

The metaphor of holey adaptive landscapes provides a pictorial representation of the process of speciation as a consequence of genetic divergence. In this metaphor, biological populations diverge along connected clusters of well-fit genotypes in a multidimensional adaptive landscape and become reproductively isolated species when they come to be on opposite sides of a ``hole'' in the adaptive landscape. No crossing of any adaptive valleys is required. I formulate and study a series of simple models describing the dynamics of speciation on holey adaptive landscapes driven by mutation and random genetic drift. Unlike most previous models that concentrate only on some stages of speciation, the models studied here describe the complete process of speciation from initiation until completion. The evolutionary factors included are selection (reproductive isolation), random genetic drift, mutation, recombination, and migration. In these models, pre- and post-mating reproductive isolation is a consequence of cumulative genetic change. I study possibilities for speciation according to allopatric, parapatric, peripatric and vicariance scenarios. The analytic theory satisfactorily matches results of individual-based simulations reported by Gavrilets et al. (1998). It is demonstrated that rapid speciation including simultaneous emergence of several new species is a plausible outcome of the evolutionary dynamics of subdivided populations. I consider effects of population size, population subdivision, and local adaptation on the dynamics of speciation. I briefly discuss some implications of the dynamics on holey adaptive landscapes for molecular evolution.

Speciation has traditionally been considered to be one of the most important and intriguing processes of evolution. In spite of this consensus and significant advances in both experimental and theoretical studies of evolution, understanding speciation still remains a major challenge (Mayr 1982a;Coyne 1992). The main reason for such a discouraging situation is that direct experimental approaches, which are widely used for solving other problems of evolutionary biology, are not effective for studying speciation because of the time scale involved. Experimental work necessarily concentrates on distinct parts of the process of speciation intensifying and simplifying the factors under study (Rice and Hostert 1993;Templeton 1996). In situations where direct experimental studies are difficult or impossible, mathematical modeling has proved to be indispensable for providing a unifying framework.
Although numerous attempts to model parts of the process of speciation have been made, a quantitative theory of the dynamics of speciation is still missing. Currently, verbal theories of speciation are far more advanced than mathematical foundations. As often is the case with verbal theories (both scientific and otherwise), different deduced (or induced) aspects of speciation are emphasized by different workers resulting in confusion and controversy. The situation is not helped by the absence of general agreement on a species definition (e.g., Claridge et al. 1997).
Here, I attempt to develop some foundations for a general dynamical theory of speciation. One possible approach to this goal would be to begin with a species definition, then to define speciation accordingly and to develop an appropriate dynamical model. I do not think such approach would be very useful, due to a lack of generality. My models are not based on a specific "species concept". I reason that species are "different" with respect to some characteristics, and that whatever these differences, they have a genetic basis. Thus, modeling the dynamics of speciation is equivalent to modeling the dynamics of genetic divergence. I use a bottom-up approach: begin with a model incorporating a range of factors thought to lead to speciation (e.g. selection, mutation, population subdivision etc.) and then try to interpret its dynamic behavior in terms of different "species concepts". As expected, many aspects of speciation that are emphasized by different species concepts (such as reproductive isolation, separate genotypic clusters or common evolutionary trajectories) emerge from the same processes. This clearly indicates that different species concepts are not mutually exclusive.
The choice of a modeling approach depends upon the purpose of the model. A common view in (evolutionary) biology is that mathematical models are mainly useful for making predictions that can be used in experimental work. Although such a pragmatic approach is probably what should be expected in contemporary society, a model's testable predictions are not necessarily its main contribution to science. Insights provided by models, their ability to train our intuition about complex phenomena, to provide a framework for studying such phenomena and to identify key components in complex systems are at least as important as specific predictions. For these purposes the most useful tools are simple models and metaphors.
Sewall Wright's (1932) metaphor of "rugged adaptive landscapes" is a well-known and widely used metaphor in evolutionary biology. In the standard interpretation, a rugged adaptive landscape is a surface in a multidimensional space that represents the mean fitness of the population as a function of gamete (or allele) frequencies which characterize the population state (see Fig.1a). It is envisioned that this surface has many peaks and valleys corresponding to different adaptive and maladaptive population states, respectively. The population is imagined as a point on the surface which is driven by selection up hill but can get stuck on a local peak. Two general points about scientific metaphors should be kept in mind. The first is that specific metaphors (as well as mathematical models) are good for specific purposes only. The second is that accepting a specific metaphor necessarily influences and defines the questions that are considered to be important. The metaphor of "rugged adaptive landscapes" is very useful for thinking about adaptation. However, its utility for understanding speciation is questionable. From a pragmatic point of view, the process of splitting a population into two different species is impossible to describe in a framework where a population is the smallest unit. Finer resolution is necessary for describing the splitting of populations. Accepting the metaphor of rugged adaptive landscapes immediately leads to a problem to be solved: how can a population evolve from one adaptive peak to another across an adaptive valley when selection opposes any changes away from the current adaptive peak? Wright's solution to this problem, his shiftingbalance theory (Wright, 1931(Wright, , 1982, does not seem to be satisfactory (Gavrilets 1996;Coyne et al. 1997). Provine (1986), Barton and Rouhani (1987), Whitlock et al. (1995), Gavrilets (1997a) and Coyne et al. (1997) discuss other weaknesses of Wright's metaphor. I have argued elsewhere (Gavrilets 1997a) that his metaphor of rugged adaptive landscapes with its emphasis on adaptive peaks and valleys is, to a large degree, a reflection of the three-dimensional world we live in. Both genotypes and phenotypes of biological organisms differ in numerous characteristics, and, thus, the dimension of "real" adaptive landscapes is much larger than three. Properties of multidimensional adaptive landscapes are very different from those of low dimension. Consequently, it may be misleading to use three-dimensional analogies implicit in the metaphor of rugged adaptive landscapes in a multidimensional context. I believe that understanding speciation requires a different metaphor.
The metaphor of "holey" adaptive landscapes. An individual organism can be considered as a combination of genes. All possible combinations of genes form a geno- Adaptive landscapes. a. A rugged adaptive landscape. The main emphasis is on "peaks" and "valleys" corresponding to different well-adapted and maladaptive genotypes. b. A holey adaptive landscape. The main emphasis is on clusters of well-fit genotypes that extend throughout the genotype space. All other genotypes are treated as "holes".
type space (which, mathematically, can be represented by a hypercube). In discussing the evolution of populations, it is useful to visualize each individual as a point in this genotype space. Accordingly, a population will be a cloud of points, and different populations (or species) will be represented by different clouds. Selection, mutation, recombination, random drift and other factors change the size, location and structure of these clouds. To construct an adaptive landscape one assigns "fitness" to each genotype (or each pair of genotypes) in genotype space. Different forms of selection and reproductive isolation can be treated within this conceptual framework. For example, fitness can be a genotype's viability (in the case of viability selection); it can be fertility or the probability of successful mating between a pair of genotypes (in the case of fertility selection or premating isolation, respectively). A finite population subject to mutation is likely to be represented by genotypes with fitnesses within a fitness band determined by the balance of mutation, selection and random drift. A general property of adaptive landscapes with a very large number of dimensions is that genotypes with fitnesses within a specified band form connected "clusters" that extend throughout the geno-type space (Gavrilets 1997;Gavrilets and Gravner 1997). Thus, populations can evolve and diverge along bands of highly-fit genotypes without going across the states with a large number of low-fit genotypes (that is without crossing any adaptive valleys).
The metaphor of "holey" adaptive landscapes puts special emphasis on these clusters of highly-fit genotypes, disregarding fitness differences between them and treating all other genotypes as "holes" (Gavrilets 1997;Gavrilets and Gravner 1997). The justification for the latter is a belief that selection will be effective in moving the population away from these areas of genotype space on a time scale that is much faster than the time scale for speciation. Accordingly, microevolution and local adaptation can be viewed as the climbing of the population from a "hole" towards the holey adaptive landscape, whereas macroevolution can be viewed as a movement of the population along the holey landscape with speciation taking place when the diverging populations come to be on opposite sides of a "hole" in the adaptive landscape. In this scenario, there is no need to cross any "adaptive valleys"; reproductive isolation between populations evolves as an inevitable side effect of accumulating different mutations. The metaphor of holey adaptive landscapes can be traced to a verbal two-locus two-allele model of reproductive isolation proposed by Dobzhansky (1937) and similar ideas discussed by Bateson (cited in Orr, 1997), Muller (1942), Maynard Smith (1970Smith ( , 1983, Nei (1976), Barton and Charlesworth (1984), Kondrashov and Mina (1986). For more discussion of this metaphor see Gavrilets (1997ab) and Gavrilets and Gravner (1997). Orr (1995) and Gavrilets (1997a) gives a summary of relevant experimental evidence. The metaphor of "holey" adaptive landscapes is illustrated graphically in Figure 1b.
Mathematical models for holey adaptive landscapes. Here I briefly review previously published work on the evolutionary dynamics on "holey" adaptive landscapes. Nei (1976) and Wills (1977) were the first to present formal analyses of the Dobzhansky model. Nei et al. (1983) studied one-and two-locus multi-allele models with stepwise mutations and considered both postmating and premating reproductive isolation. In their models genotypes were reproductively isolated if they were different by more than 1 or 2 mutational steps. In these situations, speciation was very slow. They conjectured, however, that increasing the number of loci may significantly increase the rate of speciation. Bengtsson and Christiansen (1983) presented a deterministic analysis of mutationselection balance in the Dobzhansky model. Bengtsson (1985), Barton and Bengtsson (1986) and Gavrilets (1997b) analyzed the properties of hybrid zones arising under Dobzhansky-type epistatic selection. Wagner et al. (1994) considered a two-locus, two-allele model of stabilizing selection acting on an epistatic character. For a specific set of parameters, the interaction of epistasis in the trait and stabilizing selection on the trait resulted in a fitness "ridge". The existence of this ridge simplified stochastic transitions between alternative equilibria. Gavrilets and Hastings (1996) formulated a series of two-and three-locus Dobzhansky-type viability selection models as well as models for selection on polygenic characters. They studied these models in the context of founder effect speciation and noticed that the existence of ridges in the adaptive landscape made stochastic divergence much more plausible. Similar conclusions were reached by Gavrilets and Boake (1998) who studied the effects of premating reproductive isolation on the plausibility of founder effect speciation. Derrida (1991, 1992) proposed a model where the probability of mating between two haploid individuals is a decreasing function of the proportion of loci at which they are different. Here, any two sufficiently different genotypes can be considered as sitting on opposite sides of a hole in a holey adaptive landscape. These authors as well as Manzo and Peliti (1994) studied this model numerically assuming that the number of the loci is infinite, the loci are unlinked and highly mutable, and mating is preferential. Orr (1995) and Orr and Orr (1996) studied speciation in a series of models in which viability of a diploid organism depends on the number of heterozygous loci. All these papers postulated the existence of ridges of highlyfit genotypes. Gavrilets and Gravner (1997) studied a general class of multilocus selection models and showed the existence of ridges to be inevitable under fairly general conditions. Independently, a similar conclusion was reached in Reidys et al. (1997). Most previous studies of the dynamics of speciation on holey adaptive landscapes were numerical. To develop a dynamical theory of speciation it is desirable to have a simple model that can be treated analytically.

The Model
I consider finite populations of haploid individuals with discrete, non-overlapping generations. I assume that reproduction involves gene exchange (amphimixis) between individuals. The restriction to haploids is for algebraic simplicity. Models for diploids will be discussed later. Individuals are different with respect to L possibly linked diallelic loci. Without any loss of generality each individual's genotype can be represented as a sequence of 0's and 1's. Let l α = (l α 1 , ..., l α L ) where l α i is equal to 0 or 1, be such a sequence for an individual α. In standard population genetics models, the population state is usually described in terms of gamete frequencies. In systems with many loci such an approach is not practical. For instance, with 10 diallelic loci there are 2 10 different gametes. Thus, one would need to analyze more than 1000 coupled equations.
Another complication follows from the fact that even in very large populations with hundreds of thousands individuals each specific genotype is represented only by a small number of copies or is not represented at all. Thus, the notion of a gamete frequency in multilocus evolution might be very difficult to justify. Here, I will be interested in the levels of genetic variation within subpopulations and genetic divergence between subpopulations. Both can be characterized in terms of genetic distance d defined as the number of loci at which two individuals are different. More formally, the genetic distance d αβ between individuals α and β is Genetic distance d is the standard Hamming distance.
It is analogous to the number of segregating sites in a sample of two gametes, which is widely used in molecular evolutionary genetics (Li, 1997), and to the number of heterozygous loci in a diploid organism. Genetic distance d is also closely related to the notion of the overlap q between two sequences, d = L 2 (1 − q), which is commonly used in statistical physics (e.g., Derrida and Peliti 1991). I model the expected dynamics of average genetic distances within and between populations, using D w for the former and D b for the latter.
I assume that reproductive isolation is caused by cumulative genetic change. I will use a very simple symmetric model which is closely related to the models discussed above and which allows one to treat both pre-and postmating isolation within the same framework. I posit that an encounter of two individuals can result in viable and fecund offspring only if the individuals are different at no more than K loci. Otherwise the individuals do not mate (premating reproductive isolation) or these offspring are inviable or sterile (postmating reproductive isolation). More formally, I assign "fitness" w to each pair of individuals depending on the genetic distance d between them (see Appendix for an outline of more complicated approaches). In this formulation, any two genotypes different at more than K loci can be conceptualized as sitting on opposite sides of a hole in a holey adaptive landscape (cf. Derrida 1991, 1992). At the same time, a population can evolve to any reproductively isolated state by a chain of single locus substitutions. The adaptive landscape corresponding to this model is both "holey" and "correlated". The latter means that the probability that two genotypes are reproductively isolated correlates with the genetic distance between them. In Nei et al. (1983) and Gavrilets and Boake (1998) models individuals separated by more than one mutational step were reproductively isolated which corresponds to K = 1. The neutral case (no reproductive isolation) corresponds to K equal to the number of loci. The mathematical model presented above was interpreted as describing sexual haploid populations with fitnesses assigned to pairs of individuals depending on the genetic distance between them. However, there is an alternative interpretation in that the model describes randomly mating diploid populations. In the diploid case, the genetic distance (1) between the two gametes forming an individual is equivalent to individual's heterozygosity, and fitness function (2) specifies fitness as a function of individual heterozygosity. Therefore, most conclusions of this paper will also be applicable to situations when postmating reproductive isolation is in the form of reduced (or zero) viability of hybrids due to incompatibility of the genes they receive from their parents (Wu an Palopoli 1994).

Dynamics in the neutral case
Before developing a theory for the dynamics of speciation in the above model, it is illuminating to start with the neutral case. Here I summarize some relevant results that are presented in (or follow directly from) classical papers (Watterson 1975;Li 1976;Slatkin 1987;Strobeck 1987). Let µ be the probability of mutation per locus per generation. The approximations below assume that mating is random, the number of loci L is large, but µ is very small so that the probability of mutation per individual per generation v ≡ Lµ << 1. The migration rate m and the inverse of the population size 1/N are small as well.
Genetic variation within an isolated population. Let us consider an isolated population of size N T . The expected change in the average genetic distance within the population per generation is where the first term in the right-hand side is the contribution of mutation whereas the second term is the random drift reduction of D w . Asymptotically, a mutation-drift equilibrium is reached with Genetic divergence between isolated populations. Let us consider several isolated populations of arbitrary size. The probability that a specific mutation gets fixed in a population is one over the population size. Different mutations will get fixed in different populations resulting in their genetic divergence. The average genetic distance between any two of them increases with the rate equal twice the mutation rate per gamete The rate of neutral divergence does not depend on population sizes. In particular, it is the same independent of whether there are many small populations or a few large populations. Because the number of loci L is finite, an indefinite increase of D w , which is implied by equation (5) is impossible. This equation as well as equation (3) above and equations (6) below approximate the dynamics when genetic distances D w and D b are small relative to the number of loci L. To treat the general case, one has to substitute v for v(1 − 2D w /L) in equations (3) and (6a) and for v(1 − 2D b /L) in equations (5) and (6b). With a finite number of loci, the genetic distance D b between isolated populations approaches L/2 asymptotically. Subdivided populations. The effect of migration on the average genetic distances depends on the spatial structure of populations. Assume that a population of size N T is subdivided into n subpopulations of size N = N T /n and that a proportion m > 0 of individuals migrate to any of the other n − 1 subpopulations. The expected changes in the average genetic distances within and between subpopulations are Equations (6) assume that v, m and 1/N are small. Asymptotically, a mutation-migration-drift equilibrium is reached with where θ is given by equation (4). The average genetic distance within a subpopulation (of size N ) does not depend on the number of subpopulations n and migration rate m and is the same as is expected in a single population with size N T . The average genetic distance between sub-populations increases with the population subdivision and decreasing migration. Figure 2 illustrates the dynamics of neutral divergence in a system of two subpopulations. Peripheral population. Assume a "peripheral" population of size N is receiving migrants from a very large "main" population. Genetic variation in the main population is assumed to be constant (and not influenced by migration from the peripheral population). The expected changes in the average genetic distances within the pe- where D 0 is the average genetic distance within the main population and m is the proportion of individuals in the peripheral population replaced by migrants from the main population. Asymptotically, a mutation-migration-drift equilibrium is reached with where the former equation assumes that genetic variation in the main population is sufficiently large (D 0 >> vN ) and the number of migrants, N m, is not too small. The average genetic distance within the peripheral population is always larger than that for an isolated population of its size (D * w > 2vN ). If the number of migrants is large (N m >> 1), the average genetic distance within the peripheral population is about the same as in the "main" population.

Dynamics with reproductive isolation
The main feature of the model for reproductive isolation introduced above and other models of holey adaptive landscapes is the existence of chains of equally-fit combinations of genes separated by single substitutions that extend throughout the genotype space. These chains can be though of as "neutral paths" in the adaptive landscape. It is important to realize, however, that the exis-tence of "holes" in a holey adaptive landscape makes the actual dynamics of genetic divergence not neutral. In this section I summarize some analytical results on the evolutionary dynamics in the case of reproductive isolation described by equation (2).  have studied the possibilities for speciation in this model numerically. Details of the analytical methods used are outlined in the Appendix. To derive the dynamic equations below, I have used the same assumptions as described above at the beginning of the section on the neural case substituting the assumption of random mating for the assumption of random encounters. In addition, I have assumed that the distributions of genetic distances both within and between populations are Poisson. There are several sets of approximations resulting in Poisson distribution of genetic distances. In the present contest, the weakest set seems to be the assumption that genetic variation at each locus is small most of the time (rare-alleles approximation) and that the population is approximately at linkage equilibrium. These assumptions are standard in analyzing the dynamics of multilocus systems under the joint action of selection, mutation, and random drift (e.g., Barton 1986;Barton and Turelli 1987;Bürger et al. 1989;Gavrilets and de Jong 1993). The fit of individualbased simulations with analytic predictions is satisfactory both at qualitative and quantitative levels (see below and . Genetic variation within an isolated population. After the population becomes polymorphic at K loci, new mutations are selected against when rare because individuals carrying them have a reduced probability of producing viable and fecund offspring. Selection experienced by individual loci underlying reproductive isolation is frequencydependent (and is similar to that arising in the case of underdominant selection on a diploid locus). The change in D w per generation in an isolated population of size N is approximately where and Γ(x, y) is an incomplete gamma function (e.g., Gradshteyn and Ryzhik, 1994). The value of D * w at the mutation-drift-selection equilibrium can be found by equating the right-hand side of (10) with zero and solving for D w . Figure 3a illustrates the dependence of D * w on the parameters of the model. This figure indicates that the equilibrium values of D w are close to the corresponding neutral predictions (4) if K is larger than 2 or 3 times θ (where θ = 2N v is the average genetic distance within a finite population in the neutral case). Figure 3b gives the values of the effective selection coefficient s. With moderately large K (that is with K ≥ 10), s is very small. The effective selection coefficient s can also be thought of as the strength of induced selection on each locus underlying reproductive isolation. Figure 3b shows that very strong selection on the whole genotype (implied by the existence of complete reproductive isolation at finite values of K) results in very weak selection at the level of individual loci. The mean fitness of the population, w w , can be defined as the proportion of pairs of individuals that can mate and produce fertile and viable offspring (cf., Nei et al. 1983). For a population with an average genetic distance D w , where Γ(x + 1) is a gamma function (e.g., Gradshteyn and Ryzhik, 1994. For integer x, Γ(x + 1) = x!). Figure  3c shows that in spite of relatively high levels of genetic variation maintained in the population, the genetic load (that is the proportion of reproductively isolated pairs of individuals, 1 − w w ) is very low. This seems to be a general property of holey adaptive landscapes (cf. Wills 1977; Bengtsson and Christiansen 1983). Genetic divergence between isolated populations. Even after the genetic distance within an isolated population has reached an equilibrium level, the population keeps evolving as different mutations get fixed. As a consequence, isolated populations will continuously diverge genetically. The asymptotic rate of divergence of two isolated populations of size N each is where is the rate of divergence relative to the neutral case. Here, S = N s/2, s is defined by equation (11) with D w corresponding to the mutation-selection-drift equilibrium, and erf (x) is the error function (= 2/ √ π x 0 exp(−y 2 )dy). In the neutral case, s = 0, U = 1/N and equation (13a) reduces to equation (5). Figure 4 illustrates the dependence of the relative rate of divergence R on model parameters. In the neutral case, the rate of genetic divergence ∆D b does not depend on the population size (equation 5). In contrast, with reproductive isolation the rate of divergence decreases with increasing population size. After the population becomes polymorphic at K loci, new mutations are selected against when rare. Genetic drift operating in finite populations overcomes the effect of selection and allows genetic divergence. For example with K = 20 and v = .0384, a population of size N = 800 will accumulate about 5 substitutions per 1000 generations. A few thousand generations will be sufficient for D b to exceed K significantly. In contrast, very large randomly mating populations will diverge very slowly. Figure 4 indicates that the rate of substitutions is close to the corresponding neutral predictions if K is larger than 2-3 times θ (θ = 2N v). Note that as in the neutral case considered above, an implicit assumption in equation (13) is that genetic distance D b is small relative to the number of loci L. In the general case, At what moment can the two diverging population be considered as two different species? The answer obviously depends on what one means by a species. Let us say that the two populations are different species if the proportion, w b , of encounters between individuals from different populations that can result in mating and viable and fertile offspring is less than a small number γ. (This definition uses the biological species concept.) During initial stages of divergence, this proportion can be approximated by the right-hand side of equation (12) with D b taking the place of D w . Figure 5 shows the minimum genetic distance between populations required for speciation as a function of K for several values of γ. One can see that a The rate of difergence R relative to the neutral case in an isolated population of size N as a function of K for v = .0384. The circles, squares, diamonds and triangles give estimates from individual-based simulations for N = 100, 200, 400 and 800, respectively (thirty runs for each parameter configuration).
genetic distance between the populations on the order of 2 or 3 times K will be sufficient for the status of separate "biological" species. Note that there is very little effect of the magnitude of γ.
Speciation in a subdivided population. In the deterministic limit (that is with N T → ∞), the genetic variation of a subdivided population can be maintained by migration. This can happen if initially alternative alleles are close to fixation in different subpopulations and selection is sufficiently strong relative to migration (e.g., Karlin and McGregor 1972). Let k be the number of loci at which alternative alleles are close to fixation in different subpopulations. Respectively, L − k will be the number of loci at which the same allele is close to fixation in different subpopulations. In the deterministic limit, k does not change. In the n-island model the dynamics of D w and D b are described by equations where s is defined by equation (11), and the "effective" migration rate if k ≤ K and m e = 0 otherwise. Here, w w is given by equation (12)  k)/Γ(K + 1 − k) is the probability that two randomly chosen individuals from different populations are not reproductively isolated. The effective migration rate m e can be thought of as half the probability of mating between individuals from different subpopulations. With no reproductive isolation (with very large K) or no genetic divergence between subpopulations (with D b ≈ D w , k = 0), the effective migration rate is equal to the actual migration rate (m e = m). Comparing equations (14) with their neutral analogs (8) shows that reproductive isolation results in two effects. First, it directly reduces genetic variation within subpopulations and genetic divergence between subpopulations. These effects are described by the first terms in the right-hand side of equations (14). Also, reproductive isolation reduces the gene flow between populations. Given that D b ≥ D w , then m e ≤ m reflecting the fact that genes brought by migrants have a smaller probability of being incorporated in the resident population. In the deterministic limit, both D w and D b always evolve to finite equilibrium values.
Random genetic drift results in two effects. First, it reduces the genetic variation within sub-populations by the amount D w /N . The dynamic equation for D w becomes Second, genetic drift might change k. The expected change in k can be approximated as where R is given by equation (13b). The first term in the right-hand side of equation (17) can be thought of as the rate at which an allele that is initially rare in both sub-populations becomes close to fixation in one of the subpopulations. This rate was found by Lande (1979) using a diffusion approximation and assuming that migration is weak. The second term is the rate at which the loci with different alleles initially close to fixation in different subpopulations become fixed for the same allele in both of them. To find this term I used Barton and Rouhani's (1987) method. Depending on parameter values and initial conditions there are two different dynamic regimes. In the first regime, both D w and D b evolve to finite values, which are smaller than those in the neutral case (and which are much smaller than the number of loci L). Here, selection (reproductive isolation) reduces genetic divergence both within and between subpopulations. In the second regime, D w stays small (relative to L) whereas D b increases effectively indefinitely (to values order L/2). Here, selection (reproductive isolation) reduces the effective migration rate to zero resulting in speciation. These dynamics can be understood in the following way. Changes in D w and D b induced by selection are expected to happen on a faster time scale than changes in k induced by random genetic drift. Thus, D w and D b values should be close to the equilibrium values predicted by equations (14b,16) when k is treated as a constant. The dynamic behavior depends on whether k reaches a finite equilibrium value or keeps increasing. In the latter case, the effective migration rate m e reduces to zero, the rate of change of k approaches 2vR with D b increasing at the same rate (cf. equation 13b). Figure 6 illustrates the dynamics observed by numerically iterating the model equations. The iterations started with all N individuals identical. During the first 1000 generations there were no restrictions on migration between subpopulations and the whole population evolved as a single randomly mating unit (cf. ). The average genetic distance within the population D w evolves according to equation (10). Starting with generation 1000, restrictions on migration were introduced and the dynamics are described by equations (14b-17) afterwards. After generation 1000, each of these figures has two sets of three curves corresponding to two different values of the parameter(s) under consideration. The curves within each set represent D b , D w and k. With migration rate m = 0.01, all these variables evolve towards finite equilibrium values (see Fig.6a) whereas with a smaller migration rate (m = 0.001), D b and k increase effectively indefinitely signifying that speciation has taken place. Thus, reducing migration makes speciation more plausible. Figure 6b shows that increasing mutation rate (from v = 0.0384 to 5 times this value) has a similar ef-fect. These two figures describe the dynamics expected in a system of two subpopulations. Figure 6c compares the dynamics observed in a population subdivided into 2 and 4 subpopulations. This figure shows that increasing population subdivision makes speciation more plausible. Note that the process of genetic divergence described in Fig.6c results in a simultaneous emergence of 4 species. In the cases where speciation takes place (as signified by continuous increase in the genetic distance between subpopulations), the curves representing D b and k are parallel meaning that asymptotically the genetic divergence is due to fixation of different mutations in different subpopulations. In the cases where speciation does not take place, k is close to zero. Altogether, at the qualitative level the results presented in Fig.6 correspond to both biological intuition and the results of individual-based simulations in . At the quantitative level, there is a very good fit between simulations and analytical predictions for levels of genetic variation maintained in subpopulations and the asymptotic rate of divergence between subpopulations. However, the conditions for speciation as predicted by iterating equations (14b-17) appear to be more strict than those observed in the individual-based simulations performed by . For example, for parameter values used in Figure 6c, no speciation in a system of four subpopulations occurs if m > 0.0035. In contrast, in individual-based simulations speciation was observed for m = 0.01 (Figure 3b in ). The main reason for this disperancy seems to be an inadequacy of equation (17) at moderate levels of migration (e.g. Lande 1979; Barton and Rouhani 1987).

Speciation in a peripheral population.
Here I consider the case of a peripheral population of size N receiving migrants from a very large main population. The dynamics of the average genetic distance within the peripheral population, D w , the average genetic distance between the peripheral and main populations, D b , and the number of diverged loci k are approximated by equations ∆k =vR2 −N me − km e R(e/2) N me .
Here D 0 is the average genetic distance within the main population. Figure 7 illustrates the dynamics observed by numerically iterating equations (18). For D 0 I used mutation-selection balance values for a very large isolated population predicted by equations (10, 11). The initial values of D w and D b were equal to D 0 . The parameter values in Fig.7a and Fig.7b are the same as those that resulted in speciation in Fig.6a and 6b, respectively. The outcome of the dynamics is the same -speciation -but the rate of divergence is smaller than when all subpopulations are uniformly small. This is apparent from the level of genetic distance between subpopulations achieved after 1000 generations of divergence which are about twice as small in Figures 7a and 7b relative to those in Figures 6a and 6b.

Discussion
The theory developed above together with earlier numerical simulations (see references above) show that rapid speciation is a plausible outcome of the evolutionary dynamics in subdivided populations. Here speciation is a consequence of two fundamental factors. The first factor is the existence of various and possibly significantly different highly-fit combinations of genes underlying diverse solutions (genetical, ecological, behavioral, developmental etc.) to the problem of survival and reproduction. In multidimensional genotype space these combinations of genes tend to form connected clusters that extend throughout genotype space. At the same time these genotypes are not mutually compatible -they are separated by "holes". The second factor is mutation pressure. Because the population size is finite and the number of loci is very large whereas the probability of a specific mutation is very small, different mutations tend to appear (and increase in frequency) in different subpopulations (cf., Barton 1989;Mani and Clarke 1990). Metaphorically speaking, mutation tends to tear apart the cloud of points representing the population in genotype space. Combining genes from two different organisms in one offspring can counteract the disruptive effect of mutation, keeping a single randomly mating population together in genotype space. But restricting gene exchange as a consequence of limited migration between subpopulations gives mutation a significant advantage. Eventually the population cloud will be broken and smaller clouds representing the subpopulations will drift apart in genotype spacean event representing speciation. Given sufficient genetic divergence, restoring migration to high levels will not return the system back to the state of free gene exchange between subpopulations which now can be considered as different species. It is not necessary to invoke strong selection for local adaptation to explain speciation in a subdivided population, as studied here, or after a founder event (Gavrilets and Hastings 1996;Gavrilets and Boake 1998 Unlike most previous models that concentrate only on some stages of speciation, the model studied here describes the complete process of speciation from initiation until completion. I assumed that reproductive isolation is caused by cumulative genetic change. The model is described in terms of dynamic equations for the variables analogous to those used in molecular evolutionary biology -the average genetic distances between and within subpopulations. Average genetic distances within (sub)populations always evolve towards finite equilibrium values. Depending on parameter values and initial conditions average genetic distances between subpopulations either converge to a finite equilibrium or increase effectively indefinitely. The former regime is interpreted as no speciation. In the latter regime, three effects take place simultaneously: (1) genetic distances between subpopulations significantly exceed genetic distances within them, (2) encounters between individuals from different subpopulations do not result in viable and fertile offspring, (3) evolutionary changes in a subpopulation do not affect other subpopulations. Thus, subpopulations form separate genotypic clusters in genotype space, become reproductively isolated and undertake changes as evolutionary independent units. This regime is interpreted as speciation according to any of the species concept common in the literature (e.g., Mallet 1995; Claridge et al. 1997).
The dynamic equations derived above describe the expected changes in the average genetic distances neglecting stochastic fluctuations around the expected values. The predicted dynamics have two clearly distinct regimes: convergence towards a finite equilibrium (no speciation) or effectivly indefinite divergence (speciation). Stochastic fluctuations around the expected values, which are present in natural populations (and individual-based simulations), make the boundary between these two regimes less strict and may result in the population escaping the first regime and entering the second regime after some time (see . My analysis has been based on approximations which are standard in studying multilocus systems. I assumed that alleles are rare, that linkage disequilibrium, mutation and migration rates are small and used a theory developed by Lande (1979), Walsh (1982) and Barton and Rouhani (1987) for describing stochastic transitions driven by random genetic drift. The analytic theory presented here fits satisfactorily with the results of individual-based simulations. The model can be used to evaluate qualitative effects of different factors on the dynamics of speciation, the order of magnitude of parameters resulting in or preventing speciation, and the time scale involved. According to both biological intuition and previous numerical simulations, increasing mutation rate and decreasing migration promote speciation. Increasing the number of loci has significantly increased the plausibility of speciation relative to that in earlier models (Nei et al. 1983;Wagner et al. 1995;Gavrilets and Hastings 1996). Note that the actual number of loci influences the dynamics only through the mutation rate per gamete, v, and parameter K. For realistic parameter values the time scale for speciation can be as short as a few thousands or even hundreds of generations. This is compatible with rates observed in several cases of rapid speciation in natural populations described recently (Schluter and McPhail 1992;Yampolsky et al. 1994;Johnson et al. 1996;McCune 1996McCune , 1997 including the most spectacular case -the origin of hundreds of species of Lake Victoria cichlids in 12,000 years (Johnson et al. 1996). The model has demonstrated the plausibility of speciation with relatively low levels of both initial genetic variation and new genetic variation introduced into the population each generation (both supplied by mutation). With higher levels of the former (as in laboratory experiments on speciation, reviewed by Hostert 1993, andTempleton, 1996) or of the latter (for instance as a result of natural hybridization, reviewed by Bullini 1994, Rieseberg 1995, Arnold 1997, the rate of speciation is expected to be even higher.

Local adaptation and speciation
The model analyzed above shows that rapid speciation in a subdivided population can occur even without any differences between selection regimes operating in different subpopulations (that is without selection for local adaptation). An important question is how genetic changes brought about by selection for local adaptation would affect the dynamics of speciation (e.g., del Solar 1966;Ayala et al. 1974;Kilias et al. 1980;Dodd 1989;Schluter 1996;Givnish and Sytsma 1997). These effects will depend on whether the genes responsible for local adaptation are different from or are the same as the genes underlying reproductive isolation.
Assume first that the two sets of genes are completely different. Let the strength of selection per locus induced by reproductive isolation be very small so that these loci can be considered as effectively neutral. (For the model studied here, this seems to be the case if K is larger that 2-3 times θ, where θ is the average genetic distance maintained by mutation in a finite population in the neutral case.) Then, Birky and Walsh's (1988) results tell us that the rate of substitution in these loci will not be affected by selection on other loci independently of linkage. However, given that reproductive isolation is a result of genetic incompatibilities, the loci underlying reproductive isolation will be under frequency-dependent selection against rare alleles which is analogous to underdominant selection in diploid populations. Birky and Walsh (1988) have shown that linkage to advantageous alleles slightly increases the rate of fixation of detrimental mutations. This suggests that selection on linked loci will increase the rate of substitutions in the loci underlying reproductive isolation and, thus, will promote speciation to some degree. No results seem to be known on how linkage to advantageous alleles increases the rate of fixation of underdominant mutations or alleles experiencing frequencydependent selection. No quantitative predictions can be made here, but most likely if the two sets of loci are not extremely tightly linked, effects of selection for local adaptation on the rate of speciation will not be significant.
Assume now that the loci under consideration pleiotropically affect both survival in a given environment and reproductive isolation. For instance, this may be the case if disruptive selection acts on habitat preferences which also define mating patterns (e.g., Rice 1984;Rice and Salt 1988) or if the probability of mating between individuals depends on the difference in their morphological traits which are under direct selection. Let s LA be the average strength of selection per locus induced by selection for local adaptation. Using Walsh's (1982) results, the relative rate of fixation of new mutations in an isolated population of size N can be approximated (see Appendix) as where S = N s/2, s is the strength of selection per locus induced by reproductive isolation, and α = s LA /s. With α = 0, equation (19) reduces to (13b). Figure 8 illustrates the dependence of R on S and α. Increasing α always increases R. Thus, selection for local adaptation always increases the rate of substitutions and promotes speciation. With sufficiently strong selection for local adaptation (s LA > s), the net effect of new alleles will be advantageous and their frequencies will tend to increase even when rare. In the limit of large population size, the probability of fixation is 2(s LA − s). This is analogous to the classical results on the probability of survival of an advantageous mutant in a very large population (Haldane 1927;Walsh 1982). The rate of accumulation of genetic differences will be 2(s LA − s)N v and can be significant. Very strong artificial selection for local adaptation has been shown to result in rapid evolution of reproductive isolation (e.g., del Solar 1966; Kilias et al. 1980;Dodd 1989). However, the changes brought about by moderately strong artificial selection may not exceed those resulting from random genetic drift only (e.g. Ringo et al. 1985).

Population subdivision and speciation
In the models considered here, speciation is a by-product of fixation of different alleles in different subpopulations.
It is well known that the rate of fixation of neutral alleles does not depend on population size, that of advantageous alleles increases with population size, and that of deleterious or underdominant alleles decreases with population size (e.g., Gillespie 1991;Ohta 1992). At the level of individual loci, selection induced by reproductive isolation in the form considered here is similar to underdominant selection (or frequency-dependent selection against rare alleles). Thus, in the absence of selection for local adaptation (or with independent loci controlling traits for local adaptation) decreasing population size will increase the rate of substitutions and promote speciation (see equation 13b). Effects of the population size on the plausibility of speciation will be similar even if the same loci control both reproductive isolation and locally beneficial traits given that selection for local adaptation is not too strong (s LA < s, see equation 19 and Fig.9). In all these cases, speciation will be driven by mutation and random genetic drift and will be fastest if the population is subdivided into small subpopulations. This conclusion about the effect of population subdivision on the probability of speciation in Dobzhansky-type models differs from that of Orr and Orr (1996). They argued that the degree of population subdivision has no effect on the rate of speciation if speciation is caused by mutation and random drift. Orr and Orr did not consider the actual process of fixation of new mutations assuming that it will be a simple neutral process. However, the existence of holes in the adaptive landscape makes the process of substitution non-neutral and new mutations are selected against when rare. Such mutations are fixed more easily in smaller subpopulations. For the discussion of the existing experimental evidence regarding effects of random genetic drift on the plausibility of speciation see Rice and Hostert (1993) and Templeton (1996). The time scale for speciation is short meaning that restrictions on migration between subpopulations do not need to be long lasting; several hundreds of generations may be sufficient for a significant divergence and evolution of reproductive isolation. It is quite possible that several new species will emerge from a highly subdivided population within a short period of time (see above; ). These theoretical conclusions are consistent with a verbal "micro-allopatric" model of speciation suggested for cichlid fishes in the East African Great Lakes (e.g., Reinthal and Meyer, 1997). Hoelzer and Melnick (1994) have emphasized that the possibility of simultaneous emergence of several new species should be incorporated more explicitly in the contemporary methods for reconstructing phylogenies.
If the same loci control both reproductive isolation and locally beneficial traits and selection for local adaptation is sufficiently strong (s LA > s), increasing the population size will result in increasing the rate of substitutions (see Fig.9). In this case, speciation will be driven by selection and will be fastest if the population is subdivided into a small number (say, two) of large subpopulations (Orr and Orr 1996) as implied by the vicariance scenario (e.g. Wiley 1988).
Many species are thought to be represented by a few large populations and many smaller "peripheral" populations. Mayr (1963Mayr ( , 1982b proposed the theory of peripatric speciation arguing that speciation is typically initiated in small peripheral populations and he attributed a special role to genetic drift in this process. Gavrilets (1996) has shown that an invasion of a new adaptive combination of genes is most successful if it is initiated in a peripheral population. The results presented here bear out Mayr's argument(see Fig. 8). Small peripheral populations will rapidly diverge genetically from the "main" large population and speciate. Although differences in selection regimes between peripheral and main populations can accelerate divergence, random genetic drift will be the most important factor. On the other hand, if a peripheral population is large enough and is under a selection regime that is sufficiently different from the one operating in "main" populations, then disruption of gene flow can cause evolutionary divergence, perhaps leading to rapid speciation, in the absence of contributions from random genetic drift (Garcia-Ramos and Kirkpatrick 1997).
Summarizing, large randomly mating populations will diverge genetically and speciate only if there is strong selection for local adaptation (for instance after a change in the environment). In contrast, small populations will diverge and speciate even without differences in selection regimes between them. Possibilities for speciation strongly depend on the geographic structure of the population. Here, analysis was restricted to the island model and the continent-island model. Manzo and Peliti (1994) and  present numerical results for stepping-stone models.

Relationship to other speciation models
Using genetic distance (1) implies the equivalence of loci. A general case of non-equivalent loci can be described by introducing a (L × L) matrix G = {G ij } of weights and defining a generalized distance between individuals α and β as where l α and l β are vectors defining the corresponding genotypes and superscript T means transpose. Considering haploid populations and premating isolation only, the model assumes that individuals can mate only if they are not too different genetically. Here, the degree of reproductive isolation was controlled by cumulative genetic difference. However, using the generalized distance (20) allows one to treat models for reproductive isola-tion controlled by quantitative traits as well as models for sexual selection within the same framework (see Appendix). The close relationship between the models of speciation as a consequence of "quasi-neutral" divergence along ridges in the adaptive landscapes and as a consequence of sexual selection was already recognized by Barton and Charlesworth (1984). A fundamental reason for speciation on a holey adaptive landscape is mutation which tends to break the population into reproductively isolated pieces. Population subdivision and the resulting reduction in gene exchange facilitates this process. Here, migration rates compatible with rapid speciation were small (that is speciation was allopatric or parapatric). An interesting question is whether speciation is possible with much higher migration rates. In other words, is sympatric speciation by mutation and random genetic drift on a holey adaptive landscape possible? Numerical simulations of similar models of sympatric speciation where mutation rates were higher Derrida 1991, 1992) or the time span studied was longer (Wu 1985) or the population size was smaller (Gavrilets and Boake 1998) than here, provide an affirmative answer. Adding disruptive selection due to either abiotic factors (say, different resources) or biotic factors (competition) should create additional pressure on the population cloud which might result in rapid sympatric speciation. Gavrilets and Gravner's (1997) results have suggested that clusters of well-fit genotypes that extend throughout genotype space are plausible. If this is so, biological populations are expected to evolve mainly within these clusters and consist most of the time of well-fit genotypes with fitnesses within some band. The metaphor of "holey" adaptive landscapes neglects the fitness differences between genotypes in the cluster but these differences are supposed to exist and should be apparent on a finer scale. If one applies a finer resolution, the movement along the cluster will be accompanied by slight increases or decreases in fitness. Evolution will proceed by fixation of weakly selected alleles which can be advantageous, deleterious, over-and underdominant, or apparently neutral depending on the specific area of genotype space the population passes through. Smaller populations will pass faster through the areas of genotype space corresponding to fixation of slightly deleterious mutations whereas larger populations will pass faster through the areas corresponding to fixation of (compensatory) slightly advantageous mutations. This pattern of molecular evolution, as predicted from the general properties of multidimensional adaptive landscapes, is similar to the patterns revealed by the methods of experimental molecular biology, which form the empirical basis for the nearly neutral theory of molecular evolution (Ohta 1992). From general considerations, one should not expect complete symmetry of "real" adaptive landscapes which are supposed to have areas varying with respect to the "width" and concentration of ridges of well-fit genotypes. Sexual populations are expected to spend more time in areas of high concentration of well-fit genotypes (Peliti and Bastolla 1994). One of the biological manifestations of this effect will be apparent reduction in the probability of harmful mutations, that is, evolution of genetic canalization (cf., Wagner 1996). The metaphor of holey adaptive landscapes may be useful for thinking about these and other evolutionary problems.

Appendix
Effects of mutation, migration and drift on the dynamics of the average genetic distances within and between subpopulations have been previously studied thoroughly (e.g., Watterson 1975;Li 1976;Slatkin 1987;Strobeck 1987). What is left is to add reproductive isolation (that is selection) to the model. I will use the deterministic framework assuming that the population size N → ∞.
The distribution of D w under rare-alleles and linkage equilibrium approximation. I will use the standard notations A i and a i for alternative alleles at the i-locus (i = 1, . . . , L). Let p i be the frequency of allele A i at the i-th locus, q i = 1−p i , and ψ w,i = 2p i q i . Variable ψ w,i can be thought of as the probability that two randomly chosen individuals (sequences) from the same subpopulation are different at the i-th locus. Let d w,i = (l α i − l β i ) 2 be the genetic distance at the i-th locus between two randomly chosen individuals α and β. Note that d w,i = 1 with probability ψ w,i and d w,i = 0 with probability 1 − ψ w,i . Because d w,i is a Binomial random variable, its generating function is γ dw,i (s) = ψ w,i s + 1 − ψ w,i , which can be approximated as e ψw,i(s−1) if ψ w,i << 1 (rare-alleles approximation). Under approximate linkage equilibrium, the generating function of d w = d w,i is γ dw (s) = Π i e ψw,i(s−1) = e i ψw,i(s−1) = e Dw (s−1) (A1) where D w = i ψ w,i . This shows that random variable d w has approximately Poisson distribution with parameter D w and, thus, Selection within an isolated population. Let w(d) be the expected number of fertile and viable offspring that can be produced by a pair of individuals different in d loci. The average fitness of the population is The dynamics of the general model of fertility selection and premating isolation in a haploid population considered here are identical to that of a symmetric viability selection model for a diploid population with viabilities w(d) depending on the number of heterozygous loci d. Under approximate linkage equilibrium, changes in allele frequencies are described by Wright's equation (Wright 1969). Using the equalities ∂ ln w/∂p i = 2(q i − p i )∂ ln w/∂ψ i and D w = i ψ i , equation (A3) can be rewritten as To describe the dynamics of allele frequencies one needs to know the mean fitness of the population.
Truncation selection. This is a selection scheme analyzed in detail in the main body of the paper. Here Using the Poisson approximation (A1), the mean fitness is where the last equality follows from equation (8.352) in Gradshteyn and Ryzhik (1994), resulting in s given by equation (11). To find equation (10), one starts with (A4a) and proceeds using the fact that ∆ψ i ≈ 2(q i − p i )∆p i and that D w = ψ i . Other selection schemes can be considered in a similar way, and some of them result in relatively compact expressions for w and s.
Linear selection. Here The mean fitness is

Quadratic selection.
Here The mean fitness is where H is the hypergeometric function (Gradshteyn and Ryzhik 1994). Exponential selection. Here The mean fitness is I have not explored how assuming these selection schemes would affect the outcome of the dynamics. Stochastic transitions in an isolated population. Adding mutation results in equation where µ is the rate of mutation (assumed to be equal for forward and backward mutations). Equation (A9) is similar to the classical equation describing underdominant selection on a single locus in a diploid population. This allows one to use Lande's results (1979;see also Hedrick 1981;Walsh 1982;Barton and Rouhani 1987) to find the rate of stochastic divergence. This rate is twice the expected number, vN , of new mutations in a population times the probability that a given one will be fixed, U . Using the diffusion approximation, U is defined by equa-tions (1a) and (2) in (Lande 1979). Lande used some approximations to evaluate U . However, the integrals in his equation (1a) can be found exactly resulting in where S = N s/2 (Walsh 1983). Expanding in a Taylor series under the assumption that 1/N << 1 results in (13b) which is equivalent to Lande's (1979) formula. The difference between Lande's approximate formula and the exact equation (A10) is negligible.
The distribution of D b under rare-alleles and linkage equilibrium approximation. Let us consider two subpopulations. Let p i and P i be the frequencies of allele A i in the first and second subpopulations, respectively. The genetic distance d b,i at the i-th locus between two randomly chosen sequences from two different subpopulations is a Binomial variable taking values 1 and 0 with probabilities ψ b,i = p i Q i + q i P i and 1 − ψ b,i , respectively (q i = 1−p i , Q i = 1−P i ). I will assume that genetic variation within each subpopulation is low so that ψ b,i is close to either 0 or 1. Let The genetic distance between individuals α and β can be represented as d b = k − 1 δ i + 2 δ i where the first sum is over k loci at which ψ b,i ≈ 1 and the second sum is over L − k loci at which ψ b,i ≈ 0. Using the assumption of linkage equilibrium, the generating function of d b becomes where φ i is the expectation of δ i and D b is the expectation of d b . Using the properties of generating functions, the distribution of d b is With fitness function (A5), the probability that two randomly chosen individuals from different subpopulations are not reproductively isolated is Deterministic dynamics in a subdivided population. With no reproductive isolation and with equal forward and backward migration rates and equal population sizes, the change in p i due to migration is ∆ m p i = m(P i − p i ). The corresponding change in D w is ∆ m D w = 2m(D b − D w ). With reproductive isolation and given that D b ≥ D w , individuals migrating from other subpopulations will have reduced probability of mating. Let w w and w b be the expected numbers of fertile and viable offspring that can be produced as a result of within and between subpopulations encounters. For simplicity I will omit the index specifying the locus under consideration. With equal population sizes and migration rates, the change in the allele frequency due to migration becomes where the "effective" migration rate is (compare with the models with migration between populations of unequal size where the effective migration rate is m times the ratio of the population sizes, e.g. Gavrilets 1996). The corresponding change in D w is Changes ∆ m p i can be thought of as changes in allele frequencies brought about by selection between groups of individuals (migrants and residents) whereas the first term in the right-hand side of equation (A9) can be thought of as the change in p brought about by individual selection.
The dynamics of allele frequencies at a specific locus under the joint action of selection, mutation and migration are described by With s =s = const and m e = const and m e < s/6, dynamic system (A14) has two types of stable equilibria: mutation-selection balance equilibria with p ≈ P, pq ≈ µ/s, ψ w ≈ 2µ/s, ψ b ≈ 2µ/s and migration-selection equilibria with p ≈ Q, pq ≈ µ/s + m/s, ψ w ≈ 2µ/s + 2m/s, ψ b ≈ 1 − 2µ/s − 2m/s. I assume that in the deterministic limit, k out of L loci evolve towards migrationselection balance equilibria whereas the remaining L − k loci evolve towards mutation-selection balance equilibria.
In the latter L − k loci, the dynamics of ψ w and ψ b are approximated by equations In the former k loci, the dynamics of ψ w are described as before by (A15a) whereas the dynamics of ψ b are approximated by equation Selection always reduces ψ w whereas mutation always increases it (see equation A15a). Selection and mutation have the same effects on ψ b for the loci evolving toward mutation-selection balance equilibria (see equation A15b). However, for the loci evolving toward migration-selection balance equilibria, selection increases ψ b whereas mutation decreases it (see equation A16).
Summing up over all loci, one finds equations (14) of the main text. Equations (18) are derived in a similar way assuming that the allele frequencies in the main population do not change.
Stochastic transitions in a subdivided population. In a subdivided population, migration tends to reduce genetic differentiation. Given that migration is sufficiently strong relative to selection, the same allele will be close to fixation in both subpopulations. If, by a chance, an alternative allele approaches fixation in one of the subpopulations creating significant differentiation at a given locus, such differentiation will be quickly lost. The number of loci k at which alternative alleles are close to fixation in different subpopulations will be close to zero on average. However, if migration is relatively weak, then the differentiation created by random genetic drift will not be lost quickly and actually can even accumulate. Let us consider a locus at which initially the same allele is close to fixation in both subpopulations (that is both p ≈ 0 and P ≈ 0). Neglecting the changes in P , the deterministic change in p due to selection and migration is approximately Lande (1979; see also Barton and Rouhani 1987) has shown that the rate at which allele A becomes close to fixation in the first subpopulation while its frequency is about zero in the second population is approximately 2 −N me times the rate of fixation in the absence of immigration. Assuming that alleles A are brought about by mutation at rate µ and summing up over L − k loci, one finds the first term in the right-hand side of equation (17). Once alternative alleles are close to fixation in different subpopulations, random drift can remove genetic differentiation. Let us consider a locus at which initially p ≈ 0 but P ≈ 1. Neglecting the changes in P , the deterministic change in p due to selection and immigration is approximately Using Barton and Rouhani's (1987) method one finds that the rate at which allele A becomes close to fixation in both subpopulations is approximately (e/2) N me times the rate of fixation in the absence of immigration. Assuming that alleles A are brought about by migration at rate m e and summing up over k loci, one finds the second term in the right-hand side of equation (17).
Stochastic divergence with local adaptation. Let us assume that the allele under consideration is favorable in a given environment with selective advantage s LA . The change in this allele frequency as defined by the joint action of selection induced by reproductive isolation and selection for local adaptation is ∆ s p = spq(p − q) + s LA pq.
This equation is identical to the one describing meiotic drive in the Appendix of Walsh (1982). Following Walsh, the fixation probability is where S = N s/2 and α = s LA /s. Expanding the numerator in a Taylor series under the assumption that 1/N << 1 and multiplying the results by the expected number of mutants, vN , results in the relative rate of fixation given by equation (19).
Genetic distance (20) and some other models. Genetic distance (1) is recovered by assuming that G is an identity matrix. Assuming that G is a diagonal matrix with non-equal diagonal elements is a simple way to introduce non-equivalence of loci. The case when the probability of mating depends on the difference in a quantitative trait can be treated within the same framework. Let c i be the contribution of the i-th locus to a quantitative trait z. Neglecting microenvironmental effects, the trait values for individuals α and β are x α = c i l α i and z β = c i l β i , respectively. The square of the difference of z α and z β is recovered from equation (20) by assuming that G ij = c i c j . A common way to model sexual selection is to assume that the probability of mating between a male and a female depends on the difference in a female phenotypic trait, z f , and a male phenotypic trait, z m , which are controlled by two different sets of loci (e.g., Lande 1981;Kirkpatrick 1982;Nei at al. 1983;Wu 1985;Turner and Burrows 1995). Let z m = c m i l i and z f = c f i l i where the sums are taken over the corresponding sets of loci. The value (z m − z f ) 2 is recovered from (20) by assuming that matrix G has a block form The diagonal L m × L m and L f × L f zero matrices correspond to the interactions within the set of L m genes controlling the male trait and within the set of L f genes controlling the female trait (L = L f +L m ), and matrix G s describing the interactions between the two sets of genes has elements G s ij = c f i c m j .