## Abstract

We analyze two models of the maintenance of quantitative genetic variance in a mixed-mating system of self-fertilization and outcrossing. In both models purely additive genetic variance is maintained by mutation and recombination under stabilizing selection on the phenotype of one or more quantitative characters. The Gaussian allele model (GAM) involves a finite number of unlinked loci in an infinitely large population, with a normal distribution of allelic effects at each locus within lineages selfed for *τ* consecutive generations since their last outcross. The infinitesimal model for partial selfing (IMS) involves an infinite number of loci in a large but finite population, with a normal distribution of breeding values in lineages of selfing age *τ*. In both models a stable equilibrium genetic variance exists, the *outcrossed equilibrium*, nearly equal to that under random mating, for all selfing rates, *r*, up to critical value, , the *purging threshold*, which approximately equals the mean fitness under random mating relative to that under complete selfing. In the GAM a second stable equilibrium, the *purged equilibrium*, exists for any positive selfing rate, with genetic variance less than or equal to that under pure selfing; as *r* increases above the outcrossed equilibrium collapses sharply to the purged equilibrium genetic variance. In the IMS a single stable equilibrium genetic variance exists at each selfing rate; as *r* increases above the equilibrium genetic variance drops sharply and then declines gradually to that maintained under complete selfing. The implications for evolution of selfing rates, and for adaptive evolution and persistence of predominantly selfing species, provide a theoretical basis for the classical view of Stebbins that predominant selfing constitutes an “evolutionary dead end.”

MANY species of flowering plants, and some hermaphroditic animals, reproduce by a mixture of self-fertilization and outcrossing, often evolving to predominant selfing (Stebbins 1957, 1974; Harder and Barrett 2006; Igic and Kohn 2006; Jarne and Auld 2006). In such mixed-mating systems inbreeding depression for fitness is a critical determinant of mating system evolution (Lande and Schemske 1985; Charlesworth and Charlesworth 1987; Charlesworth *et al.* 1990; Porcher and Lande 2005; Charlesworth and Willis 2009; Devaux *et al.* 2014). Spontaneous deleterious mutations, as well as standing genetic variation, contribute to inbreeding depression and display a strongly bimodal distribution of allelic effects on fitness (Dobzhansky 1970; Fudala and Korona 2009; Bell 2010). Lethal and semilethal mutations in standing variation are on average nearly recessive, whereas mildly deleterious mutations have slightly recessive to nearly additive fitness effects (Simmons and Crow 1977; Willis 1999; Vassilieva *et al.* 2000; Eyre-Walker and Keightley 2007).

Stabilizing selection on quantitative characters is thought to be prevalent in natural populations and, although it may be weak or fluctuating on many characters (Wright 1969; Lande and Shannon 1996; Kingsolver *et al.* 2001; Lande 2007), it may create a substantial component of the total inbreeding depression for fitness. Stabilizing selection on a quantitative character produces allelic effects on fitness that are mildly deleterious and slightly recessive (Wright 1935), in agreement with observations on the mildly deleterious component of inbreeding depression (Simmons and Crow 1977; Willis 1999).

Under stabilizing selection on a quantitative character, an allele with an additive effect on the character may be either advantageous or deleterious, depending on whether the mean phenotype is above or below the optimum, and alleles at different loci with opposite effects on the character may compensate for each other in their effects on phenotype and fitness, even when the mean phenotype is at the optimum (Fisher 1930, 1958; Wright 1931, 1969). Charlesworth (2013) highlighted the difficulty of empirically distinguishing the relative contributions of unconditionally *vs.* conditionally deleterious mutations of small effect.

With respect to primary characters of morphology, physiology, and behavior that determine fitness, Wright (1921, 1969) showed for purely additive genetic variance in the absence of mutation and selection that inbreeding can increase the genetic variance up to a factor of 2. Lande (1977) modeled the maintenance of quantitative genetic variance by mutation under stabilizing selection for regular systems of nonrandom mating, including inbreeding with no variance in inbreeding coefficient among individuals. In sharp contrast to Wright’s classical results, Lande (1977) found that a regular system of nonrandom mating has no impact on the equilibrium genetic variance maintained by mutation and stabilizing selection. Charlesworth and Charlesworth (1995) modeled the maintenance of quantitative variation by mutation under random mating *vs.* complete selfing for different models of selection, concluding that complete selfing substantially reduces the genetic variance compared to random mating, but they did not model mixed mating.

Mixed-mating systems, such as combined selfing and outcrossing, introduce the serious complication of zygotic disequilibrium (nonrandom associations of diploid genotypes among loci) and variance of inbreeding coefficient among individuals (Haldane 1949; Crow and Kimura 1970). These complications are successfully encompassed only by models of unconditionally deleterious mutations (Kondrashov 1985; Charlesworth and Charlesworth 1998). A theory of the maintenance of quantitative genetic variance does not currently exist for such mixed-mating systems. To account for zygotic disequilibrium in quantitative characters under mixed mating, we employ the selfing age structure of the population (the distribution among lineages of the number of generations since the last outcrossing) introduced by Kelly (2007) into the Kondrashov (1985) model.

