Negative Coupling: The Coincidence of Premating Isolating Barriers Can Reduce Reproductive Isolation

Model

We consider a two-island population genetics model of haploids with non-overlapping generations. There are four autosomal diallelic loci, two coding for mating traits, T^A and T^B, and two coding for mate preferences, P^A and P^B (as shown in Fig. 1A). Superscripts A and B here refer to the two sets of preference/trait loci that can independently contribute to assortative mating.

The trait loci are subject to divergent viability selection in males (i.e., are “magic traits” [Gavrilets 2014]; as have been found in nature [Servedio et al. 2011]), such that alleles $T_k^A$ and $T_k^B$ are locally advantageous in population k (where k = 1 or 2) (Fig. 1C). Females with allele $P_k^A$ (resp. with allele $P_k^B$) prefer to mate with males with the trait allele $T_k^A$ (resp. allele $T_k^B$), regardless of which population they are in at the time (Fig. 1B). We assume that alleles $P_k^A$ and $P_k^B$ are initially predominant in population k, generating divergent sexual selection at the T^A and T^B loci. There are thus 16 genotypes, $P_1^A{P}_1^B{T}_1^A{T}_1^B$, $P_1^A{P}_1^B{T}_1^A{T}_2^B$, $P_1^A{P}_1^B{T}_2^A{T}_1^B$, $P_1^A{P}_1^B{T}_2^A{T}_2^B$, through $P_2^A{P}_2^B{T}_2^A{T}_2^B$, the frequencies of which are denoted by x_1,k through x_16,k respectively, where the second subscript k again denotes the population. In the following, we will specifically focus on the alleles $P_2^A$, $P_2^B$, $T_2^A$, and $T_2^B$, whose frequencies in population k will be denoted by $p_{2,k}^A$, $p_{2,k}^B$, $t_{2,k}^A$, and $t_{2,k}^B$.

The life cycle consists of migration, viability selection, mating, and the production of zygotes, without genetic drift. After migration, the genotypic frequencies are $x_{i,k}^{\ast }=(1-m_k)x_{i,k}+m_k{x}_{i,l}$, where l = 2 when k = 1 and l = 1 when k = 2. Parameter m_k denotes the proportion of the population k that consists of migrants after migration has occurred (backward migration rate). Unless stated otherwise, we assume symmetric migration such that m₁ = m₂ = m.

Males with alleles $T_k^A$ (resp. $T_k^B$) have a selective advantage due to local adaptation, with relative fitness $1+s_k^{\mathrm{A}}$ (resp. $1+s_k^B$) in population k. After viability selection, the genotypic frequencies in males are:

where

Here, ${\delta }_{i,k}^{\mathrm{A}}=1$ if genotype i has allele $T_k^A$ and ${\delta }_{i,k}^A=0$ otherwise. Likewise, ${\delta }_{i,k}^B=1$ if genotype i has allele $T_k^B$ and ${\delta }_{i,k}^B=0$ otherwise. Unless stated otherwise, we consider symmetric viability selection and the same strength of selection on both loci T^A and T^B, such that $s_1^A=s_2^A=s_1^B=s_2^B=s$.

Nonrandom mating occurs following preference-trait rules (Kopp et al. 2018). In population k, females carrying allele $P_1^A$ (resp. $P_2^A$, $P_1^B$, and $P_2^B$) are $1+{\alpha }_k^{A,1}$ (resp. $1+{\alpha }_k^{A,2}$, $1+{\alpha }_k^{B,1}$, and $1+{\alpha }_k^{B,2}$) times as likely to mate with a $T_1^A$ male (resp. a $T_2^A$, $T_1^B$, and $T_2^B$ male) than with a male of the opposite allele, upon encounter. In population k, the frequency of mated pairs of females with genotype i and males with genotype j is thus:

where

