Fluctuations around the deterministic trajectory:

We assume a discrete-generation model in which a population is chosen by random sampling from the propagules produced by the previous generation. The frequencies of the chromosomes among the propagules are given by the deterministic recursions for either a haploid population (with selection on haploids, followed by random mating and meiosis, to produce haploid propagules) or a diploid population (with selection on diploids, followed by meiosis and random mating to produce diploid propagules). From these propagules, either 2N haploid orN diploid juveniles are sampled with replacement, implying that the chromosome frequencies follow a multinomial distribution. We census the population immediately after sampling occurs.

Random sampling in any particular generation causes genotype frequencies to differ from that expected on the basis of the propagule frequencies. We measure the perturbations generated by a single round of sampling by the random variables ζ_j, ζ_k (for perturbations in the allele frequencies at locij,k) and ζ_jk (for perturbations in the linkage disequilibrium). The exact equations for change in the allele frequencies and linkage disequilibrium due to selection, recombination, and random genetic drift are given inappendix A.

The moments of the multinomial distribution can be used to determine the expected values of the perturbations generated by a single round of random sampling, as well as their variances and covariances. For instance, because random genetic drift does not cause any directional change in allele frequency,E[ζ_j] = 0, but the allele frequency will vary around its deterministic trajectory by an amount that is inversely proportional to the population size with Var[ζ_j] =p_jq_j/(2N). The first and second moments for the perturbations are given inappendix A. Higher-order moments are of the order 1/N² or smaller and are ignored in our analysis.

To describe the dynamics over multiple generations, we must keep track of the cumulative deviations from the deterministic trajectory, which we do using the vectorz = (δp_j, δp_k, δD_jk). The random variable δp_j measures the difference between the actual allele frequency at any point in time and the allele frequency predicted in the absence of random genetic drift. Similarly, δD_jk measures the difference between the actual linkage disequilibrium and the disequilibrium predicted deterministically. Throughout, we assume that the perturbations inz remain small and keep only the leading-order terms.

In the text, we assume that there is no epistasis and that the initial level of disequilibrium, if present, is of the orderN⁻¹. Under these assumptions, the disequilibrium predicted deterministically can be set to zero, so that δD_jk measures the actual amount of disequilibrium in the presence of drift. The method can be extended, however, to account for epistasis or for strong initial disequilibrium by keeping track of perturbations from the deterministic trajectory forD_jk (seeappendix A).

Recursions for cumulative effects of drift in the presence of selection can be written asz* =f(z) and can be determined from (A2). Assuming that the perturbations are small,z* can be expanded in a Taylor series,

1

