# Fixation Probability in a Two-Locus Model by the Ancestral Recombination–Selection Graph

- Sabin Lessard
^{1}and - Amir R. Kermany

- Département de Mathématiques et de Statistique, Université de Montréal, Montréal, Québec H3C 3J7, Canada

- 1Corresponding author: Département de Mathématiques et de Statistique, Université de Montréal, C.P. 6128, Succursale Centre-ville, Montréal, Québec H3C 3J7, Canada. E-mail: lessards{at}dms.umontreal.ca

## Abstract

We use the ancestral influence graph (AIG) for a two-locus, two-allele selection model in the limit of a large population size to obtain an analytic approximation for the probability of ultimate fixation of a single mutant allele *A*. We assume that this new mutant is introduced at a given locus into a finite population in which a previous mutant allele *B* is already segregating with a wild type at another linked locus. We deduce that the fixation probability increases as the recombination rate increases if allele *A* is either in positive epistatic interaction with *B* and allele *B* is beneficial or in no epistatic interaction with *B* and then allele *A* itself is beneficial. This holds at least as long as the recombination fraction and the selection intensity are small enough and the population size is large enough. In particular this confirms the Hill–Robertson effect, which predicts that recombination renders more likely the ultimate fixation of beneficial mutants at different loci in a population in the presence of random genetic drift even in the absence of epistasis. More importantly, we show that this is true from weak negative epistasis to positive epistasis, at least under weak selection. In the case of deleterious mutants, the fixation probability decreases as the recombination rate increases. This supports Muller’s ratchet mechanism to explain the accumulation of deleterious mutants in a population lacking recombination.

THE Hill–Robertson (HR) effect (Hill and Robertson 1966) is often mentioned as one of the main arguments in favor of the evolution of recombination. In short, it predicts that beneficial mutant alleles arising at different loci in a finite population are more likely to fix in the population as the recombination rate increases even when selection acts independently upon the loci.

Since the early works of Fisher (1930) and Muller (1932), it is generally believed that an evolutionary advantage of recombination is to bring together beneficial mutant alleles arising at different loci. Accordingly the effect of recombination should be to increase the rate of evolution of the population (Crow and Kimura 1965). However, it has been shown that recombination has no effect on this rate in an infinite population if there is initial linkage equilibrium and absence of epistasis so that linkage equilibrium is maintained thereafter in the population (Felsenstein 1965; Maynard Smith 1968).

If recombination can have an effect on the rate of evolution only by breaking down linkage disequilibrium in absolute value, then the effect should be to increase this rate only when linkage disequilibrium in the population is negative (NLD). In the case of a two-locus model, this happens when the frequency of the double mutant is strictly smaller than the product of the frequencies of the mutant alleles. This situation is arguably likely to happen in the view that beneficial mutations are very rare (Crow and Kimura 1969).

On the other hand, NLD could be produced by negative epistasis (NE), with the double mutant being less fit than what it would be under independent effects of the mutant alleles. Then the double mutant would die more often than is expected with mutant alleles acting independently, leaving NLD in the population. In the opposite case of positive epistasis (PE), Eshel and Feldman (1970) showed that the frequency of the double mutant in an infinite population is always larger in an asexual population than in a population with recombination. This suggests that NE rather than PE could be advantageous for the evolution of recombination.

In taking account of a finite population size, Bodmer (1970) considered the expected time until the first formation of a double mutant from two initial single mutants. He concluded that recombination would have a greater advantage in a small population than in a large one. Karlin (1973) showed that this expected time, without selection effects, was indeed a decreasing function of the recombination rate *r*. But he showed also that the expected time until the total fixation of the double mutant was an increasing function of *r*. In other words, increasing recombination might be advantageous in speeding the time until the first formation of the double mutant, but disadvantageous by breaking apart the favored gamete type once formed.

Summing first-order terms for expected changes in gene frequencies in a large finite population under weak selection with additive gene action, which corresponds to an absence of epistasis (AE), Hill and Robertson (1966) deduced that the probability of fixation of an allele *A* initially in NLD with a beneficial allele *B* at another tightly linked locus increases with the recombination fraction between the two loci. The effects of linkage disequilibrium and epistasis on the probability of fixation of gametes and alleles in a finite population under the assumption of weak selection were further studied by diffusion approximations in the limit of a large population size (Ohta 1968).

In the case of initial linkage equilibrium (LE) in a finite population, first-order approximations fail to detect the effect of linkage on the fixation probability. Higher-order effects in the absence of epistasis were first exhibited by simulations (Hill and Robertson 1966).

In a finite population initially in LE, genetic drift creates random instances of linkage disequilibrium. Although random drift can generate both positive and negative disequilibria without any *a priori* bias on average, selection dispels positive linkage disequilibrium (PLD) more efficiently than NLD even in the absence of epistasis, so that the average linkage disequilibrium becomes negative. As shown by simulations and some analytical arguments (Hill and Robertson 1966; see also Barton and Otto 2005), this leads to an average accumulation of NLD. As a consequence, responses to selection at different loci are expected to interfere with each other, even in the absence of gene interaction. This is known as the HR effect. It is by reducing the interference caused by the randomly generated linkage disequilibrium that an increase in the recombination rate raises the rate of fixation of favorable mutants.

The relationship between the HR effect and the Fisher–Muller theory for the evolutionary advantage of recombination was pointed out by Felsenstein (1974). Moreover, it was noted that Muller’s ratchet mechanism (Muller 1964) for the accumulation of deleterious mutants in the absence of recombination, which is formally equivalent to the accumulation of advantageous mutants in the presence of recombination, can be explained by the HR effect.

The theory of evolution at a selectively neutral modifier locus that controls the recombination fraction between two major loci that are under selection in an infinite population was developed by Feldman *et al.* (1980). If the major loci are in linkage disequilibrium at a balance between selection against deleterious alleles and mutation toward them, then a mutation increasing recombination succeeds if the linkage disequilibrium is negative, which occurs when epistasis is negative, and the modifier locus is sufficiently tightly linked to the major loci. If the modifier locus is loosely linked, NE has to be weak enough (Otto and Feldman 1997). A similar conclusion has been reached for sweeps of beneficial alleles (Barton 1995a). However, including spatial heterogeneity extends the range of epistasis over which recombination can be favored, from strong NE to PE depending on environmental circumstances (Lenormand and Otto 2000).