Here, ${\gamma }_{i,j}^{A,1}=1$ (resp. ${\gamma }_{i,j}^{A,2}=1$) if the female genotype i carries allele $P_1^A$ (resp. $P_2^A$) at the preference locus P^A and the male genotype j carries allele $T_1^A$ (resp. $T_2^A$) at the trait locus T^A, and ${\gamma }_{i,j}^{A,1}=0$ (resp. ${\gamma }_{i,j}^{A,2}=0$) otherwise. The values of ${\gamma }_{i,j}^{B,1}$ and ${\gamma }_{i,j}^{B,2}$ are defined analogously on the basis of allele matching between the P^B locus in females and the T^B locus in males. The normalization in the denominator ensures that all females have equal mating success (meaning that sexual selection acts only on males). We assume symmetrical preference strengths such that ${\alpha }_k^{A,1}={\alpha }_k^{A,2}={\alpha }_k^A$ and ${\alpha }_k^{B,1}={\alpha }_k^{B,2}={\alpha }_k^B$. Additionally, unless stated otherwise, we assume that mating preference strengths are the same in the two populations, such that ${\alpha }_1^A={\alpha }_2^A={\alpha }^A$ and ${\alpha }_1^B={\alpha }_2^B={\alpha }^B$.

Mating between haploids leads to the formation of diploid zygotes, which then undergo meiosis and recombination to produce the next generation of haploids. To assess the effect of recombination on trait and preference divergence, we assume that preference loci may be physically linked with each other, with recombination rate r_P between loci P^A and P^B. Likewise, we assume that trait loci may be physically linked with each other, with recombination rate r_T between loci T^A and T^B. For simplicity, however, we assume that the preference loci are not physically linked to the mating trait loci (e.g., preference and trait loci are located on different chromosomes). In a supplementary analysis, we relax this assumption and consider that there is physical linkage between all loci, so that the gene order is P^BP^AT^AT^B. In another supplementary analysis, we consider P^AT^AT^BP^B as an alternative gene order.

Simulation Experiments

We assume secondary contact (a situation that is often found in nature; Harrison 1993), where population 2 is initially fixed for the $P_2^A{P}_2^B{T}_2^A{T}_2^B$ genotype, whereas population 1 has the $P_1^A{P}_1^B{T}_1^A{T}_1^B$ genotype present at a high frequency (=0.95, to avoid strictly symmetrical starting conditions), with all other genotypes present according to their frequencies at linkage equilibrium.

To assess the long-term equilibrium of the system (e.g., for local allele frequencies and LDs), we run numerical iterations until the change in each allelic frequency is less than 10⁻⁶ per generation. We report normalized LDs spanning the range $\left[-1,+1\right]$. For instance, LD between P^A and T^B in population k is normalized relative to $\min \{p_{2,k}^A(1-t_{2,k}^B),(1-p_{2,k}^A)t_{2,k}^B\}$ when it is positive, and relative to $\max \{-p_{2,k}^A{t}_{2,k}^B,-(1-p_{2,k}^A)$ $(1-t_{2,k}^B)\}$ when it is negative (Lewontin 1964). LDs referred to in the text are always measured within each population, as opposed to across both populations.

To understand the effect of a second set B of preferences and traits on divergence at an “original” set A, we performed comparative simulations with the effect of the preference locus P^B turned off (hereafter called simulation “without P^B”, where α^A = α and α^B = 0) or turned on (hereafter called simulations “with P^B” where α^A = α^B = α), the latter representing the presence of a second potential reproductive barrier. Note that, in the simulations without P^B, we nevertheless included the T^B locus (alongside T^A and P^A), to keep the strength of divergent ecological selection the same in both sets of simulations.

Reproductive Isolation

