In finite populations, genetic drift generates interference between selected loci, causing advantageous alleles to be found more often on different chromosomes than on the same chromosome, which reduces the rate of adaptation. This "Hill–Robertson effect" generates indirect selection to increase recombination rates. We present a new method to quantify the strength of this selection. Our model represents a new beneficial allele (A) entering a population as a single copy, while another beneficial allele (B) is sweeping at another locus. A third locus affects the recombination rate between selected loci. Using a branching process model, we calculate the probability distribution of the number of copies of A on the different genetic backgrounds, after it is established but while it is still rare. Then, we use a deterministic model to express the change in frequency of the recombination modifier, due to hitchhiking, as A goes to fixation. We show that this method can give good estimates of selection for recombination. Moreover, it shows that recombination is selected through two different effects: it increases the fixation probability of new alleles, and it accelerates selective sweeps. The relative importance of these two effects depends on the relative times of occurrence of the beneficial alleles.

---

ONE of the first general hypotheses proposed to explain the maintenance of sex and recombination states that sex increases the rate of adaptation, by allowing beneficial mutations initially present on different chromosomes to be combined onto the same chromosome (MORGAN 1913; FISHER 1930; MULLER 1932). The first analytical treatment of the question was done in the 1960s and led to some debate: using a simple calculation, CROW and KIMURA (1965) showed that the rate of incorporation of beneficial mutations is indeed higher in sexual than in asexual populations, whereas a deterministic model made by MAYNARD SMITH (1968) showed no difference between sexual and asexual populations. FELSENSTEIN (1974) pointed out that this disagreement stemmed from the different assumptions made by these models: Crow and Kimura represented a finite population in which advantageous alleles first arise as single copies, while Maynard Smith considered the case of an infinite population, in which deleterious alleles at two loci are maintained at mutation–selection equilibrium and become advantageous after a change in the environment. Because the model assumes no epistasis, the two loci always stay in linkage equilibrium, and recombination has no effect on the change in genotype frequencies. By contrast, in a finite population such as the one modeled by Crow and Kimura, selection and random drift combine to generate some linkage disequilibrium even in the absence of epistasis. For example, consider the case where a new beneficial allele enters the population while another beneficial allele at a second locus is already present. The new mutation can occur either on the good background at the second locus, in which case the linkage disequilibrium (LD) is positive, or on the bad background, in which case the linkage disequilibrium is negative. A simple calculation shows that, on average, this initial linkage disequilibrium is zero. However, when both advantageous alleles arise in coupling (LD > 0) they tend to increase in frequency rapidly, while when they are in repulsion (LD < 0) they compete against each other, which slows their rate of increase. Because situations in which LD < 0 tend to last longer, the initial variance in linkage disequilibrium generates, on average, negative linkage disequilibrium. The same "Hill–Robertson effect" occurs when the variance in LD is generated by drift due to finite population size in a single population (HILL and ROBERTSON 1966; FELSENSTEIN 1974) or in an arbitrary number of interconnected populations (MARTIN et al. 2006). The negative linkage disequilibrium generated by the Hill–Robertson effect reduces the rate of incorporation of advantageous alleles and limits the speed of adaptation. BARTON (1995b) quantified this effect by calculating the reduction in fixation probability of an advantageous allele when selection occurs simultaneously at another linked locus.

Increasing the recombination rate between loci decreases the strength of the Hill–Robertson effect and thus increases the rate of fixation of beneficial alleles. Simulation and analytical models have shown that this generates indirect selection on modifier genes that increase recombination rates (FELSENSTEIN and YOKOHAMA 1976; OTTO and BARTON 1997, 2001; ILES et al. 2003; BARTON and OTTO 2005; MARTIN et al. 2006). Recombination modifiers are not mere theoretical constructs, but have been identified in a number of experimental studies. In particular, selection experiments have shown that it is possible to artificially increase recombination within genomes, either by direct selection on recombination rates or as a by-product of selection for another trait (e.g., KIDWELL 1972; CHARLESWORTH and CHARLESWORTH 1985a,b; KOROL and ILIADI 1994; KOROL et al. 1994; and references in OTTO and BARTON 2001). In some cases, genetic factors affecting recombination rates have been mapped, revealing a variety of effects: some factors have a local effect, increasing recombination only over a small region of the genome (which can be located on the same or on a different chromosome), while others affect the overall recombination rate of whole chromosomes. Recombination rates are also known to evolve in nature: for example, TRUE et al. (1996) found that the genetic map of the island endemic Drosophila mauritiana is 1.8 times longer than the map of its sister species D. melanogaster, with a higher frequency of crossovers along each chromosome. Rapid adaptation to a new environment may have selected for a higher recombination rate, but a broader survey is needed to establish whether there is a correlation between adaptation and recombination rate.