On the other hand, Feldman *et al.* (1980) also showed that, if the major loci are at a stable equilibrium in linkage disequilibrium under selection and recombination, then a mutation at the modifier locus increases in frequency when rare if and only if it decreases the recombination fraction. This is part of a general reduction principle for genetic modifiers in an infinite population in a constant environment (Feldman and Liberman 1986).

It has been argued that a modifier allele that increases the recombination rate would be promoted in a finite population due to its role in reducing the negative effect of poor genetic backgrounds on the probability of fixation of favorable mutants, at least in the absence of epistasis. This has been shown by applying a branching process to mutant lines in an infinite population with deterministic changes in the frequencies of the genetic backgrounds (Barton 1995b; Otto and Barton 1997). The same approach has been used to study the probability that both beneficial mutants fix. The analysis of this probability has been refined to deal with the troublesome case where the second mutant is more beneficial than the first (Yu and Etheridge 2010).

Simulations have indicated that this is true across a broad range of epistatic interactions, from weak negative epistasis to positive epistasis, provided that the population size is small enough (Otto and Barton 2001). This suggests that the HR effect overwhelms the influence of epistasis on LD over a wide range of epistasis values.

More recently, a perturbation method to track fluctuations in linkage disequilibrium during the spread of beneficial alleles and to measure the impact on a modifier allele of recombination has been proposed (Barton and Otto 2005). The method consists of considering only the first and second moments of random sampling effects on the deterministic dynamics for the allele frequencies and linkage disequilibrium in an infinite population.

A different perturbation technique to approximate the probability of ultimate fixation of an allele in a multilocus setting assumes small selection effects at different loci in a population of fixed finite size (Lehman and Rousset 2009). This is an extension of a direct Markov chain approach for one-locus models based on expected changes in allele frequencies in one time step or one generation (Rousset 2003; Lessard and Ladret 2007; Lessard and Lahaie 2009). Then the first-order effect of selection can be expressed in terms of expected times that lineages of sampled genes take to merge backward in time, under neutrality. The calculation of these times for one-locus models in the limit of a large population size makes use of the coalescent (Kingman 1982) and its extension to incorporate multiple mergers in the case of highly skewed reproduction schemes (Pitman 1999; Sagitov 1999; Möhle and Sagitov 2001).

In the case of multilocus selection models with a Wright–Fisher reproduction scheme allowing for recombination in a population of fixed finite size, Lehman and Rousset (2009) considered Taylor expansions of the fixation probability with respect to the intensity of selection. They deduced exact linear recurrence systems of equations for gamete frequencies in sampled individuals backward in time under neutrality to compute the coefficients. Advanced matrix theory was used to interpret these coefficients in terms of mean sojourn times in the backward neutral process. However, a first-order expansion of the fixation probability with respect to the intensity of selection is not sufficient to detect the HR effect in a two-locus model in the absence of epistasis. Actually, a third-order expansion is necessary. In this case, the coefficients of the approximation become difficult to interpret.

Our objective in this article is to consider an ancestral recombination–selection process to deduce an analytic approximation for the probability of ultimate fixation of an allele in a finite but large population under weak selection and tight linkage. The allele is assumed to be a mutant type *A* introduced at a given locus into the population in which a previous mutant type *B* is already segregating with a wild type at another linked locus. Exact conditions for a small increase in the recombination rate to increase the probability of ultimate fixation of a single *A* are addressed.

We focus on a discrete-time two-locus selection model with a Moran reproduction scheme (Moran 1958). We consider the ancestral recombination–selection graph for sampled gametes in the limit of a large population size, which is known as the ancestral influence graph (AIG) (Donnelly and Kurtz 1999). The AIG provides a supragenealogy for a sample of individuals at linked, nonneutral loci in a limiting Fleming–Viot measure-valued diffusion process with selection and recombination. It is a supragenealogy in the sense that the true genealogy of the sample is embedded into it. It combines the ancestral recombination graph (ARG) (Griffiths and Marjoram 1996, 1997) and the ancestral selection graph (ASG) (Krone and Neuhauser 1997; Neuhauser and Krone 1997), extending the coalescent (Kingman 1982) to include both recombination and selection. The ARG and ASG, given the sample composition, have been widely used in likelihood methods to estimate the recombination rate or detect recombination hotspots (*e.g.*, McVean *et al.* 2002; Stephens and Donnelly 2003; Fearnhead *et al.* 2004; Wakeley and Sargsyan 2009) and to locate disease genes from marker loci (*e.g.*, Hudson and Kaplan 1988; Fearnhead 2003; Larribe and Lessard 2008; Larribe and Fearnhead 2011).

We make use of a discrete-time Moran model for mortality selection determined at two loci in a finite haploid population to ascertain the analysis. After recalling the definitions and assumptions, the probability of ultimate fixation of an allele is expressed in terms of sums of expected sample frequencies, which correspond to expected times with given ordered random samples. Then the ancestral graphs obtained by tracing the supragenealogy of an ordered random sample through coalescence, recombination, or selection events backward in time, whose limit as the population size increases is an AIG, are described. These graphs are used to express the expected times with ordered random samples of given types. It is shown that an approximation of any order of the fixation probability with respect to the population-scaled recombination and selection parameters in the limit of a large population size can be obtained by considering ancestral graphs with enough recombination or selection events. Finally this is applied to directional selection with either beneficial mutants or deleterious mutants, in epistatic interaction or in the absence of interaction, by considering one recombination event and one or two selection events to detect the effect of recombination.

It is expected that the results are valid in the domain of attraction of the Fleming–Viot process with recombination and selection in the same way that a wide class of Cannings exchangeable models including the Moran model and the Wright–Fisher model fall in the domain of attraction of the Kingman coalescent (Möhle and Sagitov 2001).

## Definitions and Model

Suppose a population of finite size *N* distributed over *N* distinct sites, so that each site is occupied by one and only one individual. Each individual is one of four types, *AB*, *Ab*, *aB*, or *ab*, with respect to two loci with alleles *A*, *a* segregating at locus 1 and *B*, *b* at locus 2.

Reproduction is assumed to follow a discrete-time Moran model. At each time step τ ≥ 0, two individuals are sampled at random in the population and they produce an offspring. Random sampling of the parents is assumed to take place with replacement so that selfing is permitted and then occurs with probability *N*^{−1}.

