Hybrid fitness effects modify fixation probabilities of introgressed alleles

Abstract Hybridization is a common occurrence in natural populations, and introgression is a major source of genetic variation. Despite the evolutionary importance of adaptive introgression, classical population genetics theory does not take into account hybrid fitness effects. Specifically, heterosis (i.e. hybrid vigor) and Dobzhansky–Muller incompatibilities influence the fates of introgressed alleles. Here, we explicitly account for polygenic, unlinked hybrid fitness effects when tracking a rare introgressed marker allele. These hybrid fitness effects quickly decay over time due to repeated backcrossing, enabling a separation-of-timescales approach. Using diffusion and branching process theory in combination with computer simulations, we formalize the intuition behind how hybrid fitness effects affect introgressed alleles. We find that hybrid fitness effects can significantly hinder or boost the fixation probability of introgressed alleles, depending on the relative strength of heterosis and Dobzhansky–Muller incompatibilities effects. We show that the inclusion of a correction factor (α, representing the compounded effects of hybrid fitness effects over time) into classic population genetics theory yields accurate fixation probabilities. Despite having a strong impact on the probability of fixation, hybrid fitness effects only subtly change the distribution of fitness effects of introgressed alleles that reach fixation. Although strong Dobzhansky–Muller incompatibility effects may expedite the loss of introgressed alleles, fixation times are largely unchanged by hybrid fitness effects.

where g t is the fraction of possible interactions that are possible in hybrid genomes t generations after admixture, and the -1 term ensures that δ t ≤ 0. Representing the expected number of DMIs n 2 p I in the F1 hybrid with m and making use of the Taylor series, Equation S1 becomes for small c. As DMI effects are expected to be polygenic (True et al. 1996;Presgraves and Meiklejohn 2021), it is valid to assume that the effect size of an individual DMI (c) is small. Therefore, terms of O(c 2 ) and greater can be ignored, and Equation S2 is approximated by −mcg t . This is equivalent to an additive model with initial strength (1 − c) m − 1 ≈ −mc. Then,Equation 5 in the main text is obtained by representing the initial strength of DMI effects −mc with δ 1 and substituting g t with the decay rate of possible interactions.
Heterosis effects are assumed to arise from the masking of recessive deleterious alleles, i.e., the dominance theory. The same reasoning as above for DMI effects applies to the decay of heterosis effects. The number of possibly masked recessive alleles is directly proportional to the amount of introgressed DNA in hybrid genomes, i.e., n, n/2, n/4, and n/8 in F1, BC1, BC2, and BC3 hybrids, respectively. Assuming an equal probability for each introgressed allele to mask a recessive deleterious allele (p M ), an equal effect size z for each masking event, and an additive model, we can calculate the strength of heterosis effects at time t (η t ) with respect to their initial strength (np M z = η 1 ) (Equation 4 in the main text).