The first analytical study on the evolution of recombination due to the Hill–Robertson effect was by OTTO and BARTON (1997). Their model represents two loci under directional selection: at the first locus, a new beneficial allele arises as a single copy, while at the second locus another beneficial allele is already established and sweeps through the population; finally, a third locus affects the recombination rate between the selected loci. Because recombination decreases interference between the selected loci, the new beneficial allele has a higher fixation probability when it occurs on the high-recombination background. The frequency of the modifier allele coding for higher recombination is then expected to increase by hitchhiking, on average, as the new beneficial allele goes to fixation. Otto and Barton's analysis proceeds in two steps: first, they derive the probability of fixation of the new advantageous allele, arising on each of the four possible backgrounds at the two other loci; then, they use a deterministic calculation to obtain the change in frequency of the recombination modifier, given that the new beneficial allele reaches fixation. This approach, however, has several limitations. In particular, although interference between selected loci due to the Hill–Robertson effect is taken into account in the calculation of fixation probabilities, it does not enter the deterministic calculation of the hitch given by the new beneficial allele to the modifier, given that it fixes. This last calculation indeed assumes linkage equilibrium between the selected loci. There are, however, three possible sources of LD, which, as discussed by Otto and Barton, are probably the cause of the quantitative differences they observe between analytical and simulation results. The first is the initial disequilibrium generated by the mutation event producing the new beneficial allele; although this disequilibrium is zero on average, it is sometimes positive (when the new allele arises in coupling with the other beneficial allele) and sometimes negative (when it arises in repulsion), and the effect of this disequilibrium on the change in frequency of the modifier does not cancel out. This initial disequilibrium is taken into account in some of the numerical calculations of Otto and Barton (dotted lines in Figures 4–9 of OTTO and BARTON 1997), but this still leaves important departures from simulation results in several cases. The second source of disequilibrium is random drift, which can have strong effects while the new allele is still rare. The effect of this disequilibrium is analyzed in a more recent model (BARTON and OTTO 2005), by considering small fluctuations around the deterministic increase of the new allele. This analysis, however, does not apply while the allele is present in just a few copies (in which case drift can generate important departures from the deterministic trajectory), which is just when the effect on the modifier is strongest. Finally, a third source of disequilibrium comes from the fact that when the new beneficial allele happens to reach fixation, it is more likely to be in coupling with the beneficial allele at the other selected locus. Therefore, the disequilibrium on the paths that lead to fixation of the new allele tends to be more positive than the disequilibrium observed on the paths leading to loss of the allele.

In this article, we show that these restrictions can be overcome by finding the distribution of numbers of the new mutation across the four backgrounds, after it has become established but while it is still rare. Since this initial stage involves the independent replication of rare mutants, it can be described by a branching process; the solution can then be taken forward deterministically, to find the net expected increase of the modifier allele. This assumes that population size is sufficiently large that the effects of drift can be neglected once alleles are established within the population. As we will see, this method involves solving numerically a set of differential equations; therefore, one has to resort to a computer to evaluate the strength of selection on the modifier, as with simulation methods. However, our method has two main advantages over simulations: first, it can stand in place of the millions of replicate simulations needed to get accurate estimates of the change in frequency of the modifier; second (and more importantly), it provides us with an understanding of the scaling relationships among the parameters. We will see that our method provides accurate predictions over a wide range of parameter values. It also captures two different effects of recombination modifiers: they affect the fixation probability of new advantageous alleles and also the speed of selective sweeps. These two effects both generate indirect selection at the modifier locus; we will see that their relative importance depends on the relative times of occurrence of beneficial alleles. Finally, we use our method to estimate the overall strength of selection on a modifier affecting recombination over the whole genome, when recurrent beneficial mutations occur throughout the genome. This shows that the rate of adaptive substitutions has to be relatively high for interactions between selected alleles to generate a substantial force on recombination modifiers. Although the estimated rates of adaptive amino acid substitutions for the Drosophila and the human lineages do not seem compatible with such high values, it is still possible that most adaptive changes occurred in relatively short periods, during which indirect selection on recombination may have been important. Despite this requirement for a high substitution rate, stochastic forces may still appear as a more likely candidate to explain the widespread occurrence of recombination than deterministic forces generated by epistasis among loci. Indeed, epistasis can select for increased recombination rates only when it is negative and not too variable among loci (BARTON 1995a; OTTO and FELDMAN 1997), whereas available data do not show any clear tendency for epistasis to follow this pattern. Finally, we will see that our model does not give accurate results for all possible parameter values and initial conditions: it best describes the case where a new, relatively weakly selected allele occurs when a more strongly selected allele is already sweeping through the population. We have also worked on another method that better describes the case where a relatively strongly selected allele occurs while a more weakly selected one is already established and the case where two beneficial alleles occur at about the same time (N. H. BARTON and D. ROZE, unpublished results). In these cases, recombination has little influence on fixation probabilities, but still has an effect in accelerating selective sweeps. Finally, the method of BARTON and OTTO (2005) deals with the effects of small fluctuations due to drift, once alleles are abundant. Together, these methods thus give us a complete understanding of the different regimes under which recombination can be favored.

The general setting is the same as in OTTO and BARTON (1997). We consider a haploid population of size N; mating is random, and generations are discrete and nonoverlapping. We assume that selection occurs at two loci, favoring allele A over allele a at the first locus and allele B over allele b at the second locus. We call s_A and s_B the selective advantages of A and B and assume no epistasis, so that the fitnesses of the different genotypes are given by

(1)

We assume that , and , so that the fate of alleles A and B is decided while they are still rare. We then assume that a third locus, with two alleles M and m, affects the recombination rate between the two selected loci (recombination modifier). The recombination rate between loci (A, a) and (B, b) equals r_AB(1 – r), r_AB, and r_AB(1 + r) in mm, Mm, and MM individuals, respectively; r thus measures the modifier effect, and for simplicity we suppose that this effect is additive. The recombination rate between loci (M, m) and (A, a) is denoted r_MA. We present the model in the case where loci are in order (M, m)–(A, a)–(B, b); the method applies more generally, however, and we also give results for different orderings of the three loci.

We call t the time at which allele A enters the population (as a single copy). We suppose that, at this time, allele B is already established (although it can still be rare). We see later that this assumption is not crucial and that in many cases we can still obtain accurate results even if B is not established at t. We also assume that allele M is frequent enough that stochastic fluctuations can be neglected at the (M, m) locus and suppose that, at time t, the (M, m) and (B, b) loci are in linkage equilibrium (thus, we neglect any associations between M and B that may have been generated by previous sweeps at yet other loci).