With respect to the two loci, the offspring produced is either an exact copy of one of its parents, with probability 1 – *r*, or a recombinant, with probability *r*. This probability of recombination is inversely proportional to the population size, so that *r* = ρ*N*^{−1}, where ρ represents a population-scaled recombination fraction. Weak recombination is modeled by keeping ρ constant as *N* → ∞.

On the other hand, one individual is chosen at random to be replaced by the offspring. Replacement actually occurs with some probability that depends on the type of the individual, called its mortality. It is given by 1 – *c _{AB}s*, 1 –

*c*, 1 –

_{Ab}s*c*, or 1 –

_{aB}s*c*for an individual of type

_{ab}s*AB*,

*Ab*,

*aB*, or

*ab*, respectively (see Figure 1). These can be interpreted as probabilities of dying. If replacement does not occur, then the offspring is eliminated and there is no change in the population during the corresponding time step.

Here, the parameters 0 ≤ *c _{AB}*,

*c*,

_{Ab}*c*,

_{aB}*c*≤ 1 represent coefficients of selection with respect to an intensity of selection 0 <

_{ab}*s*< 1. They can be viewed as viability parameters. Neutrality corresponds to

*s*= 0.

The intensity of selection is expressed in the form *s* = σ*N*^{−1}, where *N* is the population size. The parameter σ stands for a population-scaled intensity of selection. Weak selection is modeled by keeping σ constant as *N* → ∞.

Alleles *A* and *B* are mutant types, while alleles *a* and *b* are wild types. The mutant alleles *A* and *B* are advantageous when each one reduces the mortality of its carrier compared to what it would be without these alleles. This is the case if the coefficients of selection satisfy the inequalities*A* and *B* are deleterious.

Allele *B* is a mutant that was introduced some time ago at locus 2 into a population entirely composed of *ab* individuals and its frequency has reached some value 0 < *x* < 1. Then a single mutant allele *A* is introduced at random at locus 1 into the population, so that it is linked to *B* with probability *x* and to *b* with the complementary probability 1 – *x*. In both cases its frequency is given by the inverse of the population size, that is, *N*^{−1}. In the former case, the frequency of *aB* is reduced to *x* – *N*^{−1} and in the latter the frequency of *ab* is reduced to 1 – *x* – *N*^{−1}.

Linkage disequilibrium (LD) is measured by the difference between the frequency of the double mutant, *AB*, and the product of the frequencies of the mutant alleles, *A* and *B*, which is represented by *D*. Alternatively, *D* is equal to the difference between the product of the frequencies of *AB* and *ab* and the product of the frequencies of *Ab* and *aB*. In the present case, linkage disequilibrium following the introduction of a single *A* is initially positive (PLD) with probability *x* and given by *N*^{−1}(1 – *x*), while it is initially negative (NLD) with probability 1 – *x* and given by –*N*^{−1}*x*. This yields an average LD given by

Epistasis refers to the phenomenon in which the effect of a mutant at one locus, here *B*, is masked or enhanced by a mutant at another locus, here *A*. Population geneticists extended the concept to mean nonindependent or multiplicative effects of mutants.

Epistasis is positive (PE) if interactions between *A* and *B* are such that the double mutant is more fit in comparison to the wild gamete type than what it would be if the mutant alleles have independent effects on fitness. In terms of mortality parameters, this means the inequality

Where advantageous mutations are concerned, PE enhances the fitness increase predicted from individual mutational effects, whereas NE lessens it. It is the opposite for deleterious mutations with respect to fitness decrease.

Note that, in the limit of weak selection when *s* = σ*N*^{−1} → 0 as *N* → ∞, epistasis is positive, negative, or null if

## Expected Change in Allele Frequency

Let **x**(τ) = (*x _{AB}*(τ),

*x*(τ),

_{Ab}*x*(τ),

_{aB}*x*(τ)) be the vector of the individual type frequencies at the current time step τ ≥ 0. Then the frequency of

_{ab}*AB*at the next time step will increase by

*N*

^{−1}with probability

*N*

^{−1}with probability

*(τ) – μ*

_{AB}*(τ). Here we use the notation*

_{AB}*A*and similarly

*B*. Moreover,

Therefore, the change in the frequency of *AB* from time step τ to time step τ + 1, given by Δ*x _{AB}*(τ) =

*x*(τ + 1) –

_{AB}*x*(τ), is found to have

_{AB}*Ab*. Hence the change in the frequency of allele

*A*, which can be expressed as Δ

*x*(τ) = Δ

_{A}*x*(τ) + Δ

_{AB}*x*(τ), has conditional expectation

_{Ab}*A*at time step τ. Straightforward algebraic manipulations lead to the following conclusion.

**Proposition 1*** For the discrete-time Moran model with recombination and selection described in* *Figure 1*, *the conditional expected change in the frequency of A is given by**where N is the population size and* σ *= sN is a population-scaled intensity of selection with coefficients* 0 ≤ *c _{AB}*,

*c*,

_{Ab}*c*,

_{aB}*c*≤ 1

_{ab}*for the individual types AB*,

*Ab*,

*aB*,

*and ab*,

*respectively*.

## Probability of Fixation of an Allele at One Locus

The random process **x**(τ) = (*x _{AB}*(τ),

*x*(τ),

_{Ab}*x*(τ),

_{aB}*x*(τ)) for τ ≥ 0 is a Markov chain on a finite state space

_{ab}*S*. This is the set of all four-dimensional frequency vectors whose entries are multiples of

*N*

^{−1}.

There are four absorbing states represented by*AB*, *Ab*, *aB*, and *ab*, respectively. All other states are transient.

In virtue of the ergodic theorem for Markov chains (see, *e.g.*, Karlin and Taylor 1975; Grimmett and Stirzaker 1982), the probability of transition from state **x** to state **y** in *k* time steps, namely*P*** _{xy}**(∞). This probability is 0 unless

**y**is an absorbing state. Therefore,

**y**= (

*y*,

_{AB}*y*,

_{Ab}*y*,

_{aB}*y*) in

_{ab}*S*, converges in the same limit to

*A*given an initial population state

**x**.

On the other hand, we have*k* → ∞, this leads to*A* leads to the following result.

**Proposition 2*** For the discrete-time Moran model with recombination and selection of Proposition 1*, *the probability of ultimate fixation of A is given by**where**for z*_{1} *= AB*, *Ab and z*_{2} *= aB*, *ab*.