Stochastic dilution of introgressed DNA
To simulate the stochastic dilution of introgressed DNA, we explicitly track the fraction of introgressed DNA present in an individual at a given time t relative to the amount of introgressed DNA present in the F1 hybrid. For reproduction, two parents are randomly picked with a probability equivalent to their fitness. To avoid selfing, it is ensured that the parents are different individuals. The total amount of introgressed DNA in the child is then the sum of the amount of introgressed DNA it inherited from each parent. Assuming a genome size of 3.2 Gb, the genome can be divided into 32 chunks of 100 Mb, given that one centimorgan is approximately 1 Mb. For each chunk, a diploid individual has two haplotypes -one from each of their parents -and the probability that a parent passes on a specific haplotype to their child is 0.5. For this reason, the probability that a specific set of haplotypes that trace to the same grandparent is inherited to the next generation is modeled with a binomial distribution with 32 trials and a probability of success of 0.5 (the effects of different genome sizes can be explored using the code available at https: //github.com/LachanceLab/introgression_theory). Given the parents relative amount of introgressed DNA is Y p,i , the probability distribution of the relative amount of introgressed DNA in the offspring (Y o ) is then The strength of heterosis and DMI effects in a given offspring at time t is then given by , respectively. In each generation, the entire population is replaced by their offspring, i.e., the Wright-Fisher model.
Given that we consider a rare introgression event, that is, the hybridization of a single individual from the recipient and donor population, the chances of individuals with introgressed DNA mating are negligible in large populations during the first few generations after admixture; hence, the effects on the decay of HFEs are negligible, too. However, in small populations, it is more likely that two individuals with introgressed DNA are mating during the first few generations, slowing the dilution of introgressed DNA, and hence also slowing the decay of HFEs. This would lead to an inflation of the impact of HFEs on the fixation probability of the introgressed B allele in small populations. To counter this effect and allow an assessment of the stochastic dilution of introgressed DNA, we only allow the inheritance of introgressed DNA alongside the introgressed B allele.

Fixation probabilities from a branching process
Here, we will derive an expression for fixation probability under consideration of HFEs using a branching process (Taylor and Karlin 1998). We suppose an infinitely large, randomly mating population of diploid individuals and that the tracked B allele is initially present as a single copy, i.e., the hybridization of a single individual from the donor and recipient population. Furthermore, we assume that each parent allele independently contributes k offspring alleles to the next generation, where k is Poisson distributed. Then, the probability generating function for the number of offspring alleles in generation t is: where λ t is the hybrid fitness w AB,t in generation t. Due to the time-varying hybrid fitness, compounding the probability generating function is time-heterogeneous, and the probability generating function of the number of offspring alleles after t generations is Because the only two absorbing states are extinction and fixation of an allele, they are implicit of each other, and the fixation probability of an allele is equivalent to its survival probability (u). The probability of extinction by generation t is G t (0), and thus the survival probability, that is, there is at least one copy in generation t, is

Different epistasis models leading to DMI effects
In the main text, we focused on DMI effects arising from epistatic interactions between unlinked semidominant introgressed and semidominant non-introgressed alleles (Equation 5). As pointed out in the main text, while the dynamics are identical if introgressed alleles are dominant, no DMI effects are observed in the scenario of a rare introgression event if introgressed alleles are recessive. However, the dynamics differ if non-introgressed alleles are recessive or dominant. Generally, the overall conclusions drawn in the main text regarding the impact of HFEs on the fate of introgressed alleles remain unchanged for the different models of epistatic interactions, but the dynamics slightly differ. Here, we will discuss the dynamics under these different models of epistatic interactions.
To distinguish the dynamics under the different epistasis models, we refer to the model discussed in the main text as the semidominancedominance epistasis model (δ t ), to the model of interactions between recessive non-introgressed alleles and (semi-)dominant introgressed alleles as the recessive-dominance epistasis model (δ rd,t ), and to the model of interactions between dominant non-introgressed alleles and (semi-)dominant introgressed alleles as the dominance-dominance epistasis model (δ dd,t ). These three epistasis models differ in the decay rates and onset of DMI effects.
The logic of the decay of DMI effects under the recessive-dominance and dominance-dominance epistasis model is the same as for the semidominance-dominance epistasis model, which is described in detail above. However, their decay functions (g t ) are different. As a reminder, we presume the introgression of a single set of chromosomes with n loci, and the amount of introgressed DNA is assumed to halve every generation. We also presume polygenic DMI effects as per empirical results of True et al. (1996) and Presgraves and Meiklejohn (2021). Hence, each DMI makes a small contribution to the overall DMI effects, and the strength of DMI effects decays with the dilution of introgressed DNA in hybrid genomes. Above, we showed that in this case, the multiplicative model of DMI effects is approximated by the first term of its Taylor series, which is equivalent to an additive model.
Under the recessive-dominance epistasis model, interactions between introgressed and non-introgressed loci do not occur in the F1 hybrid ( Figure S4). This is because the F1 hybrid cannot be homozygous for non-introgressed loci. Due to repeated backcrossing within the recipient population, hybrids become homozygous for non-introgressed loci so that DMI effects become possible. For this reason, DMI effects are observed first in BC1 hybrids (δ rd,2 ). In BC1 hybrids, 1 2 n introgressed loci can interact with the 1 2 n non-introgressed loci for which BC1 hybrids are homozygous, amounting to 1 2 n × 1 2 n = 1 4 n 2 possible interactions. Extending this logic to BC2 and BC3 hybrids, 3 4 n × 1 4 n = 3 16 n 2 and 7 8 n × 1 8 n = 7 64 n 2 interactions are possible. Therefore, the decay of DMI effects in the scenario of the recessive-dominance epistasis model is given by: Thus the decay rate under the recessive-dominance epistasis model is the same as under the semidominance-dominance epistasis model, but DMI effects are delayed by one generation under the recessive-dominance (compare Figure 1B & Figure S4). When the initial strengths of DMI effects under the semidominance-dominance and recessive-dominance epistasis models are the same (δ 1 = δ rd,2 ), the compounded DMI effects are identical.
Under the dominance-dominance epistasis model, the number of possible interactions between introgressed and non-introgressed alleles inherited from different parents determines the strength of DMI effects. Interactions between introgressed and non-introgressed alleles inherited from the same parent do not add to DMI effects because of the dominance of non-introgressed loci ( Figure S5). Thus, the number of possible interactions simply halves every generation (i.e., n 2 in F1 hybrids, 1 2 n 2 in BC1 hybrids, etc.). Hence, the first part in the parenthesis of the non-reduced form of Equation 5 describes the decay of DMI effects under this model (δ dd,t ): where δ dd,1 is the strength of DMI effects in the F1 hybrid. Note that the decay rate of DMI effects under the dominance-dominance epistasis model and the decay of heterosis effects is the same. Table S2 summarizes the decay of the expected strength of DMI effects for each epistasis model up to BC4 hybrids. m is the expected number of DMIs in the F1/BC1 hybrid (i.e., the number of possible interactions in the F1/BC1 hybrid times the probability that an interaction is incompatible), and c is the effect size of an individual DMI.
Hybrid fitnesses are similarly influenced by the different epistasis models. When explaining DMI effects with the recessive-dominance epistasis model, DMI effects show the same decay rate but are delayed by one generation compared to under the semidominance-dominance epistasis model ( Figure S6). Because of this, the hybrid breakdown in BC1 hybrids is more severe under the recessive-dominance epistasis model compared to the semidominance-dominance epistasis model, given the same initial strength of DMI effects (δ 1 = δ rd,2 ). In contrast to this, the dynamics are simpler when assuming the dominance-dominance epistasis model to explain DMI effects ( Figure S7). Due to the same decay rates of heterosis and DMI effects, either hybrid vigor (i.e., when η 1 > |δ dd,1 |) or hybrid breakdown (i.e., when |δ dd,1 | > η 1 ) is observed, but no hybrid vigor followed by hybrid breakdown can be observed. Heterosis and DMI effects cancel out under the dominance-dominance epistasis model when heterosis and DMI effects have the same initial strength (η 1 = |δ dd,1 |).
Fixation probabilities are similarly affected by the different epistasis models. The accuracy of Equations 7 & 8, respectively, is robust to the choice of epistasis model (compare Figures 3, S8 & S9). Under the recessive-dominance epistasis model, fixation probabilities are increased when initial heterosis effects are more than 1.62 times as strong as DMI effects (η 1 > 1.62 × |δ rd,2 |; Figure S10). This deviates from the expected condition for weak HFEs, i.e., the fixation probability is increased if heterosis effects are initially 4 3 times as strong as DMI effects (see Results in the main text). This is because initial heterosis effects need to be stronger to overcome DMI effects than under the semidominance-dominance epistasis model since hybrid breakdown in the BC1 and following generations is more severe due to the DMI effects being delayed by one generation. In contrast to this, under the dominance-dominance epistasis model, fixation probabilities are always increased when initial heterosis effects are stronger than DMI effects (η 1 > |δ dd,1 |) and vice versa. This finding is due to identical decay rates under this epistasis model ( Figure S11).
The distribution of fitness effects (DFE) of alleles that reach fixation is similarly transformed by the different epistasis models for DMI effects. Under the recessive-dominance epistasis model, the DFE is shifted to slightly greater selection coefficients when DMI effects are strong (η 1 = 0.0; δ rd,1 = −0.8), i.e., s is 0.0649 compared to 0.0634 under the semidominance-dominance epistasis model ( Figure S12). In contrast to this, under the dominance-dominance epistasis model, the shift towards greater selection coefficients is diminished when DMI effects are strong (η 1 = 0.0; δ dd,1 = −0.8), i.e., s is 0.0633 compared to 0.0634 under the semidominance-dominance epistasis model ( Figure  S13).
Sojourn times are also qualitatively similarly affected by the different epistasis models for DMI effects. Under the recessive-dominance epistasis model, the loss of an allele is accelerated due to stronger hybrid breakdown compared to the semidominance-dominance epistasis model (compare Figure 6D & S14D). Contrarily, under the dominant-dominant epistasis model, sojourn times under the classical model are matched if heterosis and DMI effects have the same strength (η 1 = |δ dd,1 |, Figure S15). Ohta T, Kojima KI. 1968

Table S1
Comparison of fixation probability if the amount of introgressed donor DNA -and hence also HFEs -decays deterministically (i.e., halves every generation) or stochastically (i.e., Equation S3). Fixation probabilities were computed based on 1, 000, 000 WF simulations with s = 0.01. Statistically significant different fixation probabilities are indicated by * (Fisher's exact test, P < 0.05 8 = 0.00625). Strong heterosis (η 1 = 0.8, δ 1 = 0.0, α = 3.67)  The joint distribution of fitness effects and fixation probabilities, i.e., 2αs f (s). HFEs (as encompassed by the parameter α) linearly transform the joint distribution independently from the original shape of the DFE. As represented by the area under the curve, heterosis effects increase the overall likelihood of introgression, while DMI effects diminish it. The original DFE of alleles was modeled with a normal distribution with a mean of −0.001 and a standard deviation of 0.05. Fixation probabilities of the introgressed B allele were estimated based on 100, 000 Wright-Fisher simulations assuming an effective population size of N e = 10, 000 and an initial allele frequency of q 0 = 1/2N e .   Temporal dynamics of different hybrid fitness scenarios where DMI effects arise from recessive-dominance epistasis: A) classical model, B) strong heterosis effects, C) strong DMI effects, and D) strong heterosis and DMI effects. HFEs (greyscale) quickly decay and asymptotically approach zero. Therefore, the overall fitness of heterozygotes that carry the introgressed B allele (blue) is timesensitive. In the absence of HFEs, our model is equivalent to the classical model with constant hybrid fitness determined by the intrinsic selection of the tracked B allele (s = 0.1; black). When η 1 > |δ rd,2 |, heterosis is initially the dominant force, resulting in hybrid vigor until intrinsic selection takes over after approximately eight generations. Conversely, when |δ rd,2 | > η 1 , DMI effects are the stronger force, leading to hybrid breakdown, until intrinsic selection takes over. When η 1 = |δ rd,2 |, heterosis effects are the stronger evolutionary force in the F1 hybrid, inducing hybrid vigor. With the onset of DMI effects in BC1 hybrids, it comes to severe hybrid breakdown. α is the compounded effect size of HFEs (Eq. 6). Temporal dynamics of different hybrid fitness scenarios where DMI effects arise from dominance-dominance epistasis: A) classical model, B) strong heterosis effects, C) strong DMI effects, and D) strong heterosis and DMI effects. HFEs (greyscale) quickly decay and asymptotically approach zero. Therefore, the overall fitness of heterozygotes that carry the introgressed B allele (blue) is timesensitive. In the absence of HFEs, our model is equivalent to the classical model with constant hybrid fitness determined by the intrinsic selection of the tracked B allele (s = 0.1; black). When η 1 > |δ dd,1 |, heterosis is initially the stronger force, resulting in hybrid vigor until intrinsic selection takes over after approximately eight generations. Conversely, when |δ dd,1 | > η 1 , DMI effects are the stronger force, leading to hybrid breakdown, until intrinsic selection takes over. When η 1 = |δ dd,1 |, heterosis effects and DMI effects cancel each other because of the same decay rates. Thus, our model resembles the classical model. α is the compounded effect size of HFEs (Eq. 6). The distribution of fitness effects of introgressed alleles conditioned on fixation when the original DFE is modeled with a normal distribution with mean −0.001 and standard deviation 0.05, and DMI effects arise from recessive-dominance epistasis. Heterosis effects subtly shift the DFE towards smaller selection coefficients, while DMI effects shift it towards slightly greater values. Under this model, the shift is slightly more prominent than under the semidominance-dominance epistasis model due to the stronger hybrid breakdown (compare Figure 5). The mean selection coefficient is shifted from 0.0609 under the classical model to 0.0636 (compared to 0.0622 if semidominance-dominance epistasis is assumed) when η 1 = 0.8 and δ rd,2 = δ sd,1 = −0.8. The DFEs of fixed alleles was inferred from 10 7 Wright-Fisher simulations with selection coefficients drawn from the original DFE (see inset). Simulation parameters: N e = 10, 000 and q 0 = 1/2N e . The distribution of fitness effects of introgressed alleles conditioned on fixation when the original DFE is modeled with a normal distribution with mean −0.001 and standard deviation 0.05, and DMI effects arise from dominance-dominance epistasis. Heterosis effects subtly shift the DFE towards smaller selection coefficients, while DMI effects shift it towards slightly greater values. Under this model, the shift is slightly less prominent than under the semidominance-dominance epistasis model due to the weaker hybrid breakdown (compare Figure 5). The mean selection coefficient is shifted from 0.0610 under the classical model to 0.0633 (compared to 0.0634 if semidominance-dominance epistasis is assumed) when η 1 = 0.0 and δ dd,1 = δ sd,1 = −0.8. Note that the curves for the cases when η 1 = δ dd,1 = 0.0 and when η 1 = |δ dd,1 | = 0.8 overlap because both resemble the classical model. The DFEs of fixed alleles was inferred from 10 7 Wright-Fisher simulations with selection coefficients drawn from the original DFE (see inset). Simulation parameters: N e = 10, 000 and q 0 = 1/2N e . Figure S14 HFEs affect extinction times but not fixation times when DMI effects arise from recessive-dominance epistasis: A) classical model, B) strong heterosis effects, C) strong DMI effects, and D) strong heterosis and DMI effects. Strong heterosis effects prolong the time until the loss of the introgressed B allele, whereas strong DMI effects expedite its loss compared to the classical model. When heterosis effects and DMI effects are of similar strength, the extinction is slightly expedited due to the slower decay of DMI effects. This is because HFEs act on a shorter timescale than on which an allele reaches fixation, illustrating the separation of timescales. The differences between our model and the classical model regarding the extinction time are slightly more distinct under this epistasis model than compared to the semidominance-dominance epistasis model (compare Figure 6D in the main text). Distributions of sojourn times were inferred based on 100, 000 Wright-Fisher simulations with N e = 10, 000, s = 0.01, and q 0 = 1/2N e . Figure S15 HFEs affect extinction times but not fixation times when DMI effects arise from dominance-dominance epistasis: A) classical model, B) strong heterosis effects, C) strong DMI effects, and D) strong heterosis and DMI effects. Strong heterosis effects prolong the time until the loss of the introgressed B allele, whereas strong DMI effects expedite its loss compared to the classical model. This is because HFEs act on a shorter timescale than on which an allele reaches fixation, illustrating the separation of timescales. When heterosis effects and DMI effects are of the same strength, heterosis and DMI effects cancel each other out so that our model resembles the classical model. Distributions of sojourn times were inferred based on 100, 000 Wright-Fisher simulations with N e = 10, 000, s = 0.01, and q 0 = 1/2N e .