We decompose the spread of allele A into two phases: during the first phase, between times t and t*, A either is lost from the population by drift or becomes established (i.e., reaches a high enough number of copies so that it will almost certainly fix). In the following we assume that A is still rare in the population at t*, so that between t and t*, the replication of each copy of A does not depend on the presence of other A alleles in the population, in which case the change in number of A's can be described by a branching process. After t* we assume that A, if not lost, increases in frequency on a deterministic trajectory. As we said before, allele B at the other selected locus is supposed to be already established at t, and we assume that the frequency of B stays on a deterministic trajectory all the time. We suppose that t* is sufficiently far that B will have reached fixation at t*. We assume that the frequency of the modifier allele M is not too small, so that M changes in frequency solely due to its effect on the spread of allele A: when M helps A to recombine onto the good B background, it then tends to be carried along during the spread of A after time t*, by hitchhiking. Before t*, the effect of A on the frequency of the modifier is very small (because A remains at low frequency). Thus, we assume that the frequency of M remains constant between t and t*.

In summary, the only stochastic process in our model concerns the change in number of copies of A between times t and t*. During this time, the frequency of M does not change, while B increases in frequency on a deterministic trajectory, until fixation. After t*, B is fixed, A increases on a deterministic trajectory if it is still present, and M changes in frequency due to the linkage disequilibrium established with A during the first phase. Note that these different assumptions all stem from the assumption that population size is large and that M and B are common enough at t that the effects of drift at these loci can be neglected.

In APPENDIX A, we show that the expected total change in frequency of the recombination modifier after time t* is given by

(2)

where = r_MA/s_A, S is the total number of A alleles at t*, and D is the linkage disequilibrium between loci (M, m) and (A, a) at t*, multiplied by population size N. To express the expected change in frequency of the modifier, we thus need to calculate the expectation of S^–1D over the probability distribution of the number of copies of A on the different possible backgrounds at t*. We now present a method to derive this distribution.

In the following, we use index j to denote a possible genetic background on which allele A may be found (in our three-locus model, the possible backgrounds are thus mb, mB, Mb, and MB). As we said before, we assume that A is present as a single copy at time t; however, let us assume for the moment that the initial number of copies of A is arbitrary and call n_j;0 the number of copies of A on background j at t. We denote by n₀ the vector of all n_j;0. We then call n_j the number of copies of A on background j at time t* and n the vector of all n_j. Finally, (n) is the probability that the number of copies of A on the different backgrounds at t* is equal to n. Although it is difficult to derive the distribution directly, we show below that it is possible to obtain its Fourier transform , defined using the vector of dummy variables ,

(3)

where is the imaginary number , and where the integral is a multiple integral over all possible values of n. Note that although n takes only discrete values, we treat it here as a continuous variable (since we integrate over n, rather than summing); this will later allow us to use a diffusion approximation. The Fourier transform is thus a function of variables _j, and is the vector of all _j's (which has the same dimension as n). Considering rather than proves useful, because we can make use of the convolution theorem, which states that the Fourier transform of the distribution of a sum of independent random variables is equal to the product of the Fourier transforms of the individual distributions of these variables. Because we assume that copies of allele A replicate independently between t and t*, the distributions of the number of copies originating from each initial A allele are independent. Since n is the sum of these numbers, must take the form

(4)

where is the Fourier transform of the distribution of the number of copies of A at t* originating from a single copy, present on background j at t (the functions _j thus do not depend on the n_j;0's). Equation 4 can also be written

(5)

where each P_j equals . Note that instead of the Fourier transform, we could also use the Laplace transform or the generating function of , which also satisfies the convolution theorem. However, in the following we need to invert the transform, and for this the Fourier transform proves more convenient.

When t = t*, the distribution is reduced to a single point at n = n₀. The Fourier transform of such a distribution is given by , giving at t = t*

(6)

for all j. We then use a diffusion approximation to express the P_j's as a function of t, for a given value of t*. A backward diffusion on is given by

(7)

The term in Equation 7 represents the expected increase in n_j at t, assuming that all n_j;0 are small relative to N (the term _jk represents the expected change in n_j accounting for possible movement from background k and can be expressed in terms of recombination rates and selection coefficients). Assuming a Poisson number of offspring per parent, the variance in the change of n_j at t is approximately n_j;0 (for ), while the covariance between the changes in n_j and n_k is negligible for all j k. Substituting Equation 5 into Equation 7 then yields

(8)

(e.g., HARRIS 1963, chap. V.15). Equation 8, together with the boundary condition (6), can be used to express the P_j's, and thus , for arbitrary t. In general, Equation 8 will have to be solved numerically. One can then recover the distribution by taking an inverse Fourier transform or derive directly moments of from its Fourier transform ; this is illustrated in the next two-locus example. Because we use the diffusion approximation, we have to assume that selection and recombination are weak (so that the n_j's change slowly over time). However, we will see that this is not a major restriction of our model. Finally, we note that Equation 8 is identical to the equation giving fixation probabilities of mutations occurring on different backgrounds in BARTON (1995b) and OTTO and BARTON (1997), except that boundary conditions differ. As we will see later, fixation probabilities can be recovered from the limit of the Fourier transform as tends to infinity (this is because the distribution has a spike at zero, corresponding to the loss of the A allele, which is described by the limit of for large ). In the next section, we illustrate our method by considering the case of two selected loci (without a recombination modifier).

We consider here the effect of selection at the (B, b) locus on the distribution of the number of copies of A at t*, in the absence of a recombination modifier. Because we assume that B is fixed at t*, allele A can be found only on background B at that time. The distribution is thus a function of a single variable (the number of A alleles at t*, which will be on background B), which we simply denote n. The Fourier transform of (n) is given by

(9)

We also know from the convolution theorem that must take the form

(10)

(Equation 5), where P_B and P_b are functions of , and where n_B;0 and n_b;0 are the number of copies of A on the B and b backgrounds at time t, respectively. At t = t*, the distribution is reduced to a single point at n = n_B;0, and its Fourier transform is thus given by . This, with Equation 10, gives the boundary conditions

(11)