We analyze two models of the maintenance of quantitative genetic variance in a mixed-mating system of self-fertilization and outcrossing. In both models purely additive genetic variance is maintained by mutation and recombination among unlinked loci under stabilizing selection on the phenotype of one or more quantitative characters. The Gaussian allele model (GAM) involves a finite number of unlinked loci in an infinitely large population, assuming a normal distribution of allelic effects at each locus within lineages selfed for *τ* consecutive generations since their last outcross. The infinitesimal model for partial selfing (IMS) involves an infinite number of loci in a very large but finite population, assuming a normal distribution of breeding values in lineages of selfing age *τ*, with no assumption on the distribution of allelic effects within loci. Aspects of the results common to both models are considered to be robust and have fundamental implications for the evolution of plant mating systems.

## Basic Assumptions of the Models

The diploid population is partially self-fertilizing such that each zygote has a probability *r* of being produced by self-fertilization and probability of being produced by outcrossing to an unrelated individual. There is no genetic variation in the selfing rate. Genetic variance in quantitative traits under stabilizing selection is assumed to be purely additive; in many cases dominance and epistasis can be largely removed by a suitable scale transformation, *e.g.*, logarithmic for size-related characters (Wright 1968; Falconer and Mackay 1996). The population is measured in each generation before selection on a (vector of) quantitative character(s), *z*. Under random mating the total additive genetic variance (or variance–covariance matrix) can be partitioned into two additive components . Diagonal elements of *V* give the genic variances of each character (twice the variance of allelic effects on each character summed over all loci), and off-diagonal elements of *V* give the genetic covariances between characters due to pleiotropy (twice the covariances of allelic effects between pairs of characters summed over all loci with pleiotropic effects on the characters). *C* is the variance–covariance matrix of twice the total covariance of allelic effects among loci within gametes due to linkage disequilibrium (nonrandom association of alleles between loci within gametes). In the absence of selection, inbreeding reduces heterozygosity, proportionally reducing additive genetic (co)variance within families and increasing additive genetic (co)variance among families (Wright 1921, 1969; Crow and Kimura 1970).

Individual environmental effects on the phenotype are assumed to be independent among individuals, normally distributed with mean 0 and variance *E*, and additive and independent of the selfing age or breeding value (total additive genetic effect summed over all loci in an individual). In any generation before phenotypic selection a cohort of selfing age *τ* has genetic and phenotypic variances, respectively, of (1a) (1b)Here is Wright’s (1921, 1969) biometrical correlation of additive effects of alleles at the same locus within individuals, rather than Malécot’s probability of identity by descent. These two measures of *f* may be rather different since small mutational changes in additive effects of alleles cause only a small decrease in the biometrical correlation between alleles at the same locus, whereas mutation completely eliminates allelic identity. Assumptions specific to each model guarantee that for all individuals within a cohort of selfing age *τ* a single value of applies to all loci affecting a given character (GAM) or to all characters (IMS). The covariance of additive effects of alleles from different gametes equals multiplied by the variance of allelic effects within gametes (or the covariance of allelic effects at different loci within gametes).

The genetic and phenotypic variances before selection in the population as a whole are denoted without subscripts as (1c) (1d)where is the frequency of selfing age *τ* in the population before selection.

Multivariate stabilizing selection on the individual phenotype, *z*, is described by a Gaussian function of the individual deviation from an optimum phenotype, *θ*, (2)where Ω is a symmetric matrix describing stabilizing and correlational selection (Lande and Arnold 1983), and superscripts ^{T} and ^{−1} respectively denote vector or matrix transpose and matrix inverse. Using the normal approximation for phenotypes and breeding values, the mean fitness of a cohort of selfing age *τ* is thenwhere is the phenotypic distribution in the cohort of selfing age *τ*. Vertical bars ❘❘ denote the determinant of a matrix and . In subsequent formulas diagonal elements of the Ω-matrix are denoted as , and for independently selected characters the off-diagonal elements are 0. The general environment is assumed to be constant such that the mean phenotype always remains at the optimum, , so the mean fitness of a selfing age cohort simplifies to (3a)and the population mean fitness is (3b)For independently selected characters, the mean fitness within each selfing age class is the product of the mean fitnesses of the characters, but this is not true for the population as a whole.

## Gaussian allele model (GAM)

### Mutation–selection balance for one character

We employ the mutation model of Kimura (1965) and Lande (1975, 1977) with *n* loci having completely additive effects on a particular character. Each allele at a given locus mutates at a constant rate with the same distribution of changes in additive effect on the character, with no directional bias. This produces a constant mutational variance for the character, , without changing the mean phenotype. Empirical estimates of the mutational variance, scaled by the environmental variance (or as here using ) are typically about – (Lande 1975, 1995; Houle *et al.* 1996; Lynch 1996). Assuming weak stabilizing selection and a high mutation rate per locus, the distribution of allelic effects at each locus is approximately Gaussian (Kimura 1965; Bürger 2000). This is consistent with empirical observations of high genomic mutation rates per character, on the order of – observed for quantitative traits in maize (Russell *et al.* 1963), and the assumption that *n* is much less than the total number of genes in the genome, on the order of –100 for the effective number of loci (Lande 1975, 1977).

