General results

The expected change in genetic diversity at locus A per generation can be written as(8)where is the change in diversity due to selection and the change in diversity due to drift, given by(9)to leading order. In the absence of selection, = 0 while at quasi-equilibrium is given by (from Equation A10)(10)with and the expected change in neutral diversity per generation thus becomes(11)This corresponds to the classical result that = under partial selfing (Pollak 1987): the increased homozygosity caused by selfing amplifies the effect of drift, since the same allele is sampled twice every time a homozygote is sampled. When selection acts at locus B, we obtain the following expression for to the first order in s:(12)An expression to the second order in s is provided in Appendix B; however, both expressions generally yield very similar quantitative results, although adding the terms in may slightly improve the predictions under loose linkage. The different terms that appear in Equation 12 may be interpreted as follows. From the definitions given in the previous section, the term may also be written as (where again stands for the average over all individuals), thus measuring a covariance between and (note that ). The quantity is higher in individuals carrying rarer alleles at the neutral locus: for example, in the case of a biallelic neutral locus, it equals and in 00, 01, and 11 individuals, where is the frequency of allele 1. Furthermore, the quantity is higher in individuals carrying more deleterious alleles at the selected locus (it is nearly 0, 1, and 2 in individuals carrying 0, 1, and 2 deleterious alleles at locus B, assuming is small). Therefore, a positive value of indicates that rarer alleles at the neutral locus tend to be found in individuals carrying higher numbers of deleterious alleles at the selected locus, while a negative value of indicates the opposite. Recursions for the moments and over one generation are given in Appendix B, to the first order in s, and In the absence of selfing (), we obtain that at quasi-equilibrium, while is generated by the variance in linkage disequilibrium and by the effect of selection and is positive when Indeed, when a given neutral allele becomes associated (by chance) to the deleterious allele at locus B, this neutral allele tends to decrease in frequency, generating a positive (the deleterious allele at locus B tends to become associated with rarer alleles at locus A). This in turns reduces genetic diversity at the neutral locus (as shown by Equation 12), since these rarer alleles will further decrease in frequency due to their association with the deleterious allele. Because partial selfing generates cross-haplotype associations (between genes present on different haplotypes of a diploid), and are given by more complicated expressions when involving moments such as or (see Appendix B).

Similarly, the quantity that appears on the second line of Equation 12 measures a covariance between and the quantity being higher in homozygotes at locus B than in heterozygotes. A positive value of thus indicates that rarer alleles at locus A tend to be found more often in homozygotes at locus B than in heterozygotes. As shown by Equation 12, this would reduce neutral diversity when the deleterious allele is partially recessive (), since in this case homozygotes at the selected locus have a lower fitness than heterozygotes. A recursion for is given in Appendix B (Equation B14); remarkably, it shows that identity disequilibrium (covariance in homozygosity) between the two loci generates negative (that is, rarer alleles at locus A tend to be found more often in heterozygotes at locus B) even in the absence of selection. Indeed, setting in Equation B14 yields at quasi-equilibrium(13)where is the identity disequilibrium between loci A and B (Weir and Cockerham 1969). It is given by where(14)is the probability of joint identity-by-descent at the two loci. This result may be interpreted as follows. Due to identity disequilibrium, the frequency of heterozygotes at locus A is higher among heterozygotes at locus B than among homozygotes. Furthermore, the frequency of rarer neutral alleles is higher among heterozygotes than among homozygotes at locus A (this is easily seen in the case of a biallelic locus): therefore, the frequency of rarer alleles at locus A should be higher among heterozygotes at locus B. As shown by Equation 12, this effect tends to increase neutral diversity, as long as the deleterious allele is partially recessive (), so that heterozygotes at locus B have a higher fitness than homozygotes. As shown in the next subsection (high effective recombination), the effect of identity disequilibrium dominates over all other effects when the effective recombination rate is sufficiently high [more precisely, ], in which case partially recessive deleterious alleles tend to increase neutral diversity (this effect usually stays rather small). However, weaker effective recombination increases the relative importance of the terms on the first line of Equation 12 that tend to decrease neutral diversity. Furthermore, weak effective recombination changes the sign of (through the terms in s in Equation B14), due to the fact that the association between rarer alleles at locus A and the deleterious allele at locus B (represented by moments , ) generates an association between those rarer alleles and homozygosity for the deleterious allele at locus B.