where the subscripts (a,b,c) are set toj ork when referring to deviations in the allele frequencies and tojk when referring to the disequilibrium. The partial derivatives in (1) measure the sensitivity of each recursion to the perturbations caused by random sampling. They were derived from the recursions (A2) using the computer algebra packageMathematica (Wolfram 1991) and are given in supplementary information (S1;http://www.genetics.org/supplemental/). Note that all partial derivatives are evaluated along the deterministic trajectory (z =0) and are therefore independent of the perturbations. Taking the expectation of (1) gives

2

Similarly, second moments are given by the expectation of the product of (1) with itself:

3

Contributions to the first- and second-order moments arise by drift and are ofO(N⁻¹) as long as the allele frequencies are not too small (appendix A). In contrast, starting from a population with a specific genotypic composition, higher moments are initially absent, are generated by drift only toO(N⁻²), and are henceforth ignored. Note that the recursionequation (3) forEz^*_az^*_b also describes the change in the covariance,, becauseEz^*_aEz^*_b isO(N⁻²) and can be ignored. Similarly,Ez^*_a² describes the change in the variance, Varz^*_a, of the perturbations around the deterministic trajectory. Finally, theE[z_cζ_a] in (3) are alsoO(N⁻²) and are dropped.

Having dropped terms ofO(N⁻²) in (2) and (3) and using the fact that many of the partial derivatives in supplementary information (S1) equal zero, fluctuations in the allele frequencyp_j around the deterministic trajectory satisfy the following recursions:

4a

4b

Fluctuations in the allele frequencyp_k are described by analogous equations, with subscriptsj andk interchanged. Of greater relevance to the evolution of recombination are fluctuations that occur in the disequilibrium, δD_jk, around the deterministic equilibrium ofD_jk = 0:

5a

5b

5c

Eδp^*_kδD^*_jk is given by (5c), with subscriptsj andk interchanged. These equations apply even if selection coefficients vary through time, provided that there is no frequency dependence. Frequency-dependent selection could be modeled in a similar fashion, taking into account that the partial derivatives will contain additional components depending on the form of selection. With either constant or variable selection coefficients, recursions (4) and (5) can be evaluated numerically (Mathematica package available upon request) to predict the means, variances, and covariances of the perturbations over time.

The recursions for the perturbations can be further analyzed under the assumption that selection is weak (appendix B). Solution (B3) for the perturbations involving the disequilibrium is derived assuming that linkage is sufficiently loose relative to selection that the terms involving disequilibria reach a steady-state level over a faster timescale than the changes in allele frequencies; this steady-state level is known as the quasi-linkage equilibrium or QLE (Kimura 1965;Nagylaki 1993;Barton 1995a). Solution (B5) is derived assuming that both selection coefficients and recombination rates are small enough that the recursions can be well approximated by differential equations, which are then solved.

Development of disequilibria under directional selection:

We now consider how the above equations can be used to infer how disequilibrium develops in finite populations. Variance in linkage disequilibrium is generated directly by random drift, contributing thep_jq_jp_kq_k/2N term to (5b). As long as the direction of selection remains the same, all of the terms in (5c) are positive ∂δp^*_j/∂δD_jk, ∂δp^*_k/∂δD_jk, ∂δD^*_jk/∂δD_jk > 0, indicating that a positive covariance between the frequency of each beneficial allele and linkage disequilibrium will develop as a result of this variance in linkage disequilibrium. This covariance arises because allele frequencies increase more rapidly in populations that happen to have positive linkage disequilibrium. The presence of both variance in disequilibrium and covariance between allele frequencies and linkage disequilibrium then causes the expected linkage disequilibrium (5a) to depart from zero. Because ∂²δD^*_jk/∂δD²_jk is negative (appendix A), variance in disequilibrium across replicate populations leads, on average, to negative disequilibrium. Mathematically, this occurs because the mean fitness,, given by (A1) is higher with positive disequilibrium than with negative disequilibrium for any given set of allele frequencies. Consequently, the amount of positive disequilibrium gets divided by a larger quantity (²) in (A2b) and becomes relatively smaller in the next generation. Another, perhaps more intuitive, explanation, is that the expected disequilibrium across replicates becomes negative because selection is more efficient when good alleles are found on good genetic backgrounds (i.e., when disequilibrium is positive by chance), causing positive disequilibrium to decay more rapidly than negative disequilibrium whenever it occurs. The positive covariance generated between each allele frequency and the linkage disequilibrium leads the expected disequilibrium to become even more negative (∂²δD^*_jk/∂δp_j∂δD_jk, ∂²D^*_jk/∂δp_k∂δD_jk < 0). This occurs because those populations that happen to have positive disequilibrium and, consequently, faster rising allele frequencies (generating the positive covariance) also have higher mean fitness values in future generations. The linkage disequilibrium in future generations will then continue to be divided by a larger quantity (²) in (A2b), compounding the reduction in positive disequilibrium relative to negative disequilibrium. As a consequence of all of these effects, the expected value of linkage disequilibrium becomes negative even though it was initially zero and even though selection acts independently on the two loci.

Provided that allele frequencies are intermediate [i.e.,p,q areO(1)], the average amount of negative disequilibrium that develops is always proportional to 1/N: if the population size is doubled, the amount of disequilibrium generated by drift and selection is halved. When linkage is loose, disequilibrium does not accumulate appreciably over time. With tight linkage, however, substantial amounts of disequilibrium can accumulate within a population, especially when allele frequencies are initially low. This is a consequence of the fact that a small amount of disequilibrium generated by drift between rare alleles becomes a large amount of disequilibrium as the alleles approach a frequency of¹/₂. If the initial allele frequency is too small [O(1/N)], as is the case with a single favorable mutation, then higher-order terms in the Taylor series approximation (1) become appreciable and can no longer be ignored.

Figure 1 shows the distribution of allele frequencies over time in a population of size 2N = 10,000 when allelesP_j andP_k are initially at frequencyp₀ = 0.01 (s_j =s_k =r_jk = 0.1). The thick curve shows the mean trajectory (which is nearly equal to the deterministic trajectory), while the thin curves indicate the range within which allele frequencies lie 95% of the time as calculated from (4). These predictions were compared to simulation results based on a Monte Carlo simulation that sampled 2N chromosomes each generation using a multinomial distribution, with parameters equal to the frequency of each chromosome given by the deterministic recursions. Simulations (dots) indicate that (4) accurately captures the effects of drift on allele frequencies. Note that drift causes substantial variation in the trajectory of allele frequencies even in a population of size 10,000. This variation arises primarily during the early stages when only a small number of beneficial alleles are present (2Np₀ = 100) and represents random accelerations or decelerations in the spread of favorable alleles. On average, beneficial alleles spread slightly less rapidly than expected in an infinite population, indicating that the Hill-Robertson effect is present, if small, in populations of size 10,000. The effects of drift become even more important when there are initially fewer alleles within the population (e.g., 2Np₀ ≤ 10). In this case, however, iterating (4) and (5) overestimates the Hill-Robertson effect (results not shown), because higher-order moments in the perturbations are ignored and yet play an important role whenever allele frequencies are near zero.

Figure 2 shows how drift affects (a) the variance in disequilibrium scaled relative to the allele frequencies,EδD²_jk/p_jq_jp_kq_k, (b) the correlation between δp_j and disequilibrium,Eδp_jδD_jk/, and (c) the average disequilibrium. The average magnitude of the linkage disequilibrium is greater when the initial frequency of the beneficial alleles is lower, because random associations generated in the early stages persist and are amplified as the beneficial alleles rise in frequency (compare dotted, dashed, and solid curves, which correspond to initial frequencies 0.001, 0.01, and 0.1). In this example, wheres_j =s_k =r_jk = 0.1, the QLE approximation (B3) fails to capture the influence of the initial allele frequency, which is pronounced when selection is strong relative to recombination [e.g., the QLE prediction (B3a) for the scaled variance is constant at 0.00053]. On the other hand, the weak selection approximation (B5) overestimates the impact of initial fluctuations (compare thin to thick curves), although calculations with weaker selection (0.01) show that the approximation then agrees closely.

Figure 3 shows that the expected disequilibrium scales with 1/(2N) and matches the disequilibrium observed in simulations. The only exception is when the initial number of copies is very small (in this example, for 2Np₀ = 10). With few initial copies of the beneficial alleles, three issues arise that are not accounted for by our analysis. First, the fluctuations around the deterministic trajectory become so large that the Taylor approximation toEquation 1 is no longer adequate. Second, with few initial copies, beneficial alleles do not necessarily fix within the population. Third, those beneficial alleles that do fix are more likely to be positively associated with other beneficial alleles. Therefore, during the first few generations, rare beneficial alleles are more likely to be lost when negative disequilibrium develops (thereby destroying the disequilibrium) and more likely to persist when positive disequilibrium develops. Therefore, the average disequilibrium among beneficial alleles that have survived loss while rare might very well become positive, at least during the initial spread of the beneficial alleles. There is some evidence for this effect within the simulations: in the fifth generation, significantly positive disequilibrium was found when 2N = 10,000 and initial frequencies were 0.001 but in no other case. These positive associations will counteract the accumulation of negative disequilibrium to some extent, explaining, in part, the smaller amount of disequilibrium observed for this case (Figure 3,p₀ = 0.001).

Figure 4 shows the most extreme value of the expected linkage disequilibrium observed during the spread of allelesP_j andP_k as a function ofr_jk. Results from Monte Carlo simulations are well approximated by iterating Equations 4 and 5 (compare dots to thick dashed curve forp₀ = 0.01 and thick solid curve forp₀ = 0.1). The QLE approximation (B3c; thin curve) is accurate only for loose linkage (r_jk >s_j,s_k = 0.1) and fails to capture the dependence of disequilibrium on initial allele frequencies. With tight linkage, random fluctuations produced when small numbers of the favorable alleles are present persist to produce large fluctuations when the alleles become common. Consequently, the most extreme disequilibrium observed is very sensitive to starting allele frequencies when linkage is tight. Whenp₀ = 0.01 andr_jk < 0.01, drift and selection are substantial enough to drive the expected linkage disequilibrium to near its minimum value (−0.25), even though disequilibrium is initially absent and even though selection acts independently on all loci (s_j =s_k = 0.1).

Including a modifier of recombination:

The negative disequilibrium generated by drift in the presence of directional selection reduces the genotypic variability within a population and slows the response of that population to selection. A modifier allele that increases recombination among the selected genes will regenerate the genetic variation in fitness hidden in linkage disequilibrium. In particular, individuals carrying such a modifier allele are more likely to produce offspring of the fittest genotype. This effect can select for increased recombination at modifier loci, an effect examined in this section. Suppose that a modifier of recombination segregates at locusi. We assume that the modifier has no direct effect on fitness and alters recombination rates by a small amount δr. If terms ofO(δr²) are neglected, then the only significant terms involve the change in recombination between the selected loci (denoted δr_j,k|i). The modifier may also alter its own rates of recombination with other loci (δr_i,j|i, etc.) and may show dominance (δr_j,k|i,i). However, these effects are negligible for weak modifiers (Barton 1995a).

The equations for the effects of drift on the modifier frequency and its associationsD_ij,D_ik,D_ijk are given by (A2c–e). Disequilibria are measured as central moments of effects (Barton 1995a). For two loci, this definition coincides with the standard two-locus measure of gametic disequilibrium. For three loci,D_ijk is equal to ∑g[X](X_i −p_i)(X_j −p_j)(X_k −p_k), where the sum is taken over all chromosomes, each at frequencyg[X], whereX is a vector indicating whether or not an allele is present at each locus (i.e.,X_i is 0 if alleleP_i is absent and 1 if it is present). A modifier allele that increases recombination between locij andk becomes associated with the fittest genotype,P_jP_k, if that genotype is underrepresented within the population (i.e., ifD_jk is negative). This association is represented by the first term in (A2e), which generatesD_ijk in proportion to −δr_j,k|iD_jk. This three-way association then generates positive two-way associations,D_ij andD_ik, between a modifier allele that increases recombination and the beneficial alleles at locij andk (A2d). All three associations,D_ij,D_ik,D_ijk, cause changes in the frequency of the modifier allele (A2c). Note that because the modifier is assumed to have a small effect, these equations are linear inD_ij,D_ik,D_ijk, and all will change in proportion to δr_j,k|i.

Equations 2 and3 can be used to find the mean and variance of the three-locus system, as before. The first step is to find the first and second differentials of (A2c–e), such as ∂δp^*_i/∂δD_jk. These differentials are calculated as in supplementary information (S1) and are available upon request. The second step is to calculate the expected values and the covariances for the perturbations involving the modifier by incorporating these differentials intoEquations 2 and3. Following this method and dropping terms ofO(N⁻²), the recursion describing the cumulative effects of drift and selection at locij andk on the frequency of a modifier allele at locusi is

6

where φ_j = 1 +s_j(p_j −q_j)/2 and φ_k = 1 +s_k(p_k −q_k)/2. The expected change in the modifier is driven both by expectation termsE[δD_ij] and by covariances such asE[δp_jδD_ij], each of which is proportional to δr_j,k|i/(2N). The recursions for these expectation terms are complicated but can be calculated symbolically in a straightforward way from (2) and (3). These recursions may be iterated numerically or may be approximated by assuming that selection is weak as described in supplementary information (S2;http://www.genetics.org/supplemental/).

Assuming weak selection and loose linkage, the system approaches quasi-linkage equilibrium, and the per generation change in the frequency of a modifier at QLE can be found (supplementary information, S2). This gives a complicated function of the recombination rates (S2.4), which is proportional to δr_j,k|ip_iq_iV_jV_k/2N, where and are the additive genetic variances in diploid fitness at the two selected loci. Modifier alleles that increase recombination (δr_j,k|i > 0) thus increase in frequency even in the absence of epistasis in finite populations.

Expression (S2.4) takes a simple form if there is no genetic interference and if recombination rates are small in absolute terms but large relative to the selection coefficients, as required for QLE (i.e.,s_j,s_k ≪r_jk,r_ij,r_ik,r_ijk ≪ 1). Assuming, without loss of generality, that locusk is the right-most locus, we can describe the relative position of the modifier by defining α =r_ik/r_jk; α varies between zero and one when the gene order isjik and is greater than one when the gene order isijk. Then, keeping only the leading-order term in (S2.4), the per generation change in the modifier becomes

7a

where

7b

7c

Figure 5 illustratesg(α) and shows how evolution of a modifier of recombination depends on the relative positions of the loci, α. The functiong(α) rises rapidly near α = 0 and α = 1, where it is approximatelyg(α) ≈ 2/α² and 2/(1 − α)², respectively, although the QLE approximation breaks down if recombination rates become too small relative to selection. For gene orderjik,g(α) reaches a minimum value of 37¹/₃, when the modifier is halfway between the selected loci. For gene orderijk,g(α) lies between 2/(1 − α)² and 4/(1 − α)² and decreases rapidly for relatively distant modifier loci.

Equation 7a demonstrates that the effect that drift and selection have on a modifier of recombination falls precipitously as recombination becomes more frequent within the genome (∼r⁻⁵_jk). If selected loci are scattered randomly over a genome, then the effect of a modifier on tightly linked loci will dominate its dynamics. As recombination rates become comparable with the selection coefficients, however, the QLE assumptions break down, making it impossible to determine the average force on the modifier when selected loci are scattered randomly across a genome. Note that as selection and recombination decrease together toward zero in (7a), selection on the modifier approaches a positive limit, inversely proportional to population size (i.e., ifs,r, and δr are all small and of the same order, ΔE[δp_i] ∼ 1/N). This limit appears puzzling, because the average frequency of the modifier should not change in the absence of selection. There is no contradiction here, because we have assumed throughout thats ≫ 1/N when dropping terms ofO(N⁻²). Thus, our approximations leading to (7) require thatN increases toward infinity ass decreases toward zero, causing ΔE[δp_i] to decrease implicitly as selection weakens.

InFigures 6–9, we describe changes in the frequency of the modifier using the standardized function, β_i = Δp_i/(δr_j,k|ip_iq_i), which we call theselection gradient acting on the modifier allele. When multiplied by the genetic variance in recombination within the population,, the selection gradient describes the expected increase in recombination (= 2Δp_iδr_j,k|i). The selection gradient can also be related to the effective selection coefficient,s_e, defined as the function that causes Δp_i to equals_ep_iq_i. For a modifier of effect δr_j,k|i, the effective selection coefficient iss_e = β_iδr_j,k|i. It is important to recall that the modifier does not directly experience selection. Rather, the selection gradient measures the indirect selection arising from genetic associations between the modifier locus and the directly selected locij andk. To describe the total effect of selection at linked loci on recombination, we use thecumulative selection gradient, ∑β_i, which represents the sum of β_i from the initial generation up until a given generation, or thenet selection gradient, β_i,net, which represents the sum of β_i over the course of selection at locij andk until fixation of the beneficial alleles.

As a consequence of the fact that disequilibrium is more likely to develop by drift when allele frequencies are initially low at selected loci (Figure 3), the selection gradient acting on a modifier is larger under these conditions, as shown inFigure 6. Simulation results for the change in frequency of the modifier are in general agreement with iteration of (2) and (3), even though the assumptions of weak selection and a weak modifier are violated in the simulations (wheres_j =s_k = 0.1 and δr_j,k|i =r_jk/2). The largest discrepancies are observed when the initial number of alleles is low (e.g., when 2N = 10,000 andp₀ = 0.001), such that the perturbations due to drift become larger than allowed by our analysis. It is worth emphasizing that the simulation results presented inFigure 6 represent the mean ± standard errors based on 10⁶ replicates. Consequently, the standard deviations are enormous (multiply the length of the bars by 1000). Thus, while our analysis accurately predicts the mean change in frequency of a modifier, the effect of one bout of selection on a modifier is highly variable with only two selected loci and one modifier locus. Only when many loci affect fitness and recombination does the frequency of a modifier of recombination increase consistently (Otto and Barton 2001;Ileset al. 2003).

The net selection gradient, β_i,net, acting on a modifier over the entire time course of substitution at two linked loci is shown inFigure 7. The net increase of the modifier can be large with tight linkage: for example, withr_ij =r_jk = 0.01, Δp_i = 38p_iq_iδr_j,k|i. For very tight linkage, the change in the modifier is underestimated by iteration of (2) and (3) because the disequilibria caused by drift become substantial and higher-order effects begin to contribute [e.g., forr_jk = 0.001 andr_ij = 0.01, Δp_i = 334p_iq_iδr_j,k|i from simulations but only 75p_iq_iδr_j,k|i from (2) and (3)]. The QLE approximation (S2.4) is fairly accurate for loose linkage (r_ij,r_jk >s_j,s_k = 0.1) but substantially overestimates the change in modifier frequency for tight linkage, as shown by the thin lines inFigure 7.

Fluctuating polymorphisms:

The analysis described in this article can be used whether selection coefficients are constant or vary over time, as long as the population size is reasonably large and the numbers of each allele are never very small. With directional selection, indirect selection on a modifier ceases with the fixation of beneficial alleles. With fluctuating selection, on the other hand, polymorphism at the selected loci can be maintained over longer periods of time, prolonging selection for increased recombination. Even in the absence of epistatic interactions among loci and even when a population is initially in linkage equilibrium, drift in a finite population subject to fluctuating selection will lead to the accumulation of disequilibrium that reduces genetic variability in fitness. We have usedEquations 2 and3 to track the change in modifier frequency in a population subject to sinusoidal fluctuations in selection at two loci withs_j[t] = α_jCos(2π ·t/τ) ands_k[t] = α_kCos(2π ·t/τ + Φ), where τ is the period (assumed to be the same for both loci) and Φ measures the shift in phase of selection at the two loci (Figure 8). Note that when the direction of selection changes, the selection gradient on the modifier weakens (Figure 8a) or becomes negative (Figure 8b). Nevertheless, recombination rates rise over successive cycles.

The QLE result (S2.4) continues to apply when selection fluctuates over time, as long as recombination is frequent relative to the strength of selection and to the changes over time in selection (s_j[t],s_k[t], 1/τ ≪r_ij,r_jk,r_ik,r_ijk). Under sinusoidal selection, the product of the additive genetic variances in fitness at the two loci (V_jV_k) is ∼α²_jp_jq_jα²_kp_kq_k if the allele frequencies do not change substantially over a cycle. This result indicates that, when linkage is loose, indirect selection on a modifier will be stronger when the allele frequencies remain near¹/₂, such that the genetic variance in fitness remains high.

For weak selection and very tight linkage (rτ ≪ 1), we can obtain an approximate solution to the recursion equations (S2.2) and (S2.3). Assuming that several cycles of fluctuating selection have passed such that the dynamical system has reached a steady state, the solution to these recursions may be approximated as described in supplementary information (S3;http://www.genetics.org/supplemental/). The change in the modifier per generation averaged over one cycle is, to leading order in the recombination rates,

8

whereH_jk is the harmonic mean value ofp_jq_jp_kq_k, Varp^′_j and Varp^′_k measure the variance in the rate of allele frequency change over the cycle, and Covp^′_j,p^′_k measures the covariance in these rates at the two loci:

9a

9b

Figure 9 comparesEquation 8 to the iterated solution of (2) and (3) as well as to the QLE solution.Equation 8 is accurate only when selection is approximately the same strength as the product of the recombination rates and the period (α ≈rτ as inFigure 9b), which reflects the order assumptions made in supplementary information (S3). In other cases, the iterated solution of (2) and (3) should be used.

As with the QLE solution (S2.4),Equation 8 indicates that the change in allele frequency at a locus that modifies recombination is proportional to the effect of the modifier on recombination (δr_j,k|ip_iq_i), the inverse of the population size (2N), and the square of the selection coefficients at each selected locus (through Varp^′_j Varp^′_k). Unlike the QLE solution for loose linkage,Equation 8 is not maximized when the allele frequencies remain near¹/₂, because the harmonic mean of the allele frequencies in the denominator is much smaller when the allele frequencies approach zero or one. Consequently, with tighter linkage, selection on a modifier is maximized when there are high-amplitude fluctuations in allele frequencies, such that beneficial alleles reach low frequencies (where drift is stronger) as well as intermediate frequencies (where selection on a modifier is more effective). Selection on the modifier is also strongest when changes in allele frequency at the two loci occur at similar times, maximizing Covp^′_j,p^′_k².

Our analysis indicates that there are two qualitatively different regimes with respect to selection on a modifier of recombination, depending on whether recombination rates are high relative to selection [the QLE regime, (S2.4)] or not [the tight linkage regime, (8)]. In the former case, the disequilibrium generated by drift dissipates rapidly, so that disequilibrium can be approximated as the consequence of recent selection (within a time frame proportional to 1/r), which is well captured by the current additive genetic variance in fitness. When recombination rates are weak relative to selection, however, disequilibrium generated while alleles are at low frequency persists as the alleles rise to higher frequency, and the current additive genetic variance in fitness no longer determines the change in frequency of the modifier. Instead, the force on the modifier is strongest whenever the harmonic mean allele frequency at the two loci,H_jk, is small yet the selected alleles spend time at intermediate frequency (with high Varp^′_j, Varp^′_k, and Covp^′_j,p^′_k). We also see a qualitatively different relationship between the recombination rates and the change in the modifier frequency in these two regimes.Equation 8 shows that a modifier of recombination changes in frequency by an amount proportional tor⁻³ for tight linkage, while the QLE analysis (see Equations 7) predicts a change that is proportional tor⁻⁵. Thus, while both analyses predict that selection for recombination is strongest under tight linkage, the selection gradient is overestimated by the QLE analysis (Figure 9).

EXACT EQUATIONS FOR TWO AND THREE LOCI

We use the methods ofBarton and Turelli (1991) to derive exact recursions for two selected loci (labeledj,k), and a modifier of small effect at locusi. We census immediately after juveniles have been randomly sampled from the propagules produced by the previous generation. After the census, haploid populations undergo selection, random mating, and meiosis to produce the next generation of haploid propagules, where we assume that the number of potential propagules is so large that these processes can be treated deterministically. The model also applies to a diploid population, as long as fitness is the product of the fitness contribution of each haplotype. For diploids, the life cycle involves a census of diploid juveniles, selection, meiosis, syngamy of gametes to produce diploid propagules, and random sampling ofN juveniles. Equations A1.2 ofBarton (1995a) differ slightly, because diploid fitness was assumed instead to be the sum of the fitness contribution of each haplotype. We assume that the modifier alleleP_i increases recombination between locij andk by a small amount δr_j,k|i. For ease of analysis, we followBarton (1995a) and set the recombination rate between locij andk tor_jk + 2q_iδr_j,k|i inP_iP_i individuals,r_jk + (q_i −p_i)δr_j,k|i inP_iQ_i individuals, andr_jk − 2p_iδr_j,k|i inQ_iQ_i individuals. (Note that theeffect of the modifier, measured by the difference between genotypes, is not frequency dependent.) Letr_ij andr_ik be the recombination rates between the modifier locus and the selected loci,j andk, respectively, and letr_ijk be the probability of recombination between any of the locii,j,k; these terms are used only to determine the offspring of heterozygousP_iQ_i individuals and so their values in other modifier genotypes are immaterial. Disequilibrium terms involving the modifier locus (D_ij,D_ik,D_ijk) will be dominated by the indirect effects of changing recombination between the selected loci and will be proportional to δr_j,k|i. In the following, we assume that the modifier is weak and keep only linear-order terms in δr_j,k|i.

Let φ_j = 1 +s_j(p_j−q_j)/2 and φ_k = 1 +s_k(p_k −q_k)/2 measure the mean fitness of each locus considered separately. Accounting also for the disequilibrium within a population, the mean fitness is equal to

A1

The three-locus recursions are

A2a

A2b

A2c

A2d

A2e

wherer_ijk = (r_ij +r_ik +r_jk)/2 for any gene order with and without interference.p^*_k andD^*_ik can be obtained from (A2a) and (A2d), respectively, by interchanging subscriptsj andk. These equations assume only that δr_j,k|i is small and are valid for strong selection (s_j,s_k) and strong linkage disequilibrium between the selected loci (D_jk). (A2) can be derived from the three-locus recursion equations presented inFeldman (1972). A supplementary Mathematica package is available that derives (A2) as well as the other main results given inappendixesa andb.

We assume that 2N haploid (orN diploid) juveniles are sampled from the propagules produced by the parental generation according to a multinomial distribution as in the standard Wright-Fisher model (Ewens 1979). The moments of the multinomial distribution can therefore be used to determine the expected values of the perturbations in (A2) and the covariances between them. Although there is no expected change in allele frequencies due to drift,E[ζ_i] =E[ζ_j] =E[ζ_k] = 0, the expected change in linkage disequilibrium isE[ζ_jk] = −D_jk/(2N). Because disequilibrium terms involving the modifier locus are proportional to δr_j,k|i, the expected amount of drift in these terms will be of the order δr_j,k|i/N, which is assumed to be very small and is ignored. Thus, we need to quantify the effects of drift only on the following variances and covariances, which are written to order (1/N):

BecauseE[ζ_a] is either 0 or order 1/N, theE[ζ_aζ_b] required inEquation 3 are, to order 1/N, equal to the cov(ζ_a, ζ_b) given here. Note that the allele frequencies and disequilibrium in these terms are evaluated after the parental generation produces propagules. When selection is weak, however, the allele frequencies and disequilibrium from the previous census can be used to obtain leading-order approximations.

These derivations apply for any amount of disequilibrium, allowing the method to be extended to models where there is a deterministic source of disequilibrium (e.g., epistasis or an initially high level of disequilibrium). In the text andappendix B, we assume that the disequilibria remainO(1/N) throughout the process. In this case,E[ζ_jk] and the covariance terms become negligible. In the numerical iterations of (2) and (3), however, no assumptions are made about the order of the disequilibrium, so that the Mathematica package (available upon request) can be used regardless of the initial level of disequilibrium.

WEAK SELECTION APPROXIMATIONS FOR TWO SELECTED LOCI

When selection is weak (s_j ands_k of orders, wheres ≪ 1), the recursions for the perturbations given by (1) can be further approximated by

B1a

B1b

Throughout thisappendix, values for alleleP_k can be obtained from values given for alleleP_j by interchanging subscriptsj andk. In (B1b), we have included the leading-order terms for each of the perturbations with respect tos even when these areO(s²). These lower-order terms are critical to the initial development of disequilibria when δD_jk is zero.

Taking expectations as in (4) and (5), assuming that the perturbations including the disequilibria areO(N⁻¹), ignoring terms that areO(N⁻²), and using the results ofappendix A for the contributions due to drift gives

B2a

B2b

B2c

B2d

B2e

Supplementary information (S2) describes similar weak selection recursions for the three-locus case.

Equations B2 can be solved explicitly. First, consider loose linkage (r_jk ≫s_j,s_k). In this case, the disequilibrium rapidly approaches a QLE distribution, withEδD²_jk andE[δD_jk] changing slowly relative to the action of recombination. One can therefore setEδD^*_jk toE[δD_jk], etc., in (B2) and solve for the QLE values. For the disequilibrium:

B3a

B3b

B3c, B3c

shows that negative disequilibrium is generated on average in a finite population subject to multiplicative selection. The magnitude of this disequilibrium is quite small with loose linkage but increases rapidly asr_jk becomes small. It also increases with the strength of selection acting on locij andk and with the amount of genetic variability at these loci. Note that the covariance between allele frequencies and linkage disequilibrium (B3b) always contributes more than the variance in linkage disequilibrium to the accumulation of negative linkage disequilibrium by a factor 2(1 −r_jk)/r_jk.

With tight linkage (1/N ≪s ≈r_jk ≪ 1), disequilibrium decays slowly over time, and we can no longer assume quasi-linkage equilibrium. We can, however, move from the discrete recursion equations to continuous-time approximations, which can be more readily solved. We introduce the method in general and then apply the solution to (B2).

In the ODE approach, the discrete recursion equations (B2) all have the form

B4a

whereC is a constant involving recombination rates,F[t] is a function of time involving terms such as −s(p −q), andG[t] is a function of time involving the expectations,E[ ]. When selection and recombination are rare, the change inx[t] over a short time period can be approximated byx[t + Δt] −x[t], wherex[t + Δt] is obtained from (B4a) under the assumption that the amount of selection and recombination in a shorter interval of time, Δt, scales with Δt assΔt andrΔt. Taking the limit of (x[t + Δt] −x[t])/Δt as Δt goes to zero yields an ordinary differential equation:

B4b

wherec,f, andg are used to denoteC,F, andG to linear order ins andr, using a Taylor series that ignores any cross products or higher-order terms. The solution to differential equations of the form (B4b) is

B4c

The integrals involvingf[t] in (B4c) can be solved by noting that, when selection is weak and frequency independent, selection can be defined according to the change in allele frequency that it causes:p′ =dp/dt =spq. Whenf[t] = −s(p −q), for example, applying this definition and integrating yields

Lettingv[t] be the product of the appropriate allele frequencies (depending on which terms are included inf[t]), (B4c) becomes

B4d

The solution (B4d) applies regardless of how selection varies over time. Barton (2000) used a similar method to find the net hitchhiking effect of fluctuating selection on linked neutral loci and gave an approximation valid for tight linkage (r ≪s; his Equation 4).

Following this procedure and assuming that terms of orderr_jks_j are negligible relative tor_jk ands_j, we obtain the following weak selection approximations:

B5a

B5b

B5c

We have focused on the effects of drift on the disequilibrium; analogous equations can be obtained for the effects on allele frequencies. Note that if allele frequencies increase from some low value faster than linkage disequilibrium is dissipated by recombination (s_j,s_k >r_jk), the initial distribution (E₀[ ]) will make a substantial contribution. Even when disequilibrium is initially absent, substantial disequilibrium can build up if the beneficial alleles start at low frequency, because periods with low values ofp_jq_jp_kq_k contribute disproportionately to the integrals in (B5).

1. Altenberg, L., and M. W. Feldman, 1987. Selection, generalized transmission and the evolution of modifier genes. I. The reduction principle. Genetics 117: 559–572. [DOI] [PMC free article] [PubMed] [Google Scholar]
2. Barton, N. H., 1995. a A general model for the evolution of recombination. Genet. Res. 65: 123–144. [DOI] [PubMed] [Google Scholar]
3. Barton, N. H., 1995. b Linkage and the limits to natural selection. Genetics 140: 821–841. [DOI] [PMC free article] [PubMed] [Google Scholar]
4. Barton, N. H., 1998. The effect of hitch-hiking on neutral genealogies. Genet. Res. 72: 123–133. [Google Scholar]
5. Barton, N. H., 2000. Genetic hitchhiking. Philos. Trans. R. Soc. Lond. B 355: 1553–1562. [DOI] [PMC free article] [PubMed] [Google Scholar]
6. Barton, N. H., and M. Turelli, 1991. Natural and sexual selection on many loci. Genetics 127: 229–255. [DOI] [PMC free article] [PubMed] [Google Scholar]
7. Charlesworth, B., 1976. Recombination modification in a fluctuating environment. Genetics 83: 181–195. [DOI] [PMC free article] [PubMed] [Google Scholar]
8. Charlesworth, B., 1990. Mutation-selection balance and the evolutionary advantage of sex and recombination. Genet. Res. 55: 199–221. [DOI] [PubMed] [Google Scholar]
9. Charlesworth, B., 1993. Directional selection and the evolution of sex and recombination. Genet. Res. 61: 205–224. [DOI] [PubMed] [Google Scholar]
10. de Visser, J. A. G. M., R. F. Hoekstra and H. van den Ende, 1996. The effect of sex and deleterious mutations on fitness inChlamydomonas. Proc. R. Soc. Lond. Ser. B 263: 193–200. [Google Scholar]
11. de Visser, J. A. G. M., R. F. Hoekstra and H. van den Ende, 1997. Test of interaction between genetic markers that affect fitness inAspergillus niger. Evolution 51: 1499–1505. [DOI] [PubMed] [Google Scholar]
12. Elena, S. F., and R. E. Lenski, 1997. Test of synergistic interactions among deleterious mutations in bacteria. Nature 390: 395–398. [DOI] [PubMed] [Google Scholar]
13. Ewens, W. J., 1979Mathematical Population Genetics. Springer-Verlag, Berlin.
14. Feldman, M. W., 1972. Selection for linkage modification. I. Random mating populations. Theor. Popul. Biol. 21: 430–439. [DOI] [PubMed] [Google Scholar]
15. Feldman, M. W., F. B. Christiansen and L. D. Brooks, 1980. Evolution of recombination in a constant environment. Proc. Natl. Acad. Sci. USA 77: 4838–4841. [DOI] [PMC free article] [PubMed] [Google Scholar]
16. Felsenstein, J., 1974. The evolutionary advantage of recombination. Genetics 78: 737–756. [DOI] [PMC free article] [PubMed] [Google Scholar]
17. Felsenstein, J., 1988 Sex and the evolution of recombination, pp. 74–86 inThe Evolution of Sex, edited by R. E. Michod and B. R. Levin. Sinauer Press, Sunderland, MA.
18. Felsenstein, J., and S. Yokoyama, 1976. The evolutionary advantage of recombination. II. Individual selection for recombination. Genetics 83: 845–859. [DOI] [PMC free article] [PubMed] [Google Scholar]
19. Fisher, R. A., 1930The Genetical Theory of Natural Selection. Oxford University Press, Oxford.
20. Hill, W. G., and A. Robertson, 1966. The effect of linkage on the limits to artificial selection. Genet. Res. 8: 269–294. [PubMed] [Google Scholar]
21. Iles, M. M., K. Walters and C. Cannings, 2003. Recombination can evolve in large finite populations given selection on sufficient loci. Genetics 165: 333–337. [DOI] [PMC free article] [PubMed] [Google Scholar]
22. Keightley, P. D., 1996. Nature of deleterious mutation load inDrosophila. Genetics 144: 1993–1999. [DOI] [PMC free article] [PubMed] [Google Scholar]
23. Keightley, P. D., and A. Caballero, 1997. Genomic mutation rates for lifetime reproductive output and lifespan inCaenorhabditis elegans. Proc. Natl. Acad. Sci. USA 94: 3823–3827. [DOI] [PMC free article] [PubMed] [Google Scholar]
24. Keightley, P. D., and A. Eyre-Walker, 2000. Deleterious mutations and the evolution of sex. Science 290: 331–333. [DOI] [PubMed] [Google Scholar]
25. Kimura, M., 1965. Attainment of quasi-linkage equilibrium when gene frequencies are changing by natural selection. Genetics 52: 875–890. [DOI] [PMC free article] [PubMed] [Google Scholar]
26. Kimura, M., and T. Maruyama, 1966. The mutational load with epistatic gene interactions in fitness. Genetics 54: 1303–1312. [DOI] [PMC free article] [PubMed] [Google Scholar]
27. Kondrashov, A. S., 1982. Selection against harmful mutations in large sexual and asexual populations. Genet. Res. 40: 325–332. [DOI] [PubMed] [Google Scholar]
28. Kondrashov, A. S., 1984. Deleterious mutations as an evolutionary factor. I. The advantage of recombination. Genet. Res. 44: 199–217. [DOI] [PubMed] [Google Scholar]
29. Kondrashov, A. S., 1988. Deleterious mutations and the evolution of sexual reproduction. Nature 336: 435–441. [DOI] [PubMed] [Google Scholar]
30. Kondrashov, A. S., 1993. Classification of hypotheses on the advantage of amphimixis. J. Hered. 84: 372–387. [DOI] [PubMed] [Google Scholar]
31. Korol, A. B., and K. G. Iliadi, 1994. Recombination increase resulting from directional selection for geotaxis in Drosophila. Heredity 72: 64–68. [DOI] [PubMed] [Google Scholar]
32. Latter, B. D. H., 1965. The response to artificial selection due to autosomal genes of large effect. II. The effects of linkage on limits to selection in finite populations. Aust. J. Biol. Sci. 18: 1009–1023. [DOI] [PubMed] [Google Scholar]
33. Martin, G., S. P. Otto and T. Lenormand, 2005. Selection for recombination in structured populations. Genetics (in press). [DOI] [PMC free article] [PubMed]
34. Maynard Smith, J., 1968. Evolution in sexual and asexual populations. Am. Nat. 102: 469–473. [Google Scholar]
35. Maynard Smith, J., 1978The Evolution of Sex. Cambridge University Press, Cambridge, UK.
36. Maynard Smith, J., 1980. Selection for recombination in a polygenic model. Genet. Res. 35: 269–277. [DOI] [PubMed] [Google Scholar]
37. Maynard Smith, J., 1988. Selection for recombination in a polygenic model: the mechanism. Genet. Res. 51: 59–63. [DOI] [PubMed] [Google Scholar]
38. Morgan, T. H., 1913Heredity and Sex. Columbia University Press, New York.
39. Mukai, T., 1969. The genetic structure of natural populations ofDrosophila melanogaster. VII. Synergistic interactions of spontaneous mutant polygenes controlling viability. Genetics 61: 749–761. [DOI] [PMC free article] [PubMed] [Google Scholar]
40. Muller, H. J., 1932. Some genetic aspects of sex. Am. Nat. 66: 118–138. [Google Scholar]
41. Nagylaki, T., 1993. The evolution of multilocus systems under weak selection. Genetics 134: 627–647. [DOI] [PMC free article] [PubMed] [Google Scholar]
42. Otto, S. P., and N. H. Barton, 1997. The evolution of recombination: removing the limits to natural selection. Genetics 147: 879–906. [DOI] [PMC free article] [PubMed] [Google Scholar]
43. Otto, S. P., and N. H. Barton, 2001. Selection for recombination in small populations. Evolution 55: 1921–1931. [DOI] [PubMed] [Google Scholar]
44. Otto, S. P., and M. W. Feldman, 1997. Deleterious mutations, variable epistatic interactions, and the evolution of recombination. Theor. Popul. Biol. 51: 134–147. [DOI] [PubMed] [Google Scholar]
45. Otto, S. P., and T. Lenormand, 2002. Resolving the paradox of sex and recombination. Nat. Rev. Genet. 2: 252–261. [DOI] [PubMed] [Google Scholar]
46. Qureshi, A. W., 1963A Monte Carlo Evaluation of the Role of Finite Population Size and Linkage in Response to Continuous Mass Selection. Statistical Laboratory, Iowa State University, Ames, IA.
47. Rice, W. R., 2002. Experimental tests of the adaptive significance of sexual recombination. Nat. Rev. Genet. 3: 241–251. [DOI] [PubMed] [Google Scholar]
48. Seager, R. D., and F. J. Ayala, 1982. Chromosome interactions inDrosophila melanogaster. I. Viability studies. Genetics 102: 467–483. [DOI] [PMC free article] [PubMed] [Google Scholar]
49. Seager, R. D., F. J. Ayala and R. W. Marks, 1982. Chromosome interactions inDrosophila melanogaster. II. Total fitness. Genetics 102: 485–502. [DOI] [PMC free article] [PubMed] [Google Scholar]
50. Stephan, W., T. H. E. Wiehe and M. W. Lenz, 1992. The effect of strongly selected substitutions on neutral polymorphism: analytical results based on diffusion theory. Theor. Popul. Biol. 41: 237–254. [Google Scholar]
51. Szathmary, E., 1993. Do deleterious mutations act synergistically? Metabolic control theory provides a partial answer. Genetics 133: 127–132. [DOI] [PMC free article] [PubMed] [Google Scholar]
52. Turner, J. R. G., and M. H. Williamson, 1968. Population size, natural selection, and the genetic load. Nature 218: 700. [DOI] [PubMed] [Google Scholar]
53. West, S. A., A. D. Peters and N. H. Barton, 1998. Testing for epistasis between deleterious mutations. Genetics 149: 435–444. [DOI] [PMC free article] [PubMed] [Google Scholar]
54. Wolfram, S., 1991Mathematica. Addison Wesley, New York.