For simplicity, we assume *n* unlinked loci with equal mutational variance, so that a single inbreeding coefficient applies to all loci affecting a given character within a given age class. Multiple characters are assumed to be genetically and phenotypically independent and subject to independent stabilizing selection. For multiple characters that differ in their parameters (number of loci, mutational variance, and stabilizing selection) a different set of recursions across the selfing age classes is required.

### Gamete production from the selfing age classes

Summing the additive effects of alleles at exchangeable loci in Lande (1975, 1977), the genetic covariance and genic variance of gametes produced by individuals in selfing age class are (4a) (4b)For the outcrossed class , a straightforward way to derive the components of genetic variance and gametic output is through a weighted average of gametes produced by all selfing age classes. However, this approach mixes distributions with possibly rather different genetic variances, creating substantial kurtosis within the outcrossed class, particularly when large negative linkage disequilibrium builds up by selection of different compensatory mutations in long-selfed lineages. The outcrossed class then combines (1) a subclass produced by outcrossing between long-selfed individuals, with about half the total genetic variance of their parents (since ), and (2) subclasses produced by matings with at least one outcrossed parent in which recombination halves the negative linkage disequilibrium inherited from their parents, increasing the total genetic variance. We found that pooling these subclasses of outcrossed individuals into a single class, ignoring kurtosis, can create artifactual oscillations of genetic variance in the first two age classes. To eliminate this artifact, we analyze selection acting separately on all possible types of outcrosses according to the selfing age of the parents.

Outcrossing by random mating among selfing age classes implies that for pairs of parents with selfing ages *i* and *j* the frequency of the subclass of the outcrossed class 0, denoted as , is simply the product of the frequencies of parental age classes at the adult stage, after selection in the previous generation denoted by subscript , (5)The mean fitness of the outcrossed class as a whole is then

The mean fitness of each outcrossed subclass, , depends only on its total genetic variance, but selection acts differently on unequal gametic contributions to the genic variance and covariance by parents of different selfing age. Gametes produced by the outcross class, averaged over all subclasses, have genic variance and covariance (7a) (7b)where subscript denotes the grandparental generation. The relative frequencies and mean fitness of class are obtained using the components of genetic variance for the outcrossed progeny of a mating between individuals of selfing age classes *i* and *j*, with a prime denoting the next generation, (8a) (8b) (8c)These yield the total genetic and phenotypic variance (using Equations 1) and hence the mean fitness (Equations 3) of the subclass.

### Recursions for components of genetic variance

All selfing age classes after the first obey the recursions, for , (9a) (9b) (9c)In the absence of selection and mutation Equation 9c reduces to the recursion of Wright (1921, 1969) for the inbreeding coefficient under continued selfing, .

The genetic variance components of the outcrossed class are twice the weighted average of gametic outputs from all selfing age classes: (10a) (10b) (10c)Finally, for the first selfing age class ,

(11a) (11b) (11c)### Age distribution of selfing lineages

After stabilizing selection, mating, and reproduction, the distribution of selfing ages in the population is (12a) (12b)Equations 1–12 constitute the complete recursion system for the evolution of quantitative genetic variance. For numerical computation it is necessary to truncate the age distribution of selfing lineages at some upper limit. This approach is often used in analyzing the demography of age-structured populations (Caswell 2001). Recursion formulas for the frequency and genetic composition of the terminal selfing age class are given in *Appendix A*. The number of classes retained for numerical analysis should be sufficiently large that substantially increasing it does not appreciably affect the results.

## Infinitesimal model for partial selfing (IMS)

To relax the critical but controversial assumption in the GAM of a normal distribution of allelic effects within loci (Turelli 1984; Turelli and Barton 1994) and to facilitate analysis of correlated characters, we extend Fisher’s (1918) infinitesimal model to encompass mixed mating in a large but finite population. Our infinitesimal model involves an infinite number of loci, with no assumption on the distribution of allelic effects within loci. It does, however, assume a Gaussian distribution of breeding values within each cohort of a given selfing age (justified by the central limit theorem as for the classical infinitesimal model under random mating). The accuracy of this assumption is monitored numerically, using the kurtosis of breeding values in the population.

Fisher’s (1918) infinitesimal model concerns an infinite population with an infinite number of loci each having an infinitesimal effect on a quantitative character. Selection then causes no change in allele frequencies at any locus, although it can change the mean phenotype and the linkage disequilibrium among loci (Bulmer 1971). The total genic variance, *V*, in the population thus remains constant, but the total genetic variance, *G*, evolves because selection and recombination change the total linkage disequilibrium variance among loci, *C*. Fisher’s infinitesimal model for an infinite population thus does not require mutation to maintain genetic variance.

