Inbreeding depression (ID) is the reduced survival and fertility of offspring produced from inbreeding compared to offspring from mating between unrelated individuals. ID is a major genetic factor that affects the evolution of mating systems, including self-fertilization (Goodwillie et al., 2005) and biparental inbreeding (Uyenoyama, 1986). As inbreeding increases genome-wide homozygosity, there are two primary hypotheses about the genetic basis of ID: increased homozygosity at fitness loci with heterozygotes advantages (overdominance) and increased homozygosity of partially recessive deleterious mutations. While these hypotheses are not mutually exclusive, empirical evidence suggests that ID may be mainly contributed by partially recessive deleterious mutations (Charlesworth and Charlesworth, 1999; Charlesworth and Willis, 2009). Furthermore, empirical estimations suggest that deleterious mutations have a continuous distribution in fitness effects ranging from strongly to slightly deleterious, and mutations with larger-effect tend to be more recessive (Eyre-Walker and Keightley, 2007; Charlesworth and Willis, 2009; Charlesworth, 2012). Although the genomic mutation rate of small-effect deleterious mutations tends to be higher than large-effect ones (Mukai et al., 1972; Simmons and Crow, 1977; Klekowski and Godfrey, 1989), large-effect mutations may have a substantial contribution to ID due to high recessiveness (Lande et al., 1994; Porcher and Lande, 2005; Winn et al., 2011).

Selfing is a major form of inbreeding in plants (Jarne and Charlesworth, 1993), and the impacts of ID on the evolution of selfing have received long-standing attention. Compared to outcrossers, selfing individuals can contribute 50% more gametes to the next generation by providing pollen for both other individuals’ and their own ovules. This transmission advantage is weaker if selfing reduces the amount of exported pollen (pollen discounting; Harder and Wilson, 1998). By treating ID as a fixed parameter, early models show that without pollen discounting, selfing is favored when the fitness of selfed offspring relative to outcrossed offspring exceeds 0.5 (Lloyd, 1979; Lande and Schemske, 1985).

Subsequent models relaxed the assumption of fixed ID by considering the joint evolution of selfing and deleterious mutations. These studies found that genetic associations between alleles modifying the selfing rate and alleles affecting fitness can be important in mediating the evolution of selfing (Holsinger, 1988; Uyenoyama and Waller, 1991a, 1991b, 1991c). A modifier enhancing selfing can develop association with fitter alleles by promoting segregation and, thus, may invade even when ID is high (Holsinger, 1988; Uyenoyama and Waller, 1991a, 1991c). Conversely, in a highly selfing population, an outcrossing-enhancing modifier may invade under low ID (Kamran-Disfani and Agrawal, 2014; Xu, 2022) because outcrossing promotes effective recombination between loci, thus increasing the efficacy of selection (Uyenoyama and Waller, 1991a). However, that genetic associations with fitness alleles may only slightly affect the invasibility of a selfing rate modifier, unless the modifier or deleterious mutations have strong effects (Charlesworth et al., 1991, 1992; Schultz and Willis, 1995; Damgaard, 1996).

Given the potential importance of genetic associations between selfing rate and fitness alleles during mating system evolution, several studies have sought to identify this correlation in nature, as summarized in Table 1 . In general, the results are mixed. Some studies found a negative correlation between ID and selfing rate indicating that selfing-enhancing modifiers tend to be associated with fitter alleles, while other studies found positive or no relationship.

Nevertheless, theoretical models and the above empirical tests actually focus on different types of genetic associations and, thus, are not directly comparable. Theoretical studies look at the transient genetic association developed during the sweep of selfing rate modifier (Charlesworth et al., 1991; Uyenoyama and Waller, 1991a). This association is observable only in a short period since the modifiers sweep rapidly (Schultz and Willis, 1995). In contrast, what the empirical studies detected is the standing genetic association formed between fitness alleles and segregating selfing rate modifiers under mutation–selection balance. Therefore, to better understand the empirical results, we need theoretical analyses of standing genetic associations between selfing rate and ID in a population with selfing rate variation maintained under selection–mutation balance.