To assess the effect of two sets of preferences and traits (sets A and B) on reproductive isolation relative to that with only one “original” set (set A), we numerically estimate the strength of the overall barrier to gene flow at a neutral locus unlinked to all preference and trait loci (Bengtsson 1985; Barton and Bengtsson 1986; Westram et al. 2022). To do so, we calculate the effective migration rate, which is inversely proportional to the total strength of barriers to gene flow, from the rate of convergence to equilibrium at the neutral locus (Akerman and Bürger 2014). More precisely, we consider an additional autosomal diallelic locus N that is neutral and unlinked to all other loci. Starting from the equilibrium state at all other loci, we assume there is initially complete divergence at the neutral locus between populations (with alleles N₁ and N₂ fixed in population 1 and 2, respectively). We then average, over 1000 generations, the change in allele frequency at the neutral locus relative to the current divergence in allele frequency between populations at this locus. This metric corresponds to the rate of convergence to equilibrium (characterized by an allelic frequency of 0.5 in both populations), and thus corresponds to the effective migration rate, m_eff, which is inversely proportional to reproductive isolation, that is, to the strength of the overall barrier to gene flow. We thus define as a measure of reproductive isolation (Barton and Bengtsson 1986; Westram et al. 2022)

Note that gene flow at neutral loci that are physically linked to divergent loci will be lower than that at the unlinked neutral locus considered here (Barton and Bengtsson 1986; Westram et al. 2022).

The C++ and R scripts used for the simulations and the figures are available in a Zenodo repository (doi:10.5281/zenodo.10400557; zenodo.org/records/10400557).

Overview of the Key Results

Although our population genetics model is as simple as possible—with each premating isolating barrier encoded by just two diallelic loci (as shown in Fig. 1)—contrasting and complex outcomes emerge depending on the recombination rates among loci. In the following sections, we exhaustively describe six qualitatively different regimes. For the convenience of the reader, here we present the key results that we feel are most worth noting.

More importantly, our model shows that reproductive isolation can decrease following the coincidence of distinct reproductive barriers (a situation we refer to as “negative coupling”). Indeed, the presence of a second set of preference and trait loci (encoding an additional premating isolating barrier; see equilibrium state in Fig. 2) can, under some conditions, lead to decreased divergence at the first set of loci (encoding an original premating isolating barrier) and to reduced LD between the male trait loci (Figs. 3 and 4). This may induce a decrease in the strength of reproductive isolation (Fig. 5).

This result can primarily be attributed to two effects. First, the homogenizing effect of sexual selection on allelic divergence (described above; as shown in a two-locus model in Servedio and Bürger 2014) can be reinforced when there are two sets of preferences and traits. Second, under some combinations of recombination rates (with higher recombination between preference loci than between trait loci), there can be more recombinant females with mismatched preferences than recombinant males with mismatched traits. The presence of these recombinant females will then generate high mating success for these recombinant males, creating a feedback whereby LD is reduced among both preference and trait loci, and may even become negative. The combinations of these and other effects, described in more detail below, highlight that it is incorrect to assume that coincident sets of premating barriers will necessarily lead to “positive coupling” (i.e., an increase in overall reproductive isolation). Negative coupling can result instead.

A Second Set of Preference and Trait Loci Impacts Divergence According to Six Different Regimes

We first describe how an additional preference locus P^B affects divergence both at the trait locus T^B (due to sexual selection induced by P^B), and at loci P^A and T^A, which form the “original” set of preference and trait loci. Recall that both T^A and T^B are “magic traits” that are also under ecologically divergent selection.

While divergence is generally maintained at all loci in the presence of P^B (Fig. 2), it is not necessarily higher than when P^B is absent (Fig. 3). First of all, the presence of P^B reduces divergence at the trait locus T^B (Figs. 3C and 4; for an exception, see Supplemental Fig. S1), which occurs due to the homogenizing effect of sexual selection described earlier (Servedio and Bürger 2014). By contrast, the impact of P^B on divergence at the “original” set of preference and trait loci, P^A and T^A, depends on the recombination rates between the preference loci (r_P) and between the trait loci (r_T), and results from the interplay of multiple evolutionary forces. In this subsection, we focus on describing the evolutionary forces that dominate under a variety of regimes that depend on the combination of parameters (r_P, r_T). Note that, although Figures 2, 3, and 5 focus on a single choosiness value (α = 5), corresponding to moderately strong preferences, the impact of P^B on divergence at the P^A and T^A loci depends greatly on the level of choosiness (Fig. 4).