An unrealistic feature of the classical infinitesimal model of Fisher (1918) and Bulmer (1971) for an infinite population is that the equilibrium inbreeding depression due to stabilizing selection increases with the selfing rate of a population, until at sufficiently high selfing rates stabilizing selection finally creates enough negative linkage disequilibrium to purge the genetic variance (our unpublished results). This unrealistic feature of the classical infinitesimal model can be understood from the classical result of Wright (1921, 1969) for an infinite population with no selection or mutation, in which the equilibrium genetic variance increases as a linear function of the population inbreeding coefficient.

To obtain a more realistic infinitesimal model for partial selfing, we introduce the IMS in a finite population, in which the genic variance, *V*, is maintained by a balance between mutation and random genetic drift. The IMS can then incorporate the well-known influence of inbreeding in changing the effective population size and hence the genic variance maintained by mutation.

The IMS is derived from the GAM by letting the number of loci approach infinity, , and simultaneously letting the effective population size under random mating become very large and the mutational variance become very small, and , such that remains constant. Under random mating the IMS has the same dynamics in response to selection as in the classical infinitesimal model of Fisher (1918) and Bulmer (1971), but more generally the IMS also allows the genic variance to adjust to changes in the selfing rate as follows. With these assumptions the recursions for the inbreeding coefficient as a function of selfing age, Equations 9c, 10c, and 11c, reduce to Wright’s classical formula for continued selfing, (13a)The genic variance in the total population (or in any selfing age class) obeys the recursionwhere is the effective population size at selfing rate *r*. It is well known that a completely selfing population has an effective size half that under random mating (Charlesworth and Charlesworth 1995). More generally, using methods of Wright (1931, 1969), the effective size of a population with inbreeding coefficient *f* is . Substituting this into the recursion for the genic variance and using Wright’s (1921, 1969) formula for the equilibrium inbreeding coefficient in a partially selfing population with no selection or mutation, , gives, asymptotically, the equilibrium genic variance maintained by mutation–drift balance in a large partially selfing population, (13b)where in agreement with previous results on mutation–selection balance under random mating (Clayton and Robertson 1955; Lande 1980). For any given selfing rate, this equilibrium genic variance replaces the recursions (9b), (10b), and (11b). It is half as large for complete selfing as under random mating.

Equation 13b, with Equations 1a and 1c and Wright’s relation , implies that at equilibrium in the absence of selection (for which ) the genetic variance maintained in the population is actually independent of the selfing rate, . This occurs because the reduction of genic variance with larger selfing rate (and smaller effective population size) is exactly compensated by the increased inbreeding coefficient in the population. As we show below, this allows stabilizing selection to purge quantitative genetic variance from the population at high selfing rates and produces a realistic equilibrium inbreeding depression that always decreases with increased selfing rate.

A great advantage of the IMS is that it can readily model the maintenance of genetic variability in correlated characters under multivariate selection, with genetic correlations between characters due to a combination of pleiotropic mutation and correlational selection. A key ingredient of this generality is that in the IMS Wright’s recursion for the inbreeding coefficient under continued selfing (Equation 13a) applies to all loci in the genome regardless of their pleiotropic effects on different characters.

The IMS was used to investigate how pleiotropic mutation and correlational selection changed the results. With many loci of low mutability, in the limit as with very large effective population size, such that remains constant, the GAM can be converted to a general multivariate IMS by interpreting various quantities as matrices rather than scalars, *e.g.*, in Equations 4 rewriting as (Lande 1980; Lande and Arnold 1983). In the IMS, mutational and genic variance–covariance matrices can be simply transformed by rotation of axes to produce either selectively or genically independent characters [using eigenvectors of the correlational selection matrix *γ* or the genic variance–covariance matrix ]. Our numerical analysis of the IMS therefore focused on mutationally and genically independent characters under correlational selection.

In both models deviation from the assumption of a Gaussian distribution for either the allelic effects (GAM) or breeding values (IMS) within selfing age cohorts occurs due to mixing of genetic contributions among all selfing ages upon outcrossing. We considered the models to have good accuracy when the equilibrium kurtosis in the population remained close to that for a normal distribution (). For comparison to the numerical results we derived the equilibrium kurtosis of breeding values in Wright’s (1921, 1969) neutral model of partial selfing in an infinite population with no selection or mutation, assuming normality of breeding values within selfing age cohorts, (*Appendix B*). Under random mating or complete selfing the breeding value in the population is normal. The maximum kurtosis of occurs at selfing rate Thus in Wright’s neutral model of partial selfing in an infinite population the deviation from normality of breeding values at equilibrium must be rather small.

## Analytical and Numerical Results

### Purging genetic variance by stabilizing selection under continued selfing

In both the GAM and the IMS under continued selfing, stabilizing selection purges the genetic variance, but the details of how this happens and the extent of the purging differ in the two models, as shown below.

#### Continued selfing:

In the GAM, assuming weak selection ( so that ) and small mutational variance (), Equation 9c can be expanded as a Taylor series to first order in small terms and solved for a slowly changing quasi-equilibrium, which gives . This can be used with Equations 9a, 9b, and 1 to find a first-order approximation for the asymptotic dynamics of the genetic variance and its components for large selfing age, . Thus under continued selfing in the GAM the genetic variance approaches an equilibrium, (14a)This is smaller by a factor of than the weak selection approximation for the equilibrium genetic variance under random mating, (Kimura 1965; Lande 1975, 1977).