As shown in Appendix B, calculating the different terms of Equation 12 at quasi-equilibrium requires computing six two-locus moments that are generated by finite population size: and Recursions for these moments are also given in Appendix B. Although the solutions obtained are rather complicated, they are readily computed numerically using Mathematica (see File S2); furthermore, we will see that they can be approximated by simpler expressions in several cases. Importantly, all expressions obtained are in and thus vanish when N tends to infinity.

Besides its effect on selection at locus B also affects the average excess homozygosity at locus A and thus the term in Equation 8. Because is multiplied by in Equation 9, it is sufficient to compute this moment in the limit as N tends to infinity, to obtain an expression for to the first order in Using the results in Roze (2015) (also derived in File S1), we have at quasi-equilibrium, to the first order in s and (15)where again is the identity disequilibrium between loci A and B. Indeed, homozygosity at locus A is reduced when the deleterious allele at locus B is partially recessive (), due to the fact that homozygotes at locus A tend to be also homozygous at locus B, while homozygotes at locus B have a lower fitness than heterozygotes.

Putting everything together, we obtain an expression for the change in diversity at locus A over one generation of the form(16)where T is a function of s, h, r, and α. To the first order in we thus have(17)where represents the effect of background selection. Although the expression obtained for T from the equations given in Appendix B is complicated, we will now see that simple approximations can be obtained in several regimes (in particular, high effective recombination, tight linkage, and high selfing).

High effective recombination (with partial selfing)

Under partial selfing and when the effective recombination rate is high (and assuming that the dominance coefficient h of the deleterious allele is significantly different from 0.5), Equation 12 is dominated by the term on the second line, since is generated by drift even in the absence of selection and is thus proportional to (see Equation 13), while the terms and on the first line are generated by selection and drift and are thus proportional to Neglecting the term on the first line of Equation 12 and using Equations 9 and 15 yields the following expression for the change in neutral diversity, to the first order in s, and (18)Using Equation 2, we obtain for the effective population size at the neutral locus(19)independent of s. Equation 19 shows that identity disequilibrium between the neutral and the selected locus () increases the effective population size at the neutral locus when the deleterious allele is partially recessive. This is caused by two effects: (i) identity disequilibrium reduces the excess homozygosity at locus A caused by selfing (Equation 15), since homozygotes at locus A tend to be also homozygous at locus B, while homozygotes at locus B have a lower fitness than heterozygotes when the deleterious allele is partially recessive; and (ii) the higher fitness of heterozygotes at locus B increases diversity at locus A since heterozygotes at locus B tend to be also heterozygous at locus A (Equation 13).

This increase in effective population size caused by identity disequilibrium usually stays modest, however (since it is expected to occur only for high effective recombination), and is thus difficult to observe in simulations. Background selection has stronger effects when the effective recombination rate becomes low, in which case the term in the first line of Equation 12 becomes of the same order of magnitude as the term on the second line [indeed, we can show that the denominator of and is proportional to ε when both and s are of order ε], while the sign of changes due to the effect of selection. Two approximations for this regime are given below (tight linkage, high selfing).

Random mating

In the absence of selfing (), Equations B4–B9 yield the following expressions for and at quasi-equilibrium,(20)(21)while and = 0. Furthermore, Equations B10–B14 yield(22)while and = 0. Finally, the rate of decay of neutral diversity is given by(23)Under tight linkage (i.e., both r and s are of order ε, where ε is a small term), we obtain from Equations 20–22(24)(25)giving(26)in agreement with the results obtained by Hudson and Kaplan (1995) and Nordborg et al. (1996). Under loose linkage (), Equations 20–23 yield(27)Setting in Equation 27 and replacing by yields which is equivalent to Robertson’s (1961) heuristic result that selection at unlinked loci decreases the effective population size by four times the additive variance in fitness—indeed, the variance in fitness caused by selection against the deleterious allele is (see also Charlesworth 2012).

