Theoretical studies of speciation have been dominated by numerical simulations aiming to demonstrate that speciation in a certain scenario may occur. What is needed now is a shift in focus to identifying more general rules and patterns in the dynamics of speciation. The crucial step in achieving this goal is the development of simple and general dynamical models that can be studied not only numerically but analytically as well. I review some of the existing analytical results on speciation. I first show why the classical theories of speciation by peak shifts across adaptive valleys driven by random genetic drift run into trouble (and into what kind of trouble). Then I describe the Bateson-Dobzhansky-Muller (BDM) model of speciation that does not require overcoming selection. I describe exactly how the probability of speciation, the average waiting time to speciation, and the average duration of speciation depend on the mutation and migration rates, population size, and selection for local adaptation. The BDM model postulates a rather specific genetic architecture of reproductive isolation. I then show exactly why the genetic architecture required by the BDM model should be common in general. Next I consider the multilocus generalizations of the BDM model again concentrating on the qualitative characteristics of speciation such as the average waiting time to speciation and the average duration of speciation. Finally, I consider two models of sympatric speciation in which the conditions for sympatric speciation were found analytically. A number of important conclusions have emerged from analytical studies. Unless the population size is small and the adaptive valley is shallow, the waiting time to a stochastic transition between the adaptive peaks is extremely long. However, if transition does happen, it is very quick. Speciation can occur by mutation and random drift alone with no contribution from selection as different populations accumulate incompatible genes. The importance of mutations and drift in speciation is augmented by the general structure of adaptive landscapes. Speciation can be understood as the divergence along nearly neutral networks and holey adaptive landscapes (driven by mutation, drift, and selection for adaptation to a local biotic and/or abiotic environment) accompanied by the accumulation of reproductive isolation as a by-product. The waiting time to speciation driven by mutation and drift is typically very long. Selection for local adaptation (either acting directly on the loci underlying reproductive isolation via their pleiotropic effects or acting indirectly via establishing a genetic barrier to gene flow) can significantly decrease the waiting time to speciation. In the parapatric case the average actual duration of speciation is much shorter than the average waiting time to speciation. Speciation is expected to be triggered by changes in the environment. Once genetic changes underlying speciation start, they go to completion very rapidly. Sympatric speciation is possible if disruptive selection and/or assortativeness in mating are strong enough. Sympatric speciation is promoted if costs of being choosy are small (or absent) and if linkage between the loci experiencing disruptive selection and those controlling assortative mating is strong.
Speciation, that is, the origin of species, is one of the most intriguing evolutionary processes. Speciation is the process directly responsible for the diversity of life. Understanding speciation still remains a major challenge faced by evolutionary biology even now after almost 150 years since the publication of Darwin's book, On the Origin of Species. Although the argument on the appropriate way(s) to define a “species” is still unsettled (Wilson 1999), speciation is ultimately a consequence of genetic divergence. Here, I will view the dynamics of speciation as the dynamics of genetic divergence between different populations (or between parts of the same population) resulting in substantial reproductive isolation. A suite of tools for modeling the dynamics of genetic divergence is provided by theoretical population genetics and ecology. These two theories form a foundation for the quantitative/mathematical study of speciation.
Theoretical population genetics has identified a number of factors controlling evolutionary dynamics such as mutation, random genetic drift, recombination, natural and sexual selection, etc. A straight-forward approach for classifying different mechanisms and modes of speciation is according to the type and strength of the factors controlling or driving genetic divergence. In principle, any of the factors listed above can be used at any level of classification. However, traditionally in evolutionary biology the discussions of speciation mechanisms are framed in terms of a classification in which the primary division is according to the level of migration between the diverging (sub)populations (Mayr 1942). In this classification the three basic (geographic) modes of speciation are allopatric, parapatric, and sympatric.
Allopatric speciation (from Greek words allos meaning “other” and patra meaning “fatherland” or “country”) is the origin of new species from geographically isolated (sub)populations. In the allopatric case there is no migration of individuals (and gene flow) between the diverging (sub)populations. At the other extreme of the highest possible migration level lies sympatric speciation (from Greek word sym meaning “the same”). Sympatric speciation is usually defined as the origin of new species from a single local population (Mayr 1942). This definition (as well as similar definitions based on the “absence of geographic separation,” occurrence within a “cruising range,” occurrence within a “geographic range” of ancestor, etc., although intuitively appealing, is not precise enough for modeling purposes. Indeed, to use such a definition one would need to additionally define the exact meaning of “local population,” “geographic separation” or “cruising range,” etc. I will define sympatric speciation as the emergence of new species from a population where mating is random with respect to the place of birth of the mating partners (for a similar approach, see Kondrashov and Mina 1986). This definition is actually implied in most mathematical models of sympatric speciation. The intermediate cases when migration between diverging (sub)populations is neither zero nor maximum possibly fall within the domain of parapatric speciation (from Greek word para meaning “beside”, “side-by-side,” “next to”).
To illustrate this classification of speciation modes, let us consider a system of two demes that exchange a proportion m (≤1/2) of its members each generation. Assume that mating follows migration and is random within each deme. Then the allopatric case corresponds to m = 0, the sympatric case corresponds to m=1/2, and parapatric case corresponds to 0<m<1/2 (see Figure 1). Both this figure and biological intuition suggest that parapatric speciation is the most general (geographic) mode of speciation. Endler (1977) made the same point on the basis of empirical data. I conclude this discussion with a few comments some of which were made repeatedly in the past (e.g., Smith 1955, 1969; Endler 1977). The first is that a consistent use of terminology is very important. Sometimes the notion of sympatric speciation is used as an antonym of allopatric speciation, and thus includes both parapatric speciation and sympatric speciation as defined above. This lays grounds for confusion. Even more misleading is the usage of the term “sympatric speciation” as a synonym of “the origin of species in a well defined geographic area.” For example, both in the empirical and theoretical literature one can easily find statements about “sympatric speciation of Lake Victoria cichlids” which are only a little less misleading than something like “sympatric speciation of North American birds.” One should also be absolutely clear that the definitions of the geographic modes of speciation imply nothing about the forces that drive genetic divergence leading to speciation. Any given evolutionary factor can play a role within each of the three geographic modes. For example, sympatric speciation can be driven by mutation and random drift and allopatric speciation can be driven by ecological factors. Models describing both these examples exist. Finally, the traditional stress on the spatial structure of (sub)populations as the primary factor of classification rather than, say, on selection reflects both the fact that it is most easily observed (relative to the difficulties in inferring the type and/or strength of selection acting in natural populations) and the growing realization that spatial structure of populations is very important. Therefore, I do not find compelling recent suggestions to abandon this classification in favor of a classification based on types of selection (Via 2001) or make it subordinate to a “geography/prezygotic isolating mechanisms” continuum (Kirkpatrick and Ravigné 2002).
Speciation is a very complex process which is affected by many different factors (genetical, ecological, developmental, environmental, etc.) interacting in nonlinear ways. Both this complexity and the difficulties of experimental approaches coming in particular from the very long time scales that are typically involved imply that mathematical models have to play a very important role in speciation research. The ability of models to provide insights into the speciation process, to train our intuition, to provide a general framework for studying speciation, and to identify key components in its dynamics are invaluable.
For reasons that are unclear, models of speciation were slow to appear. Paradoxically, none of the four greats who are usually viewed as the founders of modern theoretical population genetics—Fisher, Wright, Haldane, Kimura—expressed much interest in developing mathematical models explicitly dealing with speciation. (Among their work most closely related to speciation are Fisher's verbal models of runaway evolution caused by sexual selection [1930], Wright's verbal theory of shifting balance [1931, 1982] and his model of assortative mating [1921], studies of stochastic peak shifts [Haldane 1931; Kimura 1985; Wright 1941], and clines [Fisher 1950; Haldane 1948]). To my knowledge, the first substantial discussions of speciation utilizing formal mathematical models are those by Maynard Smith (1962, 1966) and Bazykin (1965, 1969). Although these papers have been very important for setting the field and one of them has become a “citation classic,” as far as concrete results are concerned, they were mostly about modeling the maintenance of genetic variation rather than about speciation per se. A part of the 1966 paper by Maynard Smith did treat speciation directly, but only in a somewhat awkward numerical way for a single set of parameter values assuming complete reproductive isolation from the start. Systematic studies of speciation utilizing mathematical models started only in the early 1970s with the pioneering papers of Crosby (1970), Dickinson and Antonovics (1973), and Balkau and Feldman (1973). Crosby (1973) was the first to use an individual-based model for describing parapatric speciation. Dickinson and Antonovics (1973) were the first to use numerical iteration of dynamic equations to explore a wide range of parameter values for a specific model of parapatric and sympatric speciation. Balkau and Feldman (1973) were the first to find conditions for parapatric speciation in a certain model analytically. For comparison, according to Rice and Hostert's (1993) review (the title of which I have adapted for this paper), the earliest experimental paper on speciation is that of Koopman published in 1950, and by the early 1970s nineteen more experimental papers had been published.
However now the situation has dramatically changed and the theoretical studies well outnumber the experimental papers. This imbalance only continues to grow. For example, recent reviews of experiments on speciation (Florin and Ödeen 2002; Ödeen and Florin 2000, 2002) list only four new experimental papers published since Rice and Hostert's 1993 review, while during the same time interval dozens of new modeling papers came out. At present there have been at least a hundred papers modeling some aspects of speciation (Kirkpatrick and Ravigné 2002). Unfortunately, theoretical speciation research provides an example where the second law of dialectics has failed, at least so far. (The second law of dialectics is concerned with the transformation of quantity into quality [e.g., Engels 1940, Hegel 1975].) Instead of a sound mathematical theory of speciation, there is what I can only interpret as a (justifiable) frustration among both empiricists and theoreticians with the theoretical speciation research. For example, Via (2001), an empiricist, who is very favorable to the idea of sympatric speciation and appears to be convinced it is very probable from a theoretical point of view still asks theoreticians (who have already published dozens of papers on the subject) to “generalize, by further theoretical exploration and integration of existing models, the theoretical conditions under which sympatric speciation can occur” (p. 389). Theoreticians Kirkpatrick and Ravigné (2002) talk about the “balkanization” of the theory of speciation and the absence of clear and general results worthy to be included in the textbooks. Kirkpatrick and Ravigné devised a scheme for classifying the published models of speciation in a hope to better understand their results. Turelli et al. (2001) paint a very gloomy picture of the field of theoretical research concluding in particular that because of the complexity of the process of speciation …“progress on major issues … is more likely to emerge from empirical than from mathematical analyses” (p. 330).
I believe the current situation stems primarily from the limitations of the methods used in theoretical speciation research. Most modeling papers on speciation use numerical simulations as the basic tool. The most general feature of simulation models is that their results are very specific. Interpretation of numerical simulations, the interpolation of their results for other parameter values, and making generalizations based on simulations are notoriously difficult. For these reasons the results of simulations are typically wide open to interpretation and narrowing of interpretation is not possible without additional work. However, simulation results are usually impossible to reproduce because many technical details are not described in original publications. At the same time, the authors of the original publications (and this is not limited to theoretical research) almost never return to their work to check the generality of their conclusions under more general or reasonable parameter values or assumptions.
What is missing in the theoretical speciation research are general and transparent analytical results comparable to those in other areas of theoretical population genetics and ecology. It is simple mathematical models allowing for analytical investigation (rather than complex numerical models) that form the basis of most scientific theories, and there are no reasons why evolutionary biology, in general, and speciation research, in particular, should be an exception. The questions that can be answered depend on the tools used in theoretical research. The majority of existing models of speciation (especially models of sympatric speciation) mostly demonstrate that speciation in a certain scenario may occur. What is needed now is a shift in focus to identifying more general rules and patterns. In principle, to learn that “selection promotes speciation” or that “peak shifts by random genetic drift are unlikely” one does not need mathematical models because this is what biological intuition already tells us. Mathematical models are needed only if one wants to get much more specific conclusions and predictions about the effects of genetic architecture and different types and intensities of selection on the probability of speciation or peak shift and their dynamical characteristics. Such knowledge can be coupled with the estimates of relevant parameters from natural populations to get a much deeper understanding of the process. Mathematical models are also needed if intuition is not working because of the complexity of the process under consideration.
My goal here is to describe and illustrate some existing theoretical results on speciation that are both simple and general. The common wisdom is that a picture is worth a thousand words. In the exact sciences, an equation is worth a thousand pictures. Equations and their interpretations are the most concrete results a theoretician can come up with. Therefore, equations and their interpretations are the focus of this paper. The extensive body of purely numerical work on speciation is deemphasized. The equations to be discussed below are easily accessible to biologists who lack mathematical training. They represent a kind of a “do-it-yourself” kit which can be used by biologists to check or train their intuition about speciation by substituting specific numerical values of parameters and interpreting the results.
In the remainder of this paper, I first will show why the classical theories of speciation by peak shifts across adaptive valleys driven by random genetic drift run into troubles (and into what kind of troubles). Then I will describe the Bateson-Dobzhansky-Muller (BDM) model of speciation that does not require overcoming selection. I describe exactly how the probability of speciation, the average waiting time to speciation, and the average duration of speciation depend on the mutation and migration rates, population size, and selection for local adaptation. The BDM model postulates a rather specific genetic architecture of reproductive isolation. I then show exactly why the genetic architecture required by the BDM model should be common in general. The theoretical results to be described here have led to a major shift in focus from Wright's “rugged adaptive landscapes” to nearly neutral networks and holey adaptive landscapes. Next I consider the multilocus generalizations of the BDM model again concentrating on the qualitative characteristics of speciation such as the average waiting time to speciation and the average duration of speciation. Finally, I consider two models of sympatric speciation in which the conditions for sympatric speciation were found analytically. Results of this type, if known more widely, can clear up a lot of controversy that currently surrounds the questions about the plausibility and generality of sympatric speciation. The Appendix contains a glossary of symbols used in the paper.
Dynamics of Stochastic Peak Shifts
Here I describe two simple models of stochastic transitions between adaptive peaks driven by random genetic drift. My main goal is to show in quantitative terms that strong reproductive isolation is a very unlikely outcome of a stochastic peak shift.
One-locus two-allele model with underdominance
Following Lande (1979), let us consider a randomly mating diploid population with discrete and nonoverlapping generations. Assume that there is a single diallelic locus with alleles A and a controlling fitness (viability). Let the relative fitnesses of genotypes AA and Aa be wAA = 1, wAa = 1 – s and waa, respectively (0 ≤ s ≤ 1) . The adaptive landscape corresponding to this model has two “adaptive peaks” at genotypes AA and aa, and an “adaptive valley” at genotype Aa (see Fig. 2a).
Assume that initially the population is monomorphic for allele A. Let us allow for mutation from A to a with a small probability µ per gene per generation. Selection will tend to eliminate alleles a, whereas mutation will continuously supply them. If the population size N is very large, the population will stay at a state where allele a is maintained at a very low frequency µ/s. However, if the population size N is small, the frequency of the mutant allele a will eventually cross 1/2 as a result of random genetic drift and approach 1 resulting in a “peak shift.” Parameter s, measuring the depth of the adaptive valley that separates the peaks, also characterizes the strength of (postmating) reproductive isolation between the population states before and after the peak shift.
How long does one have to wait until a peak shift occurs in this model? If the mutation rate is very small and the product of N and s is at least moderately large (>2) , the average waiting time until a peak shift occurs can be approximated as
(Lande 1979). This equation shows that T grows exponentially with the product Ns. (Note that although the term 1/
does decrease with increasing Ns, it is dominated by the exponential term). For example, if µ = 10–5 (which is a common estimate of the mutation rate, e.g., Futuyma 1998, Griffiths et al. 1996), s=0.05 (which implies a shallow adaptive valley) and N=400 (which is a relatively small population size), then T = 1013 generations, which is a very long time, in fact, far longer than the probable age of any extant vertebrate species.
Most of this time will be spent waiting for a “lucky” mutant allele destined to be fixed to appear in the population as the overwhelming majority of mutant alleles are removed by selection and random genetic drift. One can also estimate the average time τ that it will take for the “lucky” mutant allele to take over the population. Time τ characterizes the average actual duration of stochastic transitions between the peaks. Using the diffusion approximation (e.g., Ewens 1979) one finds that this time is approximately
The term ln(
) changes very slowly with the product Ns. Therefore, the actual duration of transition is on the order of 1/s generations and, thus, is much shorter than the average waiting time to transition T. For example, with the same values of parameters as above, τ = 66 generations. The ratio of τ and T can be viewed as a measure of the likelihood to observe the actual transition. The above results show that observing the transition is very unlikely.
Additive quantitative character
Next, let us consider a diploid population of size N where individuals differ with respect to a single additive quantitative trait z. Assume that the genotypic variance of the trait is somehow maintained at a constant value G. Let the population be under viability selection defined by a fitness function
where s>0. This fitness function has two local peaks at z = –1 and z = 1 separated by a valley at z=0 (see Fig. 2b). The height of the peaks is one, and the fitness at the bottom of the valley is exp(–s), which is approximately 1 – s for small s. That is, parameter s, as before, characterizes the depth of the adaptive valley separating the adaptive peaks.
Assume that initially the average trait value is in a neighborhood of one of the peaks. Random genetic drift will eventually result in the population crossing the adaptive valley and approaching the other peak (Barton and Charlesworth 1984; Lande 1985). In this model the average waiting time to the peak shift is approximately
Thus, as in the previous model, the waiting time to transition grows exponentially with Ns. Most of time T will be spent in a neighborhood of the initial adaptive peak waiting for the “right” sequence of stochastic events for the peak shift to occur. The average actual duration of transition τ is approximately that is, it has order 1/(Gs). For example, if the genetic variance G = 0.04 (which implies that the bottom of the valley is at five standard deviations from the peak), s=0.05 and N=400, then T = 1.31 × 1020 generations and τ = 490 generations.
There are two important conclusions emerging from these simple analyses. First, unless the population size is small and the adaptive valley is shallow, the waiting time to a stochastic transition between the adaptive peaks is extremely long. This implies that it is very unlikely that a single peak shift will result in strong reproductive isolation. Second, if transition does happen, it is very quick. This implies that observing it “in action” is practically impossible. These conclusions are very general and apply to a broad class of stochastic transitions between adaptive peaks (e.g., Barton and Charlesworth 1984; Barton and Rouhani 1987). Strong reproductive isolation may follow after a chain of stochastic peak shifts each of which is across a shallow valley and, thus, results in a very small degree of reproductive isolation (Walsh 1982). If there are many adaptive peaks separated by shallow valleys, then the pattern of the evolutionary dynamics will be that of a sequence of rapid transitions between adaptive peaks separated by very long periods of no apparent change spent at each of the “intermediate” adaptive peaks. The extended stays at “intermediate” peaks implies that in this model of allopatric divergence the pattern of evolutionary dynamics is quite different from that of “punctuated equilibrium” (Eldredge 1971; Eldredge and Gould 1972; Gould 2002) in which intermediate stages are (almost) never observed. Below we will see that the pattern of “punctuated equilibrium” arises in models of parapatric speciation.
The Bateson-Dobzhansky-Muller Model
By the Bateson-Dobzhansky-Muller (BDM) model, I mean a very specific assumption about the genetic architecture of reproductive isolation, namely that there are two loci with two alleles at different loci being “incompatible” (Bateson 1909; Dobzhansky 1937; Muller 1942). The adaptive landscapes implied by this model are illustrated in Figure 3. In this Figure it is alleles a and B that are assumed to be incompatible in the sense that individuals carrying both of them have zero viability (Fig. 3a) or that females carrying allele a do not mate with males carrying allele B (Fig. 3b). Speciation occurs if two populations that have initially the same genetic composition end up at the opposite sides of the ridge of high fitness values as illustrated in Figure 3. In the BDM model, the populations are not required to cross any adaptive valleys as they simply follow a ridge of high fitness values. For example, hybrid inviability in the fish Xiphophorus behaves as if controlled by a simple BDM model (Orr and Presgraves 2000).
Assume that initially the population has genotype AABB fixed. Let us allow for mutation from A to a and from B to b at an equal rate µ per generation, and let us neglect the possibility of backward mutations. In this model the population will definitely reach the state with genotype aabb fixed (see Fig. 3b). This means that speciation is certain. How long does one have to wait until speciation and what is the actual duration of the speciation process? Can speciation happen by mutation and random genetic drift alone? What are the quantitative effects of selection for local adaptation on the dynamics of speciation? These are the questions I intend to answer in this section.
Let us define the average waiting time to speciation T as the average time that it takes to get from the ancestral state (i.e., with genotype AABB fixed) to the state of reproductive isolation (i.e., with genotype aabb fixed). Let us define the average actual duration of speciation τ as the average time that it takes to get from the ancestral state to the state of reproductive isolation without returning to the ancestral state. Time τ can also be thought of as the average duration of the intermediate stages in the actual transition to a state of reproductive isolation (speciation). In a sense, the ratio of τ and T characterizes the probability of observing incipient speciation. The dynamics of speciation can be modeled as a random walk performed by the most common genotype in the population along the ridge of high fitness values (Gavrilets 2000a). I will describe the corresponding approximations for T and τ separately for allopatric and parapatric cases.
Allopatric speciation
First assume that the two loci under consideration have no other pleiotropic effects on fitness not related to reproductive isolation. In this case, speciation will be driven by mutation and random genetic drift which represents a general null model of speciation. The average waiting time to speciation and the average duration of speciation are equal to
This last expression is very easy to understand. The average waiting time until fixation of a neutral allele is approximately the reciprocal of the mutation rate (Nei 1976). To evolve to a state with haplotype ab fixed, the population has to fix two neutral alleles—the fact reflected in the expression for T above. Note that because in the neutral case the rate of substitutions does not depend on the population size, the population size N does not enter the above equations for T and τ. Because in this model genetic divergence is irreversible, the average duration of speciation is equal to the average waiting time to speciation.
Next we allow for selection for local adaptation. I will assume that mutant alleles a and b have a small selective advantage sa over the ancestral alleles A and B in the environment experienced by the population. In this case the average waiting time to speciation and the average duration of speciation are approximately
where Sa = 4Nsa (and it is assumed that Sa is at least moderately large, i.e., >3) . The above expression strongly supports previous arguments (e.g., Schluter 2000) that selection for local adaptation can dramatically accelerate speciation. For example, increasing Sa by one order of magnitude will decrease the waiting time to speciation by one order of magnitude.
Parapatric speciation
In the parapatric case, we assume that each generation a small proportion m of the population is substituted by immigrants coming from a source population. All immigrants have ancestral genotype AABB.
First assume that speciation is driven by mutation and random genetic drift (the null model). Then the waiting time to speciation is
while the average duration of speciation is where the approximate equalities assume that the rate of mutation is much smaller than the rate of migration (µ≪m). For example, if m = 0.01 and µ = 10–5, then T≈108 generations and τ≈105 generations. Note that the mutation rate has more influence on the dynamics of parapatric speciation than the migration rate.
Next we consider two models of selection for local adaptation. The first model describes direct selection on the alleles underlying reproductive isolation. Specifically, let us assume that mutant alleles a and b have a small selective advantage sa over the ancestral alleles A and B. Then the average waiting time to speciation and the average duration of speciation are
where it is assumed that the coefficient Sa = 4Nsa is at least moderately large (>3) . For example, with the same parameters as above and with Sa=5, then T =1.73 × 105 generations and τ = 2.24 × 104 generations. That is, the average waiting time to speciation is reduced by three orders of magnitude relative to that in the model with no selection, and this brings it near a realistic range of values.
The second model considers effects of a genetic barrier to neutral gene flow. We assume that new alleles a and b have no effects on adaptation to local conditions but that the population has already diverged from the source population in some other loci controlling adaptation to local conditions. Now immigrants will have a reduced fitness, which will affect the fate of the ancestral alleles A and B that they carry. The effects of selection “induced” on neutral alleles via their association with some locally deleterious alleles can be characterized in terms of the gene flow factor (Bengtsson 1985; Barton and Bengtsson 1986; Gavrilets 1997a; Gavrilets and Cruzan 1998) defined as the probability that a neutral gene brought by immigrants makes it to the local genetic background. For example, assume that immigrating adults differ from the residents in two genes: a gene reducing the viability of F1 hybrids to 1–s (relative to viability 1 of the residents) and a neutral gene. Assume also the average number of offspring produced by matings between the immigrants and residents, α, and by matings between F1 hybrids and residents, β, is reduced relative to that of matings among residents (which is normalized to be 1) . This can happen if there is fertility and/or sexual selection against immigrants and/or hybrids. Then the gene flow factor is
where r is the rate of recombination between the selected and neutral loci (Gavrilets and Cruzan 1998). For example, let s = α = β = = 0.5. Then if r = 0.5, then γ = 0.07; if r = 0.05, then γ = 0.016. That is, joint action of viability selection and fertility selection and assortative mating can significantly decrease the effective rate of immigration of neutral alleles even if they are unlinked to the locus under selection. Returning back to modeling speciation, let us assume that the gene flow factor resulting from selection for local adaptation is γ. Then approximately These equations show that a strong genetic barrier to the neutral gene flow (i.e., low γ) will significantly decrease the waiting time to speciation but will not affect its average duration. For example, if γ = 0.1, the waiting time to speciation reduces to one tenth of that in the case of divergence driven by mutation and drift.
There are several conclusions emerging from this analysis. Unless there is selection for local adaptation, the waiting time to speciation in the BDM model is very long (but not as long as under stochastic peak shift). Migration can significantly delay speciation. It is generally agreed upon that spatial heterogeneity in selection is common and that this can have profound effects of the possibility of speciation even in the presence of substantial gene flow (Endler 1977; Ogden and Thorpe 2002; Schneider et al. 1999; Smith et al. 1997). Our model makes this intuition more precise. Selection for local adaptation (either acting directly on the loci underlying reproductive isolation via their pleiotropic effects or acting indirectly via establishing a genetic barrier to gene flow) can significantly decrease the waiting time to speciation. Direct selection is much more effective than indirect selection. In the parapatric case the average actual duration of speciation is much shorter than the average waiting time to speciation.
Nearly Neutral Networks and Holey Adaptive Landscapes
We have seen that evolution of strong reproductive isolation via stochastic transitions between isolated adaptive peaks is unlikely. At the same time, the BDM model shows that strong reproductive isolation can be achieved if the population evolves along a ridge of high fitness values in the genotype space. How likely are such ridges in general? My goal in this section is to answer this question quantitatively. I will follow the general approach of Gavrilets and Gravner (1997) and Gavrilets (1997b, 2003).
To get some intuition about this question let us first consider the set of all possible haploid genotypes that have only two genes each with a very large number of alleles. The alleles at the first locus will be denoted as Ai and the alleles at the second locus will be denoted as Bj. Assume a step-wise mutation pattern (Nei et al. 1983) allowing allele Ai mutate only to alleles Ai + 1 and Ai – 1 and in a similar way, allowing allele Bi mutate only to alleles Bi + 1 and Bi – 1. The genotype space (i.e., the set of all possible genotypes) in this model can be represented by a two-dimensional lattice of square sites where the x- and y- coordinates specify the alleles at the first and second locus, respectively.
Following Gavrilets and Gravner (1997), assume that genotype fitnesses (viabilities) are generated randomly and independently and are only equal to one (viable genotype) or zero (inviable genotype) with probabilities 𝒫 and 1 – 𝒫, respectively. Here, one might think of the set of all possible genotypes playing one round of Russian roulette with 𝒫 being the probability to get a blank. In this model, any change in genotype no matter how large or small, results in a fitness value completely independent of the initial value. Note that this assumption is, of course, a great oversimplification which however does not affect our conclusions.
Let us paint viable genotypes (sites) in black and inviable genotypes (sites) in white (see Fig. 4). For each site (genotype), its four adjacent sites (directly above, below, on the left, and on the right) represent other genotypes that can be obtained by a single mutation. In this model, viable genotypes tend to form connected networks (or ridges) that populations can evolve along by fixing single mutations without the need to go across any inviable states. The number and the structure of networks depend on the probability 𝒫. For small values of 𝒫 there are many networks of small size (see Fig. 4a). As 𝒫 increases, the size of the largest network increases (see Fig. 4a,b,c). As 𝒫 exceeds a certain threshold 𝒫c, known as the “percolation threshold,” there emerges the largest network (known as the “giant component”) that extends (“percolates”) through the whole system and includes a significant proportion of all viable genotypes (see Fig. 4d). In this model, describing a so-called “site percolation” on an infinite two-dimensional lattice, the percolation threshold is 𝒫c = 0.593 (e.g., Grimmett 1989).
The percolation threshold drops dramatically with more loci and alleles. For example, if there are ℒ loci each with 𝒜 alleles, then the percolation threshold is
For example, if ℒ = 10,000 and 𝒜 = 10, then values of 𝒫 larger than 10–5 will typically result in the existence of an extensive network of high-fitness ridges expanding throughout the genotype space. Here, the genotypes that belong to this network have the same fitness. In this sense, the network is neutral. In more general situations where a range of variation of fitness values is allowed, the above results can be generalized to demonstrate the existence of connected nearly neutral networks in which genotypes have approximately the same fitness (Gavrilets 1997b, 2003; Gavrilets and Gravner 1997).
Starting with Wright (1932) typical adaptive landscapes are usually imagined as very rough surfaces with many different peaks and valleys (see Fig. 5a). Continuous evolution on such “rugged adaptive landscapes” requires crossing adaptive valleys (as modeled above). The results on neutral and nearly neutral networks show that Wright's picture can be very misleading if adaptive landscapes have a very high dimensionality because such landscapes are characterized by the existence of extensive nearly neutral networks. Among different nearly neutral networks, those with sufficiently high fitnesses are of particular evolutionary importance because they allow for continuous evolutionary innovations without any significant loss in fitness. An important notion describing such networks is that of “holey adaptive landscape.” A holey adaptive landscape is defined as an adaptive landscape where relatively infrequent high-fitness genotypes form a contiguous set that expands throughout the genotype space. An appropriate three-dimensional image of such an adaptive landscape that focuses exclusively on the percolating network is a nearly flat surface with many holes representing genotypes that do not belong to the network (see Fig. 5b). The smoothness of the surface in this figure reflects close similarity between the fitnesses of the genotypes forming the corresponding nearly-neutral network. The “holes” include both lower fitness genotypes (“valleys” and “slopes”) and very high fitness genotypes (the “tips” of the adaptive peaks). The BDM model considered above provides one of the simplest examples of holey adaptive landscapes. Many more examples are known (Gavrilets 1997b, 2003; Gavrilets and Gravner 1997).
This section has dealt exclusively with the structure of adaptive landscapes. A number of mechanisms, such as random drift, pleiotropic selection, or random fluctuations in fitness can drive evolutionary divergence on the landscapes. The next section considers the first two mechanisms.
Multilocus Generalizations of the BDM Model
The BDM model was formulated in terms of only two loci. However, existing data on the genetics of reproductive isolation show that typically there are many different loci underlying reproductive isolation even at very early stages of divergence (Naveira and Masida 1998; Wu 2001; Wu and Palopoli 1994). If genetic architecture of reproductive isolation is known, the dynamics of speciation can be modeled, at least in principle. First, however, the genetic architecture of reproductive isolation is never known completely and, second, the mathematical treatment becomes extremely complicated as the number of genes involved increases. Therefore, rather than studying the dynamics of speciation given a specific genetic architecture, it becomes much more fruitful to look at the dynamics of speciation expected “on average.” (We have already used a similar approach in our discussion of the properties of typical adaptive landscapes in the previous section.) One such approach (Orr 1995; Orr and Orr 1996; Orr and Turelli 2001) reduces the complexity of adaptive landscapes underlying reproductive isolation to simpler effects of genetic “incompatibilities” of certain types that arise with certain probabilities. In what follows, two (or more) genes are called incompatible if their joint presence in an individual's genotype or in the genotypes of a (potential) mating pair results in reduction of a fitness component. The goal of this section is to find the average waiting time to speciation and the average duration of speciation.
Accumulation of genetic incompatibilities
Let us assume that there is a large number of diallelic loci potentially affecting reproductive isolation. Let any combination of k alleles at different loci be incompatible with probability q. (Note that in the BDM model k=2). Consider two populations that have diverged in d loci. Variable d is a measure of genetic distance separating the populations. As shown by Orr (1995), the expected number of incompatibilities between them is
where
is the binomial coefficient (i.e., the number of combinations of k objects chosen from a set of d objects). The value of J increases very rapidly (“snowballs”) with genetic distance d (approximately as the k-th order of d) . Notice that the snowball effect is much more pronounced with larger values of k.
How do the incompatibilities translate into reproductive isolation between the populations? Let us assume that complete reproductive isolation occurs when C incompatibilities separate the populations. (Note that in the BDM model, C=1). Then the probability that two genotypes at distance d are not reproductively isolated is approximately
(SG, unpubl. data). Here Γ(·, ·) and Γ(·) are the incomplete gamma function and gamma function, respectively (Gradshteyn and Ryzhik 1994) and J is given by equation (1). That the probability w decreases with d is compatible with a general empirical pattern that reproductive compatibility decreases with genetic distance between parents (Edmands 2002). The probability that two genotypes at distance d are reproductively isolated is 1 – w(d). One can also approximate the average K and the coefficient of variation CVK of the number of substitutions required for speciation (i.e., complete reproductive isolation): (S. Gavrilets, unpubl. data). As expected, K increases with the ratio C/q. The dependence of K on k is nonmonotonic with K minimized at some intermediate values of k. Figure 6 shows that as genetic divergence exceeds value K, the probability that no complete reproductive isolation occurs undergoes a rapid transition from 1 to 0. This “threshold effect” is especially strong when many complex incompatibilities are required for complete reproductive isolation. The latter feature is also apparent from the fact that the coefficient of variation CVK quickly goes to zero as both C and k become large.
It should be intuitively clear that the adaptive landscapes implied by this model are “holey.” However, in contrast to the Russian Roulette model described above, here the adaptive landscape is “correlated” in the sense that fitnesses of similar genotypes are similar which is more realistic.
Allopatric speciation
The simplest approach to model allopatric speciation is to assume that the populations accumulate substitutions at a constant rate, say ω substitutions per generation (Orr 1995; Orr and Orr 1996; Gavrilets 2000a; Orr and Turelli 2001). If “on average” K substitutions are required for complete reproductive isolation between two populations, the average waiting time to speciation is approximately
generations, where the coefficient 2 is because there are two populations accumulating substitutions. Thus, the average waiting time to allopatric speciation can be predicted on the basis of K alone without the need to specify the fitness function w(d) precisely. For example, if reproductive isolation results from C pairwise incompatibilities (k=2) , then equation (3) predicts that K =
, leading to the average waiting time to speciation
Orr (1995) and Orr and Turelli (2001) arrived to a similar expression using a different approach. How can one approximate the rate of substitutions ω? One way is to assume that genes underlying reproductive isolation have no other pleiotropic effects on fitness. In this case, genetic divergence will be driven exclusively by mutation and random drift (the null model of speciation). Let ν be the rate of neutral mutations per gamete per generation. If the within-population genetic variation is absent (or very low), each new mutation that is compatible with the current genetic state can be treated as selectively neutral. Then the rate of accumulation of substitutions can be written as
This approach predicts that the average waiting time to speciation does not depend on the population size either. However, in general the process of accumulation of incompatibilities cannot be treated as neutral even if the relevant genes do not control any other fitness components. This is because the alleles producing reproductive isolation are weakly selected against when rare (Nei et al. 1983; Gavrilets et al. 1998; Gavrilets 1999). A more careful analysis of the incompatibility model used here shows that large populations will actually evolve slower than small populations (Nei et al. 1983; Gavrilets et al. 1998, 2000; Gavrilets 1999). If alleles underlying reproductive isolation also increase fitness in the local conditions through some pleiotropic effects by a small amount sa, then
where as before Sa = 4Nsa. Strong selection for local adaptation will speed up both the genetic divergence and speciation.
One can also find the variances of these waiting times. If substitutions accumulate at a constant rate, then the coefficient of variation CVT of the time to speciation is the same as the coefficient of variation of the number of substitutions necessary for complete reproductive isolation (see eq. 3b). That is,
Orr and Turelli (2001) arrived to a similar expression for k = 2 using a different approach. Interestingly, equations (3a, 4, 5) show that speciation times of the organisms characterized by large C are expected to have both larger averages and narrower relative ranges (Orr and Turelli 2001).
Parapatric speciation
To predict the dynamics of parapatric speciation one needs to specify function w(d) precisely. A simple choice is the threshold function of reproductive compatibility
(Gavrilets et al. 1998; Gavrilets 1999). This function should be viewed as a limiting case of function w(d) given by equation (2) when complete reproductive isolation requires many complex incompatibilities (that is, when k or C are large; see the brief discussion of the “threshold effect” below eq. 3). The neutral case (no reproductive isolation) corresponds to K equal to the number of loci. Some other choices of w(d) were considered in Gavrilets (1999, 2000b).
Suppose there is a large number of loci and let ν be the probability of mutation per gamete per generation. As before, we assume that the population is subject to immigration at a constant rate m and that all immigrants have an “ancestral” genotype. Now speciation can be modeled as a random walk (Gavrilets 2000a) on the integers 0, 1, 2, … performed by the average pairwise distance d between the immigrants and the genotype most common in the population. Each new mutation fixed in the population will increase d by one. The fixation of an ancestral allele brought by the immigrants in a locus that has already diverged will decrease d by one. Speciation occurs when genetic distance d reaches the value of K + 1.
Neglecting within-population genetic variation, the average waiting time to speciation is approximately
(Gavrilets 2000a). The average duration of speciation is approximately where the number 0.577 is Euler's constant and Ψ(·) is the psi (digamma) function (Gradshteyn and Ryzhik 1994). [Function Ψ(K + 1) slowly increases with K and is equal to 0.42 at K = 1, to 2.35 at K = 10, and to 4.61 at K = 100.] For example, if m = 0.01, ν = 0.001 and K = 5, then the waiting time to speciation is very long: T = 1.35 × 1010 generations, but if speciation does happen, its duration is relatively short: τ = 1236 generations. Figure 7 illustrates the dependence of T and τ on model parameters in more detail. Notice that τ is of the order 1/ν across a wide range of parameter values.
Direct selection for local adaptation
Assume that each “new” allele improves adaptation to the local conditions. Let sa be the average selective advantage of a new allele over the corresponding ancestral allele and Sa = 4Nsa. The waiting time to speciation is approximately
where T is given by equation (6). The average duration of speciation is approximately For example, with the same values of parameters as above and Sa = 2, Ts = 2.74 × 104 generations and τs = 2170 generations. Thus, selection for local adaptation dramatically decreases T (in the numerical example, by the factor ≈ 50,000). Selection for local adaptation also somewhat increases τ relative to the case of speciation driven by mutation and genetic drift.
Effects of a genetic barrier
Assume that alleles underlying reproductive isolation have no other pleiotropic effects on fitness but that the population has already diverged from the source population in some other loci underlying adaptation to local conditions. Let γ be the corresponding gene flow factor resulting from selection for local adaptation. In this case, the migration rate m has to be replaced by the effective migration rate me = γm. The average waiting time to speciation can now be approximated as
where T is given by equation (6). Even a weak genetic barrier can significantly decrease the waiting time to speciation if there are many loci underlying reproductive isolation. For example, if γ = 0.1 and K = 5, then Tγ reduces by five orders of magnitude relative to the neutral case. The average duration of speciation is approximately and is increased by the presence of the barrier.
It would be very important to know the variances of the waiting time to speciation and of the duration of speciation. Unfortunately, no analytical results comparable to expression (5) are known for the parapatric case.
Most species consist of geographically structured populations, some of which experience little genetic contact for long periods of time (Avise 2000). Different mutations are expected to appear first and increase in frequency in different populations necessarily resulting in some geographic differentiation even without any variation in local selection regimes. This genetic differentiation by mutation and drift alone can result in allopatric and parapatric speciation. However, in the parapatric case the waiting time to speciation is relatively short only if a very small number of genetic changes is sufficient for complete reproductive isolation. If more significant genetic change is necessary, then parapatric speciation without selection for local adaptation is hardly possible. The results presented above show that even relatively weak selection (acting directly on the loci underlying reproductive isolation or indirectly via genetic barrier) reduces the waiting time to speciation by orders of magnitude.
The waiting time to speciation, T, is very sensitive to parameters: changing a parameter by a small factor can increase or decrease T by several orders of magnitude. Most of the parameters of the model (such as the migration rate, intensity of selection for local adaptation, the population size, and, probably, the mutation rate) directly depend on the state of the environment (biotic and abiotic) the population experiences. This suggests that speciation is triggered by changes in the environment (Eldredge 2003). If it is a significant environmental change that initiates speciation, the populations of many different species inhabiting the same geographic area should all be affected in a similar way. In this case, one expects more or less synchronized bursts of speciation in a geographic area, that is, a “turnover pulse” (Vrba 1985).
The results about the duration of speciation τ lead to two important generalizations. The first is that the average duration of parapatric speciation, τ, is much smaller than the average waiting time to speciation, T. This feature of the models studied here is compatible with the patterns observed in the fossil record which form the empirical basis of the theory of punctuated equilibrium (Eldredge 1971; Eldredge and Gould 1972; Gould 2002). The second generalization concerns the absolute value of τ which is on the order of one over the mutation rate for a subset of the loci affecting reproductive isolation for a wide range of migration rates, population sizes, intensities of selection for local adaptation, and the number of genetic changes required for reproductive isolation. Given a “typical” mutation rate on the order of 10–5 – 10–6 per locus per generation (Griffiths et al. 1996; Futuyma 1998) and assuming that there are at least on the order of 10–100 genes involved in the initial stages of the evolution of reproductive isolation (Singh 1990; Wu and Palopoli 1994; Coyne and Orr 1998; Naveira and Masida 1998; Wu 2001), the duration of speciation is predicted to range between 103 and 105 generations.
Sympatric Speciation
The greatest share of modeling work has been on sympatric speciation. In spite of this and for the reasons discussed in the introduction, no clear and simple quantitative results are generally known. Here I consider two simple models of sympatric speciation for which conditions for sympatric speciation can be found analytically.
The first model, developed by Udovic (1980) describes a diploid population with alleles A and a and the first locus and alleles B and b at the second locus. There are two common indices measuring the progress toward sympatric speciation in models of this kind. One is the heterozygote deficiency index
where y is the frequency of heterozygotes at a locus under consideration, p and 1 – p are the frequencies of the two alleles at the locus. The range of possible values of I is between 0 and 1. If the population is in Hardy-Weinberg proportions (which is expected if mating is random), y = 2pq and I = 0. If the population has split into two homozygous groups that do not mate and hybrids (i.e., heterozygotes) are completely absent, then y = 0 and I = 1. Another index is the normalized linkage disequilibrium, D′, between the two loci which is defined as Here D is the gametic linkage disequilibrium (defined as D = xABxaB – xAbxaB, where xAB, xab, xAb, and xaB are the frequencies of the corresponding gametes) and pA, pa, pB and pb are the corresponding allele frequencies. The range of possible values of D′ is from –1 to +1. If alleles are distributed randomly between different gametes, then D′ = 0. If the population has split into two genetic clusters so that the only gametes present are AB and ab, then D′ = 1. If the only gametes present in the population are Ab and aB, then D′ = –1.
In the second model, which deals with sexual conflict, individuals are haploid, and there are multiple alleles. The progress toward speciation in this model will be apparent from the formation of discrete genotypic clusters and from the degree of reproductive isolation between members of different clusters.
The Udovic model
Udovic's paper published in 1980 was the only one in almost a 40-year-long time interval that got very close to finding conditions for sympatric speciation analytically in a nontrivial model. However, both his approach and results were mainly ignored, being overshadowed by other papers describing numerical simulations. In this section, I describe and extend Udovic's ground-breaking results.
Consider a very large population with discrete nonoverlapping generations. There are two diallelic possibly linked loci with r being the recombination rate. The first locus with alleles A and a will be called the assortative mating locus or the AM locus. The AM locus controls assortative mating according to the symmetric version of the O'Donald model (O'Donald 1960). That is, each individual mates with probability α with another individual that has the same genotype at the AM-locus. With probability 1 – α the individual is engaged in random mating. If α = 0, the population is randomly mating. If α = 1, the population is split into three reproductively isolated genotypic clusters corresponding to the three genotypes AA, Aa, and aa, of which the cluster of heterozygotes will gradually disappear. In the O'Donald model, if there are no other loci, the allele frequencies at the AM locus do not change whereas the genotype frequencies approach an equilibrium at which the index of heterozygote deficiency is
(O'Donald 1960) The second locus with alleles B and b will be called the disruptive selection locus or the DS locus. The DS locus is subject to symmetric frequency-dependent disruptive selection favoring rare genotypes. This implies that if mating is random (that is, if α = 0), the population evolves to a stable polymorphic equilibrium where the allele frequencies are pB = pb = 1/2 and the heterozygotes experience a relative fitness loss which we will denote as S (0 ≤ S ≤ 1). If S = 1, heterozygotes are inviable and (postmating) reproductive isolation between the two homozygotes is complete.
In this model there always exists a line of equilibria where the allele frequency at the AM locus can be arbitrary, the deficiency of heterozygotes at the AM locus is given by equation (7), the allele frequencies at the DS locus are at 1/2, the genotype frequencies at the DS locus are in Hardy-Weinberg proportions (IB = 0) , and the linkage disequilibrium between the loci is absent (D′ = 0). In words, on this line of equilibria the loci behave as completely independent. Generically, this line of equilibria is locally stable if disruptive selection and assortative mating are weak and linkage is loose.
Under certain conditions the line of equilibria becomes unstable. If this happens, the population evolves to one of the two alternative polymorphic equilibria at which the frequencies of all four alleles are equal to 1/2, heterozygotes are in deficiency and there is a statistical association between the alleles in different loci. The two equilibria differ in the sign of linkage disequilibrium D′ and which one is approached by the population depends on initial conditions. At these equilibria, the population splits into two genetic clusters which arises within the population sympatrically and the loci mutually reinforce their effects leading to stronger reproductive isolation between the clusters. Evolution toward such an equilibrium represents (a step toward) sympatric speciation.
No linkage
If the loci are unlinked (r = 1/2), sympatric speciation occurs if
This condition has a very simple biological interpretation. Considering each locus in isolation, conditions S = 1 and α = 1 imply complete postmating isolation and complete premating isolation, respectively. Thus, equation (8) tells us that sympatric speciation occurs only if the cumulative strength of selection against hybrids and the strength of assortative mating, S + α, is larger than the threshold strength for each of them to cause complete reproductive isolation when acting in isolation.If condition (8) is satisfied, then at equilibrium the index of heterozygotes deficiency in the DS locus is
Notice that the condition for sympatric speciation (8) coincides with the condition for IB to be positive.
The index of heterozygote deficiency at the AM locus is
The normalized linkage disequilibrium is Note that the maximum possible value for each of these three measures is observed at S = 1 and is equal to α.
Linkage
With linkage (i.e., if r < 1/2), sympatric speciation occurs if
The equilibrium values of IA, IB and D′ are too cumbersome to be given here. Figure 8 illustrates the effects of parameters on the possibility of sympatric speciation and the degree of resulting genetic differentiation. Close linkage (i.e., small r) and strong selection are known to both make sympatric speciation easier and result in stronger genetic differentiation (Dickinson and Antonovics 1973; Udovic 1980; Felsenstein 1981). The inequality (9) and Figure 8 provide a quantitative description of these effects. I note that the equation for IA was given in the appendix of Udovic's original paper. The fact that positivity of IB is required for the stability of the equilibrium is also apparent from the numerical simulations reported by him.
As an example, let us consider the case of unlinked loci (r = 1/2) assuming linear frequency-dependent selection with fitnesses
where s > 0. This is one of the simplest cases of frequency-dependent selection (e.g., Gavrilets and Hastings 1995, 1998) which implies that rare alleles have a fitness advantage over common alleles. For example, genotypes using an under-explored resource can get a competitive advantage over other genotypes. In this model the equilibrium with allele frequencies pB = pb = 1/2 exists and is stable for any s > 0. At this equilibrium the relative fitness loss of heterozygotes is S = s/(2 + s) and condition (8) for sympatric speciation is For example, if s = 1 (that is, selection is weak), then α has to be larger than 0.67 (that is, strong assortativeness in mating is required) whereas if s = 8 (that is, selection is strong), then α has to be larger than 0.20 (that is, weak assortativeness in mating is sufficient). Alternatively, if α = 0.2 (that is, weak assortativeness in mating), then s has to be larger than 8 (that is, strong selection is required) and if α = 0.67 (strong assortativeness in mating), then s has to be larger than 1 (i.e., weak selection is sufficient). Figure 9 illustrates the distributions of genotype frequencies for three different combinations of parameters obtained by numerical iterations of dynamic equations.
A somewhat inconspicuous feature of the Udovic model inherited from the O'Donald model of assortative mating is that organisms pay no costs for being “choosy” no matter how rare their preferred mates are. Under this condition the Udovic model shows that sympatric speciation is possible if disruptive selection and/or assortativeness in mating are strong enough. Closer linkage between the AM locus and the DS locus promote sympatric speciation.
Sexual conflict
In the model considered in the previous section sympatric speciation resulted as an outcome of disruptive natural selection emerging from ecological interactions between individuals. In this section I discuss a model where sympatric speciation is a consequence of sexual selection.
Consider a large sexual haploid population with distinct nonoverlapping generations. We concentrate on two possibly linked multiallelic loci assuming step-wise mutation. The alleles Ai at the first locus are only expressed in females (or eggs), and the alleles Bj at the second locus are only expressed in males (or sperm). The genotype space in this model is identical to that in the Russian Roulette model considered above.
I will say that two individuals are “compatible” if mating, fertilization, and offspring development are not prevented by isolating mechanisms. Assume that the probability Ψij that a female (or egg) carrying an Ai allele is compatible with a male (or sperm) carrying a Bj allele is
This choice of function Ψij formalizes a general observation that successful reproduction requires certain “matching” of male and female genes (Tregenza and Wedell 2000). Specifically, function Ψij is maximized for all pairs of female and male genes with i = j. The coefficient σ controls the tolerance of matching: small σ implies that matching has to be pretty good for successful mating, and large σ implies that even when the value of |i – j| is large, mating can happen easily.
Let Pi be the proportion of the males in the population that are compatible with females carrying allele Ai. For each female genotype the value of Pi can be found if one knows the frequencies of male genotypes in the population. Pi can also be thought of as a proxy of a female mating rate. We assume that males can be involved in multiple mating and that they compete for fertilization opportunities. Although males generally benefit from high mating rates, in many situations multiple matings reduce female viability and total fitness, that is, there is a conflict between the sexes with regard to the optimum mating rate (Arnqvist and Rowe 1995; Chapman and Partridge 1996; Rice 1996; Holland and Rice 1998; Parker and Partridge 1998). To formalize this fact we assume that the overall probability that an Ai female leaves offspring is a unimodal function of Pi that reaches a maximum at a certain value Popt < 1 (Gavrilets 2000b; Gavrilets et al. 2001; Gavrilets and Waxman 2002). For example, in sea urchins, egg fitness is maximized at a level of sperm density which is much smaller than levels common under natural conditions (Franke et al. 2002). In our model, female mating rate is directly proportional to the proportion of compatible males and the assumption Popt < 1 formalizes the idea of sexual conflict over mating rate because for the males, it is optimal to have Popt = 1 since then all females are susceptible to fertilization by any male. To clarify the implications of the above assumptions, assume that the population is monomorphic for male allele Bj. Then Pi = Ψij and the females that have the optimum mating rate and the highest overall fitness are those for which Ψij = Popt. Using the definition of Ψ it is easy to see that there are two such female alleles: Aj+δ and Aj–δ both at distance δ from the male allele, where
This simple model exhibits three general dynamic regimes (Gavrilets and Waxman 2002). The first regime is an endless coevolutionary chase between the sexes in which females continuously evolve to decrease the mating rate while males continuously evolve to increase it (Holland and Rice 1998; Gavrilets 2000b; Gavrilets et al. 2001; Gavrilets and Waxman 2002). In this regime, there is a dynamic compromise between the sexes, and the proportion of compatible pairs is intermediate between Popt and 1. The coevolutionary chase is generically observed if the level of genetic variation is not too large.
Here we are concerned with two other regimes observed when the population size or mutation rates are sufficiently large. In the Buridan's Ass regime (Gavrilets and Waxman 2002) there is very low variation in male alleles maintained by mutation whereas female alleles split into two clusters both at the optimum distance δ from the male allele (see Fig. 10a). In this regime, males get trapped between the two female subclusters and have relatively low mating success.
In the sympatric speciation regime (Gavrilets and Waxman 2002) males answer the diversification of females by diversifying themselves and splitting into two clusters that start evolving toward the corresponding female clusters. As a result, the initial population splits into different genetic clusters (species) that are reproductively isolated and which have emerged sympatrically (see Fig. 10b). The regime of coevolutionary chase within-species ends after increasing genetic variation in female alleles leads to the splitting of female alleles into two subclusters within each species. By contrast, genetic variation in male alleles remains very low within each species. At equilibrium, female Pi values are close to Popt whereas males get trapped between two female subclusters and have low mating success.
The probability Ψij can be written as a function Ψ(d) of the “genetic distance” d = i – j between the female and male alleles. In the limit of very low mutation rates, sympatric speciation occurs if
If Ψ is a Gaussian function with zero mean and variance σ2, this inequality can be rewritten as (Gavrilets and Waxman 2002). If the above conditions are not satisfied, the population stays in the Buridan's Ass regime. Sympatric speciation requires small values of Popt implying that sexual conflict over mating rates must be strong. For example, if Popt = 0.67, then δ = 3 and condition (10a) is not satisfied (see Fig. 10a). If Popt = 0.45, then δ = 4 and condition (10a) is satisfied (see Fig. 10b). Sufficiently small values of Popt can result in more than two species emerging sympatrically (Gavrilets and Waxman 2002). The above results are not affected by the recombination rate between the loci.In this model sympatric speciation requires sufficiently strong selection (i.e., small Popt) and sufficiently strong assortativeness in mating (i.e., small σ). In contrast to the Udovic model, costs of being choosy are explicitly included in the sexual conflict scenario. That sympatric speciation still occurs is explained by the fact that the loci underlying reproductive isolation also experience direct selection for diversification induced by sexual conflict. In this model selection does not have to overcome the homogenizing effect of recombination that otherwise can prevent sympatric speciation (Udovic 1980; Felsenstein 1981; Rice 1984).
Conclusions
The theory of speciation needs many more simple and general analytical results allowing for transparent biological interpretation. In spite of the complexity of the processes leading to speciation, analytical approaches can be successful as demonstrated above. The major shift in the focus of theoretical studies needs to be from the demonstration that speciation is possible in a specific scenario (which has so far been the major goal of speciation models) to answering much more detailed questions about the probability of speciation, the waiting time to speciation, the duration of speciation, the degree of genetic and phenotypic divergence between sister species, the way different resources (including space) are partitioned between the sister species, etc. Several important generalizations about speciation have already emerged from analytical models.
(1) Unless the population size is small and the adaptive valley is shallow, the waiting time to a stochastic transition between the adaptive peaks is extremely long. This implies that it is very unlikely that a single peak shift will result in strong reproductive isolation. If transition does happen, it is very quick. This implies that observing it “in action” is practically impossible.
(2) Most species consist of geographically structured populations. Different mutations are expected to appear first and increase in frequency in different populations necessarily resulting in some geographic differentiation even without any variation in local selection regimes. Speciation can occur by mutation and random drift alone with no contribution from selection as different populations accumulate incompatible genes. This claim is based on models describing all three geographic modes (allopatric, parapatric, and sympatric). In the same way as divergence by random drift and mutation represents a null model of molecular evolution, speciation by random drift and mutation represents a null model of speciation. (Note that Wright himself believed that “the principal evolutionary mechanism in the origin of species must … be an essentially nonadaptive one” [Wright 1932, p. 364]).
(3) The importance of mutations and drift in speciation is augmented by the general structure of adaptive landscapes. Organisms have thousands of genes and millions of nucleotides. Although divergence in a small number of genes leading to strong (or complete) reproductive isolation is possible, in most cases divergence will include multiple genes. Theoretical considerations of multidimensional genotype spaces have led to the realization that their properties are quite different from those of spaces of low dimensionality. In particular, the multidimensional genotype spaces are characterized by the existence of nearly neutral networks and holey adaptive landscapes. Speciation can be understood as the divergence along these networks (driven by mutation, drift, and selection for adaptation to a local biotic and/or abiotic environment) accompanied by the accumulation of reproductive isolation as a by-product.
(4) The waiting time to speciation driven by mutation and drift is typically very long. Migration can significantly delay speciation. However, selection for local adaptation (either acting directly on the loci underlying reproductive isolation via their pleiotropic effects or acting indirectly via establishing a genetic barrier to gene flow) can significantly decrease the waiting time to speciation. Direct selection is much more effective than indirect selection. In the parapatric case the average actual duration of speciation is much shorter than the average waiting time to speciation.
(5) Theoretical studies predict extreme sensitivity of the probability of speciation and the waiting time to speciation on model parameters which in turn strongly depend on the environmental conditions. This suggests that in general speciation is triggered by changes in the environment. Theoretical studies also show that once genetic changes underlying speciation start, they go to completion very rapidly. This is so both for changes driven by strong selection and for changes driven by weak stochastic factors. Thus, the quantitative theory predicts the short duration of intermediate stages in speciation and the difficulties of observing their traces in the fossil record.
(6) Sympatric speciation is possible if disruptive selection and/or assortativeness in mating are strong enough. Sympatric speciation is promoted if costs of being choosy are small (or absent) and if linkage between the loci controlling assortative mating and those experiencing disruptive selection is close.
These ′ are but a few initial steps on a long road. A major future challenge for theoretical speciation research is to develop a comprehensive dynamical theory of speciation and to link microevolutionary processes with macroevolutionary patterns observed in the fossil record.
Acknowledgments
I am really grateful to N. H. Barton, C. R. B. Boake, M. Saum, and an anonymous reviewer for extremely helpful comments of the manuscript. Supported by National Institutes of Health grant GM56693 and by National Science Foundation grant DEB-0111613.
Literature Cited
Fig. 2.
Simple adaptive landscapes with two equal peaks. (a) One-locus two-allele diploid model with underdominance. (b) Disruptive selection on an additive quantitative character
Fig. 3.
Adaptive landscapes in the diploid BDM model. The height of a bar gives the fitness of the corresponding combination of genes. (a) Fitness of an individual. The presence of alleles a and B results in zero viability. (b) Fitness of a mating pair. Females carrying allele a do not mate with males carrying allele B. It is assumed that alleles at locus B are not expressed in females whereas alleles at locus A are not expressed in males. The arrows show two possible routes to the evolution of reproductive isolation
Fig. 4.
Formation of clusters of viable genotypes (painted in black) in two dimensions for different values of the probability of being viable p. Inviable genotypes are painted in white. (a) P = 0.15, (b) P = 0.30, (c) P = 0.45, (d) P = 0.60
Fig. 5.
Adaptive landscapes. (a) A rugged adaptive landscape formed by genotypes with fitnesses assigned randomly from a uniform distribution. (b) A holey adaptive landscape formed by genotypes with fitnesses within a narrow fitness band. Evolution along a holey landscape is nearly neutral
Fig. 6.
The probability of no reproductive isolation with C = 10. Different lines correspond to k = 2 (the shallowest), 4, 6, 8 (the steepest)
Fig. 7.
The average waiting time to speciation, T, and the average duration of speciation, τ, in the case of speciation driven by mutation and random drift. The horizontal axes give the ratio of migration and mutation rates. In (a), the vertical axis gives the product of T and the mutation rate ν on the logarithmic scale. In (b), the vertical axis gives the product of τ and the mutation rate ν on the linear scale. In (b), the lines correspond to K = 1,2,4,6,8, and 10 (from bottom to top). [Using composite variables for the vertical axes allows one to represent relevant dependencies in two dimensions]
Fig. 8.
The Udovic model. (a) Conditions for sympatric speciation. Speciation occurs for parameter values above the corresponding curve. (b) Isolation index IB for α = 0.50. The recombination rates are given next to the corresponding curves
Fig. 9.
The distributions of genotype frequencies after 200 generations starting with a random distribution. The loci are unlinked (r = 1/2) and s = 2, S = 0.5. The frequencies of double heterozygotes AB/ab and Ab/aB are pulled together. (a) α = 0.4 so that condition (8) is not satisfied and no association of A and B evolves (IA = 0.25, IB = 0, D′ = 0). (b) α = 0.6 so that condition (8) is satisfied but the association of A and B is weak (IA = 0.46, IB = 0.20, D′ = 0.50). (c) α = 0.8 so that condition (8) is satisfied and the association of A and B is strong (IA = 0.72, IB = 0.60, D′ = 0.82)
Fig. 10.
Population genetic states in the sexual conflict model. Parameters: mutation rate µ = 10–5, recombination rate r = 0.5, σ2 = 10, and female overall fitness wf = exp(–(P – Popt)2). (a) The Buridan's Ass regime (Popt = 0.6; the average of value of Pi is 0.64; there is a single male allele B0 and two female alleles A3 and A–3); 5000 generations of selection. (b) The sympatric speciation regime (Popt = 0.4; the average of values of Pi is 0.42); 8000 generations of selection