For the GAM the above results with Equation 9a show that with continued selfing, although approaches a constant, both and continue to increase indefinitely at the asymptotic rate (14b)In view of the dynamics of above Equation 14a, this confirms that with increasing selfing age the inbreeding coefficient approaches 1.

For the IMS, a similar analysis produces the asymptotic recursion for the genetic variance under continued selfing, , the only solution of which is(15)Under continued selfing, purging of genetic variance in the IMS occurs solely by the accumulation of negative linkage from stabilizing selection; as , increasing homozygosity reduces the effective recombination in proportion to (Lande 1977) so that . By comparison, under random mating in the IMS, assuming weak stabilizing selection (), the equilibrium genetic variance is approximately .

These analytical results concerning the genetic variance and its components as functions of selfing age were confirmed numerically. Figure 1A shows for the GAM the indefinite increase of the genic variance and the negative linkage disequilibrium variance with increasing selfing age, due to the accumulation of compensatory mutations (Equation 14b) while the genetic variance approaches a constant (Equation 14a). For populations with a high selfing rate a large negative linkage disequilibrium in the GAM leads to recombination and segregation in the second generation after outcrossing (), producing a high genetic variance and low mean fitness (Figure 1, B and D). Similar effects occur in the IMS, but with restricted magnitude (Figure 2, B and D).

#### Mixed mating and the distribution of selfing ages:

With partial selfing, the distribution of selfing ages in the population plays a crucial role in the dynamics of genetic variance. At a stable equilibrium the distribution of selfing ages (Equations 12) is , where is the probability of selfing *τ* generations in a row in the absence of selection, and for , with , is the relative fitness of lineages surviving to selfing age *τ*.

Despite continued selfing and stabilizing selection eventually purging the genetic variance in both models, at any given selfing rate in the population, the distribution of selfing ages determines the extent of purging of the genetic variance in the population as a whole. This occurs because shifts in the selfing age distribution change the balance of contributions of young and old selfing ages to the genetic variance at outcrossing (Equation 9), which is then transmitted through the selfing ages (Equations 8 and 10). A low or moderate *r* shifts the selfing age distribution to the left, toward younger ages; a high *r* shifts the selfing age distribution right, toward older ages. A crucial property of the relative survival function of selfed lineages is that the mean fitness at each selfing age depends on selection on the whole organism rather than just a single character (Equation 11); this explains why selection on multiple characters affects the purging of genetic variance in each character, as illustrated in the numerical results.

Numerical examples of the mean fitness as a function of selfing age, and the selfing age distribution, are illustrated in Figure 1, C–F, and Figure 2, C–F, for the GAM and the IMS, respectively. At high selfing rate in the GAM the long-selfed lineages become reproductively isolated from the outcrossed lineages, as shown by the very low fitness of intermediate selfing ages and the nearly disjunct bimodal distribution of selfing ages (Figure 1, D and F). Similar effects occur in the IMS at high selfing rate, but to a lesser extent (Figure 2, D and F).

### Equilibrium genetic variance in the population as a function of *r*

For intermediate selfing rates the complexity of the models precludes analytical solution. Recursions for the GAM and IMS were iterated numerically to characterize the equilibrium genetic variance in the population as a function of population selfing rate for a wide range of parameter values. Below we describe the salient similarities and differences between the GAM and IMS, emphasizing the robust results common to both models. The accuracy of the assumption of a normal distribution of breeding values within selfing age cohorts in both models was assessed using the kurtosis of breeding values for the population as a whole, *κ*. To elucidate patterns in the equilibrium genetic variance, we computed for the population as a whole the mean fitness (Equation 3b); the inbreeding depression, *δ* (loss of mean fitness upon selfing *vs.* outcrossing); and the genetic variance after selection within a generation, (*Appendix B*).

In the GAM, two types of stable equilibrium exist for the total genetic variance in the population. At an *outcrossed equilibrium*, over a wide range of selfing rates from 0 up to a critical selfing rate, , termed the *purging threshold*, the genic variance, *V*, and linkage disequilibrium, *C*, evolve to nearly compensate for nonrandom mating, maintaining nearly the same total genetic variance, *G*, as under random mating. For selfing rates above the stable outcrossed equilibrium collapses to the purged equilibrium. At a *purged equilibrium*, selfing rates even slightly above the purging threshold cause a collapse of the equilibrium genetic variance after selection within a generation, , to values equal to or less than those under pure selfing (Equation 14a). A stable purged equilibrium exists for all selfing rates.

The purged equilibrium exists and is stable at any selfing rate (), but its stability becomes weaker and its domain of attraction smaller for lower selfing rates. At selfing rates below the purging threshold, , the initial condition of genetic monomorphism always leads to the outcrossed equilibrium; the purged equilibrium is attained from initial conditions with large negative linkage disequilibrium and small total genetic variance.