Tight linkage

When r is small, Nordborg’s (1997) separation of timescales argument can be used to express associations between genes present on different haplotypes of a diploid in terms of associations between genes present on the same haplotype (see also Nordborg 2000; Padhukasahasram et al. 2008). Consider for example the association between two genes sampled at loci A and B from different haplotypes of an individual. Going backward in time, two different events may happen to the ancestral lineages of these genes: they may find themselves on the same haplotype (which may take only a few generations if both lineages stay in the same individual due to selfing) or move to different individuals due to an outcrossing event, in which case it will take a long time before they find themselves again in the same individual (assuming N is large). To leading order, the probability that these lineages join on the same haplotype before moving to different individuals is F. Considering now all possible pairs of genes at loci A and B on different haplotypes of the same individual (in all individuals of the population) and going backward in time, we may assume that a proportion F of such pairs find themselves on the same haplotype after a small number of generations, while the remaining have moved to different individual lineages (and thus become independent): therefore, —note that this approximation assumes tight linkage, as it neglects recombination events separating genes that have joined on the same haplotype, over the small number of generations considered. Using this approximation, we have(28)Similarly, (from Equation A9), which yields (assuming that is small)(29)Finally, from which we obtain (using the fact that the moment is negligible, as shown in File S1)(30)Equations 28–30 can also be obtained from Equations B4–B14, under the assumption that r is small (see File S1). Plugging Equations 28–30 into Equation 12 yields(31)Furthermore, plugging Equations 29 and 30 into Equations B4 and B10 and assuming that r is small yields the same expressions as Equations 24 and 25 (obtained for random mating) for and replacing h by r by and N by Therefore, the present model converges to Nordborg’s (1997) result when loci are tightly linked.

Figure 1A shows the decrease in effective population size at the neutral locus caused by the deleterious allele ( see Equation 17), as a function of the recombination rate between the two loci (on a log scale). When r tends to zero, tends to and selfing increases (as it decreases the frequency of the deleterious allele). When r is sufficiently high, however (∼ >0.01 for the parameter values used in Figure 1), increased selfing causes stronger background selection. Under complete selfing (), becomes independent of r and is As can be seen in Figure 1A, the tight linkage approximation accurately predicts the solution obtained from Equations B4–B14, except when r and α are high (in which case it underestimates the strength of background selection). As we will see, this discrepancy for high values of r may lead to substantial differences when integrating over a large genetic map, as the majority of deleterious alleles are only loosely linked to the neutral locus, yet significantly affect at this locus when selfing is high (this is more visible in Figure 1C, showing the strength of background selection as a function of the position of the deleterious allele along the genetic map). Note, however, that the tight linkage approximation yields the same prediction as the more general model when ().

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Background selection effect generated by a single deleterious allele ( see Equation 17) on at the neutral locus, as a function of the recombination rate r between the two loci (A and B) and of the position of the deleterious allele along a linear chromosome of total map length 10 M (C and D) (the neutral locus is located at position 0). Different colors correspond to different values of the selfing rate as indicated in A. Solid curves correspond to the results obtained from the equations in Appendix B, dashed curves in A and C correspond to the tight linkage approximation [Equation 26, replacing h by and r by ], and dotted curves in B and D correspond to the high selfing approximation (Equation 34). Parameter values:

High selfing

Simple approximations can also be obtained when the selfing rate is high, for any value of the recombination rate r. From Equations B4–B9 and assuming that α is close to 1, we obtain(32)Furthermore, Equations B10–B14 yield(33)which, using generates the following approximation for the rate of decay of neutral diversity:(34)Note that Equation 34 again yields when As shown in Figure 1D, this approximation better matches the general model than the tight linkage approximation when the selfing rate is high (but <1) and when linkage is loose. However, Figure 1B shows that it differs substantially from the general model under tight linkage and is thus expected to perform poorly when the selfing rate is not high (since in this case background selection is mainly caused by deleterious alleles that are tightly linked to the neutral locus).

Multilocus extrapolation and simulations

