# Mating system and speciation I: Accumulation of genetic incompatibilities in allopatry

Lucas Marie-Orleach, Christian Brochmann, Sylvain Glémin

► To cite this version:

Lucas Marie-Orleach, Christian Brochmann, Sylvain Glémin. Mating system and speciation I: Accumulation of genetic incompatibilities in allopatry. PLoS Genetics, 2022, 18 (12), pp.e1010353. 10.1371/journal.pgen.1010353. hal-03927196

HAL Id: hal-03927196

https://hal.science/hal-03927196

Submitted on 6 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
RESEARCH ARTICLE

# Mating system and speciation I: Accumulation of genetic incompatibilities in allopatry

Lucas Marie-Orleach¹,²*, Christian Brochmann¹, Sylvain Glémin²,³

¹ Natural History Museum, University of Oslo, Oslo, Norway, ² CNRS, Université Rennes 1, ECOBIO (écosystèmes, biodiversité, évolution) – UMR 6553, Rennes, France, ³ Department of Ecology and Genetics, Evolutionary Biology Centre, Uppsala University, Uppsala, Sweden

* lucas.marie-orleach@univ-rennes1.fr

<figure><img src="image_2.png" /></figure>

## Abstract

Self-fertilisation is widespread among hermaphroditic species across the tree of life. Selfing has many consequences on the genetic diversity and the evolutionary dynamics of populations, which may in turn affect macroevolutionary processes such as speciation. On the one hand, because selfing increases genetic drift and reduces migration rate among populations, it may be expected to promote speciation. On the other hand, because selfing reduces the efficacy of selection, it may be expected to hamper ecological speciation. To better understand under which conditions and in which direction selfing affects the build-up of reproductive isolation, an explicit population genetics model is required. Here, we focus on the interplay between genetic drift, selection and genetic linkage by studying speciation without gene flow. We test how fast populations with different rates of selfing accumulate mutations leading to genetic incompatibilities. When speciation requires populations to pass through a fitness valley caused by underdominant and compensatory mutations, selfing reduces the depth and/or breadth of the valley, and thus overall facilitates the fixation of incompatibilities. When speciation does not require populations to pass through a fitness valley, as for Bateson-Dobzhanzky-Muller incompatibilities (BDMi), the lower effective population size and higher genetic linkage in selfing populations both facilitate the fixation of incompatibilities. Interestingly, and contrary to intuitive expectations, local adaptation does not always accelerate the fixation of incompatibilities in outcrossing relative to selfing populations. Our work helps to clarify how incompatibilities accumulate in selfing vs. outcrossing lineages, and has repercussions on the pace of speciation as well as on the genetic architecture of reproductive isolation.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Marie-Orleach L, Brochmann C, Glémin S (2022) Mating system and speciation I: Accumulation of genetic incompatibilities in allopatry. PLoS Genet 18(12): e1010353. https://doi.org/10.1371/journal.pgen.1010353

**Editor:** Lindi Wahl, University of Western Ontario, CANADA

**Received:** July 22, 2022

**Accepted:** November 25, 2022

**Published:** December 15, 2022