The quantity *N*^{2}/2 time steps and over all time steps that two individuals chosen at random *with replacement* in the population at the same time step τ ≥ 0 will be of types *z*_{1} and *z*_{2} in this order.

## Ancestral Recombination–Selection Graph

An ancestral recombination–selection graph is a Markov chain on ordered samples obtained by tracing backward in time the ancestors, real or virtual, of a given number of individuals chosen at random without replacement in the population at a given time step. It is characterized by a sequence of changes in the ancestry of the sample and times between these events.

As in Krone and Neuhauser (1997), this process is considered in the framework of a Moran model, but in discrete time and with recombination allowed, so that a change in the ancestry can involve simultaneous events of coalescence, recombination, or selection. In the limit of a large population size, however, with time and parameters for recombination and selection appropriately scaled, only one event of coalescence, recombination, or selection can occur at a time with probability one. The limiting process corresponds to the AIG introduced by Donnelly and Kurtz (1999), as described in Fearnhead (2003). An exact description of the ancestral graph incorporating recombination and selection in a discrete-time Moran model could not be found in the literature, although it might exist. Such a description is actually necessary to establish rigorous approximation results for the probability of fixation in the presence of recombination and selection.

Consider the model of the previous section with 1 – *c _{AB}s*, 1 –

*c*, 1 –

_{Ab}s*c*, and 1 –

_{aB}s*c*as mortalities associated to the individual types

_{ab}s*AB*,

*Ab*,

*aB*, and

*ab*, respectively, under the conditions 0 ≤

*c*,

_{AB}*c*,

_{Ab}*c*,

_{aB}*c*≤ 1 and 0 <

_{ab}*s*< 1.

The replacement rule for an individual chosen at random can be described as follows. Replacement is inevitable irrespective of the type of the individual with probability 1 – *s*, which corresponds to the lowest possible mortality. On the other hand, replacement is type specific with probability *s*. In this case, replacement occurs with conditional probability*Ab*, for instance, replacement will occur with probability 1 – *s* + *s*(1 – *c _{Ab}*), which is the same as 1 –

*c*. With the complementary probability, there is no replacement.

_{Ab}sA type-specific replacement is considered to be a selection event. Its probability in one time step is expressed in the form *s* = σ*N*^{−1}. Recall that the probability of a recombination event in one time step is expressed in a similar form, namely *r* = *ρN*^{−1}.

The scaling used for the probabilities of selection or recombination events, along with *N*^{2}/2 time steps as unit of time, will simplify the ancestral process in the limit of a large population size. This timescale is standard for a discrete-time Moran model (see, *e.g.*, Ewens 1990).

Consider a sample of *n* distinct individuals in the population at a given time step and label them arbitrarily with the integers *i* = 1, … , *n*. Label arbitrarily the other *N* – *n* individuals in the population at the same time step with the integers *i* = *n* + 1, … , *N*.

Following the lineages of the sampled individuals in one time step back, there will be pure coalescence of *i* and *j*, for *i*, *j* = 1, … , *n* with *i* ≠ *j*, if the offspring produced was an exact copy of *j* [probability *N*^{−1}(1 – *ρN*^{−1})] and the individual replaced irrespective of its type was the individual that occupied the site of *i* [probability *N*^{−1}(1 – σ*N*^{−1})] or vice versa. We conclude that*n* in one time step back. Then the sample size is reduced by one by merging the lineages of two sampled individuals.

On the other hand, there will be pure recombination of *i* in one time step back, for *i* = 1, … , *n*, if the offspring produced was a recombinant of *k* and *l* not in the sample and different from each other, that is, for *k*, *l* = *n* + 1, … , *N* with *k* ≠ *l* [probability ρ*N*^{−1}(*N* – *n*)(*N* – *n* – 1)*N*^{−2}], and the individual replaced irrespective of its type was the individual that occupied the site of *i* [probability *N*^{–1}(1 – *σN*^{−1})]. Therefore, we find that*n* in one time step back. In this case the sample size is increased by one by splitting the lineage of one sampled individual into two, each one being actually ancestral to the sampled individual at only one of the loci.

Finally there will be pure selection of *i* in one time step back, for *i* = 1, … , *n*, if the offspring produced was an exact copy of *k* not in the sample, that is, for *k* = *n* + 1, … , *N* [probability (*N* – *n*)*N*^{−1}(1 – ρ*N*^{−1})], and the individual chosen to be replaced according to its type is the one that occupied the site of *i* (probability σ*N*^{−2}). We conclude that*n* in one time step back.

In the case of a pure selection event, the sample size is increased by one by branching the lineage of one sampled individual into two, each one being potentially ancestral to the sampled individual at both loci. The incoming lineage is the lineage of the offspring produced one time step back, while the continuing lineage is the lineage of the individual chosen to be replaced by the offspring. One of these lineages is real and the other virtual, but both lineages must be traced back until ancestors of known types are reached. Then the conditional probability of replacement can be determined.

Note that the probabilities of pure coalescence, recombination, or selection events in one time step back for a sample of fixed size *n* are all functions of order *N*^{−2}, denoted by *O*(*N*^{−2}). On the other hand, the probabilities of multiple events involving simultaneous coalescence, recombination, or selection events that would affect the lineages of the sampled individuals in one time step back are all functions of order *O*(*N*^{−3}). In all cases the sample size can decrease by at most one, when a pure coalescence event occurs, and increase by at most two, when a selection event and a recombination event occur simultaneously but without any coalescence event occurring.

Given a sample of size *n*, the total number of pure coalescence events to consider is *n*(*n* – 1)/2, while this number is *n* for pure selection events and for pure recombination events. Therefore, the total probability of change in one time step back for the whole sample is given by*N*^{2}/2 time steps as unit of time. Moreover, given a change in one time step back, the conditional probability of each pure coalescence, recombination, or selection event is

or

Let the time back, in number of time steps, for a sample of size *n* to be affected by any coalescence, recombination, or selection event be represented by τ* _{n}*. This sojourn time is a geometric random variable independent of all previous transition events and sojourn times, whose expected value is given by

*N*

^{2}/2 time steps, namely

*in the limit of a large population size. As a matter of fact,*

_{n}*t*> 0. Moreover,

*N*large enough. Therefore,

Let us summarize.

**Proposition 3*** Consider the discrete-time Moran model of Proposition 1 with population-scaled recombination fraction* ρ *= rN and population-scaled intensity of selection* σ *= sN in the case of coefficients of selection* 0 ≤ *c _{AB}*,