Figure 3 illustrates for the GAM that at selfing rates below the purging threshold, under weak stabilizing selection the equilibrium genetic variance after selection within a generation, , is only slightly smaller than that before selection, *G*. For selfing rates above the purging threshold, the equilibrium genetic variance before selection, *G*, blows up to infinity. This occurs because of the unlimited negative linkage disequilibrium (compensatory mutations) and genic variance built up under continued selfing, which segregates in the second generation after outcrossing and recombination. However, the recombinants among long-selfed lineages are strongly selected against because of their large phenotypic deviations from the optimum. For this reason, at selfing rates above the purging threshold, the equilibrium genetic variance after selection within a generation, , collapses to , the purged equilibrium after selection. For the same reason, at selfing rates above the purging threshold the equilibrium kurtosis of breeding value in the population, *κ*, blows up, but the kurtosis after selection nearly equals that for Wright’s neutral model. Thus, under the basic assumptions of the GAM, we consider the numerical results to be reasonably accurate at all selfing rates.

Figure 4 shows that for the IMS the equilibrium genetic variance in the population, *G*, also nearly equals that under random mating for selfing rates up to a purging threshold. Just above the purging threshold the equilibrium genetic variance dips sharply and then declines gradually with increasing *r*. The IMS displays only a single stable equilibrium genetic variance at every selfing rate. The kurtosis of breeding value in the population, both before and after selection, becomes large at selfing rates substantially above the purging threshold. Because excess kurtosis increases the strength of selection on the variance (Turelli and Barton 1987), the IMS overestimates the genetic variance maintained at selfing rates above the purging threshold. Qualitative differences between the models at high selfing rates arise from constancy of the genic variance at any given selfing rate in the IMS, which limits accumulation of negative linkage disequilibrium under continued selfing (compare Figure 1A to Figure 2A).

Figure 5 and Figure 6 depict, respectively for the GAM and IMS, how stabilizing selection on multiple independent characters acts synergistically to reduce and sharpen the purging threshold in comparison to that for a single character under the same intensity of stabilizing selection per character. The sharpness of the purging thresholds has low sensitivity to substantial differences in strength of stabilizing selection among characters (not shown) due to their joint influence on and by the selfing age distribution. Figure 5 and Figure 6 also display the equilibrium inbreeding depression, *δ*, and the mean fitness in the population, as functions of the selfing rate. For both models the mean fitness of the population at high selfing rates exceeds that at low selfing rates, consistent with the purging of genetic variance at high selfing rates.

Figure 7 shows analogous results for the IMS with multiple characters with no pleiotropy but under correlational selection on independent pairs of characters. The genetic covariance between pairs of mutually selected characters, *B*, is then caused solely by linkage disequilibrium, and at high selfing rates their equilibrium genetic correlation, , makes the multivariate distribution of breeding values conform closely to the shape of the fitness surface.

### General approximation for the purging threshold

For selfing rates below the purging threshold, in both the GAM and the IMS the equilibrium mean fitness in the population at the outcrossed equilibrium remains nearly the same as under random mating, . More remarkably, in both models at selfing rates above the purging threshold the mean fitness nearly equals the product of the selfing rate, *r* and the equilibrium mean fitness under completely selfing, . This can be seen most explicitly from the dotted lines for the purged equilibrium in Figure 5C and by extrapolation of the corresponding lines in Figure 6C and Figure 7D from their values at complete selfing back to the origin. The simplicity of these results indicates that in both models the purging threshold, can be accurately located by the intersection of these two lines, . The purging threshold can thus be accurately approximated as the ratio of mean fitnesses at equilibrium under random mating *vs.* pure selfing, (16)With many characters, the numerical analysis indicates a sharp purging threshold, and the analytical approximation is fairly accurate. A small inaccuracy arises, most notably in the GAM for or 10, because as a function of *r* dips slightly near the intersection with .

## Discussion

Wright (1921, 1969) showed that for a neutral model of additive genetic variance, inbreeding increases the equilibrium genetic variance in proportion to the average inbreeding coefficient in the population, so that complete selfing doubles the equilibrium genetic variance in comparison to random mating. Lande (1977) found for regular systems of mating, in which every individual performs the same type of nonrandom mating, that the equilibrium genetic variance maintained by mutation and stabilizing selection is independent of the mating system. In contrast, Charlesworth and Charlesworth (1995) analyzed the maintenance of quantitative genetic variance by mutations under stabilizing or purifying selection, concluding that complete selfing substantially reduces the equilibrium genetic variance compared to random mating; despite neglecting linkage disequilibrium, their findings agree qualitatively with Equations 14a and 15. These disparate results are reconciled in our models of mixed mating, showing that in both the GAM and the IMS the equilibrium genetic variance remains nearly the same as under random mating for selfing rates up to the purging threshold, above which a substantial purging of the genetic variance occurs.