**Copyright:** © 2022 Marie-Orleach et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All raw datasets used to generate our graphics are available on the zenodo open repository (https://doi.org/10.5281/zenodo.7384200).

**Funding:** Funding was obtained from the Research Council of Norway (https://www.forskningsradet.no) to CB (project 274607: SpeciationClock—How fast does the 'speciation clock' tick in selfing versus outcrossing lineages?). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## Author summary

Hermaphroditic organisms may use their male gametes to fertilise their own female gametes, but species vary greatly in how much they self-fertilise. Self-fertilisation (selfing) induces many genetic modifications in populations, which may ultimately affect the rates at which lineages diversify. Here we aim to build predictions on how selfing affects the rate at which hybrid incompatibility arises between geographically isolated populations.
**Competing interests:** The authors have declared that no competing interests exist.

We developed theoretical models in which populations with variable selfing rates may fixate mutations leading to hybrid incompatibility. First, we explored scenarios wherein hybrid incompatibility is caused by mutations whose fixations necessitates a population to experience temporally deleterious effects (i.e., a fitness valley), and found that selfing reduces the breadth and depth of the fitness valley and thereby overall facilitates the accumulation of such mutations. Second, we explored scenarios wherein hybrid incompatibility is caused by interactions between derived alleles of different genes (i.e., BDMi). We found that selfing reduces the manifestation of BDMi within population, and thereby facilitates their fixation. This effect prevails even in the face of local adaptation. Thus, our study clarifies how fast hybrid incompatibilities should accumulate in selfing populations, and helps to understand how new species are expected to arise in selfing lineages.

**Introduction**

*Species belonging to the same section or species group usually cross freely in woody plants and perennial herbs, but are usually separated by incompatibility barriers in annual herbs. [...] It is possible, therefore, that the correlation [...] between sterility and life form is a reflection of a more fundamental relationship between type of breeding system and the formation of sterility barriers* [1].

The wide variety of mating systems observed in animals, plants, fungi and algae has multiple ecological and evolutionary consequences that might impact higher-level evolutionary processes such as species extinction and speciation. For instance, hermaphroditic species vary in their rate of self-fertilisation—spanning from obligate outcrossing via all degrees of mixed mating to predominant selfing species [2–4]—and this has long been argued to affect macroevolutionary processes [5–9]. Because selfing tends to reduce both the genetic diversity and the adaptive potential of populations, selfing lineages have been argued to be ‘evolutionary dead-ends’ as they can be expected to go extinct at faster rates than outcrossing lineages [5, 8, 10, 11]. However, the study of the macroevolutionary effects of selfing has mostly focused on extinction, whereas the effects of selfing on speciation have received less attention, both empirically and conceptually (but see [7, 9]).

Given that fully selfing species are extremely rare, perhaps even non-existent [2, 3], we consider the biological species concept to be applicable to selfing species too. Thus, the process of speciation in both selfing and outcrossing lineages involves the gradual build-up of genetic incompatibilities between two diverging entities. The mechanisms and pace of reproductive isolation (RI) may however differ across mating systems.

Phylogenetic studies of plants that include taxa with different mating systems can be used to assess whether rates of species diversification, and potentially of speciation and extinction, differ between selfing and outcrossing lineages (e.g., [12–15]). In the Solanaceae, outcrossing is enforced by a self-incompatibility mechanism, but this has broken down several times, resulting in independently derived self-compatible lineages. Compared to the self-incompatible lineages (i.e., obligate outcrossers), the self-compatible ones (i.e., potential selfers) show lower net diversification rates [16] which, interestingly, are likely to be caused by higher rates of both speciation and species extinction in the selfing lineages [12]. Phylogenetic studies in the Primulaceae [13] and in the Onagraceae [15] suggest that young selfing taxa experience a burst of speciation that fades away with time (i.e., a ‘senescing diversification rate’ [17]). In contrast,
evidence in the Polemoniaceae is mixed: alternative phylogenetic methods provide either posi-
tive or no associations between selfing and speciation rates [14].

The effects of selfing on speciation have also been tested using experimental crosses, to address whether RI between populations or species evolves at different rates in selfing vs. out-crossing species. To our knowledge however, there are only a few such studies. In the Arctic flora, intraspecific crosses between geographically isolated populations of eight predominantly selfing species resulted in F1 hybrids with low pollen fertility and seed set, whereas no reduced fertility was observed in the single outcrossing species for which successful crosses could be made [18, 19]. In the selfing species, it was estimated that RI may have developed over just a few millennia, and these reproductively isolated populations cannot be distinguished morphologically. Such ‘cryptic species’ are found in selfing and outcrossing species, but raise the issue of species delimitation, which has been argued to be uneven between mating systems [20]. For instance, convergent evolution of floral traits in the context of selfing syndrome [21] may lead to selfing species that are morphologically similar but reproductively isolated. In contrast, because morphological variation among populations may be more discrete in selfing species compared to outcrossing species, populations of selfing species may be more likely considered as different taxonomic species even if they belong to a single biological species [20]. Therefore, using experimental crosses to test the effects of selfing on speciation requires both measurements of RI and estimates of genetic divergence between populations or species. Such data would allow to test whether RI evolves at different paces within selfing vs. outcrossing lineages.

Experimental crosses also allow to address the genetic architecture of RI. Quantitative trait loci analyses of F2 populations of the predominantly selfing species *Draba nivalis* showed that post-zygotic incompatibilities in this Arctic species are due to single-locus underdominance, a putative chromosomal translocation, and nuclear-nuclear and cyto-nuclear epistatic incompatibilities [22, 23].

Theoretical expectations on the effects of selfing on speciation are poorly studied and not straightforward. Selfing may affect speciation by excluding evolutionary forces such as sexual conflicts (e.g., conflict between parents for resource allocation to the developing embryos) and sexual selection (e.g., pollen and sperm competition for fertilisation) [24], which often are found to promote speciation [25–27]. Importantly, selfing may affect speciation by modifying key population genetics parameters. First, selfing decreases gene flow within and among populations, enhancing their isolation and thus possibly facilitating speciation [7]. Second, the non-random sampling of gametes used for reproduction by selfing individuals reduces the effective population size. For instance, the effective population size is expected to be halved in purely selfing species compared with a randomly mating outcrossing species with the same population size [28, 29]. A reduction of the effective population size has cascading effects. It elevates genetic drift, reduces genetic polymorphism, and overall weakens selection. Third, selfing increases homozygosity, which makes recombination less efficient because homologous chromosomes tend to be identical [29]. Thus, recombination breaks down linkage disequilibrium less efficiently in strongly selfing populations, reducing the evolutionary advantages of recombination [30], and overexposing the populations to the deleterious effects of linked selection, such as background selection [31], further reducing the effective population size [32, 33].

Here, we developed analytical and simulation models of population genetics to better understand the effects of selfing on speciation. We did not consider the effects of gene flow and potential heterosis, but focused on understanding how the interplay between genetic drift, genetic linkage, and selection efficacy affects the accumulation of hybrid incompatibility. We studied how mutations leading to genetic incompatibility between populations accumulate within populations that differ in selfing rates, and provide predictions on (i) the pace of
accumulation of genetic incompatibilities and (ii) the genetic architecture of RI in selfing vs. outcrossing lineages.

We sequentially explored three types of mutations: underdominant mutations, compensa-
tory mutations, and Bateson-Dobzhansky-Muller incompatibility mutations. Underdominant
mutations have deleterious effects in heterozygotes, but have no deleterious effects in homozy-
gotes, and may for instance be due to structural variants [34]. Compensatory mutations are a
pair of mutations that are both deleterious when they occur alone in a genome, but neutral
when they occur together [35] (e.g., compensatory evolution of *cis*- and *trans*-regulation of
gene expression [36–38]). Bateson-Dobzhansky-Muller incompatibility (BDMi) mutations are
a pair of mutations that have no deleterious effects when they occur alone in a genome, but
cause genetic incompatibilities when they occur together [39–42]. Importantly, fixation of
underdominant and compensatory mutations requires the population to pass through a fitness
valley, where the mutations may be counter-selected. In contrast, fixation of BDMi mutations
may be neutral or even positively selected [43]. Because homozygotes are formed more readily
in selfing species, selfing has been shown to facilitate the fixation of underdominant mutations
[44]. It is however unknown if and how selfing modulates the accumulation of mutations with
epistatic effects, such as compensatory and BDMi mutations.

Overall, we hypothesised that the effects of selfing on the accumulation of genetic incom-
patibilities depend on the mode of speciation. If genetic incompatibilities arise through genetic
drift, selfing should promote their accumulation because underdominant and compensatory
mutations are more likely to get fixed through genetic drift. In contrast, if genetic incompati-
bilities arise as a by-product of selection (e.g., ecological speciation), we expect that they arise
most readily in outcrossing lineages because selection is less efficient in selfing populations.

**Methods**

Accumulation of BDMi in allopatry has mainly been modelled as a combinatorial process of
substitutions, each predicted by single-locus theory [42, 43]. Extending such models to selfing
populations would be straightforward but partly misleading, as they do not explicitly consider
the underlying multi-locus population genetics dynamics and the possible interactions among
alleles that can be affected by selfing. Instead, we studied the effects of selfing on speciation by
modelling—in a single population—the fate of different types of mutations that create genetic
incompatibilities among populations. We determined the probability of and the time to fixa-
tion. Fixation of a single incompatibility mutation is, in most cases, not sufficient to complete
speciation but determines the overall pace at which hybrid incompatibility builds up.

We built analytical models of single-locus incompatibilities (*i.e.*, underdominant muta-
tions) and multi-locus incompatibilities (*i.e.*, compensatory and BDMi mutations), aiming to
predict the effects of selfing on the time to fixation. We then tested the analytical predictions
by building two types of simulation models. First, we performed single-locus simulations (for
the underdominant mutations) and two-locus simulations (for the compensatory and BDMi
mutations) to characterize the underlying mechanisms by which mutations go to fixation.
These models allowed us to further test the effects of background selection (see below). Second,
to extend this to a genome scale, we performed multi-locus simulations in which mutations
can occur recurrently throughout the genome.

**Analytical models**

For all models, we considered a population of *N* hermaphroditic individuals reproducing by
selfing at rate 0 ≤ *σ* ≤ 1 (Table 1). The effective size of partially selfing populations is given by
<table><caption>Table 1. Glossary of the main notations.</caption><thead><tr><td>Symbol</td><td>Biological meaning</td></tr></thead><tbody><tr><td>N</td><td>Population size</td></tr><tr><td>N<sub>e</sub></td><td>Effective population size</td></tr><tr><td>σ</td><td>Selfing rate</td></tr><tr><td>F</td><td>Wright's fixation index</td></tr><tr><td>μ</td><td>Mutation rate</td></tr><tr><td>r</td><td>Recombination rate</td></tr><tr><td>L</td><td>Genome length (used in multi-locus simulation models only)</td></tr><tr><td>h</td><td>Dominance coefficient (besides genetic incompatibilities)</td></tr><tr><td>s</td><td>Strength of selection (besides genetic incompatibilities)</td></tr><tr><td>h<sub>c</sub>, h<sub>b</sub></td><td>Dominance coefficient of the compensatory and BDMi mutations</td></tr><tr><td>s<sub>u</sub>, s<sub>c</sub>, s<sub>b</sub></td><td>Strength of selection of the underdominant, compensatory and BDMi mutations</td></tr><tr><td>k<sub>c</sub>, k<sub>b</sub></td><td>Dominance coefficient in double heterozygotes for compensatory and BDMi mutations</td></tr></tbody></table>

https://doi.org/10.1371/journal.pgen.1010353.t001

Pollak [28]:

$$N_e = \frac{N}{1 + F} \qquad (1)$$

where F, the Wright's fixation index, is:

$$F = \frac{\sigma}{2 - \sigma} \qquad (2)$$

Note that the effective population size, $N_e$, can be further reduced in selfing populations due to background selection, which we included in the single-locus and two-locus simulations (see Simulations).

We focused on the mean time to fix the first incompatibility allele or haplotype under recurrent mutations, which can be decomposed into the mean waiting time to occurrence of the first mutation destined to be fixed ($T_{wait}$) and the mean fixation time conditioned on fixation ($T_{fix}$):

$$T = T_{wait} + T_{fix} \qquad (3)$$

The different forms of incompatibilities can be encapsulated within a single general two-locus model with epistatic interactions as summarized in Fig 1. This illustrates the unity and connections among the different models but to avoid manipulating too many parameters, which may be not relevant for each specific case, we re-parametrised the models with some composite parameters that simplify equations and have clearer biological interpretations (Fig 1 and details below).

**Single-locus incompatibility.** The model of underdominant mutations has already been studied by Charlesworth [44] but we summarize it for the sake of completeness and provide additional results. We considered a single bi-allelic locus, with the ancestral allele $A_1$ that can mutate to the derived allele $A_2$ at rate $\mu$. Following Fig 1, the fitness of genotypes $A_1A_1$, $A_1A_2$, and $A_2A_2$ are, $1, 1 + h_{AS_A}$, and $1 + s_A$ with $h_A < 0$, respectively, and the frequency of allele $A_2$ is noted $x$. To re-parametrised fitness as in standard underdominance models we set $s_u = -h_{AS_A}$ and remove the $A$ superscript, so the three fitness are: $1, 1 - s_u$, and $1 + s$.
<figure><img src="image_4.png" /><figcaption>Fig 1. Fitness matrix of a symmetrical general case made of two biallelic loci with dominance coefficients ($h_A$ and $h_B$), selection coefficients ($s_A$ and $s_B$), and epistasis coefficients ($\epsilon_{aa}$, $\epsilon_{ad}$ and $\epsilon_{dd}$). The combinations of value ranges for these parameters explain our three models of genetic incompatibility. The 2D or 3D fitness landscapes represent unique combinations of parameters.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g001

The change in allele frequencies in one generation is given by:

$$ \begin{aligned} \Delta x &= x(1-x)((1-F)(2(2s_u + s)x - 2s_u) + Fs)/\bar{W} \\ &\approx x(1-x)((1-F)(2(2s_u + s)x - 2s_u) + Fs) \end{aligned} \quad (4) $$

where $\bar{W}$ is the mean fitness of the population, approximately equal to 1 when selection is weak. Eq (4) can also be written as:

$$ \begin{aligned} \Delta x &\approx 2(1-F)(2s_u + s)x(1-x)(x-x_{eq}) \quad \text{if } F < 1 \\ &\approx sx(1-x) \quad \text{if } F = 1 \end{aligned} \quad (5) $$
with

$$x_{eq} = \frac{2s_u - F(2s_u + s)}{2(1 - F)(2s_u + s)}$$

According to diffusion theory, the probability of fixation of a single A₂ mutant is given by:

$$P_{fix} = \frac{\int_{0}^{1/2N} \exp\left(-\frac{2M_{\delta x}}{V_{\delta x}}\right) dx}{\int_{0}^{1} \exp\left(-\frac{2M_{\delta x}}{V_{\delta x}}\right) dx} \quad (6)$$

where M$_{δx}$ = Δx is the expected infinitesimal change in allele frequency, and V$_{δx}$ = x(1-x)/2N$_{e}$ is the expected infinitesimal variance.

**Two-locus incompatibilities.** We also considered models with two bi-allelic loci. A₁ and A₂, B₁ and B₂ denote the ancestral and derived alleles at the first and second locus, respectively, with frequencies x and 1-x, y and 1-y. We assumed the same mutation rate, μ, at the two loci. The recombination rate between the two loci is 0 ≤ r ≤ 0.5, where 0 corresponds to fully linked loci and 0.5 corresponds to independent loci, for example located on different chromosomes. The frequency of the four haplotypes {A₁B₁, A₁B₂, A₂B₁, A₂B₂} are denoted as {X₁, X₂, X₃, X₄}, and the frequency of the ten genotypes, the combination of haplotypes Xᵢ and Xⱼ, are denoted as Gᵢⱼ with i ≤ j (for example G₁₂ = [A₁A₁:B₁B₂]). Note that we must distinguish G₁₄ from G₂₃, which correspond to identical genotypes [A₁A₂:B₁B₂] but may differ in the haplotypes produced through gametogenesis. Changes in genotype frequencies can be obtained as follows [44]:

After meiosis, haplotype frequencies are given by:

$$X'_{1} = G_{11} + \frac{1}{2}(G_{12} + G_{13} + (1 - r)G_{14} + rG_{23})$$

$$X'_{2} = G_{22} + \frac{1}{2}(G_{12} + G_{24} + (1 - r)G_{23} + rG_{14})$$

$$X'_{3} = G_{33} + \frac{1}{2}(G_{13} + G_{34} + (1 - r)G_{23} + rG_{14})$$

$$X'_{4} = G_{44} + \frac{1}{2}(G_{24} + G_{34} + (1 - r)G_{14} + rG_{23})$$

After syngamy, genotypic frequencies are given by:

$$G'_{11} = (1 - \sigma)X'_{1} + \sigma\left(G_{11} + \frac{1}{4}(G_{12} + G_{13} + r^2G_{23} + (1 - r)^2G_{14})\right)$$

$$G'_{22} = (1 - \sigma)X'_{2} + \sigma\left(G_{22} + \frac{1}{4}(G_{12} + G_{24} + r^2G_{14} + (1 - r)^2G_{23})\right)$$

$$G'_{33} = (1 - \sigma)X'_{3} + \sigma\left(G_{33} + \frac{1}{4}(G_{13} + G_{34} + r^2G_{14} + (1 - r)^2G_{23})\right)$$

$$G'_{44} = (1 - \sigma)X'_{4} + \sigma\left(G_{44} + \frac{1}{4}(G_{24} + G_{34} + r^2G_{23} + (1 - r)^2G_{14})\right)$$

$$G'_{14} = (1 - \sigma)2X'_{1}X'_{4} + \sigma\frac{1}{2}(r^2G_{23} + (1 - r)^2G_{14})$$

$$G'_{23} = (1 - \sigma)2X'_{2}X'_{3} + \sigma\frac{1}{2}(r^2G_{14} + (1 - r)^2G_{23})$$

$$G'_{ij} = (1 - \sigma)2X'_{2}X'_{3} + \sigma\frac{1}{2}(G_{ij} + r(1 - r)(G_{14} + G_{23})) \text{ for other cases}$$
Finally, after selection:

$$G_{ij}^{sel} = \frac{w_{ij} G'_{ij}}{\sum_{i \le j} w_{ij} G'_{ij}} \qquad (7)$$

where $w_{ij}$ is the fitness of genotype $G_{ij}$, which depends on the type of two-locus incompatibility studied (see compensatory mutations and BDMi mutations below). Note that in either scenario, we do not distinguish cis and trans effects on fitness (i.e., $w_{14} = w_{23}$).

We extended previous models of compensatory mutations [45, 46] by including the effects of partial selfing. Compensatory mutations at two loci can be viewed as a generalization of the one-locus underdominant model presented above where two haplotypes are equally fit, $A_1B_1$ and $A_2B_2$, but the intermediate paths, $A_1B_2$ and $A_2B_1$ are deleterious. Thus, as for underdominant mutations, the evolution of pairs of compensatory mutations requires crossing of a fitness valley. We consider a symmetrical fitness landscape where initial mutation are deleterious so $s_A = s_B = -s_c$ and $h_A = h_B = h_c$ (the c stands for compensatory). To simplify the notation, we summarized epistatic effects such that genotypes have the following fitness (see Fig 1):

$$\begin{aligned}
w_{11} &= w_{44} &= 1 \\
w_{22} &= w_{33} &= 1 - s_c \\
w_{12} &= w_{13} &= w_{24} &= w_{34} &= 1 - h_c s_c \\
w_{14} &= w_{23} &= 1 - h_c k_c s_c
\end{aligned}$$

where $s_c \ge 0$ and $0 \le h_c \le 1$ are the strength and the dominance coefficient of the deleterious effects of each mutation respectively, and $k_c$ is the dominance coefficient for the double heterozygous genotype $A_1A_2B_1B_2$.

We developed a model of BDMi mutations by incorporating the effects of partial selfing into the model of Kimura and King [47]. BDMi mutations generate genetic incompatibility only when an individual carries both derived alleles. Thus, individuals carrying either derived allele do not experience deleterious effects, meaning that the fixation of BDMi mutations does not require crossing of a fitness valley. We assumed either a purely neutral BDMi model, that is $s_A = s_B = 0$, or a BDMi model with one of the locus involved in local adaptation, say $s_A > 0$. Implicitly, positive selection at the A locus occurs in the focal population that is modelled but not in the other diverging population, where selection could occur on the other locus (so $s_A > 0$ and $s_B = 0$ in one population and $s_B > 0$ and $s_A = 0$ in the other). To lighten notations, we set $s_A = s$ and $h_A = h$ and wrote epistasis coefficients as a function of the selection coefficient against the two derived alleles $A_2$ and $B_2$, $s_b$ (b standing for BDMi) and two dominance coefficients, $h_b$ for single heterozygotes, and $k_b h_b$ for double heterozygotes. The fitness of the genotypes are as follows:

$$\begin{aligned}
w_{11} &= w_{12} = w_{22} &= 1 \\
w_{13} &= 1 + hs \\
w_{33} &= 1 + s \\
w_{24} &= (1 + hs)(1 - h_b s_b) \\
w_{34} &= (1 + s)(1 - h_b s_b) \\
w_{14} &= w_{23} &= (1 + hs)(1 - h_b k_b s_b) \\
w_{44} &= (1 + s)(1 - s_b)
\end{aligned}$$
Finally, we also considered convergent (and symmetrical) selection where both loci can contribute to adaptation: $s_A = s_B = s > 0$ (see [43] for a classification of the different cases presented here). Under this scenario, the fitnesses are:

$$w_{11} = 1$$

$$w_{12} = w_{13} = 1 + hs$$

$$w_{22} = w_{33} = 1 + s$$

$$w_{24} = w_{34} = (1 + hs)(1 + s)(1 - h_b s_b)$$

$$w_{14} = w_{23} = (1 + hs)^2(1 - h_b k_b s_b)$$

$$w_{44} = (1 + s)^2(1 - s_b)$$

## Simulations

We simulated all three models using single-locus, two-locus, and multi-locus models to check our analytical and numerical predictions.

**Single-locus and two-locus simulations.** We developed programs to simulate populations of N hermaphroditic individuals producing gametes with mutations from ancestral to derived alleles at a rate $\mu$ and, for the two-locus simulations, with recombination events between the two loci at a rate $0 \le r \le 0.5$. Self-fertilisation may occur at a rate $0 \le \sigma \le 1$. Then, selection on offspring occurs differently in the underdominant, compensatory, and BDMi mutations simulations (as described above). Finally, genetic drift was included by sampling offspring from the genotype frequencies using a multinomial distribution using the `gsl_ran_-` multinomial function from the GNU Scientific Library [48]. For each iteration, we measured the number of generations required to fixate the derived allele (underdominant mutation), a pair of derived alleles (compensatory mutations), or either derived allele (BDMi mutations).

We included the effects of background selection in our single-locus and two-locus simulation models. Selfing is known to increase the strength of background selection, which further reduces effective population size [32, 33, 49]. The strength of this effect critically depends on genomic recombination rate. When the genomic recombination rate is high, only highly selfing populations suffer from background selection. When the genomic recombination rate is low, the effect of background selection rises linearly with selfing [32, 33]. To account for this effect, we modelled background selection as a simple additional reduction of $N_e$. We used a Dirichlet-multinomial distribution in the genetic drift function (instead of a multinomial distribution) that allowed us to tailor the effective size of the population to its selfing rate, independently of census population size, N. The Dirichlet-Multinomial distribution is a multinomial distribution where the vector of probability is itself drawn in a Dirichlet distribution of parameter $\alpha$. Compared to a multinomial distribution, the variance of genotypic frequency changes is inflated by a factor $\frac{N+\alpha}{1+\alpha}$, so that:

$$N_e = \frac{1 + \alpha}{N + \alpha} N \qquad (8)$$

We used analytical approximations in Roze [33] to determine $N_e$ under two background selection scenarios, corresponding to two levels of genomic recombination. The first scenario corresponds to a low rate of genomic recombination and leads to a linear decrease in $N_e$ with selfing rate (hereafter 'linear BG effects'). The second scenario corresponds to a high rate of genomic recombination and leads to an accelerating decrease in $N_e$ with selfing rate (hereafter 'curved BG effects'). Specifically, using the BS1 function from the supplemental material of
Roze [33] (File S2), we modelled the background selection effects assuming that deleterious alleles with selection and dominance coefficients of $s = 0.05$ and $h = 0.2$ occurred with a geno-mic mutation rate of $U = 0.1$ in a genome map length of $R = 0.5$ ('linear BG effects') or $R = 40$ ('curved BG effects'). This informed us about by how much $N_e$ was reduced due to background selection for different selfing rates. Then, we used the predicted $N_e$ to compute the correspond-ing Dirichlet parameter $\alpha$ to use in the simulation, simply by inverting Eq (8).

**Multi-locus simulations.** We used individual based forward simulations in which the individuals have diploid genomes of length $L$ elements representing loci on which mutations can occur at rate $\mu$ and between which recombination can occur at rate $r$. We modelled all the types of genetic incompatibility outlined above.

The single-locus and two-locus simulations were performed using C++ scripts, and the multi-locus simulations were performed with the software SLiM [50], both using GNU parallel [51] (see S1 Scripts).

# Results

## Underdominant mutations

Our analytical and simulation models allowed us (i) to study the effects of $s_u$ and $s$ on the accumulation of underdominant mutations; (ii) to test how background selection alters this dynamic; and (iii) to develop simple analytical predictions when the strength of the selection coefficient ($s_u$) is not fixed but drawn from a distribution of fitness effects.

**Effects of $s_u$ and $s$ on underdominant mutations.** Underdominant mutations arising in a population are first counter-selected, and it is well known that either genetic drift or selfing facilitates crossing the fitness valley [44]. For the sake of completeness, we first summarize pre-vious results (S1 Fig and S1 Appendix). With higher selfing and lower heterozygote frequen-cies, selection against underdominant mutations becomes less important, increasing their fixation probability and decreasing their time to fixation. A new mutation may be positively selected if the population's selfing rate is high enough, that is if:

$$ \sigma > \sigma_{lim} = \frac{2s_u}{s + 2s_u} \qquad (9) $$

**Effects of background selection on underdominant mutations.** Including the effects of background selection increases the time to fixation when the homozygotes $A_2A_2$ are posi-tively selected ($s > 0$) in highly selfing populations (S1(E) Fig). This is because background selection reduces the effective population size, $N_e$, and thereby the efficacy of selection. Such a low efficacy of selection reduces both the fitness valley in the heterozygotes, $A_1A_2$ and the fitness peak in the homozygotes, $A_2A_2$. However, because background selection mostly impacts highly selfing populations, in which selection on heterozygotes is irrelevant, the effect of background selection mostly manifests as a lower probability of fixation in highly selfing populations.

**Distribution of the underdominant effect on fitness.** Previous studies have only considered a single mutation with constant effect [44]. However, if we consider a flow of underdominant mutations, it is more realistic that their fitness effects vary. For example, small chromosome inversions (which can generate underdominance) are more likely to have a smaller effect than bigger ones [52]. As is often assumed for deleterious mutations [53] and to gain sufficient flexibility, we assumed that $4Ns_u$ follows a Gamma distribution with mean $\gamma$ and shape $\beta$. We found that the proportion of underdominant mutations that go to fixation in
partially selfing populations compared to outcrossing populations is well approximated by:

$$ (1 - \sigma)^{-\beta} \qquad (10) $$

which is independent of $\gamma$. Hence, because the time to fixation is mostly determined by the waiting time for the mutation that will be fixed to arise, the time to fixation in partially selfing compared to outcrossing populations is given by:

$$ (1 - \sigma)^{\beta} \qquad (11) $$

which our multi-locus simulation results confirm as long as $\gamma$ is not too small (Fig 2). It shows that when $\beta < 1$ there is a non-linear accelerating effect of selfing on the accumulation of underdominant mutations, and the lower $\beta$ is the stronger the effect.

## Compensatory mutations

We found that the strength of the deleterious effects of compensatory mutations ($s_c$) affects how fixation takes place in outcrossing populations. We derived analytical predictions, and confirmed them with two-locus and multi-locus simulations, when selection is (i) weak or (ii) strong, before (iii) incorporating the effects of background selection.

<figure><img src="image_5.png" /><figcaption>Fig 2. Underdominant mutations accumulate more rapidly in selfing populations. Analytical predictions (lines) and outcomes of multi-locus simulations (symbols) of the averaged number of generations required to fixate an underdominant mutation compared to an outcrossing population ($\sigma = 0$). $4N_{s_u}$ follows a Gamma distribution with a shape parameter $\beta$ of 0.1 (blue), 0.5 (yellow), or 1 (red), and mean parameter $\gamma$ of 10 (circles), 100 (triangles), or 1,000 (squares). $L = 100$, $N = 1,000$, $\mu = 10^{-6}$, $r = 0.01$. The analytical predictions are $(1 - \sigma)^{\beta}$, 1,000 iterations.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g002
**Compensatory mutations with weak effects.** How selfing affects the fixation of mutations depends first on the strength of the deleterious effect relative to neutrality ($N_e s_c$ of the order of 1 or lower). When the deleterious effect is small, fixation occurs in two steps. A first weakly deleterious mutation goes to fixation and then a second, compensatory, mutation restoring fitness goes to fixation, behaving as weakly beneficial. The total fixation time can thus be approximated by:

$$T_{\text{weak}} \approx T_{\text{deleterious}} + T_{\text{beneficial}} \qquad (12)$$

It is important to note that dominance of the first and second mutations is reverse. The compensatory mutation arises in a genetic background of fitness $1 - s_c$, so rescaling the fitness of heterozygote and homozygote genotypes for the compensatory mutation leads to the following relative fitness: $1 - h_c s_c / 1 - s_c \approx 1 + (1 - h_c) s_c$ and $1 / 1 - s_c \approx 1 + s_c$, respectively. The effect of selfing thus depends on the dominance coefficient ($h_c$), and the fixation is either faster or slower in selfing populations than in outcrossing populations. Selfing speeds up the fixation of recessive beneficial and dominant deleterious mutations [54]. So, the fixation of compensatory mutations is faster under selfing when $h_c > 1/2$ whereas it is faster under outcrossing when $h_c < 1/2$ (S2 Fig).

**Compensatory mutations with strong effects.** When the deleterious effect is too high to allow the fixation of a single mutation ($N_e s_c > 1$), the two compensatory mutations must segregate together in the population and go to fixation together. In this case, the key parameters are the recombination rate between loci ($r$), and the dominance coefficient in the double heterozygote ($k_c$). When $r = 0$ and $k_c = 0$, we found that the fixation time of a pair of compensatory mutations ($T_{0,0}$) may be approximated by (see S1 Appendix):

$$T_{0,0} \approx \frac{(h_c(1-F) + F)s_c}{2\mu^2} \qquad (13)$$

which shows that selfing always increases the fixation time (Figs 3B and 4, and S3 Fig). This is because evolution can follow the diagonal path of the fitness landscape (i.e., through the double heterozygous genotype $A_1A_2B_1B_2$), and thus avoid the corner path made of genotypes with low fitness (i.e., $A_1A_1B_2B_2$ or $A_2A_2B_1B_1$), and this is easier in outcrossing than in selfing populations (Fig 3 and S4 Fig). In comparison to the corner path, the diagonal path represents a lower fitness valley, which flattens out as $k$ approaches 0. Thus, when the diagonal path is neutral (i.e., $k_c = 0$) and recombination cannot occur between the two loci (i.e., $r = 0$), outcrossing populations fix compensatory mutations by forming the double mutated haplotype ($A_2B_2$). This is stable because it is not broken down by recombination and can spread in the population because it is not counter-selected.

When $k_c > 0$ and/or $r > 0$, the dynamics are very different. Either selection against double heterozygotes or recombination that breaks down double-mutated haplotypes considerably reduces the probability of crossing the fitness valley. When $r$ is small, we can show that:

$$T_{r,k} \approx T_{0,0} \frac{\sqrt{\pi e^R Erf(\sqrt{R})}}{2\sqrt{R}} \qquad (14)$$

where Erf is the error function and $R = N(1 - \sigma)(r + hks)$. The right-hand function increases very rapidly with $R$. These approximations and our simulations show that both recombination and selection against heterozygotes strongly impede the fixation of compensatory mutations, and that both effects are reduced by selfing (Figs 3B and 4, and S3 Fig). When $k_c > 0$, the diagonal path represents a fitness valley that selects against the double mutated haplotype ($A_2B_2$). And, when $r > 0$, the double mutated haplotype ($A_2B_2$) may break, forming genotypes that are
<figure><img src="image_6.png" /><figcaption>**Fig 3.** The effect of selfing on the fixation of compensatory mutations depends on the dominance coefficient in the double heterozygote ($k_c$), and the recombination rate($r$). (A) Schematic representations of the fitness landscapes for compensatory mutations with $k_c = 0$ (left) or $k_c = 1$ (right). $A_1$ and $B_1$ are the ancestral alleles. $A_2$ and $B_2$ are the derived alleles. (B) Time to fixation of a pair compensatory mutations with $k_c = 0$ (left) or $k_c = 1$ (right). The lines correspond to the analytical approximations (Eq 14) for $r = 0$ (black), $r = 0.001$ (dark grey), $r = 0.01$ (grey) and $r = 0.1$ (light grey). Dots show the corrected mean times to fixation of a pair of compensatory mutations from our two-locus simulations. We corrected the raw means because the threshold made our data right censored (i.e., missing estimates above $10^9$ generations), which we accounted for by first estimating the full distribution by fitting gamma distributions on our simulation outputs using the *fitdistriplus* R package [55], from which we estimated the corrected mean times to fixation. The dashed horizontal lines indicate the generation threshold after which simulations stop. $N = 1, 000$, $\mu = 10^{-5}$, $h_c = 0.5$, $s_c = 0.01$. 1,000 iterations. (C) Population fitness (black dots—right y axis) and the frequencies of the 10 possible genotypes on the two loci fitness landscapes (solid lines—left y axis) over the last 4200 generations preceding the fixation of the pair of compensatory mutations. The line colours match the genotype colours on the fitness landscapes. $N = 1, 000$, $\mu = 10^{-5}$, $h_c = 0.5$, $s_c = 0.01$. 100 iterations. See S3 and S4 Figs for the visualisation of additional parameter combinations.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g003

off the diagonal path and eliminating the derived alleles ($A_2$ and $B_2$) (Figs 3 and 4). In contrast, the fixation time in highly selfing populations is less influenced by $r$ and $k_c$ because their high level of homozygosity (i) makes selfing populations not following the diagonal path of the fitness landscape, and (ii) makes the double mutated haplotype ($A_2B_2$) more stable over time due to the inefficient recombination in selfers [29].
<figure><img src="image_7.png" /><figcaption>Fig 4. Fixation regime of compensatory mutations as a function of the selfing rate ($\sigma$), the strength of deleterious effect (deleteriousness), the recombination rate between the two loci ($r$), and the dominance coefficient in double heterozygotes ($k_c$). Populations with high selfing rates always take the corner path on the fitness landscape (i.e., the deepest point of the fitness valley represented by genotypes that are homozygous for the ancestral allele at one locus, and homozygous for the derived allele at the other locus), while populations with low selfing rates can take the diagonal path (i.e., double heterozygous genotype). The diagonal path may however be counter selected (when $k_c$ is relatively high) and unstable (when the recombination rate is relatively high). Note that, for the sake of simplicity, all effects are here presented as discrete categories while they are in fact gradual.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g004

Our multi-locus simulations confirmed that $r$ and $k_c$ critically affect the fixation time of compensatory mutations with strongly deleterious effects (S6 Fig).

**Effects of background selection on compensatory mutations.** Regardless of the dominance coefficient in the double heterozygote ($k_c$) and the recombination rate ($r$), background selection shortens fixation time in highly selfing populations (S5 Fig). By increasing drift, background selection reduces the efficacy of selection against the deleterious mutation, which may thus segregate at higher frequency and be more likely to be associated with the compensatory mutation. Because of background selection, $N_{es}$ can be small enough under selfing for fixation to occur in two steps (which is rapid see (13)) whereas, for the same population size, $N_{es}$ can be too high under outcrossing so that fixation can only occur in a single step (which takes much longer (14)).

## BDMi mutations

We were not able to find an analytical solution or approximation for the fixation time of mutations with BDMi mutations. Therefore, (i) we compared our simulations with neutral expectations to understand how genetic incompatibilities affect the fixation of mutations within populations (Fig 5), and how this dynamic differs in selfing populations. Moreover, we (ii) studied the effects background selection, and (iii) explored scenarios where selection decreases with selfing.

**Analytical predictions for neutral mutations ($s = 0$) without incompatibilities.** For neutral BDMi mutations ($s = 0$), when the mutation rate is low compared to neutrality ($4N_c\mu < 1$), the two incompatible alleles ($A_2$ and $B_2$) rarely segregate at the same time in the population, and one can consider that they go to fixation independently. According to basic results of the neutral theory, Eq 3 reduces to:

$$T = \frac{1}{2\mu} + 4N_e \qquad (15)$$
<figure><img src="image_8.png" /><figcaption>Fig 5. Fixation time of BDMi mutations (red) compared to neutral mutations (blue) in outcrossing populations (σ = 0) (two-locus model). (A) Higher population sizes (N) decelerate the fixation of BDMi mutations ($s_B = 2.5 \cdot 10^{-2}$, $h_B = 0.5$; red) compared to neutral mutations ($s_b = 0$; blue). $r = 0.5$, $\mu = 2.5 \cdot 10^{-5}$. 1,000 iterations. (B) Higher coefficients of dominance ($h_b$) decelerate the fixation of BDMi mutations ($s_b = 2.5 \cdot 10^{-4}$; red) compared to neutral mutations ($s_b = 0$; blue). $N = 10,000$, $r = 0.5$, $\mu = 2.5 \cdot 10^{-5}$. 10,000 iterations. (C) Higher rates of recombination (r) decelerate the fixation of BDMi mutations ($s_b = 2.5 \cdot 10^{-4}$, $h_b = 0.5$; red) compared to neutral mutations ($s_b = 0$; blue). $N = 10,000$, $\mu = 2.5 \cdot 10^{-5}$. 100,000 iterations. (D) Higher rates of mutation (μ) decelerate the fixation of BDMi mutations ($s_b = 2.5 \cdot 10^{-3}$, $h_b = 0.5$; red) compared to neutral mutations ($s_b = 0$; blue). $N = 10,000$, $r = 0.5$. 10,000 iterations.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g005

The $2\mu$ term comes from the fact that mutation can occur at both loci. In this case, selfing only shortens the time a mutation needs to spread through the population to fixation ($4N_e$), but it does not affect the waiting time ($1/2\mu$). Because the former may be negligible compared to the latter (as implicitly assumed in [43]), selfing barely impacts the fixation time as confirmed by simulations (Fig 6A).

However, when the mutation rate increases ($4N_e\mu > 1$), the mean time a mutation needs to go to fixation is shorter than $4N_e$ because mutations of the same type may arise in different individuals in the population, so that multiple mutations can get fixed simultaneously. Taking this effect into account, Kimura [56] showed that, at a single locus, the mean time to fixation under continuous mutation pressure was:

$$T_{\text{Kimura}} = \frac{4N_e}{4N_e\mu - 1} \times (\gamma + \psi(4N_e\mu)) \quad (16)$$

where $\gamma$ is Euler's constant and $\psi$ is the digamma function. Note that (16) converges to (15) for
<figure><img src="image_9.png" /><figcaption>**Fig 6.** Selfing rate accelerates the fixation of BDMi mutations for high mutation rates. The fixation time estimated from the two loci models with mutation rates (μ) for A₂ and B₂ of either 2.5·10⁻⁷ (top), 2.5·10⁻⁶ (middle), or 2.5·10⁻⁵ (bottom), and the strength of selection on the derived alleles (s) is either 0 (left), or 2.10⁻⁴ (right). When there is selection, the dominance coefficient (h) is either recessive (h = 0.1, red), co-dominant (h = 0.5, yellow), or dominant (h = 0.9, blue). The solid lines in the neutral scenario correspond to analytical approximations: 1/2μ (light grey), and Eq (17) (dark grey) (see BDMi Results section for details on the approximations). Each phase portrait shows, for a single simulation, the change in allele frequencies of A₂ and B₂ plotted from the beginning, and then every 100 generations until the fixation of a derived allele. The isoclines represent the expected benefits (warm colours) and costs (cold colours) on population fitness (multiplied by Nₑ, and with an increment of 1). The direction of the arrow indicates the expected allele change (which is the balance between mutation rates and selection), and their size indicates the strength of the change. N = 10, 000, hₐ = kₐ = 0.5, sₐ = 10⁻³, r = 0.5. 10, 000 iterations.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g006

small Nₑμ. Because our BDMi models have two mutation types (for the A and B loci), we can-
not use Eq (16) as such because mutations of different types cannot get fixed simultaneously in
the populations. As a heuristic argument, we can decompose T<sub>Kimura</sub> = T<sub>wait</sub> + T<sub>fix</sub> as in Eq (3),
where T<sub>wait</sub> = 1/μ. Thus, by using T<sub>fix</sub> = T<sub>Kimura</sub> - 1/μ instead of 4N<sub>e</sub> in (15), a more accurate
expression under high mutation rate is:

$$
T = T_{\text{Kimura}} - \frac{1}{2\mu} \tag{17}
$$

Importantly, mutations can segregate at the two loci at the same time and be jointly selected against, which is not taken into account in (17). Eq 17 thus serves as a neutral reference to assess whether incompatibilities between mutations segregating within populations affect the
<figure><img src="image_10.png" /><figcaption>**Fig 7.** Within-population interferences between BDMi mutations slow down their fixation, and occur less in highly selfing populations (two-locus model). The A panel shows the fixation time of co-dominant BDMi mutations ($h_b = 0.5$), and with a strength of selection ($s_b$) of either 0 (blue), $2.5 \cdot 10^{-4}$ (yellow), or $2.5 \cdot 10^{-3}$ (red). The B panel shows the fixation time of BDMi mutations with a strength of selection ($s_b$) of $2.5 \cdot 10^{-4}$ that are either recessive ($h_b = 0.1$; light blue), co-dominant ($h_b = 0.5$; green), or dominant ($h_b = 0.9$; light green). The solid lines in the neutral scenario correspond to analytical approximations: $1/2\mu$ (light grey), and Eq (17) (dark grey) (see BDMi Results section for details on the approximations). Each phase portrait shows, for a single simulation, the change in allele frequencies of $A_2$ and $B_2$ plotted from the beginning, and then every 100 generations until the fixation of a derived allele. The isoclines represent the expected benefits (warm colours) and costs (cold colours) on population fitness (multiplied by $N_e$, and with an increment of 1). The direction of the arrow indicates the expected allele change (which is the balance between mutation rates and selection), and their size indicates the strength of the change. $N = 10, 000$, $k_b = h_b$, $s = 0$, $r = 0.5$. 10, 000 iterations.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g007

fixation time of BDMi mutations, and whether this occurs independently of the population selfing rate.

Our simulations showed that several BDMi mutations could often segregate together (S9 Fig), and cause incompatibilities within populations (Fig 7A and S7 Fig). We showed that, in outcrossing populations, a high mutation rate and/or a large effective population size slowed down the fixation of BDMi mutations (Fig 5A and 5D). Counter-selection of such segregating incompatibilities hampers the fixation of the derived alleles (Fig 7A), and more dominant incompatibilities ($h_b > 0.5$; Fig 5B) increase fixation time. Finally, recombination helps forming incompatible $A_2B_2$ haplotypes and also reduces the accumulation of BDMi mutations (Fig 5C).

Thus, selfing has opposing effects on these dynamics. On the one hand, it increases selection by exposing the incompatible haplotype in homozygous state (especially when $h_b$ is low). On the other hand, it increases drift and reduces genetic shuffling, which reduces the occurrence of the incompatible haplotype. Overall, the second effect dominates and selfing globally reduces the time to fixation of BDMi mutations (Fig 6A, S7 and S9 Figs).

**Positively selected BDMi mutations ($s > 0$).** When there is positive selection ($s > 0$), the effect of selfing on the accumulation of BDMi mutations depends on the interaction between the mutation rate ($\mu$) and the dominance coefficient ($h$) (Fig 6B, and S7 Fig). The results are very similar between local adaptation (positive selection on one locus only) and convergent selection (same positive selection on both loci). We thus only present the case of local adaptation in the main text (see S8 Fig for convergent selection).

When the mutation rate is low, BDMi mutations get fixed as beneficial mutations, whose probability of fixation depends on selfing rate and the dominance coefficient ($h$) [57]. When the mutation is recessive (i.e., $h < 0.5$), the probability of fixation increases with selfing. When the mutation is dominant (i.e., $h > 0.5$), the probability of fixation decreases with selfing. And, when the mutation is codominant (i.e., $h = 0.5$), selfing does not affect the probability of fixation [57]. Therefore, when mutation rate is low, selfing speeds up fixation of recessive BDMi
mutations, and slows down fixation of dominant BDMi mutations and can be approximated by equation (13) in [57] (without standing variation, $P_{sv} = 0$, their notations) (Fig 6B and S7 Fig). In contrast, when mutation rate increases, selfing speeds up fixation of BDMi mutations regardless of the dominance coefficient ($h$). This is because multiple BDMi mutations segregate in the population (S9 Fig) and cause genetic incompatibilities, hampering more the fixation of BDMi mutations in outcrossing population than in selfing populations, as described above.

**Effects of background selection on BDMi mutations.** Background selection has two main effects on the fixation time of BDMi mutations: higher drift makes selection less efficient when selfing increases and the occurrence of segregating mutations at both loci less likely. When there is no selection ($s = 0$), BDMis are less counter selected within populations and accumulate faster under selfing than under outcrossing. Background selection thus reinforces the effect of selfing. When there is selection ($s > 0$), selfing reduces both selection on the beneficial allele and against the incompatible haplotype. When mutation rate is low and selection against the incompatible haplotype is limited, reducing selection on the beneficial allele predominates and selfing slows down the accumulation of BDMi mutations (Fig 8 and S10 Fig).

<figure><img src="image_11.png" /><figcaption>**Fig 8.** Effects of background selection on the accumulation of BDMi mutations (two-locus model). The graph displays the fixation time, either without (left) or with (right) selection on the derived alleles, and under different scenarios of background selection, 'linear BG effects' (yellow) or 'curved BG effects' (red), which we compared to a scenario without background selection (blue) (see methods for details on the implementation of background selection effects in our simulations). $N = 10,000$, $h = 0.5$, $s = 2.5 \cdot 10^{-4}$, $h_b = k_b = 0.5$, $s_b = 10^{-3}$, $r = 0.5$. 10,000 iterations.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g008
<figure><img src="image_12.png" /><figcaption>**Fig 9.** Effects of selfing on the fixation time of BDMi mutations when the strength of selection (s) varies with selfing rate (σ) (two-locus model). The strength of selection (s) decreases with selfing rate as follows: s = s₀(1 - σ), with s₀ = 0 (red), s₀ = 2.10⁻⁴ (yellow), s₀ = 4.10⁻⁴ (green), and s₀ = 8.10⁻⁴ (blue). N = 10, 000, μ = 2.5·10⁻⁵, h = 0.5, hB = kB = 0.5, sB = 10⁻³, r = 0.5. 10, 000 iterations.</figcaption></figure>

https://doi.org/10.1371/journal.pgen.1010353.g009

However, for higher mutation rates, reduced selection against the incompatible haplotype predominates and selfing speeds up fixation. Overall, under a wide range of conditions, even with selection, selfing facilitates the accumulation of BDMi mutations.

**BDMi mutations when selection decreases with selfing rate.** In the previous analyses we considered the selection on the A₂ or B₂ alleles to be independent of selfing (s was a constant), as for the adaptation to a new environment. Selfing only modulated the efficacy of selection through its effect on N_c and homozygosity. However, selfing can also directly affect selection, in particular in case of genetic or sexual conflicts [58]. If such conflicts play an important role in the build-up of RI, it has been proposed that selfing may slow down speciation [9]. We did not explore an explicit model involving conflicts, but to mimic such a situation we considered the simple case where selfing directly reduced the selection coefficient: s = s₀(1 − σ). This is a strong effect as selection vanishes under full selfing. We considered the case of high mutation rate with background selection, which corresponds to the optimal conditions under which selfing promotes speciation. If “conflict” selection is weak, selfing still facilitates the accumulation of BDMi mutations. However, for stronger “conflict” selection (i.e., N.s > 4), despite negative interaction among segregating BDMi, fixation of BDMis is faster under outcrossing than selfing (Fig 9).

Discussion

The role of mating systems in speciation is an old question, in particular among plant evolutionary biologists [1, 20, 59, 60]. Depending on the underlying mechanisms, selfing has been proposed to either promote or hamper speciation [8, 9]. Surprisingly, despite this long-standing interest, specific models on the role of selfing in the evolution of RI are scarce and mainly concern RI caused by underdominant mutations [44]. Here, we partly filled this gap by
expanding previous theoretical work on underdominant mutations, and by considering RI
caused by epistatic mutations, such as compensatory mutations and Bateson-Dobzhansky-
Muller incompatibility mutations. Overall, we showed that selfing facilitates the accumulation
of genetic incompatibilities among populations for a wide range of parameters, and may there-
fore promote allopatric speciation. Our results further predict that the genomic architecture of
reproductive isolation varies with mating system.

Selfing helps cross fitness valleys

The Bateson-Dobzhansky-Muller model was initially proposed as a possible solution of the puzzling question of the evolution of hybrid incompatibilities as it does not require crossing fitness valleys [61]. However, some mechanisms can facilitate crossing fitness valleys, and it is well known that selfing facilitates the fixation of underdominant mutations [44]. We extended this model to a two-locus fitness landscape (i.e., compensatory mutations) for which selfing also helps crossing the valley under most conditions. Selfing has two main effects. First, it increases drift (in particular if background selection is strong). Second, because of reduced recombination and heterozygosity, it limits the breakdown of the new fittest genotype, once it has been created. At the genome scale, assuming several loci for which underdominant fitness landscapes may occur, we showed that the effect of selfing is stronger when there is a highly skewed distribution of the depth of the valleys to be crossed (11). A similar conclusion is also likely for the compensatory mutation model although we did not obtain an equivalent analytical result.

The role of interference among mutations segregating within populations

In the simplest form of the Bateson-Dobzhansky-Muller model of speciation, genetic incom-
patibilities occur between derived alleles that are supposed to arise and become fixed indepen-
dently in different populations [7, 62]. The phase during which BDMi mutations emerge,
spread through the population and eventually get fixed is often dismissed (e.g., [42]). This nec-
essary phase has recently been argued to have important implications in speciation genetics
[63, 64]. Considering that BDMi alleles may segregate in natural populations at polymorphic
frequencies allows for instance to better explain (i) why hybrid incompatibility may be variable
between different pairs of individuals originating from the same two populations (reviewed in
[63]), and (ii) why genetic incompatibilities are widespread within species, as found in *Dro-
sophila melanogaster* [65], *Caenorhabditis elegans* [66], *Arabidopsis thaliana* [67], or the genus
*Draba* [18, 19].

To mimic genetic incompatibilities arising from mutations at multiple sites in the genome, we used elevated mutation rates (4N<sub>e</sub>μ > 1) in the two-locus model, which allowed us to dis-sect the dynamics of multiple incompatible alleles segregating within populations. Under these conditions, we showed that epistatic interactions among segregating BDMis delay their fixa-tion, especially when they are unlinked and not too recessive. These results were validated by multi-locus simulations. This effect was not predicted by previous approximations which showed that, when mutations are rare enough, RI only depends on mutation rate and not on other population parameters and reproductive mode [43]. However, our results are in agree-ment with some phenomenological models, which assumed that incompatibility was a func-tion of the genetic distance between individuals, and thus can be counter-selected in large polymorphic populations [68].

Our results show that mating system impacts the purging of segregating BDMi mutations,
and thereby their fixation time. In selfing populations, the likelihood for two BDMI mutations
to arise and be gathered in the same genotype is reduced because of higher genetic drift,
reduced polymorphism (including for BDMi mutations), and reduced genetic shuffling. The accumulation of BDMi mutations is thus favoured in selfing populations compared to outcrossing ones—even under scenarios of ecological speciation where the BDMi mutations are positively selected.

## Mating systems and the pace of speciation

Although it has been discussed for decades, whether and how selfing affects the pace of speciation remains unclear. Our results show that selfing overall facilitates the accumulation of underdominant, compensatory, and BDMi mutations in allopatry. Remarkably, this effect may even persist in the face of local adaptation. Overall, selfing broadens the spectrum of genetic incompatibilities that can fix and shortens the time to fixation under many conditions.

However, the pace of accumulation of genetic incompatibilities may not be sufficient to predict the pace of speciation. Indeed, hybrid fitness is indeed not only determined by hybrid incompatibilities, but also by heterosis effects [69–71]. In addition to accelerating the accumulation of hybrid incompatibilities, selfing is also more permissive to the fixation of weakly deleterious mutations degrading the mean population fitness (drift load) [72, 73]. Higher fitness can then be restored by outbreeding between diverged populations, which may promote gene flow, counteracting the effects of incompatibilities and hampering speciation [74]. However, this effect is expected to manifest mainly in the F1 generation [69, 70] and be rather transient [75]. In addition, heterosis generated by the accumulation of deleterious mutations has been shown to have only a weak to moderate effect on hybrid zones maintained by selection [76]. This suggests that selfing should globally facilitate speciation under a wide range of conditions, but this remains to be firmly confirmed by modelling explicitly the balance between hybrid incompatibilities and heterosis over population divergence.

The clear condition under which outcrossing should promote speciation is when speciation is driven by genomic or sexual conflicts, which may vanish under selfing. It is for instance known that sexual conflicts over maternal provisioning during seed development usually are stronger in outcrossers than in selfers (e.g., [24]), i.e. sexually antagonistic co-evolution of male and female traits is expected to go faster in outcrossing populations, promoting speciation [26]. Our basic model, in which selection decreases linearly with selfing rate, confirms this prediction as soon as antagonistic selection is strong enough (i.e., of the order of $N_e s > 5$).

Empirical results are so far limited but tend to support our predictions, in particular the accumulation of numerous incompatibilities between recently diverged populations of several selfing Arctic plant species [18, 19], and the phylogeny-based analyses suggesting self-compatible lineages in the Solanaceae have higher speciation rates than self-incompatible [12, 16]. However, the process of speciation in these cases remains poorly known. In Arctic species, divergence in allopatry seems most likely, but in the Solanaceae, selfing may have promoted speciation through limitation of gene flow, which we did not study here.

## Mating systems and the genetic architecture of reproductive isolation

Beyond the effect of mating systems on the pace of speciation, our outcomes clearly suggest that mating system also affects the genetic architecture of speciation. In particular, underdominant and compensatory mutations are expected to be found relatively more often as reproductive barriers among selfers than among outcrossers. Classical examples of genetic modifications leading to underdominant effects include chromosomal rearrangements, which can generate and maintain RI between populations or species [34, 77]. To our knowledge, there are no studies specifically comparing the occurrence of underdominant chromosomal rearrangement in selfing vs. outcrossing species. Reproductive isolation due to underdominant chromosomal
rearrangements is however more often found in plants than in animals [78], which possibly is
due to the higher frequency of selfing among plant species. Our results also suggest that repro-
ductive barriers caused by a few strongly underdominant mutations are more likely to differ
between selfers and outcrossers than reproductive barriers caused by many weakly underdomi-
nant mutations.

Compensatory effects are often discussed in the context of the evolution of gene expression for which stabilising selection may lead to co-evolution of *cis*- and *trans*-regulatory mutations (*e.g.*, a *cis*-regulatory mutation increasing gene expression may be compensated by a *trans*-regulatory mutation decreasing gene expression, or *vice versa*) [79]. Although compensatory mutations are expected to take a long time to become fixed [35, 46], co-evolution of *cis*- and *trans*-regulatory mutations have been found to contribute to RI between outcrossing species of *Drosophila* [36] and mice [38], and also between selfing species of nematode [37].

Finally, our models predict that in outcrossing species, BDMi mutations are more likely to go to fixation when they are clustered (but in repulsion) than when they are widespread. In contrast, we found no specific constraint on genomic location in selfing species such that pairs of incompatible alleles could arise everywhere in a genome. This conclusion is in line with the prediction that genes involved in local adaptation [80] or in domestication [81] should be less clustered under selfing than under outcrossing.

**Conclusions**

Our analytical and simulation models show that selfing overall fosters the accumulation of
underdominant, compensatory, and BDMi mutations in allopatry. This outcome helps us pre-
dict the speciation rates as well as the architecture of RI of selfing vs. outcrossing species. Our
results provide a theoretical basis for long-standing ideas [1, 5] and are tentatively supported
by both phylogenetic studies and crossing experiments—though additional empirical work is
needed. Future theoretical work will need to account for the effect of selfing on the stability of
RI in the face of gene flow.

**Supporting information**

S1 Fig. Effects of selfing on the accumulation of underdominant mutations.

(PDF)

S2 Fig. Effects of selfing on the accumulation of a pair of compensatory mutations when the strength of the deleterious effect ($s_c$) is low (two-locus model).

(PDF)

S3 Fig. Effects of selfing on the accumulation of a pair of compensatory mutations (two-locus model).

(PDF)

S4 Fig. Effects of selfing rate on the fitness of a population, and the path taken on the fitness landscape, over the 4200 generations preceding the fixation of the pair of compensatory mutation (two-locus model).

(PDF)

S5 Fig. Effects of background selection on the accumulation of compensatory mutations (two-locus model).

(PDF)
S6 Fig. Effects of selfing on the accumulation of compensatory mutations (multi-locus model).

(PDF)

S7 Fig. Effects of selfing and selection on the accumulation of BDMi mutations (multi-locus model).

(PDF)

S8 Fig. Effects of selfing and convergent selection on the fixation time of BDMi mutations (two-locus model).

(PDF)

S9 Fig. Effects of selfing and selection on the fixation time and the number of BDMi mutations segregating in a population (multi-locus model).

(PDF)

S10 Fig. Effects of selfing, selection and background selection on the time to fixation fixation of BDMi mutations (multi-locus model).

(PDF)

S1 Appendix. Mathematical derivations.

(PDF)

S1 Scripts. Mathematica notebook, C++ scripts, SLIM scripts.

(ZIP)

## Acknowledgments

We are grateful to Ehouarn Le Faou, Steve Ramm for providing thoughtful comments on an earlier version of the manuscript. We thank Charles Mullon, Vitor Sudbrack and Thomas Lesaffre for the careful reading of the manuscript that helped clarifying the explanation of the Dirichlet-multinomial sampling. We acknowledge the GenOuest bioinformatics core facility (https://www.genouest.org) for providing the computing infrastructure.

## Author Contributions

**Conceptualization:** Lucas Marie-Orleach, Sylvain Glémin.

**Data curation:** Lucas Marie-Orleach, Sylvain Glémin.

**Formal analysis:** Lucas Marie-Orleach, Sylvain Glémin.

**Funding acquisition:** Christian Brochmann.

**Methodology:** Lucas Marie-Orleach, Sylvain Glémin.

**Visualization:** Lucas Marie-Orleach.

**Writing – original draft:** Lucas Marie-Orleach.

**Writing – review & editing:** Lucas Marie-Orleach, Christian Brochmann, Sylvain Glémin.

## References

1. Grant V. Plant Speciation. Columbia University Press; 1981.

2. Goodwillie C, Kalisz S, Eckert CG. The evolutionary enigma of mixed mating systems in plants: occurrence, theoretical explanations, and empirical evidence. Annu Rev Ecol Evol Syst. 2005; 36:47–79.
https://doi.org/10.1146/annurev.ecolsys.36.091704.175539
3. Jarne P, Auld JR. Animals mix it up too: the distribution of self-fertilization among hermaphroditic animals. Evolution. 2006; 60(9):1816-1824. https://doi.org/10.1554/06-246.1 PMID: 17089966

4. Igic B, Kohn JR. The distribution of plant mating systems: study bias against obligately outcrossing species. Evolution. 2006; 60(5):1098-1103. https://doi.org/10.1111/j.0014-3820.2006.tb01186.x PMID: 16817548

5. Stebbins GL. Self fertilization and population variability in the higher plants. The American Naturalist. 1957; 91(861):337-354. https://doi.org/10.1086/281999

6. Lande R, Schemske DW. The evolution of self-fertilization and inbreeding depression in plants. I. Genetic models. Evolution. 1985; 39(1):24-40. https://doi.org/10.2307/2408514 PMID: 28563655

7. Coyne JA, Orr HA. Speciation. vol. 37. Sinauer associates Sunderland, MA; 2004.

8. Wright SI, Kalisz S, Slotte T. Evolutionary consequences of self-fertilization in plants. Proceedings of the Royal Society B: Biological Sciences. 2013; 280(1760):20130133. https://doi.org/10.1098/rspb. 2013.0133 PMID: 23595268

9. Cutter AD. Reproductive transitions in plants and animals: selfing syndrome, sexual selection and speciation. New Phytologist. 2019; 224(3):1080-1094. https://doi.org/10.1111/nph.16075 PMID: 31336389

10. Takebayashi N, Morrell PL. Is self-fertilization an evolutionary dead end? Revisiting an old hypothesis with genetic theories and a macroevolutionary approach. American Journal of Botany. 2001; 88 (7):1143-1150. https://doi.org/10.2307/3558325 PMID: 11454614

11. Igic B, Busch JW. Is self-fertilization an evolutionary dead end? New Phytologist. 2013; 198(2):386-397. https://doi.org/10.1111/nph.12182 PMID: 23421594

12. Goldberg EE, Kohn JR, Lande R, Robertson KA, Smith SA, Igić B. Species selection maintains self-incompatibility. Science. 2010; 330(6003):493-495. https://doi.org/10.1126/science.1194513 PMID: 20966249

13. de Vos JM, Hughes CE, Schneeweiss GM, Moore BR, Conti E. Heterostyly accelerates diversification via reduced extinction in primroses. Proceedings of the Royal Society B: Biological Sciences. 2014; 281 (1784):20140075. https://doi.org/10.1098/rspb.2014.0075 PMID: 24759859

14. Landis JB, Bell CD, Hernandez M, Zenil-Ferguson R, McCarthy EW, Soltis DE, et al. Evolution of floral traits and impact of reproductive mode on diversification in the phlox family (Polemoniaceae). Molecular Phylogenetics and Evolution. 2018; 127:878-890. https://doi.org/10.1016/j.ympev.2018.06.035 PMID: 29958983

15. Freyman WA, Höhna S. Stochastic character mapping of state-dependent diversification reveals the tempo of evolutionary decline in self-compatible Onagraceae lineages. Systematic Biology. 2019; 68 (3):505-519. https://doi.org/10.1093/sysbio/syy078 PMID: 30476308

16. Zenil-Ferguson R, Burleigh JG, Freyman WA, Igić B, Mayrose I, Goldberg EE. Interaction among ploidy, breeding system and lineage diversification. New Phytologist. 2019; 224(3):1252-1265. https://doi.org/ 10.1111/nph.16184 PMID: 31617595

17. Ho EK, Agrawal AF. Aging asexual lineages and the evolutionary maintenance of sex. Evolution. 2017; 71(7):1865-1875. https://doi.org/10.1111/evo.13260 PMID: 28444897

18. Grundt HH, Kjølner S, Borgen L, Rieseberg LH, Brochmann C. High biological species diversity in the arctic flora. Proceedings of the National Academy of Sciences. 2006; 103(4):972-975. https://doi.org/ 10.1073/pnas.0510270103 PMID: 16418291

19. Gustafsson ALS, Gussarova G, Borgen L, Ikeda H, Antonelli A, Marie-Orleach L, et al. Rapid evolution of post-zygotic reproductive isolation is widespread in Arctic plant lineages. Annals of Botany. 2022; 129(2):171-184. https://doi.org/10.1093/aob/mcab128 PMID: 34643673

20. Baker HG. Reproductive methods as factors in speciation in flowering plants. In: Cold Spring Harbor Symposia on Quantitative Biology. vol. 24. Cold Spring Harbor Laboratory Press; 1959. p. 177-191.

21. Woźniak NJ, Kappel C, Marona C, Altschmied L, Neuffer B, Sicard A. A similar genetic architecture underlies the convergent evolution of the selfing syndrome in Capsella. The Plant Cell. 2020; 32 (4):935-949. https://doi.org/10.1105/tpc.19.00551 PMID: 31964802

22. Skrede I, Brochmann C, Borgen L, Rieseberg LH. Genetics of intrinsic postzygotic isolation in a circum-polar plant species, *Draba nivalis* (Brassicaceae). Evolution. 2008; 62(8):1840-1851. https://doi.org/10.1111/j.1558-5646.2008.00418.x PMID: 18485112

23. Gustafsson ALS, Skrede I, Rowe HC, Gussarova G, Borgen L, Rieseberg LH, et al. Genetics of cryptic speciation within an arctic mustard, *Draba nivalis*. PloS one. 2014; 9(4):e93834. https://doi.org/10.1371/journal.pone.0093834 PMID: 24691072

24. Raunsgard A, Opedal ØH, Ekrem RK, Wright J, Bolstad GH, Armbruster WS, et al. Intersexual conflict over seed size is stronger in more outcrossed populations of a mixed-mating plant. Proceedings of the National Academy of Sciences. 2018; 115(45):11561-11566. https://doi.org/10.1073/pnas.1810979115 PMID: 30282740
25. Baack E, Melo MC, Rieseberg LH, Ortiz-Barrientos D. The origins of reproductive isolation in plants. *New Phytologist*. 2015; 207(4):968–984. https://doi.org/10.1111/nph.13424 PMID: 25944305

26. Lafon-Placette C, Köhler C. Endosperm-based postzygotic hybridization barriers: developmental mechanisms and evolutionary drivers. *Molecular Ecology*. 2016; 25(11):2620–2629. https://doi.org/10.1111/mec.13552 PMID: 26818717

27. Janicke T, Ritchie MG, Morrow EH, Marie-Orleach L. Sexual selection predicts species richness across the animal kingdom. *Proceedings of the Royal Society B: Biological Sciences*. 2018; 285 (1878):20180173. https://doi.org/10.1098/rspb.2018.0173 PMID: 29720412

28. Pollak E. On the theory of partially inbreeding finite populations. I. Partial selfing. *Genetics*. 1987; 117 (2):353–60. https://doi.org/10.1093/genetics/117.2.353 PMID: 3666446

29. Nordborg M. Linkage disequilibrium, gene trees and selfing: an ancestral recombination graph with partial self-fertilization. *Genetics*. 2000; 154(2):923–929. https://doi.org/10.1093/genetics/154.2.923 PMID: 10655241

30. Felsenstein J. The evolutionary advantage of recombination. *Genetics*. 1974; 78(2):737–756. https://doi.org/10.1093/genetics/78.2.737 PMID: 4448362

31. Barrett SC, Arunkumar R, Wright SI. The demography and population genomics of evolutionary transitions to self-fertilization in plants. *Philosophical Transactions of the Royal Society B: Biological Sciences*. 2014; 369(1648):20130344. https://doi.org/10.1098/rstb.2013.0344 PMID: 24958918

32. Glémin S. Mating systems and the efficacy of selection at the molecular level. *Genetics*. 2007; 177 (2):905–916. https://doi.org/10.1534/genetics.107.073601 PMID: 17954924

33. Roze D. Background selection in partially selfing populations. *Genetics*. 2016; 203(2):937–957. https://doi.org/10.1534/genetics.116.187955 PMID: 27075726

34. Zhang L, Reifová R, Halenková Z, Gompert Z. How important are structural variants for speciation? *Genes*. 2021; 12(7):1084. https://doi.org/10.3390/genes12071084 PMID: 34356100

35. Phillips PC. Waiting for a compensatory mutation: phase zero of the shifting-balance process. *Genetics Research*. 1996; 67(3):271–283. https://doi.org/10.1017/S0016672300033759 PMID: 8690275

36. Landry CR, Wittkopp PJ, Taubes CH, Ranz JM, Clark AG, Hartl DL. Compensatory cis-trans evolution and the dysregulation of gene expression in interspecific hybrids of *Drosophila*. *Genetics*. 2005; 171 (4):1813–1822. https://doi.org/10.1534/genetics.105.047449 PMID: 16143608

37. Barrière A, Gordon KL, Ruvinsky I. Coevolution within and between regulatory loci can preserve promoter function despite evolutionary rate acceleration. *PloS Genetics*. 2012. https://doi.org/10.1371/journal.pgen.1002961 PMID: 23028368

38. Mack KL, Nachman MW. Gene regulation and speciation. *Trends in Genetics*. 2017; 33(1):68–80. https://doi.org/10.1016/j.tig.2016.11.003 PMID: 27914620

39. Bateson W. Heredity and variation in modern lights. *Darwin and modern science*. 1909.

40. Dobzhansky T. Studies on hybrid sterility. II. Localization of sterility factors in *Drosophila pseudoobscura* hybrids. *Genetics*. 1936; 21(2):113. https://doi.org/10.1093/genetics/21.2.113 PMID: 17246786

41. Muller H. Isolating mechanisms, evolution, and temperature. In: *Biol. Symp. vol. 6*; 1942. p. 71–125.

42. Orr HA. The population genetics of speciation: the evolution of hybrid incompatibilities. *Genetics*. 1995; 139(4):1805–1813. https://doi.org/10.1093/genetics/139.4.1805 PMID: 7789779

43. Gavrilets S. Fitness Landscapes and the Origin of Species (MPB-41). Princeton University Press; 2004.

44. Charlesworth B. Evolutionary rates in partially self-fertilizing species. *The American Naturalist*. 1992; 140(1):126–148. https://doi.org/10.1086/285406 PMID: 19426068

45. Kimura M. The role of compensatory neutral mutations in molecular evolution. *Journal of Genetics*. 1985; 64(1):7–19. https://doi.org/10.1007/BF02923549

46. Stephan W. The rate of compensatory evolution. *Genetics*. 1996; 144(1):419–426. https://doi.org/10.1093/genetics/144.1.419 PMID: 8878705

47. Kimura M, King JL. Fixation of a deleterious allele at one of two “duplicate” loci by mutation pressure and random drift. *Proceedings of the National Academy of Sciences*. 1979; 76(6):2858–2861. https://doi.org/10.1073/pnas.76.6.2858 PMID: 288072

48. Galassi M, Davies J, Theiler J, Gough B, Jungman G, Alken P, et al. GNU scientific library. *Reference Manual Edition 14, for GSL Version 14*. 2003.

49. Charlesworth B, Morgan M, Charlesworth D. The effect of deleterious mutations on neutral molecular variation. *Genetics*. 1993; 134(4):1289–1303. https://doi.org/10.1093/genetics/134.4.1289 PMID: 8375663

50. Haller BC, Messer PW. SLIM 3: Forward Genetic Simulations Beyond the Wright-Fisher Model. *Molecular Biology and Evolution*. 2019; 36(3):632–637. https://doi.org/10.1093/molbev/msy228 PMID: 30517680
51. Tange O, et al. Gnu parallel-the command-line power tool. The USENIX Magazine. 2011; 36(1):42-47.

52. Mérot C, Oomen R, Tigano A, Wellenreuther M. A roadmap for understanding the evolutionary significance of structural genomic variation. Trends in Ecology and Evolution. 2022; 35:561-572.

53. Eyre-Walker A, Keightley PD. The distribution of fitness effects of new mutations. Nature Reviews Genetics. 2007; 8(88):610-618. https://doi.org/10.1038/nrg2146 PMID: 17637733

54. Glémin S. Extinction and fixation times with dominance and inbreeding. Theoretical Population Biology. 2012; 81:310-316. https://doi.org/10.1016/j.tpb.2012.02.006 PMID: 22425679

55. Delignette-Muller ML, Dutang C, Pouillot R, Denis JB, Delignette-Muller MML. Package 'fitdistrplus'; 2015.

56. Kimura M. Average time until fixation of a mutant allele in a finite population under continued mutation pressure: Studies by analytical, numerical, and pseudo-sampling methods. Proceedings of the National Academy of Sciences. 1980; 77(1):522-526. https://doi.org/10.1073/pnas.77.1.522

57. Glémin S, Ronfort J. Adaptation and maladaptation in selfing and outcrossing species: new mutations versus standing variation. Evolution. 2013; 67(1):225-240. https://doi.org/10.1111/j.1558-5646.2012.01778.x PMID: 23289574

58. Glémin S. Balancing selection in self-fertilizing populations. Evolution. 2021; 75(5):1011-1029. https://doi.org/10.1111/evo.14194 PMID: 33675041

59. Stebbins GLS. Variation and Evolution in Plants. Columbia University Press; 1950.

60. Stebbins GLS. Evolution above the Species Level. Harvard University Press; 1974.

61. Orr HA. Dobzhansky, Bateson, and the Genetics of Speciation. Genetics. 1996; 144:1331-1335. https://doi.org/10.1093/genetics/144.4.1331 PMID: 8978022

62. Presgraves DC. Darwin and the origin of interspecific genetic incompatibilities. The American Naturalist. 2010; 176(S1):S45-S60. https://doi.org/10.1086/657058 PMID: 21043780

63. Cutter AD. The polymorphic prelude to Bateson-Dobzhansky-Muller incompatibilities. Trends in ecol- ogy & evolution. 2012; 27(4):209-218. https://doi.org/10.1016/j.tree.2011.11.004 PMID: 22154508

64. Maya-Lastra CA, Eaton DA. Genetic incompatibilities do not snowball in a demographic model of speci- ation. bioRxiv. 2021.

65. Corbett-Detig RB, Zhou J, Clark AG, Hartl DL, Ayroles JF. Genetic incompatibilities are widespread within species. Nature. 2013; 504(7478):135-137. https://doi.org/10.1038/nature12678 PMID: 24196712

66. Snoek LB, Orbidans HE, Stastna JJ, Aartse A, Rodriguez M, Riksen JA, et al. Widespread genomic incompatibilities in *Caenorhabditis elegans*. G3: Genes, Genomes, Genetics. 2014; 4(10):1813-1823. https://doi.org/10.1534/g3.114.013151 PMID: 25128438

67. Chae E, Bomblies K, Kim ST, Karelina D, Zaidem M, Ossowski S, et al. Species-wide genetic incompatibility analysis identifies immune genes as hot spots of deleterious epistasis. Cell. 2014; 159(6):1341-1351. https://doi.org/10.1016/j.cell.2014.10.049 PMID: 25467443

68. Gavrilets S. A dynamical theory of speciation on holey adaptive landscapes. The American Naturalist. 1999; 154(1):1-22. https://doi.org/10.1086/303217 PMID: 29587497

69. Lynch M. The genetic interpretation of inbreeding depression and outbreeding depression. Evolution. 1991; 45(3):622-629. https://doi.org/10.2307/2409915 PMID: 28568822

70. Simon A, Bierne N, Welch JJ. Coadapted genomes and selection on hybrids: Fisher's geometric model explains a variety of empirical patterns. Evolution Letters. 2018; 2(5):472-498. https://doi.org/10.1002/evl3.66 PMID: 30283696

71. Dagilis AJ, Kirkpatrick M, Bolnick DI. The evolution of hybrid fitness during speciation. PLoS genetics. 2019; 15(5):e1008125. https://doi.org/10.1371/journal.pgen.1008125 PMID: 31059513

72. Whitlock MC, Ingvarsson PK, Hatfield T. Local drift load and the heterosis of interconnected populations. Heredity. 2000; 84(44):452-457. https://doi.org/10.1046/j.1365-2540.2000.00693.x PMID: 10849069

73. Glémin S, Ronfort J, Bataillon T. Patterns of Inbreeding Depression and Architecture of the Load in Sub-divided Populations. Genetics. 2003;(4):2193-2212. https://doi.org/10.1093/genetics/165.4.2193 PMID: 14704197

74. Ingvarsson PK, Whitlock MC. Heterosis increases the effective migration rate. Proceedings of the Royal Society of London Series B: Biological Sciences. 2000; 267(1450):1321-1326. https://doi.org/10.1098/rspb.2000.1145 PMID: 10972127

75. Harkness A, Brandvain Y, Goldberg EE. The evolutionary response of mating system to heterosis. Journal of Evolutionary Biology. 2019; 32:476-490. https://doi.org/10.1111/jeb.13430 PMID: 30779390

76. Bierne N, Lenormand T, Bonhomme F, David P. Deleterious mutations in a hybrid zone: can mutational load decrease the barrier to gene flow? Genetics Research. 2002; 80(3):197-204. https://doi.org/10.1017/S001667230200592X
77. Stathos A, Fishman L. Chromosomal rearrangements directly cause underdominant F1 pollen sterility in *Mimulus lewisii*Mimulus cardinalis hybrids. Evolution. 2014; 68(11):3109–3119. https://doi.org/10.1111/evo.12503 PMID: 25125144

78. Rieseberg LH. Chromosomal rearrangements and speciation. Trends in ecology & evolution. 2001; 16 (7):351–358. https://doi.org/10.1016/S0169-5347(01)02187-5 PMID: 11403867

79. Signor SA, Nuzhdin SV. The evolution of gene expression in cis and trans. Trends in Genetics. 2018; 34(7):532–544. https://doi.org/10.1016/j.tig.2018.03.007 PMID: 29680748

80. Hodgins KA, Yeaman S. Mating system impacts the genetic architecture of adaptation to heteroge- neous environments. New Phytologist. 2019; 224(3):1201–1214. https://doi.org/10.1111/nph.16186 PMID: 31505030

81. Le Thierry d'Ennequin M, Toupane B, Robert T, Godelle B, Gouyon PH. Plant domestication: a model for studying the selection of linkage. Journal of Evolutionary Biology. 1999; 12(6):1138–1147. https:// doi.org/10.1046/j.1420-9101.1999.00115.x