*c*,

_{Ab}*c*,

_{aB}*c*≤ 1

_{ab}*. In addition*,

*consider an ordered sample without replacement of size n in a population of large size N. Backward in time*,

*each pair of lineages merges as a result of a coalescence event with approximate probability*2

*N*

^{–2},

*while each lineage splits into two as a result of a recombination event with approximate probability*ρ

*N*

^{–2}

*or branches into two as a result of a selection event with approximate probability*σ

*N*

^{–2},

*for an approximate total probability of change*2λ

_{n}N^{–2}

*= n*(

*n*– 1 + ρ + σ)

*N*

^{–2}

*. In number of N*

^{2}/2

*time steps in the limit of a large population size*,

*the expected time for a change is*

*Moreover*,

*in the case of a change caused by a*

*selection event*,

*the incoming lineage is real with conditional probability*1 –

*c*, 1 –

_{AB}*c*, 1 –

_{Ab}*c*,

_{aB}*or*1 –

*c*,

_{ab}if the type of the individual on the continuing lineage is AB*Ab*,

*aB*,

*or*

*ab*,

*respectively*.

## Calculation for Fixation Probability

An ordered sample of *n* individuals is represented by an *n*-dimensional vector **z** = (*z*_{1}, … , *z _{n}*), where

*z*=

_{i}*AB*,

*Ab*,

*aB*, or

*ab*, for

*i*= 1, … ,

*n*. The sample configuration is given by the vector

**n**= (

*n*,

_{AB}*n*,

_{Ab}*n*,

_{aB}*n*) with

_{ab}*n*+

_{AB}*n*+

_{Ab}*n*+

_{aB}*n*=

_{ab}*n*.

Let **z**(τ) be an ordered sample of *n* individuals chosen at random *without replacement* at time step τ ≥ 0. The probability distribution of this sample will depend on the ancestral recombination-selection graph from time step τ to time step 0, represented by *G*(τ), and the type frequencies at time step 0, given by **x**(0). What will matter most is the supragenealogy from time step τ to time step 0. It is represented by sequence events backward in time written in the form*m*_{τ} is the total number of events. These are events of coalescence, recombination, or selection in one step back, either pure or multiple.

Let *n _{G}* be the number of ancestors, real or virtual, after the occurrence of the last event of

*G*backward in time. This last event is assumed to take place at time back τ

*. On the other hand, the time with*

_{G}*n*ancestors is represented by

_{G}*G*to be an admissible supragenealogy from time step τ to time step 0, it is necessary that

*and*

_{G}*is a sum of independent geometric random variables.*

_{G}The probability of the event **z**(τ) = **z**, given **x**(0), can be expressed in the form*P*(*G _{k}*) defined by (29)–(31) for pure events of coalescence, recombination, or selection and by (32) for multiple events. Moreover, we have

**Proposition 4*** Let* **z**(τ) *be an ordered sample of n individuals chosen at random without replacement at time step* τ ≥ 0 *and* **x**(0) *be the vector of the individual type frequencies at time step* 0*. Then we have**where G is a sequence of pure or multiple events of coalescence*, *recombination*, *or selection from time step* τ *to time step* 0, *n _{G} is the number of ancestors at time step* 0,

*and*

*is a time back with this number of ancestors*.

Actually the conditional probability of **z**, given *G* and **x**(0), in Proposition 4 depends on the types of the *n _{G}* ordered ancestors at time step 0, represented by

**z**(0), so that

**z**is incompatible with

*G*and

**z**(0). Otherwise, this conditional probability is 1 times a product of conditional probabilities of replacement, which is different from 1 only in the case of selection events in

*G*. On the other hand,

*N*)

*=*

_{n}*N*× (

*N*− 1) × … × (

*N*−

*n*+ 1) denotes a falling factorial, while

**z**(0), respectively. Moreover, this sample satisfies

**z**(0) to be compatible with

*G*and

**x**(0).

Consider, for instance, a sequence of events backward in time, *G* = (*R*_{2}, *S*_{3}, *C*_{4}) for an ordered sample of size *n* = 2, as illustrated in Figure 2. Here we have a pure recombination event, a pure selection event, and a pure coalescence event, in this order backward in time. In the case of recombination, one lineage splits into two, a left lineage assumed to be ancestral at locus 1 and a right lineage assumed to be ancestral at locus 2. In the case of selection, one lineage branches into two, a continuing lineage and a new incoming lineage, both potentially ancestral. And last, in the case of coalescence, two lineages merge. The probability of the whole sequence of events is*n _{G}* = 3. The time back to the last event,

*τ*, can be expressed in the form

_{G}_{2}, τ

_{3}, and τ

_{4}are independent geometric random variables with parameters

*p*

_{2},

*p*

_{3}, and

*p*

_{4}, respectively. On the other hand, the time back spent with

*n*ancestors,

_{G}*p*

_{3}.

Finally, given that *G* = (*R*_{2}, *S*_{3}, *C*_{4}) and *τ _{G}* ≤ τ <

*τ*+

_{G}*τ*, the ordered sample

_{nG}**z**(τ) = (

*AB*,

*ab*) occurs with probability 1 –

*c*, if the ancestral state at time step 0 is

_{aB}**z**(0) = (

*Ab*,

*aB*,

*ab*). The probability of this initial ancestral state given the initial type frequencies is

*G*and then all possible

*G*for this particular ordered sample.

## Approximation Results

We are now ready to approximate the probability of ultimate fixation of *A* under the assumptions that the population size is large and the population-scaled recombination and selection parameters are small.

Consider **z** = (*z*_{1}, *z*_{2}), where *z*_{1} = *AB* or *Ab* and *z*_{2} = *aB* or *ab*. First note that

Let |*G*| denote the minimum number of ancestors along a sequence of events *G* for the ordered sample **z** = (*z*_{1}, *z*_{2}). Note that*G*| = 1, since alleles *A* in *z*_{1} and *a* in *z*_{2} cannot have the same ancestor, while*G*| ≥ 2, since allele *A* in *z*_{1} must be present in at least one ancestor. This leads to the inequality*W*_{2} is the time back in number of *N*^{2}/2 time steps for the number of ancestors in the ancestral graph starting from a sample of size 2 to reach one for the first time. This occurs when the most recent ultimate ancestor (MRUA) is found. It can be shown that *E*(*W*_{2}) is finite and bounded by a constant that does not depend on *N*. This is also true for *E*(*W _{n}*) ≥