at t = t* (Equation 11 is a special case of Equation 6, constrained so that allele A is present only on the B background at t*). From now on, we assume that A enters the population as a single copy at t. This copy will be on background B with probability p_B (which is the frequency of B at t) or on background b with probability p_b (which is the frequency of b). Equation 10 can thus be written

(12)

Indeed, is the Fourier transform () conditional on the initial copy of A being on background B (in which case n_B;0 = 1 and n_b;0 = 0 in Equation 10), while is () conditional on the initial copy of A being on background b. Assuming that selection is relatively weak at the (B, b) locus, p_B and p_b are given by the logistic equations

(13)

where t = 0 is at the midpoint of the sweep of allele B.

The -coefficients that appear in Equation 8 are easily obtained. The expected number of offspring copies produced by an A allele present on a b background at time t is (1 + s_A)/(1 + s_B p_B) 1 + s_A – s_B p_B. Of these, a proportion r_AB p_B moves to the B background, while a proportion 1 – r_AB p_B stays on the b background. The expected number of copies produced by an A allele present on a B background at t is (1 + s_A)(1 + s_B)/(1 + s_Bp_B) 1 + s_A + s_Bp_b; of these, a proportion r_ABp_b moves to the b background, while a proportion 1 – r_ABp_b stays on the B background. This gives (using Equation 8)

(14)

with

(15)

These are the same as Equations 5a and 5b in BARTON (1995b). Equation 14, with boundary condition (11), can be solved numerically using the NDSolve routine of Mathematica (notebook available on request), giving (Equation 12).

