## Abstract

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.

## GENERAL METHOD

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*the selective advantages of

_{B}*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*(1 − δ

_{AB}*r*),

*r*, and

_{AB}*r*(1 + δ

_{AB}*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*. We present the model in the case where loci are in order (

_{MA}*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*

^{ρ−1}

*D*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**

_{0}the vector of all

*n*

_{j}_{;0}. We then call

*n*the number of copies of

_{j}*A*on background

*j*at time

*t** and

**n**the vector of all

*n*. Finally, ϕ(

_{j}**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*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.

_{j}When *t* = *t**, the distribution ϕ is reduced to a single point at **n** = **n**_{0}. 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*at

_{j}*t*, assuming that all

*n*

_{j}_{;0}are small relative to

*N*(the term σ

_{jk}represents the expected change in

*n*accounting for possible movement from background

_{j}*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*at

_{j}*t*is approximately

*n*

_{j}_{;0}(for ), while the covariance between the changes in

*n*and

_{j}*n*is negligible for all

_{k}*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*'s, and thus ψ, for arbitrary

_{j}*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*'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

_{j}*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).

## TWO-LOCUS CASE

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*are functions of ω, and where

_{b}*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*(which is the frequency of

_{B}*B*at

*t*) or on background

*b*with probability

*p*(which is the frequency of

_{b}*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*and

_{B}*p*are given by the logistic equations(13)where

_{b}*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*) ≈ 1 +

_{B}p_{B}*s*−

_{A}*s*. Of these, a proportion

_{B}p_{B}*r*moves to the

_{AB}p_{B}*B*background, while a proportion 1 −

*r*stays on the

_{AB}p_{B}*b*background. The expected number of copies produced by an

*A*allele present on a

*B*background at

*t*is (1 +

*s*)(1 +

_{A}*s*)/(1 +

_{B}*s*) ≈ 1 +

_{B}p_{B}*s*+

_{A}*s*; of these, a proportion

_{B}p_{b}*r*moves to the

_{AB}p_{b}*b*background, while a proportion 1 −

*r*stays on the

_{AB}p_{b}*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 ψ_{0}) is a constant: ψ_{0}(ω) = ϕ(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 ψ_{0} 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*= 0.1,

_{B}*r*= 0.02,

_{AB}*t*= −100,

*t** = 100, and

*N*= 10

^{6}(and where

*p*is given by Equation 13). Note that the population size

_{B}*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*).

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*= 0.1,

_{B}*r*= 0.02 (solid line), and

_{AB}*r*= 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 Δ

_{AB}*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*= 0.04 than with

_{AB}*r*= 0.02); this generates an indirect selective pressure to increase recombination rates, which we quantify in the next section.

_{AB}## CHANGE IN FREQUENCY AT THE MODIFIER LOCUS

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*is the frequency of allele

_{m}*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*and

_{MB}*n*the numbers of copies of

_{mB}*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*'s are functions of ω

_{j}_{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*and

_{B}*p*are given by Equation 13. For example, the expected number of copies produced by an allele

_{b}*A*present on background

*MB*at

*t*is 1 +

*s*+

_{A}*s*(to the first order in selection coefficients). Of these, a proportion

_{B}p_{b}*r*(1 + δ

_{AB}*r*×

*p*)

_{M}*p*will move to the

_{b}*Mb*background—since the expected recombination rate between selected loci, experienced by a chromosome carrying the

*M*allele, is

*r*(1 + δ

_{AB}*r*×

*p*)—a proportion

_{M}*r*will move to the

_{MA}p_{m}*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*would not affect the equations above. Specifically, we may denote the recombination rates between loci (

_{MA}*M*,

*m*) and (

*A*,

*a*) in

*mm*,

*Mm*, and

*MM*individuals by

*r*(1 − δ

_{MA}*r*),

_{MA}*r*, and

_{MA}*r*(1 + δ

_{MA}*r*), respectively. However, because recombination between (

_{MA}*M*,

*m*) and (

*A*,

*a*) would affect genotype frequencies only when it occurs in

*Mm*heterozygotes, δ

*r*would not enter into the equations.

_{MA}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*for all backgrounds

_{A}*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*is independent of the initial frequency of the modifier (

_{M}*p*) and can be written in the form , where is a function of

_{M}*T*, θ, ρ

_{AB}, ρ

_{MA}, and

*Ns*. 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).

_{A}Figure 3 shows a test of our method against simulations, for θ = 0.1, ρ_{MA} = 0.01, ρ_{AB} = 0.2, and *Ns _{A}* = 10

^{4}(in the simulations

*s*= 0.1, δ

_{B}*r*= 0.5,

*N*= 10

^{6}, and the initial frequency of the modifier is

*p*= 0.5). In the simulation program, we start from a population where allele

_{M}*B*is at frequency

*p*= 1/(1 +

_{B}*e*

^{−T}) (the deterministic frequency of

*B*at time

*T*), allele

*M*is at frequency

*p*, and

_{M}*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*and

_{A}*s*. Finally, Figure 5 compares results obtained using our method and Otto and Barton's method for

_{B}*s*=

_{A}*s*= 0.1,

_{B}*r*= 0.001 (Figure 4A), and

_{MA}*r*= 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.

_{MA}## NET SELECTION GRADIENT OVER RECURRENT SWEEPS

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 Λ^{2}, 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*. Averaging over time gives the net selection gradient, for given positions of alleles

_{A}*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*. Moreover, we expect that should scale roughly with . Indeed, from Equation C23, takes the form , multiplied by a function of ρ

_{A}_{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*= 1/

_{B}*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.

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*is to reduce the genetic variance at the modifier locus,

_{MA}*p*. Indeed, in the absence of selection at locus (

_{M}p_{m}*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*tends to reduce the rate of increase of the modifier. Increasing the recombination rate ρ

_{MA}_{MA}decreases the variance in

*D*and thus decreases this effect.

_{MA}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*). Indeed, simulations indicate that the method gives accurate predictions as long as

_{B}*A*occurs after

*B*(

*e.g*., Figure 5). However, in the opposite case (when

*A*occurs before

*B*), and when

*s*=

_{A}*s*, we can simply swap

_{B}*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*, while Equation 31 depends only on the product

_{A}*Ns*. We performed numerically the triple integration of Equation 31 for the case where

_{A}*s*=

_{A}*s*= 0.05 and obtained

_{B}*s*

_{M}R^{2}/Λ

^{2}≈ 0.151, 0.126, and 0.102 for

*N*= 10

^{5},

*N*= 10

^{6}, and

*N*= 10

^{7}, 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

^{5}, 10

^{6}, and 10

^{7}, 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 Λ

^{2}/

*R*

^{2}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.

## DISCUSSION

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 this stochastic approach to a deterministic calculation of the hitch given to the modifier by the new selected allele, after establishment. By comparing our results with simulations, we found that this method provides very accurate predictions as long as the basic assumptions of the model are satisfied, that is, when allele *A* appears once allele *B* is already established.

As illustrated by Figure 6, our method also provides good results when the weakly selected allele *A* enters the population before the strongly selected allele *B* (both alleles being initially present as single copies), as long as *A* occurs not too long before *B*. Furthermore, Figure 6 shows that the modifier increases in frequency through two different mechanisms, whose relative importance depends on the relative times of occurrence of the two favored alleles: the modifier increases the probability of fixation of allele *A* (this effect being marked mostly when *A* occurs around the time of appearance of *B*), and it also accelerates the sweep of *A* (this effect being marked mostly when allele *A* is already established when allele *B* occurs). These two effects generate an indirect pressure selecting for the modifier. We have also seen that in some cases (when allele *A* occurs before allele *B*) increasing the recombination rate *r _{MA}* can lead to a greater expected increase in frequency of the modifier. This is due to the fact that increasing

*r*decreases the variance in linkage disequilibrium

_{MA}*D*. Indeed, the variance in

_{MA}*D*tends to reduce the genetic variance at the modifier locus, which in turn reduces the expected change in frequency of the modifier.

_{MA}Finally, we have seen that the net selection gradient for a modifier affecting recombination over a whole genome, during recurrent selective sweeps, takes the form *s _{B}*Λ

_{A}Λ

_{B}/

*R*

^{2}(where Λ

_{A}and Λ

_{B}are rates of substitutions and

*R*is the map length of the genome), multiplied by a function of θ =

*s*/

_{A}*s*and of

_{B}*Ns*; this function has to be obtained by integrating numerically over time and over the possible locations of selected loci. Because our method does not describe accurately the case where allele

_{A}*A*occurs a long time before allele

*B*, however, we could perform this integration only for the case where

*s*=

_{A}*s*. For

_{B}*s*=

_{A}*s*= 0.05, we obtained that the net selection gradient for the modifier is given by

_{B}*s*

_{M}R^{2}/Λ

^{2}≈ 0.151, 0.126, and 0.102 for

*N*= 10

^{5}, 10

^{6}, and 10

^{7}, respectively. Combining the present method with the method of N. H. Barton and D. Roze (unpublished data), however, we should be able to obtain more results for different values of θ and of

*Ns*.

_{A}According to current estimates of rates of amino acid substitutions for Drosophila and humans, the opportunity for indirect selection on recombination rates is rather small. A much higher estimate for the rate of adaptive substitution (of about one every 10 generations) has been advanced for Drosophila from a study of differences in genetic variability across loci (Schlötterer *et al*. 1997; Nurminsky 2001). While other factors than selective sweeps may also account for these differences, it remains possible that a proportion of advantageous mutations do not cause any change in protein composition (mutations in regulatory sequences, for example). In a multilocus simulation study, Otto and Barton (2001) obtained substantial changes in recombination rates due to Hill–Robertson effects between selected loci (in particular, Figure 4 in Otto and Barton 2001); similar results were obtained in much larger populations that are spatially structured (Martin *et al*. 2006). These results are not incompatible with the present results, however, as sweeps were occurring at several loci at the same time (high Λ), in a small-sized genome (small *R*), while population size was relatively small (*N* ≤ 1000). Such conditions (high Λ, small *N*) may explain increases in recombination rates obtained after selection experiments. Similarly, Iles *et al*. (2003) found substantial increases in recombination when selective sweeps occur simultaneously at more than two loci. However, such a situation may appear unlikely given the current estimates of rates of amino acid substitutions. Still, this does not mean that Hill–Robertson interactions have little effect on the evolution of sex and recombination. First, it is possible that adaptation occurs mostly during relatively short time periods, during which substitution rates are higher. Also, the direction of selection may fluctuate over time, in which case rates of substitutions may not be a good indicator of the frequency of sweeps occurring within populations. Second, population structure increases the effects of Hill–Robertson interference among loci (Martin *et al*. 2006). Third, the effect of deleterious alleles at mutation–selection balance may be far greater than the effect of beneficial alleles sweeping through populations: indeed, given the high rates of deleterious mutation that have been estimated for several species, populations should be polymorphic for deleterious mutations at many loci. The Hill–Robertson effect will also generate negative linkage disequilibria among these loci, reducing the efficiency of selection against deleterious alleles. More analytical and simulation work is thus needed before we can assess the plausibility of the hypothesis that stochastic interactions among loci play an important role in the evolution of sex and recombination.

## APPENDIX A

We express here the total change in frequency of the modifier after *t**, for given numbers of *A* alleles present on the *M* and *m* backgrounds at *t**. We assume that the change in genotype frequencies is deterministic after *t**; in this case, the change in *p _{M}* and

*p*(the frequencies of alleles

_{A}*M*and

*A*) over one generation can be written in the form(A1)(

*e.g*., Barton and Turelli 1991), where

*D*is the linkage disequilibrium between the two loci,

_{MA}*a*is a function of

_{A}*s*and

_{A}*p*, and

_{A}*q*= 1 −

_{A}*p*. We thus have(A2)Furthermore, it can be shown that

_{A}*D*/(

_{MA}*p*) decays at a rate 1 −

_{A}q_{A}*r*per generation, so that(A3)where

_{MA}*D*is the linkage disequilibrium

*D*at time

_{MA}*t** multiplied by population size

*N*, and

*S*is the number of copies of allele

*A*in the population at time

*t**. Furthermore, the dynamics of selection give(A4)giving(A5)with ρ =

*r*/

_{MA}*s*. We thus have(A6)The total change in frequency of allele

_{A}*M*is then obtained by integrating over

*p*(going from 0 to 1):(A7)(Barton 1998). For small ρ, is close to 1. The average change in frequency of the modifier over all possible realizations at

_{A}*t** is thus obtained by averaging the quantity

*S*

^{ρ−1}

*D*over the probability distribution ϕ at

*t**.

## APPENDIX B

#### Inverting the Fourier transform:

We explain here how to derive the probability distribution of the number of *A* alleles (*n*) at *t** from the Fourier transform ψ(ω), in the two-locus model. Because *n* can take only discrete values, the true probability distribution of *n* is discrete. We explained in the text how one can obtain the probability of loss of allele *A* (*n* = 0) from the value of ψ for large ω. Here we denote the true (and discrete) probability distribution of *n* at time *t**, for *n* > 0. Because we assume that allele *A* is still rare at time *t**, we can choose an integer *n*_{max} such that becomes vanishingly small for *n* > *n*_{max} (this allows us to reduce computation time). In the following we assume that *n*_{max} is an even number.

The discrete Fourier transform (DFT) of is given by(B1)The DFT takes *n*_{max} distinct values (for *k* = 0,…, *n*_{max} − 1) and is periodic: for any integer *j*, we have (periodicity comes from the fact that *e*^{ιx} = cos *x* + ι sin *x*). The sum in Equation B1 can also be written(B2)and can be approximated by the integral(B3)The function ψ_{>0} is the continuous transform of the distribution ϕ(*n*) for *n* > 0 and is given by ψ(ω) − ψ(∞), where ψ(∞) is the limit of ψ(ω) as ω tends to infinity (see text). The distribution can be recovered from by the inverse DFT:(B4)Because is periodic, this is also(B5)The number of terms to be evaluated can then be reduced by noting that the real parts of and are equal, while their imaginary parts cancel out (this can be seen from Equation B3). Thus, to invert the transform it is sufficient to calculate the real part of for *n*_{max}/2 + 1 values of ω_{k} between 0 and π,(B6)where stands for the real part of *x*. Computations could in principle be made more efficient by using Fast Fourier algorithms, although we did not explore this possibility here.

#### Calculating moments of ϕ:

Moments of the probability distribution ϕ(*n*) can be obtained directly from its Fourier transform ψ. For example, the expected number of copies of allele *A* at time *t** is given by(B7)Using Equation 12, and noting that *P _{B}* and

*P*= 0 when ω = 0, we obtain(B8)the derivatives being evaluated at ω = 0. Differential equations of the derivatives of

_{b}*P*and

_{B}*P*with respect to ω are obtained by differentiating Equations 14, giving(B9)From this, and noting that ∂

_{b}*=*

_{t}p_{B}*s*, one obtains(B10)and thus , corresponding to the deterministic increase of allele

_{B}p_{B}p_{b}*A*. Our method thus does not detect the effect of interference on the expected number of copies of

*A*at

*t*(because changes in frequencies are expressed only to first order); however, it shows an effect on higher moments of the distribution ϕ. For example, the variance of the number of copies of

*A*at

*t** is given by(B11)the derivatives being evaluated at ω = 0. After differentiating Equations 14 twice with respect to ω and simplifying, one obtains(B12)Equations B9 and B12 can be solved numerically to express

*V*as a function of

*t*. Figure 2 illustrates that selection at the (

*B*,

*b*) locus increases the variance of the number of copies of

*A*at

*t**.

## APPENDIX C

Here we develop a method to obtain the expected change in frequency of the recombination modifier, in the three-locus model. Equation 2 expresses this change in frequency in terms of the expectation of *S*^{ρ−1}*D* at time *t**, where *S* and *D* are given by(C1)Changing variables in Equation 17 (by expressing *n _{MB}* and

*n*in terms of

_{mB}*S*and

*D*), we can express the Fourier transform of the probability distribution ϕ(

*S*,

*D*) as(C2)with(C3)Using these new variables, the boundary conditions (19) become(C4)at

*t*=

*t**. To evaluate , we then proceed as follows. The derivative of ψ with respect to ω

_{D}, evaluated at ω

_{D}= 0, is given by(C5)where(C6)Equation C5 thus shows that the Fourier transform of is given by

*i*(∂ψ/∂ω

_{D}). Inversion of this Fourier transform gives(C7)We then have (for ρ < 1)(C8)Combining Equation C8 with Equation 2 leads to an expression for the expected increase of the modifier, in terms of ψ. This allows us to calculate the selection gradient,

*s*, acting on the modifier, as defined by (26),(C9)where the derivative of ψ is evaluated at ω

_{M}_{D}= 0, δ

*r*= 0.

The derivative of ψ that appears in Equation C9 can then be expressed in terms of partial derivatives of the *P _{j}*'s, using Equation 21. Using the fact that

*P*variables are of order

_{j}*s*, one obtains to first order in

_{A}*s*(C10)with(C11)(C12)where .

_{A}Using transformed variables, it can be shown that and do not depend on the initial frequency of the modifier, *p _{M}*. From Equation 24, one can show that the pairs , remain in the ratio −

*p*/

_{M}*p*for all

_{m}*T*(the derivatives being evaluated at , δ

*r*= 0). We can thus define the quantities(C13)We then define(C14)Using these variables, and after differentiating Equations 24 with respect to and δ

*r*, one obtains the following differential equations for and ,(C15)(C16)with boundary conditions at

*T*=

*T**. Equations C15 and C16 use the fact that , and when and δ

*r*= 0. Differential equations for the -variables are also obtained from Equation 24, giving(C17)(C18)with boundary conditions at

*T*=

*T**, and(C19)(C20)with boundary conditions at

*T*=

*T**. Finally, differential equations for and are given by(C21)(C22)with boundary conditions at

*T*=

*T**, with .

In practice, to evaluate the strength of selection on the modifier (*s _{M}*) when allele

*A*appears in the population at a given time

*T*, we solve numerically the system of differential Equations C15–C22 for a large array of values of . The selection gradient

*s*is then given by (from Equation C9)(C23)where the integral is calculated using the NIntegrateInterpolatingFunction routine of Mathematica (notebook available on request). Finally, it is important to note that as Equations C15–C22 do not depend on the initial frequency of the modifier (

_{M}*p*),

_{M}*s*is also independent of

_{M}*p*.

_{M}In the case where loci are in order *A*–*M*–*B*, Equations C17–C20 are changed to(C24)(C25)(C26)(C27)while in the case where loci are in order *M*–*B*–*A*, they become(C28)(C29)(C30)(C31)with ρ_{MB} = *r _{MB}*/

*s*. Differential equations on , , , remain the same.

_{B}## Acknowledgments

We thank Alex Kondrashov, Thomas Lenormand, Sally Otto, François Rousset, and Jay Taylor for helpful comments and discussions. N.B. is funded by the Royal Society and by Engineering and Physical Sciences Research Council (GR/T19537/01); D.R. is funded by the European Molecular Biology Organization (long-term fellowship ALTF 280-2004).

## Footnotes

Communicating editor: R. Nielsen

- Received March 25, 2006.
- Accepted May 4, 2006.

- Copyright © 2006 by the Genetics Society of America