*E*(

*W*

_{n}_{–1}), which is defined analogously for a sample of any size

*n*≥ 3. (See

*Appendix*.)

Now, suppose that the population-scaled recombination and selection parameters, ρ and σ, are small and of the same order of magnitude, so that ρ *= d*σ << 1 for some constant *d* > 0. Let *G* for the ordered sample **z** = (*z*_{1}, *z*_{2}). If *k*, then*G*^{(1)} is a sequence of events such that *G*^{(2)} just after the first increase in the number of ancestors that brings this number above *k*. Note that there is a finite number of *G*^{(1)} satisfying these conditions and that the probability of each one can be neglected compared to σ* ^{k}*; that is,

*P*(

*G*) =

*P*(

*G*

^{(1)})

*P*(

*G*

^{(2)}) and

**Proposition 5*** Ignoring terms of order x*_{A}(0)*O*(σ* ^{l}*ρ

*)*

^{m}*for l + m*≥

*k +*1

*where*ρ

*= rN and*σ

*= sN are the population-scaled parameters for recombination and selection*,

*respectively*,

*the expected times in the probability of ultimate fixation of A given in Proposition*2

*are approximated by*

*for*

**z**

*=*(

*z*

_{1},

*z*

_{2}),

*with z*

_{1}

*= AB*,

*Ab and z*

_{2}

*= aB*,

*ab*,

*where all terms in the summation*,

*given by*(40), (46), and (50),

*are approximated by their leading terms in the case of a large population size. Here*,

*the summation is over all sequences G of pure coalescence*,

*recombination*,

*or selection events backward in time with at most k pure recombination or selection events and a number of ancestors always larger than two with final value n*.

_{G}Note that the coefficient of σ* ^{l}*ρ

*for*

^{m}*l*+

*m*≤

*k*in Proposition 5 is obtained by considering all sequences of events

*G*involving up to

*l*pure selection events and

*m*pure recombination events.

Using MATHEMATICA and (69), a polynomial of degree *k* with respect to σ and ρ approximating the quantity (23) in Proposition 2 for σ and ρ small enough can be calculated. This approach leads to the main results of this article.

**Proposition 6*** Consider the discrete-time Moran model with small population-scaled recombination fraction* ρ *= rN and small population-scaled intensity of selection* σ *= sN with coefficients of selection* 0 ≤ *c _{AB}*,

*c*,

_{Ab}*c*,

_{aB}*c*≤ 1

_{ab}*such that*ε =

*c*≠ 0

_{AB}– c_{Ab}– c_{aB}+ c_{ab}*. Given the initial conditions in*

*Figure*1

*with x*(0) =

_{A}*N*

^{–1}

*and ignoring terms of order*

*N*

^{–2}or

*N*

^{–1}

*O*(σ

*ρ*

^{l}*)*

^{m}*for l + m*≥

*4*,

*the probability of ultimate fixation of A is approximated by*

**Proposition 7**

*Under the conditions of Proposition*6

*but in the case where the coefficients of selection satisfy*ε

*= c*0

_{AB}− c_{Ab}– c_{aB}+ c_{ab}=*and terms of order*

*N*

^{–2}or

*N*

^{−1}

*O*(σ

*ρ*

^{l}*)*

^{m}*for l + m*≥ 5

*are ignored*,

*the probability of ultimate fixation of A is approximated by*

## Discussion

### Effect of selection at linked loci on the probability of fixation of a single mutant allele

The conditional expected change in the frequency of an allele *A* in a two-locus, two-allele Moran model from one time step to the next as expressed in (15) can be interpreted [see (13)] as the current frequency of *A* in the proportion of the population chosen to be replaced (here, *N*^{−1}) times its average excess in fitness (Fisher 1930). In this average excess, the differences between the coefficients of selection of *A*-bearing and *a*-bearing individuals have relative weights given by the products of the frequencies of the individual types. Proposition 2 says that the difference between the probability of ultimate fixation of *A* and its initial frequency takes the same form with weights given by expected times with given pairs of individual types over all time steps. This corresponds to a projected average excess in fitness (Lessard and Lahaie 2009).

The expected times in Proposition 2 depend on the population-scaled parameters, σ = *sN* and ρ = *rN*, for the intensity of selection and the recombination fraction, respectively. As shown in Proposition 5, expansions in σ* ^{l}*ρ

*for*

^{m}*l*+

*m*≤

*k*and a population size

*N*large enough are obtained by considering up to

*k*pure recombination or selection events in the ancestral recombination-selection graph of pairs of individuals chosen at random.

This has been applied to get analytical approximations for the probability of ultimate fixation of a single mutant *A* introduced at random into a population in which a previous mutant *B* is segregating at another locus and has reached some frequency *x*. Here we make the assumptions of weak selection (σ << 1), tight linkage (ρ << 1), and large population size (*N* >> 1). In the case of positive or negative epistasis, the first-order effect of recombination on the fixation probability is of order ρσ^{2}. It is detected as soon as one pure recombination event and one pure selection event are considered. This is best understood from the leading recombination terms when the single mutant *A* is initially linked to *B* and LD is positive, which occurs with probability *x*, and when it is the opposite, which occurs with the complementary probability 1 – *x*. In the case of coefficients of selection given by *c _{ab}* = 0,

*c*=

_{aB}*c*=

_{Ab}*c*,

*c*= 1 with 0 ≤

_{AB}*c*≤ 1, so that

*A*and

*B*are equally advantageous, and epistasis given by ε = 1 – 2

*c*, these leading terms are approximated by

*c*= 0), tertiary effects obtained by taking into account a second pure selection event have to be considered, and their weighted average is approximated by

The following explanation for the HR effect has been given in Barton and Otto (2005, p. 2354): “…beneficial alleles that arise in coupling rise rapidly and fix within the population, leading to the disappearance of the positive disequilibrium. The negative disequilibrium persists for a much longer period of time, until one or the other allele becomes fixed. Therefore, with multiplicative selection, the variance in disequilibrium present in the first generation ultimately leads to negative disequilibrium, on average.” It is the asymmetric action of selection upon the initial positive and negative disequilibria and further disequilibria generated by random genetic drift that is responsible for the accumulation of negative disequilibrium, on average. This in turn provides an evolutionary advantage to recombination. If this is true in the case of AE, then it should also be true in the case of PE, which enhances the effect of selection. Actually, the average leading recombination term in the fixation probability for *A* that takes into account two pure selection events in the presence of epistasis is given by