Following previous work (e.g., Hudson and Kaplan 1995; Nordborg 1997; Glémin 2007; Kamran-Disfani and Agrawal 2014), results from the two-locus model may be extrapolated to the case of deleterious alleles segregating at multiple loci by assuming multiplicative effects of the different deleterious alleles on neutral diversity. When the neutral locus is located at the midpoint of a linear genome with total map length R, and when all deleterious alleles have the same selection and dominance coefficients (as in the simulation program), this yields (using Equation 2)(35)where U is the deleterious mutation rate per haploid genome and corresponds to the term T in Equation 16, in which the recombination rate r is expressed in terms of the genetic distance x (in morgans) between the neutral locus and the deleterious allele, using Haldane’s mapping function (Haldane 1919). This integral can be computed numerically as shown in File S2.

Figure 2 shows the effective size of a population of census size 20,000 as a function of the selfing rate, for a haploid genomic deleterious mutation rate and genome map length M. As can be seen in Figure 2, the analytical model (solid curves) fits well with the simulation results for all values of α. Note that in Figure 2, Figure 3, Figure 4, and Figure 5, solid curves have been obtained by calculating numerically the integral in Equation 35, using a system of recursions expressed to the second order in s for the different moments shown in Appendix B (see File S1 and File S2); however, the results obtained using the recursions given in Appendix B (expressed to the first order in s) are undistinguishable in most cases (results not shown). Explicit forms can be obtained for the integral in Equation 35 (as a function of R and the different parameters of the model), when using either the tight linkage approximation or the high selfing approximation, and are given in File S2. The results shown in Figure 2 confirm that the tight linkage approximation underestimates the effect of background selection when the selfing rate is high, as loosely linked deleterious alleles significantly affect neutral diversity. In this case, the high selfing approximation is more accurate, closely matching the simulation results when

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Effective population size as a function of the selfing rate α, for different values of the dominance coefficient of deleterious alleles h. Solid circles show simulation results (in this and subsequent figures, error bars are smaller than the size of the solid circles), solid curves show analytical predictions from the complete model (see File S2), dashed curves show the tight linkage approximation [Equation 26, replacing h by and r by ], and dotted curves show the high selfing approximation (Equation 34). Parameter values:

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Effective population size as a function of the selfing rate α, for different values of genome map length R: 0.1 (red), 1 (orange), 10 (green), and 100 (blue). Solid circles show simulation results, solid curves show analytical predictions from the complete model, dashed curves (left) show the tight linkage approximation, and dotted curves (right) show the high selfing approximation. Parameter values:

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Effective population size as a function of the selfing rate α, for different values of the strength of selection against deleterious alleles s: 0.005 (blue), 0.05 (green), and 0.5 (red). Solid circles show simulation results, solid curves show analytical predictions from the complete model, dashed curves (left) show the tight linkage approximation, and dotted curves (right) show the high selfing approximation. Parameter values:

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Effective population size as a function of the selfing rate α, for different values of the dominance coefficient of deleterious alleles h. Circles and curves are the same as in Figure 2 with

Figure 3 and Figure 4 show that the full model provides accurate predictions for for different values of map length R (from 0.1 to 100) and strength of selection against deleterious alleles s (from 0.005 to 0.5). As expected, the tight linkage approximation works better at lower values of R. When the selfing rate is low, increasing s magnifies the strength of background selection, while it has the opposite effect under high selfing (due to the fact that mutation-free individuals are more abundant when selection against deleterious alleles is stronger). Finally, Figure 5 shows that increasing the deleterious mutation rate U (to 0.5 per haploid genome) increases the strength of background selection, leading to very low values of at high selfing rates ( estimated from the simulations is close to 45 when ). When discrepancies between the full model and the simulations appear at low values of h ( in particular) and intermediate values of α: these are probably caused by genetic associations between selected loci (identity disequilibrium in particular), which are neglected in the analysis. Using an expression for the frequency of deleterious allele that takes into account the effect of identity disequilibria (Roze 2015) does not significantly improve the results (not shown). Discrepancies also appear at higher values of h and high selfing rate ( 0.8–0.9), the analytical model underestimating the strength of background selection. Finally, the model overestimates the effect of background selection when the selfing rate is very high (), for all values of h: although this is not visible in Figure 5 (as is very small when α is high), it becomes apparent when is plotted on a log scale, as shown in Figure S1. For these parameter values (), selection against deleterious alleles is no longer efficient and these accumulate over time in the population (results not shown).