To recover the distribution (n) from the Fourier transform (), it is useful to decompose into two parts: a distribution that is different from zero only when n = 0 (loss of allele A) and that is equal to (0) at this point [this distribution is thus given by (0)(n), where is Dirac's "delta function") and the rest of the distribution for n > 0 that we approximate by a continuous function on the range n 1. The Fourier transform of the first part (which we call ₀) is a constant: ₀() = (0) for all [because the Fourier transform of (n) is 1], while the Fourier transform of the second part (which we call _>0) tends to zero as tends to infinity (because the Fourier transform of a piecewise continuous function tends to zero as tends to infinity, e.g., LIGHTHILL 1958). Because is the sum of ₀ and _>0, its limit as tends to infinity is thus (0): therefore, we can obtain the probability that allele A is extinct at time t* from the value of for large . The rest of the distribution (for n > 0) can be recovered by inverting the Fourier transform _>0. In APPENDIX B, we show that this can be done for each value of n by calculating a sum of n_max terms, where n_max is an arbitrary value such that the probability that allele A is present in a number of copies greater than n_max at time t* can safely be ignored (remember that we assume that A is still rare at t*). One obtains

(16)

and where () is the limit of as tends to infinity (which can be determined by evaluating for a very large value of or from Equation 6 in BARTON 1995b, giving fixation probabilities). The number of terms to be evaluated can then be reduced by using some symmetries of (see APPENDIX B). In Figure 1 we compare this solution with simulation results, for s_A = 0.01, s_B = 0.1, r_AB = 0.02, t = –100, t* = 100, and N = 10⁶ (and where p_B is given by Equation 13). Note that the population size N does not appear in the derivation, but appears in the simulation program, which represents a Wright–Fisher model with multinomial sampling in each generation (our simulation program is written in C++ and is available on request). Figure 1 shows that Equation 16 accurately describes the probability distribution (n).

View larger version (12K): In this window In a new window Download PPT slide   FIGURE 1.— Two-locus model: probability distribution of the number of copies of allele A at time t* (n), in the absence of recombination modifier, from the analytical model (using Equation B6 with nmax = 5000, solid curve) and from simulations (shaded points). Parameter values are sA = 0.01, sB = 0.1, rAB = 0.02, t = –100, t* = 100, and N = 106. In the simulations, the distribution is obtained from 108 replications.

 
Moments of the distribution can be obtained directly from , without having to invert the Fourier transform (in essence, the method is similar to Ohta and Kimura's method for deriving moments using the diffusion approximation, e.g., OHTA and KIMURA 1969, 1971). For example, we show in APPENDIX B that the variance of is obtained by solving numerically a system of three differential equations, derived from Equation 14. Figure 2 shows the logarithm of the variance at time t* = 100, as a function of the time t at which allele A first entered the population for s_A = 0.01, s_B = 0.1, r_AB = 0.02 (solid line), and r_AB = 0.04 (dashed line). The dotted line shows the logarithm of the variance in the absence of interference (in the one-locus model), given by , with t = t* – t (e.g., FELLER 1951). Figure 2 illustrates the fact that increasing recombination between the selected loci decreases stochastic effects at the (A, a) locus (the variance is lower with r_AB = 0.04 than with r_AB = 0.02); this generates an indirect selective pressure to increase recombination rates, which we quantify in the next section.

View larger version (8K): In this window In a new window Download PPT slide   FIGURE 2.— Logarithm of the variance of the number of copies of allele A at time t* = 100, as a function of the time t when allele A first appeared. Curves correspond to the model predictions, circles and squares to simulation results (variance observed among 107 replications). Solid curve, filled circles: two-locus model, sA = 0.01, sB = 0.1, rAB = 0.02, N = 106. Dashed curve, squares: two-locus model, sA = 0.01, sB = 0.1, rAB = 0.04, N = 106. Dotted curve: single selected locus (no interference), sA = 0.01, N = 106.

 

We now introduce the recombination modifier locus (M, m). As we said before, because N is large we can assume that the frequency of allele M does not change between times t and t*; we denote this frequency p_M, while p_m is the frequency of allele m. Furthermore, since we suppose that allele B is fixed at t*, allele A can be found only on backgrounds MB and mB at t*. Calling n_MB and n_mB the numbers of copies of A on these two backgrounds at t*, the Fourier transform of the probability distribution is given by

(17)

From the convolution theorem, we also know that must take the form

(18)

(Equation 5), where the sum is over all possible initial backgrounds j = mb, mB, Mb, and MB, and where the P_j's are functions of _MB and _mB. At t = t*, we have boundary conditions

(19)

From Equation 8, one then obtains the system

(20)

where s_– and s₊ are given by Equation 15, while p_B and p_b are given by Equation 13. For example, the expected number of copies produced by an allele A present on background MB at t is 1 + s_A + s_Bp_b (to the first order in selection coefficients). Of these, a proportion r_AB(1 + r x p_M)p_b will move to the Mb background—since the expected recombination rate between selected loci, experienced by a chromosome carrying the M allele, is r_AB(1 + r x p_M)—a proportion r_MAp_m will move to the mB background, and the rest will stay on the MB background (note that because we need to express the -coefficients of Equation 8 only to first order in selection coefficients and recombination rates, we can neglect double-recombination events). Allowing the modifier to affect the recombination rate r_MA would not affect the equations above. Specifically, we may denote the recombination rates between loci (M, m) and (A, a) in mm, Mm, and MM individuals by r_MA(1 – r_MA), r_MA, and r_MA(1 + r_MA), respectively. However, because recombination between (M, m) and (A, a) would affect genotype frequencies only when it occurs in Mm heterozygotes, r_MA would not enter into the equations.

The system (20), with boundary conditions (19), can be solved numerically to obtain the P_j's for arbitrary t. In the case where A is present as a single copy at t, and assuming linkage equilibrium between the (M, m) and (B, b) loci at t, the expression for the Fourier transform becomes (from Equation 18)

(21)

Using Equations 19–21, one can express for any time t of occurrence of the A mutation.

In the following, we use scaled parameters and variables; this allows us to factor out selection coefficients and obtain more general results, in terms of the ratio of selection coefficients at both loci. In the diffusion limit, _j and P_j are of order s_A for all backgrounds j. We define the scaled variables and as

(22)

We then use the scaled parameters

(23)

The frequency of allele B as a function of time thus becomes . Using these scaled parameters and variables, Equations 20 become

(24)

From Equation 2, the change in frequency of the modifier after T* can be obtained from , where = r_MA/s_A = _MA/, S is the total number of A alleles at T*, and D is the linkage disequilibrium between loci (M, m) and (A, a) at T*, multiplied by population size N (see APPENDIX A). The variables S and D are given by

(25)

In APPENDIX C, we present a method to derive from the Fourier transform , in the case of a modifier of small effect (r small), and when linkage between the modifier and the (A, a) locus is sufficiently tight ( < 1). We then calculate the selection gradient on the modifier allele M, defined as

(26)

The expected change in frequency of the modifier over the whole process is given by

(27)

Interestingly, APPENDIX C shows that s_M is independent of the initial frequency of the modifier (p_M) and can be written in the form , where is a function of T, , _AB, _MA, and Ns_A. A measure of the strength of selection at the modifier locus, generated by the Hill–Robertson effect at the selected loci, is thus given by . As shown in APPENDIX C, obtaining for a given set of parameter values involves solving numerically a system of differential equations (Equations C15–C22) and performing a numerical integration (Equation C23). These operations have been implemented in a Mathematica notebook, which is available on request. Modifying our model to represent different orderings of the three loci is straightforward and leads to minor changes to our equations (see end of APPENDIX C).

Figure 3 shows a test of our method against simulations, for = 0.1, _MA = 0.01, _AB = 0.2, and Ns_A = 10⁴ (in the simulations s_B = 0.1, r = 0.5, N = 10⁶, and the initial frequency of the modifier is p_M = 0.5). In the simulation program, we start from a population where allele B is at frequency p_B = 1/(1 + e^–T) (the deterministic frequency of B at time T), allele M is at frequency p_M, and B and M are in linkage equilibrium. Allele A is introduced as a single copy in a random background; then for each generation we calculate the expected change in genotype frequencies and sample the new frequencies according to a multinomial distribution. We stop when allele A has disappeared or has reached fixation and measure the change in frequency of M over the whole process. Figure 3 shows that our method provides a good prediction of selection on the modifier, as long as T is not too negative (bottom, dashed curve). The discrepancy for T < –10 comes from the fact that the model assumes that allele B is established at T. However, when T is too negative, the initial frequency of B is too low for the allele to be established; this is shown in Figure 3, top, which plots the deterministic trajectory of allele B (solid curve) and its probability of fixation as a function of T (dashed curve), where we calculate the fixation probability from the standard one-locus diffusion formula

(28)

where _A = Ns_A (e.g., CROW and KIMURA 1970). When T < –15, Figure 3 shows that allele B has a high chance of being lost from the population, in which case selection on the modifier disappears. If we multiply the expression for the change in frequency of the modifier given by Equation C23 with the probability of fixation of B (given by Equation 28), we obtain a result that matches the simulations very well (Figure 3, bottom, solid curve). To generate the curves of Figure 3, we set the value of time T* to 50; however, other values of T* between 30 and 350 lead to indistinguishable curves. When T* < 30, allele A can be still present in the population but not established (which leads to an overestimation of the change in frequency of the modifier), while when T* > 350, allele A, when still present in the population, has already reached a substantial frequency (which may lead either to an underestimation or to an overestimation of the change in frequency of the modifier, depending on the parameter values, as the model ignores changes in modifier frequency before T*, but tends to overestimate the change after T*). The fact that the results are insensitive to T* over a large range of values shows that the method is consistent. In Figure 4, we show that our model also gives very good results for different values of the selection coefficients s_A and s_B. Finally, Figure 5 compares results obtained using our method and Otto and Barton's method for s_A = s_B = 0.1, r_MA = 0.001 (Figure 4A), and r_MA = 0.01 (Figure 4B), corresponding to the values used in Figures 7 and 9 of OTTO and BARTON (1997). Figure 5 shows that the present method provides a better prediction of the change in frequency of the modifier. In Figure 5A, it can be seen that both the model and simulations give a bimodal curve; this effect is explained in the next section.

View larger version (9K): In this window In a new window Download PPT slide   FIGURE 3.— (Top) Deterministic increase in frequency of allele B (solid curve) and probability of fixation of B (dashed curve, from Equation 28) from initial frequency pB (given by Equation 13) as a function of time T. (Bottom) Scaled selection gradient on the modifier , as a function of T. The dashed curve shows the results obtained from Equation C23, and the solid curve shows the results obtained from (C23) multiplied by the probability of fixation of B (from Equation 28). Dots show simulation results (average over 107 replications); error bars measure ±1.96 SE. Parameter values are: = 0.1, MA = 0.01, AB = 0.2, and NsA = 104; in the simulations sB = 0.1, r = 0.5, with N = 106 and pM = 0.5, while in the analytical model the splice point T* is set to 50.

 

View larger version (8K): In this window In a new window Download PPT slide   FIGURE 4.— Scaled selection gradient on the modifier times the probability of fixation of B as a function of T. Curves: results obtained from Equation C23, multiplied by the probability of fixation of B (from Equation 28). Circles and squares: simulation results (average over 107 replications); error bars measure ±1.96 SE (when the error bars do not appear, they are smaller than the dots). (A) sA = 0.01 and sB = 0.1 (solid curve, solid circles), 0.05 (dashed curve, squares), 0.02 (dotted curve, open circles). (B) sB = 0.1 and sA = 0.01 (solid curve, solid circles), 0.02 (dashed curve, squares), 0.05 (dotted curve, open circles). In each case rMA = 0.001, rAB = 0.02, N = 106. In the simulations r = 0.5 and pM = 0.5, while in the analytical model the splice point T* is set to 50.

 

View larger version (6K): In this window In a new window Download PPT slide   FIGURE 5.— Scaled selection gradient on the modifier times the probability of fixation of B as a function of T. Solid curves, results obtained from Equation C23; dotted curves, results obtained from Equations 7–10 and 19 in OTTO and BARTON (1997); dashed lines, results obtained by the method of OTTO and BARTON (1997), taking into account the initial linkage disequilibrium between selected loci. Dots: simulations results (average over 2 x 107 replications), error bars measure ±1.96 SE (in B, the error bars are smaller than the size of dots). Parameter values: sA = 0.1, sB = 0.1, rAB = 0.02, N = 106, and rMA = 0.001 (A) and 0.01 (B). In the simulations, r = 0.5 and pM = 0.5.

 

What will be the net selection gradient for a modifier increasing the total map length of a genome, under a constant input of beneficial mutations? We assume here that the modifier increases the recombination rate uniformly over the genome and that adaptive sweeps are sufficiently rare that we can neglect interactions among more than two beneficial alleles. More specifically, we derive an expression for the net selection gradient to the order ², where is the rate of adaptive substitutions, and neglect terms of higher order in . Our model allows us to consider two classes of mutations. One can assume, for example, that some strongly selected alleles enter the population at a given rate (sufficiently low that the effect of interactions among these mutations can be neglected), while weakly selected mutations occur at a higher rate and have a reduced fixation probability due to the recurrent sweeps of strongly selected alleles. The interaction between strongly and weakly selected alleles will select for higher recombination rates. Interactions among weakly selected alleles will also select for more recombination, but as the effect of these interactions is weaker (e.g., BARTON 1995b), they may be neglected (this should be checked, however, as this depends on how much more numerous weakly selected alleles are). To obtain the average effect of strong substitutions on weak ones, we can consider the case where allele B is strongly selected and allele A weakly selected. Equation C23 then gives the rate of increase of the modifier, given that B fixes, and for given positions of the (A, a) and (B, b) loci. The selection gradient, conditional on fixing both alleles, is then given by s_M/P_fix(A), where P_fix(A) is the fixation probability of allele A, which can be obtained from Equation 6 in BARTON (1995b). Because P_fix(A) can be written as , where is a function of T, , and _AB, while , the selection gradient given the fixation of both alleles is , which depends on T, , _MA, _AB, and Ns_A. Averaging over time gives the net selection gradient, for given positions of alleles A and B:

(29)

The net selection gradient, under constant rates of substitutions _A and _B for weakly and strongly selected alleles, and for all possible positions of selected loci, is then obtained by integrating over the genetic map,

(30)

where R is the total map length (because only selected loci that are sufficiently tightly linked to the modifier contribute significantly to its increase in frequency, we can take + as the integration limit when integrating over the genetic map). Equations 29 and 30 finally give

(31)

where the triple integral depends on and on Ns_A. Moreover, we expect that should scale roughly with . Indeed, from Equation C23, takes the form , multiplied by a function of _MA, _AB, and , say . Expanding f using a Taylor series and integrating over = _MA/ gives

(32)

which is dominated by the first term when is large. This slow decrease in the strength of selection on the modifier as population size increases probably reflects the fact that what matters most in a large population are stochastic effects occurring when selected alleles are present in small numbers of copies, which are not affected much by population size. Finally, it is important to note that the calculation of given above considers only pairwise interactions among selected loci and thus cannot be applied to the case where more than two selective sweeps overlap in time. However, as is discussed later, rates of amino acid substitutions in the human and Drosophila lineages are rather low, suggesting that interactions between more than two sweeps may be quite rare.

Equation 31 involves the selection gradient for the modifier, conditional upon the fixation of allele B (). Moreover, the integral over T includes cases where allele A enters the population before allele B. In the previous section, we compared simulation results with the model's predictions for the unconditional selection gradient for the modifier, . If we now condition upon the fixation of B, selection on the modifier can be substantial even when allele A occurs before allele B, but is still segregating at the time where B occurs. This is shown, for example, in Figure 6, which plots the selection gradient for the modifier, conditional upon the fixation of allele B (). The solid and dashed curves give the model predictions for _MA = 0.01 and _MA = 0.005 (respectively), obtained from Equation C23; dots represent simulation results, for s_B = 0.1. In the simulation program, allele B appears as a single copy at the time where p_B = 1/N according to the deterministic expression (13); this time is –13.8 for the parameter values used in Figure 6. For values of T < –13.8, allele A enters the population before allele B and appears on a random background at the modifier locus. Then, at time –13.8, allele B enters the population on a random background. The program stops when both selected loci are fixed, and the change in frequency of the modifier is measured among those cases where B has reached fixation. For values of T > –13.8, allele B enters the population before allele A; again, the program runs until both loci are fixed, and the change in frequency of the modifier is measured only among cases where B is fixed.

View larger version (10K): In this window In a new window Download PPT slide   FIGURE 6.— Scaled selection gradient for the modifier () as a function of the time of appearance of allele A (T), conditional on the fixation of allele B. Curves: results obtained from Equation C23, with MA = 0.01 (solid line) and MA = 0.005 (dashed line). Dotted curve: result obtained from Equations 7–10 and 19 in OTTO and BARTON (1997), with MA = 0.01. Dots: simulation results (average over 6 x 107 replications), with MA = 0.01 (solid circles) and MA = 0.005 (open circles); error bars measure ±1.96 SE. In the simulations, allele B enters the population as a single copy (on a random background) at time T = –13.8. Other parameters: = 0.1, AB = 0.2, NsA = 104; in the simulations, sB = 0.1, r = 0.5, N = 106, and pM = 0.5. In the analytical model the splice point T* is set to 50.

 
Several aspects of Figure 6 deserve attention. First, the model does not provide a good prediction of the change in frequency of the modifier for very negative values of T, that is, when allele A appears long before allele B. This is due to the fact that the model assumes that, at a time T* after the sweep of allele B, allele A is either lost or established, but still at low frequency in the population. However, for very negative values of T, A may have reached a high frequency at time T* or even be fixed. For example, for T = –125, simulations indicate that the frequency of A at time T* = 50, given that A will fix, is 0.99. In that case, the branching process approximation (independent replication of the different copies of A) has ceased to apply long before T* and probably even before the sweep of B. In this case, most of the change at the modifier locus occurs before T*, while the model assumes no change before T*. However, the model greatly overestimates the change after T*, because it overestimates at T* (this is confirmed by simulations), and also because Equation A7 assumes that allele A is still in low frequency at T*. Another potential problem is that the model allows p_B > 0 before allele B actually appears in the population (that is, before time –13.8). However, this does not seem to cause problems when A appears at time –50, for example. In any case, it appears unlikely that these problems can be solved by simple adjustments of our model: indeed, the discrepancy occurs when the sweeps of A and B overlap, which requires using the full three-locus recursions for describing the deterministic phase. We are currently working on an alternative method to treat this case, as is discussed next.

A second important result shown by Figure 6 is that the curves giving the selection gradient on the modifier as a function of T are bimodal (this can also be seen in Figure 5A and is confirmed by the simulations). This is due to the fact that the modifier has two different effects on the spread of allele A, whose relative importance depends on the time of appearance of A relative to B. When allele A appears after or shortly before allele B, its probability of fixation is reduced by the sweep of B. By reducing this interference, the recombination modifier increases the fixation probability of A and then benefits from the sweep of A (by hitchhiking). However, for values of T < –40, the sweep of allele B has very little effect on the probability of fixation of A (this can be seen, for example, in Figure 4 of OTTO and BARTON 1997): A is either lost or already established when B enters the population. In this case, the modifier does not affect the probability of fixation of A, but still has an effect in accelerating the sweep of A, by helping it to recombine onto the B background. This effect also generates an indirect selective force on the modifier and explains the leftmost maxima of the curves of Figure 6. Because the two peaks are separated by a time (which is the time needed for allele A to become established, in the absence of allele B), they will move apart as N increases (this is confirmed by simulations—not shown). It is possible in principle to calculate the advantage of the modifier via accelerating sweeps (left peak) using a deterministic calculation; indeed, one can relatively easily describe changes in genotype frequencies when A and B are rare and then solve for the change when they become common using eight coupled differential equations, describing evolution at the three loci. Such a calculation is presently being assembled (N. H. BARTON and D. ROZE, unpublished data).

Finally, Figure 6 shows that increasing the recombination rate _MA can have different qualitative effects on the rate of increase of the modifier, depending on the time of appearance of A relative to B (this can also be seen by comparing Figure 5A and 5B). When A appears after or shortly before B, increasing _MA decreases the rate of increase of the modifier, by decreasing the effect of hitchhiking. When A appears long before B, however, increasing _MA has the opposite effect. Before the time when B enters the population, the only effect of recombination between loci (M, m) and (A, a) is to decrease the variance in linkage disequilibrium D_MA. By using the method of BARTON and OTTO (2005), we found that the only effect of an initial variance in D_MA is to reduce the genetic variance at the modifier locus, p_Mp_m. Indeed, in the absence of selection at locus (B, b), the expected change in genetic variance at the modifier locus, over one generation, is given by

(33)

This corresponds to the classical result that selection tends to reduce genetic diversity at linked loci. Because the expected change in frequency of the modifier, at any generation, is proportional to the genetic variance at the modifier locus (at this generation), the variance in D_MA tends to reduce the rate of increase of the modifier. Increasing the recombination rate _MA decreases the variance in D_MA and thus decreases this effect.

Because the present method does not deal well with the case where allele A occurs a long time before allele B (as shown by Figure 6), it has to be coupled to another method describing this case better before we can use Equation 31 to obtain the net selection gradient for the modifier under recurrent sweeps (N. H. BARTON and D. ROZE, unpublished data). Nevertheless, we can still obtain some results for the case where = 1 (s_A = s_B). Indeed, simulations indicate that the method gives accurate predictions as long as A occurs after B (e.g., Figure 5). However, in the opposite case (when A occurs before B), and when s_A = s_B, we can simply swap A and B and continue to use the same method. However, to do so we need to incorporate the fact that allele B starts from frequency 1/N; we thus need to specify the values of N and s_A, while Equation 31 depends only on the product Ns_A. We performed numerically the triple integration of Equation 31 for the case where s_A = s_B = 0.05 and obtained s_M R²/² 0.151, 0.126, and 0.102 for N = 10⁵, N = 10⁶, and N = 10⁷, respectively. These results confirm the fact that the selection gradient roughly scales with for large ; indeed, gives 1.29, 1.36, and 1.34 for N = 10⁵, 10⁶, and 10⁷, respectively. The values given above lead to exceedingly small selection gradients, using current estimates for and R for Drosophila and humans. From comparisons between D. simulans and D. yakuba, the rate of adaptive amino acid substitutions has been estimated at one every 450 generations (SMITH and EYRE-WALKER 2002). This estimate has since been decreased to one every 800 generations (BIERNE and EYRE-WALKER 2004). The total map length R is estimated to be 3.6 for D. simulans (TRUE et al. 1996). For humans, the rate of adaptive amino acid substitution is estimated to be one every 200 years since the divergence with old-world monkeys (FAY et al. 2001), which, using a generation time of 25 years (EYRE-WALKER and KEIGHTLEY 1999) gives = per generation. The total map length is R = 34 for humans. In both cases, the ratio ²/R² is very small (of order 10^–5 for humans and 10^–7 for Drosophila). Combining the present method with the method of N. H. BARTON and D. ROZE (unpublished data) that studies the effect of recombination modifiers in accelerating sweeps and the method of BARTON and OTTO (2005) that deals with the effects of small departures from deterministic trajectories once alleles are common, we should be able to obtain estimates over a range of selection coefficients and population sizes. However, it appears unlikely that the selection gradient for recombination modifiers will be large, except during periods of frequent selective sweeps or in genomes with very tight linkage.

Apart from the idea that recombination brings a mechanistic advantage to individuals in terms of DNA repair (e.g., BERNSTEIN et al. 1988), the different hypotheses proposed to explain the maintenance of high rates of recombination in higher organisms state that recombination is advantageous because it reduces LD among selected loci (KONDRASHOV 1993; BARTON and CHARLESWORTH 1998; OTTO and LENORMAND 2002). Different possible sources of LD have been identified. First, LD may be generated by epistatic interactions among loci. In this case, increased recombination is selected when epistasis is weakly negative (FELDMAN et al. 1980; KONDRASHOV 1982, 1988; BARTON 1995a) and not too variable among pairs of loci (OTTO and FELDMAN 1997) or when epistasis fluctuates over short periods of time (CHARLESWORTH 1976; BARTON 1995a). However, experimental evidence does not indicate any strong trend for epistasis to follow this pattern, and it seems likely that the variance in epistasis is high in most cases (e.g., SEAGER and AYALA 1982; SEAGER et al. 1982; DE VISSER et al. 1996, 1997a,b; ELENA and LENSKI 1997; ELENA 1999; DE LA PEÑA et al. 2000; WLOCH et al. 2001). A second possible source of LD is environmental heterogeneity, in which case recombination is selected when directional selection at different loci covaries negatively between habitats (LENORMAND and OTTO 2000). Finally, LD may be generated by drift and directional selection: in a finite population, some negative disequilibrium does develop, on average, between selected loci, even in the absence of epistasis (HILL and ROBERTSON 1966; FELSENSTEIN 1974). By decreasing the effect of these interactions among loci, recombination increases the rate of fixation of advantageous alleles (BARTON 1995b). In this article, we presented a new method to quantify the strength of selection acting on a recombination modifier, due to this last effect.

Although the idea that sex and recombination increase the rate of adaptation has a long history in population genetics, theoretical models addressing the question of the effect of interference among loci on the evolution of recombination have been developed relatively late. OTTO and BARTON (1997) calculated the expected change in frequency of a modifier allele affecting the recombination rate between two selected loci in a finite population, using a method derived from branching processes. This model has led to important insights, but it does not always provide good quantitative estimates of the strength of selection at the modifier locus. This is due to the fact that interference among selected loci is neglected in part of the analysis: the linkage disequilibrium generated by the Hill–Robertson effect is taken into account in the calculation of probabilities of fixation of selected alleles, but not in the calculation of the hitch given by a selected allele to the recombination modifier, given that it fixes. In a more recent work, BARTON and OTTO (2005) used a different method to calculate the effect of linkage disequilibrium generated by drift on the change in frequency of a recombination modifier; however, their analysis assumes small departures from deterministic trajectories and thus does not apply when selected alleles start from a small number of copies (which is the case where drift will have the strongest effect, when population size is large).

In this article, we have presented a new method that overcomes these difficulties, in the case where population size is sufficiently large. The essence of the method lies in the calculation of the probability distribution of the number of copies of a new selected allele on different genetic backgrounds, once it is established but while it is still rare. Assuming a large population size, one can then neglect further stochastic fluctuations and express the change in frequency of the recombination modifier in terms of this probability distribution. Because the new allele remains rare during the establishment phase, the change in number of its copies on the different backgrounds can be described by a branching process, which generates a system of differential equations that can be solved numerically to obtain the Fourier transform of the probability distribution. This method fully takes into account the linkage disequilibrium between selected loci generated by the initial mutation event, by genetic drift during the establishment phase, and by the fact that the allele tends to be found more often in coupling with the other selected allele on the paths leading to its fixation than on the paths leading to its loss. We used this method to calculate the full probability distribution of the number of copies of a new selected allele when selection also occurs at a linked locus and found a good match with simulation results (Figure 1). To obtain an expression for the change in frequency of the recombination modifier, we coupled