Therefore, our analytical approximations indicate that recombination increases the probability of ultimate fixation of beneficial mutants in the case of positive epistasis or weak negative epistasis, at least under weak enough selection. This is supported by simulations (see Figure 4A for ε = 1). Even in the case of strong negative epistasis, an increase in the recombination fraction may increase the fixation probability in a range of values of the recombination fraction (see Figures 3 and 4B for ε = –1). This effect, though small, should provide some evolutionary advantage to recombination. This suggests that random drift may be an important factor in the evolution of recombination in the presence of selection.

Since negative epistasis is generally considered to be favorable to recombination in an infinite population, our conclusion that positive epistasis or weak negative epistasis, at least under weak selection, is a condition for recombination to help ultimate fixation of beneficial mutants may appear surprising. Some numerical results for approximations based on branching processes (Barton 1995b) and simulations for modifiers of the recombination fraction (Otto and Barton 2001) obtained under the assumption of stronger selection (*sN* >> 1) suggested that this could be the case. However, this seems to have been little noticed up to now.

Note also that AE in an exact finite population is modeled by multiplicative gene action on fitness. If advantageous mutants act additively and selection is weak, then epistasis is negative and weak. This is the situation that was actually simulated in Hill and Robertson (1966). However, comparisons with results obtained under multiplicative selection showed no significant differences.

### Conditions on the selection coefficients for recombination to be favored

In the presence of epistasis and assuming that allele *B* is advantageous, so that the coefficients of selection satisfy *c _{AB}* –

*c*> 0 and

_{Ab}*c*–

_{aB}*c*> 0, the coefficient of ρσ

_{ab}^{2}in the probability of ultimate fixation of

*A*given in Proposition 6, which is the leading recombination term under weak selection and low recombination, is positive if and only if ε =

*c*–

_{AB}*c*–

_{Ab}*c*+

_{aB}*c*> 0; that is,

_{ab}*c*–

_{AB}*c*>

_{Ab}*c*–

_{aB}*c*> 0. This means that

_{ab}*A*enhances the beneficial effect of

*B*. Under these circumstances, the probability of ultimate fixation of

*A*increases as the recombination rate increases. Note, however, that this does not require that allele

*A*itself is beneficial. It is favored to go to fixation as a result of a hitchhiking effect.

In the absence of epistasis, that is when ε = 0, an approximation to the next order is necessary to detect the effect of recombination on the probability of ultimate fixation of *A*. Then the leading recombination term is given by the term in ρσ^{3} in Proposition 7, whose coefficient is positive if *A* is advantageous, so that *c _{Ab}* –

*c*=

_{ab}*c*–

_{AB}*c*> 0. Under this condition, an increase of the recombination rate results in an increase of the fixation probability. This is in agreement with the HR effect when both mutants are advantageous (Hill and Robertson 1966). Note, however, that allele

_{aB}*B*does not have to be beneficial, since the sign of

*c*–

_{aB}*c*=

_{ab}*c*–

_{AB}*c*does not matter.

_{Ab}In both cases, the effect of recombination is more pronounced when the differences in the coefficients of selection are larger and when the frequency of *B* is closer to *x* =

The situation with deleterious mutants is comprised in the above discussion, since an allele is deleterious if and only if the alternate allele is advantageous. Therefore, the probability of ultimate fixation of an allele *A* decreases as the recombination rate increases if this allele either is in positive epistatic interaction with a deleterious allele *B* or is itself deleterious in the case of no epistasis with *B*. This explains how recombination can reduce the rate at which the ratchet-like mechanism suggested by Muller (1964) slows down the rate of evolution.

### Comparisons with previous results

In the absence of selection, the expected allele frequencies do not change and the probability of fixation of a gamete can be obtained from a transformation of a transition matrix (Karlin and McGregor 1968).

In the case where alleles *B* and *b* are neutral at locus 2 while alleles *A* and *a* are under selection at locus 1, so that *c _{AB}* –

*c*=

_{Ab}*c*–

_{aB}*c*= 0 and

_{ab}*c*–

_{AB}*c*=

_{aB}*c*–

_{Ab}*c*= 1, there is no epistasis (ε = 0). It follows from (71) that the probability of ultimate fixation of

_{ab}*A*when its initial frequency is the inverse of the population size, that is

*N*

^{−1}, is approximated by

When selection acts on both loci, it is possible to use the Kolmogorov forward or backward diffusion equation for two-locus gamete frequencies (Kimura 1955) to approximate the probability of ultimate fixation of a gamete and then that of an allele. The approximations obtained in this way up to the first-order effect of selection (Hill and Robertson 1966; Ohta 1968) are consistent with Propositions 6 and 7.

Lehman and Rousset (2009) expressed the probability of ultimate fixation as a sum of expected changes in frequencies in an exact Wright–Fisher model and then considered a Taylor expansion of this expression with respect to the intensity of selection. The coefficients in this expansion were obtained from a backward approach in a neutral model and symbolic calculation. An interpretation makes a clever use of matrix theory. Their approximation (A.34) for the reduction in the probability of ultimate fixation of *A* due to interference in the case of no epistasis is consistent with Proposition 7 with the correspondences 2*s _{A}* = (

*c*–

_{Ab}*c*)σ

_{ab}*N*

^{–1}, 2

*s*= (

_{B}*c*–

_{aB}*c*)

_{ab}*σN*

^{–1}, 2

*r*=

*ρN*

^{–1}, and

*p*(0) =

_{B}*x*, ignoring terms of order

*O*(

*N*

^{–2}). Therefore, the Moran and Wright–Fisher models lead to the same results after appropriate rescaling in the limit of a large population size.

One of the main interests in the approach based on the ancestral recombination–selection graph presented in this article is that the term in σ* ^{l}*ρ

*in the fixation probability is known to come up when considering at most*

^{m}*l*– 1 selection events and

*m*recombination events in the supragenealogy of ordered sampled gametes backward in time to compute expected times in given states. The approach has been applied to the probability of ultimate fixation of a single mutant allele at one locus in initial average linkage equilibrium with another segregating locus. It could be applied as well to an allele or a gamete given any initial conditions,