Overall, Figures 2 and 3 allow us to distinguish six different regimes, indicated by roman numerals in the figures and in the last column of Table 1. The six regimes depend on five primary effects of the presence of P^B, which are given in the first column of Table 1 (note that regimes IIIa and IIIb depend on the same primary effect).

In regime I (low r_P and high r_T), the presence of the preference locus P^B leads to increased divergence at both the T^A and P^A loci (Fig. 3A,B). This is because the presence of P^B increases LD between these loci (Figs. 3D and 4A) (primary effect [1] in Table 1). Indeed, there is positive LD both between T^A and T^B (Fig. 2G; due to sexual selection induced by $P_1^A{P}_1^B$ and $P_2^A{P}_2^B$ females) and between P^A and P^B (Fig. 2H) (as explained in Box 1). These positive LDs imply that $P_1^A$ females (resp. $P_2^A$ females) are more likely to mate with $T_1^A$ males (resp. $T_2^A$ males) when P^B is present. Increased LD between P^A and T^A, in turn, favors the maintenance of high divergence at the P^A locus through stronger indirect divergent selection on P^A. This thereby reduces the homogenizing effect of sexual selection induced by P^A on T^A, increasing divergence at the T^A locus as well.

In regime II (low r_P and lower r_T than in regime I), the presence of P^B still increases divergence at the P^A locus, but decreases divergence at the T^A locus (Figs. 3A,B and 4B). Strong LD is again maintained both between the trait loci (Fig. 2G) and between the preference loci (Fig. 2H). In this case, though, these strong LDs increase the maximum level of nonrandom mating and therefore the strength of homogenizing sexual selection acting on the T^AT^B gene complex (primary effect [2] in Table 1; through sexual selection acting on the $T_1^A{T}_1^B$ and $T_2^A{T}_2^B$ genotypes), which explains why divergence at the T^A locus decreases (Fig. 3A). The key difference with regime I is that, in regime II, r_T is lower and LD between the trait loci is thus higher (darker red in Fig. 2G; with $T_1^A{T}_1^B$ and $T_2^A{T}_2^B$ genotypes in high frequency); therefore, homogenizing sexual selection acts more effectively on the T^AT^B gene complex. At the same time, however, the presence of P^B actually increases divergence at the P^A locus (Fig. 3B), because the increased LD between P^A and T^A causes stronger indirect divergent selection on P^A (as in regime I).

In regimes IIIa, IIIb, IV, and V, the recombination rate r_T becomes low enough relative to r_P that “recombinant” $T_1^A{T}_2^B$ and $T_2^A{T}_1^B$ males become less frequent than “recombinant” $P_1^A{P}_2^B$ and $P_2^A{P}_1^B$ females that prefer them, and therefore benefit from high mating success. As explained in Box 1, this effect can decrease LD between the trait loci, and ultimately also between the preference loci. With this background, we can explore the distinction between these four regimes, in which the presence of P^B either decreases divergence at both the T^A and P^A loci (regimes IIIa, IIIb, and IV) or slightly increases divergence at both the T^A and P^A loci (regime V).

The hallmark of regimes IIIa and IIIb (low r_P and low r_T) is negative LD between the preference loci (Figs. 2H and 4C,D) (primary effect [3] in Table 1), which arises once sexual selection favoring “recombinant” males lowers the initially positive LD between the trait loci (as explained in Box 1). In turn, this negative LD between the preference loci ultimately generates negative LD between the trait loci T^A and T^B, and also between P^A and T^B due to sexual selection and mate choice (Fig. 2F,G, and hatched areas in Fig. 3E,F). Through these negative LDs, viability selection on T^B indirectly reduces divergence at both the T^A and P^A loci (Fig. 3A,B). This effect is particularly strong in regime IIIa, where divergence eventually vanishes at all loci (Fig. 2A–D; allele frequencies reach 0.5 in both populations). As a result, in regime IIIa, even weak asymmetrical viability selection on T^A or weak asymmetrical preference strength induced by P^A between populations can lead to the loss of genetic variation at the T^A locus (e.g., for $s_1^A=1.01s_2^A$ or ${\alpha }_1^A=1.01{\alpha }_2^A$ in Supplemental Figs. S2 and S3, respectively). In regime IIIb, this process does not have as strong an effect as in regime IIIa and reduces divergence at all loci without making it vanish.