In the GAM two possible stable equilibria exist for the genetic variance as a function of the population selfing rate. The outcrossed equilibrium has genetic variance nearly equal to that under random mating and exists only for selfing rates below the purging threshold. The purged equilibrium has lower genetic variance after selection, approaching that in long-selfed lineages, and exists for all selfing rates; in long-selfed lineages the genetic variance remains constant but in the whole population the genetic variance before selection blows up (Figure 3A). Initial outcrossing between long-selfed lineages decreases the genetic variance by half (since ), resulting in outcrossed hybrid vigor (“heterosis”); subsequent selfing or outcrossing of these allows recombination to express the accumulated genic variance previously hidden by negative linkage disequilibrium, producing breakdown in fitness. A similar process produces increasing segregation variance and reproductive isolation between geographically isolated populations with identical optimum phenotypes (Slatkin and Lande 1994; Chevin *et al.* 2014).

In the IMS only a single stable equilibrium genetic variance exists at any selfing rate; a purged equilibrium with reduced genetic variance exists only at selfing rates above the purging threshold (Figure 4A). Properties of the population at a purged equilibrium are similar in both models, but more extreme in the GAM with genic variance and negative linkage disequilibrium (compensatory mutations) increasing indefinitely under continued selfing and stabilizing selection, whereas the constant genic variance in the IMS limits the negative linkage disequilibrium and the breakdown after outcrossing of long-selfed lineages. In any model of inheritance, finite population size must eventually limit the accumulation of compensatory mutations. The IMS is therefore likely to be more realistic for most quantitative characters, especially if effective population sizes are not very large. But if the GAM is accurate for even a single character, a stable purged equilibrium will exist at all positive selfing rates (Figure 3A).

In both models the purging threshold approximately equals the ratio of the equilibrium mean fitness under random mating relative to that under pure selfing. Stabilizing selection of quantitative characters in a constant environment with the mean phenotype at the optimum produces the highest mean fitness for the population with the lowest genetic variance, which occurs under pure selfing. Selfing also increases the mean fitness by facilitating evolution of adaptive genetic correlations between characters by correlational selection and linkage disequilibrium as shown in Figure 7B (and by Lande 1984 for a regular system of inbreeding). The complexity of the phenotype and the overall strength of stabilizing selection on it are thus of paramount importance in determining the purging threshold. For a single character under moderately strong stabilizing selection the purging threshold occurs at a very high selfing rate (Figure 5), but for a complex phenotype with many sets of characters under independent stabilizing selection, the purging threshold is considerably reduced (Figure 7A).

### Genetic variance and inbreeding depression

Inbreeding depression is generally thought to be the main genetic factor counteracting Fisher’s automatic advantage of selfing and preventing the evolution of complete selfing in self-compatible species (Fisher 1941; Lande and Schemske 1985; Porcher and Lande 2005; Devaux *et al.* 2014). Previous theory indicates that inbreeding can cause rapid purging of inbreeding depression due to nearly recessive lethals and that purging of inbreeding depression due to nearly additive mildly deleterious mutations is more difficult and occurs to a lesser extent (Lande and Schemske 1985; Charlesworth *et al.* 1990; Lande *et al.* 1994; Charlesworth and Willis 2009; Porcher and Lande 2013). Consistent with this is the observation that highly selfing species show substantially reduced inbreeding depression (Husband and Schemske 1996; Winn *et al.* 2011).

Winn *et al.* (2011) found that species with intermediate selfing rates maintain a substantial inbreeding depression, comparable to that for predominantly outcrossing species. The constancy of average inbreeding depression maintained across a wide range of selfing rates below the purging threshold in the present theory is qualitatively consistent with these observations. Winn *et al.* (2011) suggest this may be due in part to selective interference with the purging of lethals at a high total inbreeding depression (Lande *et al.* 1994). This seems especially likely if the observed average inbreeding depression for species with intermediate selfing rates, , underestimates the actual total inbreeding depression, *e.g.*, due to stronger inbreeding depression in wild *vs.* experimental environments (Armbruster and Reed 2005). Porcher and Lande (2013) showed that a background inbreeding depression due to mildly deleterious mutations augments selective interference among lethals.

Charlesworth and Charlesworth (1995) reviewed limited available evidence indicating that predominantly selfing populations maintain lower quantitative genetic variance than predominant outcrossers on average. This agrees with our finding of reduced genetic variance after selection maintained at selfing rates above the purging threshold.

### Evolution of selfing rate near outcrossed and purged equilibria

Inbreeding depression is an important genetic constraint that, in combination with ecological constraints, controls evolution of the selfing rate by small genetic steps (Lande and Schemske 1985; Johnston *et al.* 2009; Porcher and Lande 2013; Devaux *et al.* 2014). In the absence of ecological constraints, selection favors gradual evolution of decreased selfing rate when the total inbreeding depression upon selfing, *δ*, exceeds 0.5, which typifies species with low or intermediate selfing rates (Husband and Schemske 1996; Winn *et al.* 2011).