*e.g.*, a double mutant given an initial single double mutant or two initial single mutants with mutants being deleterious when alone but beneficial when coupled as in simulations by Michalakis and Slatkin (1996). The results, however, would involve high-order polynomials with respect to the initial frequencies in the case of general initial conditions.

### Implementation of the approach

The approximations given in Propositions 6 and 7 for the fixation probability rely on expected times in given ordered-sample states. These expected times are computed by conditioning on events of coalescence, recombination, or selection in the ancestry of the sample. Note that the calculation time can be shortened by considering only the material that is potentially ancestral to the original sample at either of the loci at every state change backward in time. Therefore, if the current number of ancestors at only one locus is *n*_{1} and that of ancestors at both loci is *n*_{2}, then the total rate of change is**n** = (*n*_{1}, *n*_{2}) with *n* = *n*_{1} + *n*_{2}. This is the case since a recombination event on an ancestral sequence will change the state of the ancestral material only if the sequence is ancestral at both loci.

### Extensions to other models

The model considered in this article is a particular discrete-time viability model of the Moran type (Moran 1958) expressed in terms of probabilities of mortality that are linear functions with respect to some intensity of selection. This model was used to make as clear as possible the main ideas of the approach and to simplify as far as possible the rigorous justifications of the approximations. In the case of relative mortalities, for instance, the probabilities of replacement would generally be functions of any order with respect to the intensity of selection. Then Propositions 1 and 2 for the expected change in frequency and fixation probability for an allele would give only approximations up to terms of order *O*(*N*^{–2}). This would introduce technical details in the limit of a large population size. Note, however, that a fertility model with the probability of replacement depending on the type of the offspring produced instead of the type of the individual chosen to be replaced could be analogously treated.

The Wright–Fisher model would introduce other kinds of difficulties. These are related to the possible number of ancestors in the exact ancestral recombination–selection graph and the expected time to reach the most recent ultimate ancestor. The analysis presented in the *Appendix* to justify the approach would have to be refined.

It is relatively easy to understand that the fixation probability in a finite population can be expressed in terms of expected times. It is less obvious to establish the corresponding result in the limiting process for an infinite population. Actually, (69) and (40) show that*z*_{1} = *AB*, *Ab* and *z*_{2} = *aB*, *ab*, with time *t* measured in number of *N*^{2}/2 time steps, for a large enough population size *N* and small enough population-scaled parameters ρ and σ for recombination and selection, respectively. This suggests that the fixation probability corresponding to Proposition 2 in the limiting AIG is given by the same expression but with integrals instead of sums.

Like those for the coalescent (Kingman 1982), the results obtained for the Moran model in the limit of a large population size are expected to be valid for a wide range of exchangeable models. This could be established by considering an exact model extending the Cannings neutral model (Cannings 1974) to take into account not only selection (Lessard and Ladret 2007) but also recombination in the limit of a large population size. It could also be done by showing that the fixation probability obtained from a Kolmogorov backward or forward diffusion equation for gamete frequencies (Kimura 1955; Hill and Robertson 1966; Ohta 1968) takes the form given in Proposition 2 with expected values computed according to (75) in the corresponding AIG.

## Acknowledgments

We are most grateful to one of the reviewers for the suggestion that there must be some small minimum negative value of epistasis above which recombination should help fixation of beneficial mutants. This research was supported in part by the Natural Sciences and Engineering Research Council of Canada.

## Appendix: Bound for the Expected Time to Reach the Most Recent Ultimate Ancestor

Consider the Markov chain describing the number of ancestors in the ancestral recombination–selection graph backward in time for an ordered sample in a population of size *N*. For σ = *sN* and ρ = *rN* small enough, the probabilities of transition *p _{i}*

_{,}

*, from*

_{j}*i*to

*j*for

*i*,

*j*= 1, … ,

*N*, satisfy the inequalities

and*p _{i}*

_{,}

*= 0 if*

_{j}*j*<

*i*– 1 or

*j*>

*i*+ 2.

The sojourn time in state *i* ≥ 2 in number of time steps is a geometric random variable of parameter 1 – *p _{i}*

_{,}

*, whose expected value is*

_{i}*N*

^{2}/2 time steps, the upper bound is 2. Therefore, the time in number of

*N*

^{2}/2 time steps to reach state 1 for the first time from state

*i*≥ 2, denoted by

*W*, has an expected value satisfying

_{i}*M*is the number of state changes before reaching state 1 for the first time from state

_{i}*i*.

Given that the chain leaves state *i* ≥ 2, the conditional transition probabilities satisfy*X _{n}*}

_{n}_{≥0}be a Markov chain on the integers

*i*≥ 1 with

*q*

_{i}_{,}

*above as transition probabilities for*

_{j}*i*≥ 2 and state 1 absorbing. Such a Markov chain can be constructed from a sequence of independent random variables {

*U*}

_{n}

_{n}_{≥1}, each one being uniformly distributed on (0, 1]. Given

*X*

_{n}_{–1}=

*i*≥ 2, the

*n*th increment for

*n*≥ 1 is defined as

*I*

_{(}

_{a}_{,}

_{b}_{]}(

*u*) = 1 if

*a*<

*u*≤

*b*and 0 otherwise. Analogously, given

*Y*

_{n}_{–1}=

*i*≥ 2, define

*X*

_{0}=

*Y*

_{0}=

*i*≥ 2, we have

*n*≥ 1. In particular,

*N*represents the number of transitions for the Markov chain {

_{i}*Y*}

_{n}

_{n}_{≥0}to reach state 1 for the first time from state

*i*≥ 2.

The Markov chain {*Y _{n}*}

_{n}_{≥0}has transition probabilities

*i*to states

*i*+ 2,

*i*+ 1, and

*i*– 1, respectively, for

*i*≥ 2. State 1 is an absorbing state. Reaching state 1 from state

*i*≥ 2 in {

*Y*}

_{n}

_{n}_{≥0}corresponds to extinction in a Galton–Watson process {

*Z*}

_{n}

_{n}_{≥0}starting with

*Z*

_{0}=

*i*– 1 individuals, each individual leaving independently of all others 2, 1, or 0 offspring with probability

*E*(

*W*) ≤ 2

_{i}*E*(

*N*) < ∞.

_{i}## Footnotes

*Communicating editor: L. M. Wahl*

- Received June 10, 2011.
- Accepted November 14, 2011.

- Copyright © 2012 by the Genetics Society of America