Equation 35 can be extended to the more realistic case where deleterious alleles at different loci have different fitness effects (assuming that selection remains sufficiently strong at most loci so that deleterious alleles stay near mutation–selection balance), by replacing s and h by functions of map position x (e.g., Charlesworth 2012). If the number of deleterious loci is large and if the distribution of fitness effects of mutations does not depend on genomic location, the integral in Equation 35 can be replaced by a double integral, over map position x and over the joint distribution of s and h. To test this, the simulation program was modified to include a log-normal distribution of selection coefficients s across loci, assuming a constant heterozygous effect of mutations (), generating a negative covariance between s and h (see figure 5 in Roze 2015). Figure 6 shows results for and for different values of the variance of deleterious effects of mutations across loci, setting the average of the log-normal distribution and the heterozygous effect of mutations so that and in all cases. As can be seen in Figure 6, increasing the variance of s slightly increases and this effect is captured by integrating Equation 34 (high selfing approximation) over the distribution of s.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6

Left: effective population size as a function of the standard deviation σ of across loci, assuming a log-normal distribution of s and setting the average of so that for all values of σ. Solid curve shows prediction obtained by integrating Equation 34 over the distribution of s and over the genetic map; solid circles show simulation results. Parameter values: In the simulations the heterozygous fitness effect is the same for all mutations and adjusted so that for all values of σ. Right: the distribution of s for (solid curve), 0.5 (dashed curve), and 1 (dotted curve).

Temporal change in population size

It is worth emphasizing that the moment-based method presented in this article yields an expression for the “inbreeding effective size” (as it is based on the rate of decay of neutral diversity or the instantaneous rate of coalescence), while coalescent-based methods (e.g., Hudson and Kaplan 1995; Nordborg 1997) yield an expression for the “coalescence effective size,” corresponding to half the average coalescence time between two sequences randomly sampled from the population. After a single change in population size, and if selection and recombination rate are sufficiently large relative to the different moments computed in Appendix B should equilibrate quickly relative to the change in neutral diversity, and the inbreeding effective population size should thus rapidly converge to its equilibrium value for the new census population size. By contrast, expected coalescence times will take longer to equilibrate. Assuming that the inbreeding effective size instantaneously reaches its new equilibrium value, the expected coalescence time t generations after the change is given by(36)where is the new inbreeding effective size and the average coalescence time at the time of the change in population size. To test this prediction, the simulation program was modified to include a neutral sequence with an infinite number of sites at the midpoint of the chromosome. Indeed, the expected coalescence time between two sequences is simply given by where μ is the mutation rate within the sequence and D the average number of differences between two sequences, given by (where the sum is over all segregating neutral sites within the sequence). Figure 7 shows an example in which during the first generations (starting from a monomorphic population, which is equivalent to the coalescence of all lineages at generation zero), while during the last generations. As can be seen in Figure 7, the quasi-equilibrium argument leading to Equation 36 correctly predicts the dynamics of both during the initial phase and after the change in population size. The quasi-equilibrium argument could also possibly be used in situations where population size changes continuously over time (i.e., exponential population growth), although this was not explored here.

• Download figure
• Open in new tab
• Download powerpoint

Figure 7

Average coalescence time between two randomly sampled sequences at a neutral locus located at the midpoint of the chromosome, as a function of time. The population is monomorphic at time zero (all lineages coalesce), during the first generations, and during the last generations. Parameter values: Dashed curve shows prediction from Equation 36, where the inbreeding effective size is obtained from Equation 35 and assumed to reach instantaneously its equilibrium value for a given N. Solid curve shows simulation results (average over 4000 replicate simulations), where is estimated from the diversity of a neutral sequence with an infinite number of sites (see text) and mutation rate per generation.