At a purged equilibrium the inbreeding depression is reduced or even negative as in the GAM (Figure 5B and Figure 6B). heterosis and breakdown of fitness in crosses between long-selfed lineages strongly select against outcrossing and augment the impact of reduced (or negative) inbreeding depression favoring the evolution of increased selfing. This *delayed outbreeding depression* among long-selfed lineages renders them reproductively isolated not only from their outcrossed relatives but also from each other. This mechanism for *speciation by predominant selfing* accords with the finding of Goldberg and Igic (2012) that phylogenetic transitions to predominant selfing often coincide with species originations. The models therefore indicate that in a large population the evolution of predominant selfing, at a purged equilibrium of quantitative genetic variance, is an irreversible evolutionary absorbing state (Bull and Charnov 1985; Charlesworth and Charlesworth 1998; Goldberg and Igic 2008; Galis *et al.* 2010), supporting the view of Stebbins (1957, 1974) that predominantly selfing plant species generally occupy terminal branches in plant phylogenies.

Stebbins also inferred that highly selfing species have an increased extinction rate and do not persist long in phylogenetic time. The present theory shows that a population at a purged equilibrium maintains less genetic variance than a population at an outcrossed equilibrium; hence in a constant environment predominantly selfing populations have a higher mean fitness at equilibrium than under random mating. Long-term environmental trends or cycles, or extreme environments such as the edge of a species range [where new selfing species often arise from outcrossers (Wright *et al.* 2013)], may exert strong directional selection. Under this scenario, in comparison to populations at an outcrossed equilibrium, highly selfing populations with reduced genetic variance will lag farther behind changes in the optimum phenotype, with consequently a lower mean fitness through time and a higher extinction rate (Lande and Shannon 1996; Lande *et al.* 2003). Our model thus provides a theoretical foundation for the classical view of Stebbins (1957, 1974) confirmed by recent empirical findings (Takebayashi and Morrell 2001; Goldberg *et al.* 2010; Goldberg and Igic 2012; Igic and Busch 2013; Wright *et al.* 2013) that predominant selfing constitutes an “evolutionary dead end.”

## Acknowledgments

We thank B. and D. Charlesworth, W. G. Hill, and J. Pannell for discussions. Support for this work was provided by a Royal Society Research Professorship and a grant from the Balzan Foundation to R.L., the French Centre National de la Recherche Scientifique program PICS grant 5273 to E.P., and time on the computing cluster (UMS 2700 OMSI) at MNHN.

## Appendix A

### Truncation of the selfing age distribution

For selfing ages larger than some value *L*, the genetic parameters within age classes should be nearly equivalent, so that all selfing age classes can be lumped into a single terminal class of age ≥*L*, which then becomes the upper limit instead of ∞ in numerical summations. Additional recursion formulas are then required for the frequency and genetic composition in the terminal age class.

Recursions for the genic variance and covariance (due to linkage disequilibrium) and the inbreeding coefficient in the terminal selfing age class areThe recursion for the terminal age class frequency iswhere and .

## Appendix B

### Population Statistics

#### Kurtosis

By assumption each selfing age cohort, and the population as a whole, has its mean breeding value at the optimum, so there is no skew (asymmetry) in the population. However, because the genetic variance within selfing age cohorts varies among selfing ages, *τ*, the breeding value in the population will be leptokurtic. The standardized kurtosis of breeding value in the population is the weighted average fourth central moment within cohorts, divided by the square of the population variance in breeding value. For a normal distribution with variance the fourth central moment is , with a standardized kurtosis of Assuming that the distribution of breeding value within each selfing age (and within subclasses of the outcrossed class) is normal, the standardized kurtosis of breeding value in the population isFor comparison, the kurtosis of breeding value in the population can be derived for Wright’s (1921, 1969) classical model of partial selfing with no selection and no mutation. The recursion for the inbreeding coefficient (Equation 13a) with has the solution . In the absence of selection the equilibrium selfing age distribution is . The average inbreeding coefficient in the population is then in agreement with Wright (1921, 1969). With no selection under any amount of outcrossing (), the population eventually approaches linkage equilibrium, The equilibrium kurtosis of breeding value in the population is . In Wright’s neutral model, as in the GAM, at or 1 the breeding value in the population is normal (). But in the IMS at , (Equation 15) and so *κ* is not defined.

#### Inbreeding depression

The total inbreeding depression in the population caused by selfing, *δ*, is one minus the ratio of mean fitness of selfed individuals divided by the mean fitness of outcrossed individuals,This formula can be evaluated from the model only for partial selfing, For random mating and complete selfing, and , the equilibrium inbreeding depression can be obtained by calculating the phenotypic distributions of offspring that would be produced by experimental selfing and outcrossing.

#### Genetic variance after selection

With multiple loci and a high selfing rate, the average genetic variance in the population *G* blows up due to a bimodal distribution of selfing ages, with long-selfed lineages dominating. This blowup is caused by accumulation of linkage disequilibrium in the long-selfed lineages. The large genetic variance expressed by recombination and the decay of linkage disequilibrium in the early selfing age classes reduce their mean fitness. With sufficiently high selfing rates the first few selfing age classes are rare, because of the low outcrossing rate and their reduced fitness. The genetic variance in the population as a whole after selection iswhere the mean fitness in the population is the same as in Equation 3b.

## Footnotes

*Communicating editor: N. H. Barton*

- Received March 24, 2015.
- Accepted May 1, 2015.

- Copyright © 2015 by the Genetics Society of America