In regime IV (for high r_P and r_T), the presence of P^B decreases divergence at the T^A and P^A loci (Fig. 3A,B). LD between all loci remains positive, but the presence of the preference locus P^B decreases LD between loci T^A and T^B (Figs. 2G,H, 3F, and 4E) (primary effect [4] in Table 1). In turn, the strength of indirect viability selection on T^A, and consequently also on P^A, is reduced and thus leads to a decrease in divergence at these loci.

Finally, in regime V (for high r_P, and low r_T relative to r_P), the presence of P^B very slightly increases divergence at the T^A and P^A loci (Fig. 3A,B). The same decrease in LD between loci T^A and T^B occurs as in regime IV, but another primary effect prevails: the presence of the preference locus P^B weakens the strength of homogenizing sexual selection on the T^A locus, and hence on the P^A locus indirectly (Figs. 3A,B and 4F) (primary effect [5] in Table 1). This effect occurs because homogenizing sexual selection acting on the $T_1^A{T}_1^B$ and $T_2^A{T}_2^B$ genotypes is particularly weak (the variation in male mating success is low because $P_1^A{P}_1^B$ and $P_2^A{P}_2^B$ females are not overrepresented as shown in Fig. 2H). As a result, divergence at the T^A and P^A loci increases very slightly when P^B is present (Fig. 3A,B). This effect is best considered negligible.

We obtain qualitatively the same results when we consider that viability selection on T^A, preference strength induced by P^A, or migration rate is weakly asymmetric between populations (e.g., for $s_1^A=1.01s_2^A$, ${\alpha }_1^A=1.01{\alpha }_2^A$, and m₁ = 1.01m₂) (Supplemental Figs. S2–S4), except for the loss of genetic variation in regime IIIa mentioned above. Likewise, the same primary effects occur when viability selection is weaker than that implemented in the figures (e.g., for s = 0.05 instead of s = 0.5) (Supplemental Fig. S5), when viability selection acts in both sexes (Supplemental Fig. S6), or when the preference and trait loci are physically linked (with recombination occurring at rate 0.1, instead of free recombination) (Supplemental Fig. S1). Note that in simulations where viability selection is weak (s = 0.01), we also assumed weak choosiness (α = 1), so that frequency-dependent sexual selection does not lead to a loss of allelic divergence (Supplemental Fig. S5). Under these conditions, we obtain qualitatively the same results.

A Second Set of Preference and Trait Loci Can Decrease Reproductive Isolation

The way the preference locus P^B influences divergence at the T^A and P^A loci has a notable effect on how P^B affects reproductive isolation, as assessed from the effective migration rate at a neutral unlinked locus (Fig. 5). When the presence of P^B increases divergence at both T^A and P^A, reproductive isolation tends to increase, constituting a case of “positive coupling” (traditionally referred to simply as “coupling”) (Butlin and Smadja 2018). We find, however, that this is the case only in regime I, that is, if preference loci are tightly linked but trait loci are essentially unlinked. Even under these conditions, increased divergence at the T^A and P^A loci does not necessarily increase reproductive isolation (as highlighted by the comparison between Figs. 3A,B and 5B) (this is also the case in regime V where divergence at the T^A and P^A loci slightly increases and reproductive isolation is reduced). This is because reproductive isolation can be reduced due to the decrease in divergence at the T^B locus caused by homogenizing sexual selection induced by P^B (Fig. 3C).

When, in contrast, the presence of P^B decreases divergence at both the T^A and P^A loci (in regimes II, IIIa, IIIb, and IV), reproductive isolation can be seen to decrease (Fig. 5). We thus, surprisingly, find broad ranges of parameter values for which the coincidence of reproductive barriers reduces overall reproductive isolation. This situation can be referred to as “negative coupling.”