More importantly, the association between selfing rate and ID raises the question on how population-level ID should be defined and measured. Empirical studies usually estimate the population-level ID either by averaging the fitness of selfed and outcrossed offspring, or by averaging the family-level ID, referred to as population and family ID, respectively, by Johnston and Schoen (1994). Nevertheless, it is questionable whether this average-level metric truly reflects the strength of selection on selfing caused by deleterious mutations, especially when there is association between selfing rate and ID.

Here, I first show that standing genetic association between selfing rate and ID at equilibrium depends on the selection coefficient of deleterious mutations, which decreases from positive to negative as the selection coefficient becomes larger. I then show how genetic association between selfing rate and fitness alleles can impact the evolution of selfing under two scenarios: the invasion of a rare selfing rate modifier, and the evolution of the mean selfing rate from standing variation when treating selfing rate as a quantitative trait. I discuss how the population-level ID can be better measured in future empirical studies.

I use individual-based simulations to investigate the genetic association between the selfing rate and deleterious mutations in a population at selection–mutation–drift balance (C++ code is available at https://doi.org/10.5281/zenodo.10702627). Briefly, the simulation considers a hermaphroditic diploid population with non-overlapping generations and a constant population size N. Selfing rate is determined by a quantitative trait (e.g., herkogamy) controlled by multiple loci. Deleterious mutations occur at an infinite number of loci, which are linked with selfing-related loci. The level of linkage is determined by the total number of crossovers in the genome. Each generation starts with the adult population, followed by reproduction, mutation, and viability selection of the offspring.

Each individual carries two chromosomes with the length scaled to be 1. I assume the individual selfing rate α depends on the phenotype of a quantitative trait z as α = 1 1 + e − k ( z − z c ) (Degottex-Féry and Cheptou, 2023). The parameter k determines sensitivity of the selfing rate to the phenotype z. The parameter z_c is the phenotype at which the selfing rate is α = 0.5, which is altered to change the mean selfing rate of the population. I assume phenotype z is controlled by n_z identical loci, with their positions on the chromosome being 0,1/n_z ,…, ( n z − 1 ) / n z . The positions of selfing-related loci should not qualitatively affect the results, as will be discussed in the Discussion section. Each locus has two alleles A and a with additive effects. Allele A and a change the phenotypic value by – 1 / n z and 1 / n z , respectively. Here the effect of each locus is scaled by 1 / n z to eliminate the dependency of genetic variance on the number of loci n_z . The mean genomic mutation rate of the trait is U_z per generation.

Since the trait is subject to selfing-related selective forces exerting directional selection (e.g., transmission advantage and ID; Lande and Schemske, 1985) to maintain genetic variation, I assume the trait is also under stabilizing viability selection. Specifically, the viability of individuals with phenotype z is w 1 ( z ) = e − λ z 2 , where λ determines the strength of selection. Large values of λ are chosen to make viability selection strong enough to outweigh selective forces related to selfing. This ensures that the phenotypic distribution of trait z remains nearly the same when we investigate the effects of different parameter values of key factors, such as selection coefficient of deleterious mutations.

I assume deleterious mutations occur at an infinite number of loci (Kondrashov, 1985) and have identical selection coefficient s and dominance h. Each generation, the number of new mutations on each chromosome, is drawn from a Poisson distribution with parameter U/2, with their positions being drawn from a uniform distribution ~U(0,1) by excluding the positions of loci controlling the trait. The fitness component contributed by deleterious mutations is w 2 ( 1 − h s ) n h e t ( 1 − s ) n h o m , where n_het and n_hom are the number of deleterious mutations in heterozygous and homozygous states. The overall individual fitness is thus w = w 1 w 2 .

To generate the adult population after viability selection, I sample N individuals from the juvenile population with replacement, with the probability being proportional to individual fitness. To generate the offspring population, for each offspring, I first randomly sample a parent individual i from the adult population and obtain its selfing rate α_i . To determine whether the offspring is produced from self-fertilization or outcrossing, a random number ϵ is generated from uniform distribution U(0,1). When ϵ > α_i , the offspring is produced by selfing. When ϵ< α_i , the offspring is produced by outcrossing, and a second parent j ≠ i is sampled. To obtain gametes generated by meiosis, the number of crossovers between two chromosomes is drawn from a Poisson distribution with parameter L, and the position of each crossover is drawn from U(0,1).

During the simulation, family-level ID is estimated by generating an outcrossed and a selfed offspring from each adult individual, and calculating the fitness of the selfed relative to the outcrossed offspring contributed by deleterious mutations (the calculation is described two paragraphs prior). The population-level ID is calculated as the reduction of the average fitness of selfed offspring relative to the average fitness of outcrossed offspring. For each parameter, I run the population for 5,000 generations, which is sufficiently long for the population to reach the selection–mutation–drift balance. I then calculate the average of each tracked metric over the last 2,000 generations.

This study investigates the following two questions: (1) how key genetic factors affect the standing genetic associations between selfing rate and deleterious mutations at equilibrium and (2) how to measure population-level inbreeding depression to incorporate the impacts of this genetic association on the evolution of selfing.

Although individuals with a higher selfing rate will carry fewer deleterious mutations, it does not mean they will have lower ID, and it is more informative to measure the association between selfing rate and the fitness of selfed and outcrossed offspring. In fact, the association between selfing rate and family-level ID will decrease from positive to negative as the selection coefficient of deleterious mutations s increases. This is because when s is small, due to higher genome-wide homozygosity, individuals with a higher selfing rate will have less fit selfed offspring but more fit outcrossed offspring. When s is large, both selfed and outcrossed offspring of individuals with higher selfing rate are fitter.

The model thus predicts that the correlation at equilibrium between selfing rate and family-level ID tends to be negative in populations with a low mean selfing rate, but positive in highly selfing population. Results from previous studies are equivocal in supporting this prediction (see Table 1 ) perhaps because the basis of ID differs across species. Specifically, highly recessive, large-effect mutations can be effectively purged when the selfing rate is high enough, but partially recessive, weak-effect mutations are hard to be purged (Lande and Schemske, 1985; Lande et al., 1994). Therefore, large-effect mutations, which cause a negative association between selfing rate and ID, should contribute more to ID in predominantly outcrossing populations than highly selfing populations (Winn et al., 2011). The correlation may be weakest in populations with intermediate selfing rates, since small- and large-effect mutations may have comparable contributions to ID.

More importantly, the current model reveals that standing genetic associations between selfing rate modifiers and fitness alleles will impact the evolution of selfing. The average-level inbreeding depression, commonly computed in empirical studies (Johnston and Schoen, 1994), is valid in predicting whether a selfing rate modifier can invade or not, but it can be often greatly biased in determining the evolution of the mean selfing rate.

To better measure population-level inbreeding depression, it is useful to first estimate the correlation coefficient between individual selfing rate (or a phenotype related to the selfing rate) and the fitness of selfed and outcrossed offspring. Two metrics of population-level inbreeding depression can then be calculated based on Equations (3) and (5), and the values can be compared with the average-level inbreeding depression. Nevertheless, individual-based simulations suggest that there can be large variation in family-level inbreeding depression among individuals with the same selfing rate, due to variation in the number of new mutations and inbreeding history (Kelly, 2005). Large variation in family-level inbreeding depression can render the genetic correlations statistically insignificant, as found in several previous studies (see Table 1 ).

In addition, although the current simulation assumes that selfing-related loci are distributed evenly along the chromosome, the results should be qualitatively robust to the position of these loci. For intuition, note that the results are independent of the positions when all loci are nearly free recombining or completely linked. The positions of loci may make a difference only when the number of crossovers along the chromosome is intermediate, and the results should be at the intermediate between the results when the number of crossovers along the chromosome is low and high. Nevertheless, it is found that results from individual-based simulations are qualitatively similar when the number of crossovers is low and high ( Supplementary Figures 1, 3 ). Therefore, the position of loci should not qualitatively affect the conclusions above.

The original contributions presented in the study are included in the article/ Supplementary Material . Further inquiries can be directed to the corresponding author.

KX: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing.