Coevolution of male and female mate choice can destabilize reproductive isolation

Sexual interactions play an important role in the evolution of reproductive isolation, with important consequences for speciation. Theoretical studies have focused on the evolution of mate preferences in each sex separately. However, mounting empirical evidence suggests that premating isolation often involves mutual mate choice. Here, using a population genetic model, we investigate how female and male mate choice coevolve under a phenotype matching rule and how this affects reproductive isolation. We show that the evolution of female preferences increases the mating success of males with reciprocal preferences, favouring mutual mate choice. However, the evolution of male preferences weakens indirect selection on female preferences and, with weak genetic drift, the coevolution of female and male mate choice leads to periodic episodes of random mating with increased hybridization (deterministic ‘preference cycling’ triggered by stochasticity). Thus, counterintuitively, the process of establishing premating isolation proves rather fragile if both male and female mate choice contribute to assortative mating.

R eproductive isolation among taxa can be caused by different isolating barriers, such as ecological divergence, sexual isolation, hybrid sterility, and microspatial partitioning 1,2 . While those barriers have been well described, their temporal stability is little studied either theoretically or empirically. Yet, if the barriers that isolate diverged populations are not stable, the consequent recurring hybridization may have profound consequences for species evolution (e.g., speciation 3 , transgressive segregation 4 , adaptive introgression 5 , and genetic swamping 6 ).
Assortative mating-the tendency of individuals of similar phenotype to mate together more often than expected by chanceis widespread in animals 7 and plays a key role in generating premating reproductive isolation 8 . Assortative mating can arise as a by-product of adaptive divergence via temporal or spatial isolation 9 , or can be driven by various behavioral processes 10 . Of particular interest is the case of homotypic mate preferences ("matching mating rule" 11 ), where individuals preferentially choose mates with which they share phenotypic traits such as colors 12,13 or acoustic signals 14 . When these traits are simultaneously under divergent or disruptive ecological selection (socalled "magic traits" 15,16 ), choosy individuals with homotypic mate preferences are less likely to produce offspring with unfit intermediate phenotypes. Theoretically, this effect creates indirect selection favoring a further increase in choosiness (i.e., stronger homotypic preferences), establishing, or strengthening premating isolation between diverging populations 11,[17][18][19][20][21][22] . In addition, choosiness often induces positive frequency-dependent sexual selection that favors the most common phenotype; if mating success varies among individuals and assortment is not perfect, individuals having the most common phenotype have the highest mating success. Once diverging populations are sufficiently differentiated with respect to a magic trait, disruptive ecological selection on that trait is therefore complemented by disruptive sexual selection due to choosiness, which in turn can drive the evolution of even stronger choosiness [23][24][25] (but note that this sexual selection pressure can also inhibit the initial evolution of strong choosiness in sympatry [24][25][26] or parapatry 16,27,28 ). Empirically, homotypic preferences based on "magic traits" subject to disruptive selection are, indeed, often associated with premating reproductive isolation during speciation in the presence of gene flow 16,29,30 . Here, we are interested in whether this isolating barrier is stable when mate choice can be expressed by both sexes.
Previous theoretical developments have focused on the evolution of choosiness in each sex separately [17][18][19][20][21][22]31,32 . In principle, the indirect selection effect explained above (i.e., reduced production of unfit intermediate offspring) favors increased choosiness in both females and males. Nevertheless, evolutionary pressures acting on the two sexes are profoundly different. Under idealized conditions of polygyny and unlimited male mating potential, males can be thought of as an unlimited resource for females and all females have therefore equal mating success 33 . Females, on the other hand, often represent a limited resource for males, with male-male competition for access to females generating differential male mating success. Consequently, having a preference directly affects how much competition a male faces for gaining a mate. Typically, males place themselves in a disadvantageous competitive setting if they preferentially court "popular" females-that can be phrased as sexual selection acting directly against male preferences 31,34,35 .
The extent of courtship effort (resources 36,37 or time 38 ) allocated toward preferred females is another key, but an underappreciated, factor in the evolution of male mate choice. Male preferences can only evolve if the lack of mating attempts with unpreferred females improves mating opportunities with preferred ones. Typically, this is conceptualized as reallocation of courtship effort, and male mate choice is hampered if complete reallocation is difficult to achieve. Male preferences can nevertheless evolve if direct or indirect benefits (e.g., increased probability to mate 34,39 , fertility 34,40 , or offspring quality 31,41 ) outweigh these costs. Previous theoretical studies have shown that male choosiness can evolve by reinforcement 31 , and that strategic male courtship allocation can generate polymorphic male preferences 32 , which may ultimately lead to reproductive isolation between populations. Overall, however, female choosiness is considered more likely to evolve than male choosiness, at least as long as it does not associate with competitive and opportunity costs 31 .
In some interspecific sexual interactions, both males and females discriminate against heterospecifics, and therefore engage in mutual mate choice with respect to species identity [42][43][44][45] . In cichlid fishes 46,47 and Heliconius butterflies 2,13,48-50 , which are textbook examples of potential speciation via premating isolation, both males and females can display homotypic preferences based on color. However, the consequences of mutual mate preferences for reproductive isolation remain to be explored. Preferences have been shown to evolve independently if female and male choices are based on distinct traits 34 . However, females and males with mutual homotypic preferences often use the same trait to evaluate potential mates. Choosiness in one sex therefore influences the evolution of choosiness in the other through genetic linkage disequilibrium 34,[51][52][53][54] . Female choosiness may also strongly favor the evolution of male choosiness by directly increasing the mating success of choosy males focusing their courtship effort on females that are likely to accept them as mates.
Here, by analyzing a population genetic model, we characterize the coevolutionary dynamics of female and male mate choice based on the same phenotypic trait under disruptive selection. We then assess its effects on the stability of reproductive isolation. We show that female choosiness favors the evolution of male choosiness, and that selection for mutual mate choice should be common. In turn, because female and male choosiness are based on the same phenotypic trait, male choosiness weakens indirect selection on female choosiness. In finite populations, this causes coevolutionary dynamics of "preference cycling", initiated by drift and completed by selection, with strong potential to destabilize reproductive isolation.

Results
Model overview. We model the evolution of assortative mating in sympatry, based on three diploid biallelic loci that segregate independently (no physical linkage). Disruptive viability selection acts on an ecological locus A, but without mate choice, ecological divergence is hampered by random mating that brings divergent ecotypes (AA and aa) together to hybridize. In addition, we implement two distinct choosiness loci F and M, which are independently expressed in females and in males, respectively. Both sexes can therefore use the trait under disruptive viability selection as a basis for mate choice (by using a matching rule). Female and male choosiness are ecologically neutral, but they can experience indirect selection via linkage disequilibrium with the ecological locus. Hybridization rates between ecotypes may decline due to assortative mating caused by female, male, or mutual preferences. Unless stated otherwise, we assume that the alleles coding for choosiness are recessive (only FF females and MM males are choosy).
Each generation first undergoes disruptive viability selection with strength s, such that heterozygotes at the ecological locus (Aa) suffer increased mortality. Males then court females and are "visible" to them (i.e., available as potential mates) proportionally to the courtship effort they invest. Choosy males (MM) prefer to court females that match their own ecological trait. In the case of a mismatch, they reduce their courtship effort to a very small fraction ϵ m << 1 of what nonchoosy males would invest. The courtship effort thus saved can be reallocated toward preferred females, where the extent of this reallocation is described by the parameter α. In particular, if choosy males reallocate all saved courtship effort toward preferred females (α ¼ 1), they enjoy a strong mating advantage over nonchoosy males with these particular females. Females likewise express different propensities to accept courting males. Choosy females (FF) prefer males that match their own ecological trait. We assume that in the case of a mismatch, they reduce the probability of mating to a very small value ϵ f << 1. Small ϵ f and ϵ m therefore reflect strong choosiness. Unlike males, all females have the same mating success as is the case in many polygynous mating systems. The expected genotype frequencies in the next generation depend on the probabilities of mating between different genotypes, with the new generation being formed by assuming independent Mendelian inheritance at all loci.
This three-locus diploid model of mutual mate choice is too complex to produce analytical solutions 34 . The behavior of the model can be assessed, however, by numerical analyses and computer simulations. We first analyze the deterministic behavior of the model, assuming an infinite population. Subsequently, we perform stochastic simulations in populations of finite, yet appreciable, size to account for genetic drift affecting traits under weak selection. In each generation, stochasticity is introduced by sampling K offspring individuals following the distribution of genotype frequencies predicted by the deterministic model (just as in the Wright-Fisher model of genetic drift). In addition, in these stochastic simulations, we allow for mutation of alleles in offspring.
Viability and sexual selection on female and male choosiness. Viability and sexual selection may act directly on ecological or choosiness loci (through differential viability and male mating success among genotypes, respectively; Supplementary Fig. 1), and despite free recombination, also generate linkage disequilibrium between loci by creating statistical associations between alleles (in diploid individuals and in haploid gametes; Supplementary Figs. 2 and 3). In particular, the majority of choosy females and choosy males are homozygous at the ecological locus. In addition, choosy males often carry alleles coding for female choosiness (which are neutral during courtship). This linkage disequilibrium between choosiness loci arises because both choosy females and choosy males use the same ecological trait as the basis of mate choice, and therefore, tend to mate with each other. Consequently, selection that changes frequencies at a given locus will also change frequencies at other loci. As detailed below, such indirect selection plays a key role in the evolution of male and female choosiness.
To understand the intricate interplay of selective forces, we first consider cases where choosiness can evolve in only one sex. Disruptive viability selection directly acts on the ecological locus (black arrow, Fig. 1) with homozygotes having a high viability (hereafter, "homozygous" and "heterozygous" refer to the genotype at the ecological locus). In addition, female choosiness can induce positive frequency-dependent sexual selection on the ecological locus (green arrow, Fig. 1), such that homozygous males have the highest mating success. Consequently, in the case where only female choosiness can evolve, female choosiness is favored by indirect viability and sexual selection due to linkage disequilibrium with the ecological locus ( Fig. 1a; Supplementary Fig. 4).
The situation is rather different if only male choosiness can evolve. Unlike female choosiness, male choosiness does not induce direct sexual selection on the ecological locus, because females in our idealized polygyny scenario do not differ in their mating success, even if some receive less courtship than others. However, male choosiness intensifies male-male competition for the preferred female type, which induces sexual selection on the male choosiness locus itself (pink arrow, Fig. 1). Because the majority of choosy males are homozygous at the ecological locus (linkage disequilibrium), males courting "popular" homozygous females face strong competition for mating opportunities ( Supplementary Fig. 1). Choosy homozygous males therefore place themselves in a disadvantageous competitive setting and have low mating success. In addition, if reallocation of courtship effort is only partial (α < 1), choosy males lose courtship opportunities, which further lowers their mating success (blue arrow, Fig. 1). Consequently, in the case where only male choosiness can evolve, male choosiness is favored only if negative sexual selection is offset by indirect viability selection due to linkage disequilibrium with the ecological locus ( Fig. 1b; Supplementary Fig. 4).
We now turn to our main case, where choosiness can evolve in both sexes. In addition to selection acting on female and male choosiness separately (Fig. 1a, b), choosiness in each sex now induces selective forces on choosiness in the opposite sex. First, because choosy females mainly reject nonchoosy (nonmatching) males, female choosiness directly increases the mating success of choosy males (red arrow, Fig. 1c). Second, sexual selection induced by female choosiness on the ecological locus also indirectly favors male choosiness (dashed green arrow, Fig. 1c). Finally, due to the linkage disequilibrium between choosiness loci, all selective forces acting on male choosiness also indirectly affect the evolution of female choosiness (dashed pink, red, and blue arrows in Fig. 1c). Indirect sexual selection on female choosiness might, in the simplest settings, reflect direct sexual selection on male choosiness, transmitted via linkage disequilibrium between the M and F loci. However, the causalities can be more complex than that, as evidenced by cases where male choosiness increases in frequency, while female choosiness decreases in frequency in the life cycle stage that involves mating and reproduction ( Fig. 2d-f, colors indicate changes during that life cycle stage). The reason is that net indirect selection results from linkage disequilibrium between all three loci. Choosy males that also carry alleles for female choosiness (genotypes FF-MM) are particularly likely to be homozygous at the ecological locus (AA or aa), whereas other types of choosy males (genotypes ff-MM and Ff-MM) are somewhat more likely to be Aa heterozygotes ( Supplementary Fig. 2). Under these particular conditions, choosy homozygous males (unlike their heterozygous competitors) suffer intensified male-male competition, which generates indirect selection against female choosiness via linkage disequilibrium. As a whole, sexual selection can therefore favor male choosiness and simultaneously inhibit female choosiness.
Finally, if males are choosy, females are at low risk of producing unfit hybrids, regardless of their own choosiness, as males focus their efforts on matching females. This weakens the linkage disequilibrium between the female choosiness locus and the ecological locus ( Supplementary Fig. 2), and reduces indirect selection favoring female choosiness (cf. selection gradients in Fig. 4), with the consequence that female choosiness may easily drift in finite populations. The implication of such weak selection for choosiness coevolution will be assessed below by using stochastic simulations.  (Fig. 2c). Recall that while choosy individuals avoid courting/mating across ecotype boundaries, premating isolation is not complete (ϵ m ≠ 0 and ϵ f ≠ 0). Consequently, with mutual mate choice, female and male choosiness synergistically reduce hybridization rate between ecotypes ( Supplementary Fig. 5). Finally, changing the dominance hierarchy (i.e., making alleles F or M dominant) does not change our results qualitatively ( Supplementary Fig. 6).
The consequences of drift on coevolutionary dynamics. We next ran stochastic simulations to investigate the coevolutionary dynamics of male and female choosiness in populations of finite, yet appreciable, size (K ¼ 500). Unless stated otherwise, we here consider scenarios with strong disruptive selection (s ¼ 0:2), for which the deterministic outcome is female mate choice (for α 0:01) or mutual mate choice (for α > 0:01). We define a frequency threshold (¼ 0:85) above which female or male populations are considered to be mainly choosy (i.e., mainly express a choosy behavior before mating). We thereby characterize four regimes of choosiness: female choice only (F ), male choice only (M), mutual choice (F M), and partial choice (i.e., both female and male populations are at most partly choosy, P) (Fig. 3a). Note that regime P includes the regime of complete random  . Diamond arrows and classic arrows represent viability and sexual selection, respectively. Viability selection acts through differential survival (i.e., change in frequencies during the disruptive viability selection process), whereas sexual selection acts through differential male mating success (detailed in Supplementary Fig. 1). Bold and dashed arrows represent direct and indirect selection, respectively. AA and aa refer to homozygotes at the ecological locus and MM to choosy males. Selective forces represented with a low opacity in panel c are the ones already shown in panels a and b mating and is therefore different from the deterministic equilibrium of partial male choice described above (Fig. 2a).
For α > 0:01, our deterministic analysis predicts a stable equilibrium of mutual mate choice, yet with drift, choosiness traits can temporarily evolve away from this equilibrium (Fig. 3b), entering the regimes of male choice only (regime M) or partial choice (regime P). This is caused by drift-induced and selectiondriven coevolutionary dynamics of female and male choosiness that we describe below.
Despite selection favoring mutual mate choice, assortative mating is often based solely on male choosiness in stochastic simulations (regime M, Fig. 3b). When females are choosy, male choosiness is strongly favored, with drift playing an insignificant role. However, when males are choosy, selection favoring female choosiness is weak (as explained above); the frequency of choosy females may then decrease through drift (Fig. 4). Nonchoosy females can persist for significant periods of time, during which assortative mating is maintained by male choosiness only (regime M). the population at a given point in time has a significant chance of observing the regime P.
Hereafter, we refer to this coevolutionary dynamics as "preference cycling", since female and male choosiness go through deterministic cycles triggered by stochasticity, involving departure from regime F M into regimes M, P, and sometimes F (before returning to regime F M). Preference cycling also occurs if we assume that females do not all have the same mating success (i.e., if we add a weak cost of female choosiness; Supplementary Note 1) or if choosiness is allowed to vary as a continuous trait (Supplementary Note 2).
Preference cycling strongly increases hybridization rate. Since we assume that neither sex can ever achieve perfect choosiness (ϵ m ≠ 0 and ϵ f ≠ 0), it is not surprising that hybridization rate is the lowest in the regime of mutual mate choice (F M) and becomes somewhat higher during drift-induced excursions into the male choice regime (M) (blue area in Fig. 5d). More importantly, hybridization strongly increases during deterministic excursions into the partial choice regime (P) (gray area in Fig. 5d). The greatest increase in hybridization occurs if the population stays in this regime P for extended periods of time. As shown in Fig. 5, while the regime P is reached most frequently for α ' 0:1, the average time spent in this regime is maximal for α ' 0:8, and during these episodes, maladapted heterozygotes reach frequencies of up to 35%. Overall, preference cycling leads to temporary peaks of hybridization, which periodically homogenize populations (as shown by fluctuations of F ST measured at neutral loci between genotypes AA and aa in Supplementary  Fig. 7). In other words, even though mutual mate choice, whenever it occurs, achieves the strongest degree of assortative mating, it is also particularly prone to periodic breakdowns, which as a whole hamper the maintenance of premating isolation. Importantly, these periodic breakdowns of reproductive isolation do not occur if choosiness is allowed to evolve in one sex only ( Supplementary Figs. 8 and 9).
If disruptive viability selection is weak (low s) or population size is small (low K), drift remains strong relative to indirect selection acting on female choosiness, and preference cycling occurs more frequently, increasing the overall hybridization rate (Supplementary Figs. 10 and 11). If alleles coding for choosiness are dominant instead of recessive, the final approach to fixation of the male choosiness allele is less rapid, and more nonchoosy males remain in the system. As a result, preference cycling induced by drift of female choosiness occurs rarely (Supplementary Fig. 12). Likewise, if male choosiness is less perfect, preference cycling occurs less frequently, and the overall hybridization rate is decreased (high ϵ m ¼ 0:03 instead of 0.01 in Supplementary Figs. 13 and 14). This somewhat counterintuitive result is explained by the fact that with imperfect male choosiness, selection on female choosiness never becomes so weak as to be overwhelmed by drift. More generally, preference cycling may occur if the evolution of choosiness leads to nearly perfect reproductive isolation between ecotypes. This is the case if choosiness per se is nearly perfect (Supplementary Figs. 13 and 14), or if reproductive isolation is strengthened by additional barriers to gene flow (but reproductive isolation cannot drop below a certain value in that case; Supplementary Fig. 15). Under these conditions, which intuitively seem conducive to speciation, we show that coevolution of male and female choosiness has the potential to strongly destabilize reproductive isolation.

Discussion
Surprising coevolutionary dynamics of male and female mate choice occur when the preferences of both sexes are based on the same phenotypic trait under disruptive selection. We showed that choosiness (i.e., the strength of preference) in one sex influences the evolution of choosiness in the other, a factor that has not been considered in previous models of speciation with gene flow [17][18][19][20][21][22]31,32 . Based on the predictions of our model, genetic drift, incomplete reallocation of courtship effort, and strength of choosiness all prove to be important when examining the outcome of this coevolutionary dynamics in terms of reproductive isolation.
In our model, male and female preferences are based on the same phenotypic trait, but are themselves governed by different loci. This genetic basis gives scope for preference coevolution. Selection generates linkage disequilibrium between choosiness and ecological loci (despite free recombination), and indirect viability and sexual selection resulting from this linkage disequilibrium has profound consequences for the evolution of Female choosiness therefore favors the evolution of male choosiness if at least some of the courtship effort saved by refraining from courting unpreferred females can be reallocated to gain a mating advantage with preferred females. As a consequence, male choosiness evolves more easily in our model than in models investigating the evolution of choosiness in each sex separately 31 . Indeed, our results show that mutual mate choice should often be favored by selection, and can induce very strong reproductive isolation in infinite populations.
Counterintuitively, however, in finite populations, this regime of mutual mate choice is particularly unstable. In the presence of even weak genetic drift, the coevolutionary dynamics of female and male mate preferences can lead to transient but periodic breakdowns of premating isolation, strongly increasing the hybridization rate. The fact that either sex can cause assortative mating makes it difficult for mutual mate choice to be maintained; more precisely, if male preferences are sufficiently strong to establish assortativeness, female choosiness has little effect on the mating outcome and is therefore free to drift. When female choosiness is reduced, male choosiness becomes disfavored by selection, temporarily leading to a regime of random mating. This coevolutionary dynamics of "preference cycling", initiated by drift and completed by selection, strongly destabilizes reproductive isolation and leads to periods of increased hybridization, which homogenizes populations.
The establishment of premating isolation is often considered to be the first step toward speciation with gene flow, and its stability is therefore a key prerequisite for other isolating barriers to evolve. For instance, only stable premating isolation should allow the establishment of neutral genetic incompatibilities leading to subsequent postzygotic isolation [55][56][57] . If selection favoring mutual mate choice induces preference cycling and dynamic instability of premating isolation, our results suggest that it can become an obstacle to speciation. This scenario contrasts with the traditional view of speciation as a gradual process characterized by a constant accumulation of barriers to gene flow (so-called "speciation continuum") 2,58-60 . Speciation can also be "undone"; like assortative mating in our model, barriers to gene flow can dissolve, and genetic discontinuities may vanish, thereby merging two taxa into a single population by hybridization 61,62 . Our model predicts such cycles of divergence, and gene flow may actually characterize the process of diversification in nature.
To track the "speciation continuum", empirical research often estimates isolating barriers between pairs of populations varying in their level of differentiation 2,58 . Yet, it is important to ARTICLE remember that such measures, being snapshots in time, do not yield information on the long-term stability of reproductive isolation. Our predictions regarding the coevolutionary dynamics of male and female mate preferences mean that premating isolation caused by mutual mate choice should be interpreted cautiously: whether gene flow is reduced over long periods of time is an open question. If a species' range is partially fragmented, preference cycling could also cause variation in the degree of reproductive isolation among local populations, as observed, for instance, in Catostomus fish species 63,64 . Yet, variation in hybridization is rarely quantified across several natural populations, and our study highlights the value of studies characterizing the strength of isolating barriers at a broader spatial and temporal scale. The coevolutionary dynamics of female and male mate preferences (and the resulting reproductive isolation) depends crucially on how much courtship effort males can reallocate toward preferred females, as a result of foregoing courting unpreferred females. By treating the reallocation of courtship effort as a parameter (α), our model covers a wide variety of courtship and mating systems in animals. In particular, "courtship effort" may refer to time (e.g., for mate searching or performing complex displays 38 ) or energy (e.g., for resource-demanding spermatophores 36 or nuptial gifts 37 ). In addition, foregoing certain courtship opportunities and searching for more preferred mates might entail mortality costs, which also affect the extent of reallocation of courtship effort. Our results suggest that obtaining estimates of courtship reallocation in nature is important to increase our understanding of divergence in the presence of gene flow. Based on our predictions, selection should favor mutual mate choice with even little reallocation; without reallocation, in contrast, male choice should be deleterious. Importantly, if reallocation is partial, preference cycling may occur, possibly limiting divergence.
More generally, while female choosiness is abundant in nature, male choosiness usually associates with particular mating systems (e.g., monogamy) and often evolves in addition to female choosiness [42][43][44][45] (in agreement with our predictions). Yet, we showed that the evolution of male choosiness can destabilize reproductive isolation. Therefore, the mating system characterizing each taxon should determine the stability of reproductive isolation and therefore the likelihood of speciation. Further empirical studies could usefully add this consideration when testing for links between the mating system and speciation; our model suggests that taxa with mating systems that are prone to the evolution of male mate choice may not necessarily associate with high speciation rates.
The extent and stability of reproductive isolation also depends on how accurately existing preferences can be expressed. Imperfect preference may occur for many reasons. For instance, in the butterflies Heliconius melpomene and Heliconius cydno, female and male preference loci are physically linked with different color pattern loci 50 . Therefore, choosy individuals may not completely stop courting/mating across ecotype boundaries because each sex relies on different aspects of the phenotype. Our model predicts that this error-proneness could strengthen selection favoring mutual mate choice, which could in turn inhibit preference cycling and stabilize reproductive isolation. Thus, perhaps counterintuitively, our model suggests that imperfect male preferences lead to strong reproductive isolation in the longterm by maintaining selection favoring female preferences, and preventing drift-induced preference cycling. Such counterintuitive results are not unheard of: in the context of local adaptation, imperfect female choice has been shown to be more strongly favored by selection than perfect choice because it Reallocation of courtship effort ( ) Fig. 5 Coevolutionary dynamics of choosiness ("preference cycling") and the resulting hybridization rate in stochastic simulations (s ¼ 0:2, K ¼ 500). To describe the coevolutionary dynamics of female and male choosiness, we record the mean probability of reaching regime P from regime F M or M at each time step (a), the mean number of consecutive time steps in regime P (b), and the mean frequency of maladapted heterozygotes (Aa at locus A) in regime P before viability selection (c) as a function of the reallocation of courtship effort (α). To assess the resulting hybridization rate, we record the mean frequency of maladapted heterozygotes over evolutionary time (d). In panel d, we also represent the mean contribution of each regime of choosiness to hybridization. In stochastic simulations, hybridization rate is increased by the coevolutionary dynamics of female and male choosiness despite selection favoring mutual mate choice (when α > 0:01) (d). In particular, the coevolution of female and male choosiness leads to periodic episodes of random mating, strongly increasing the hybridization rate (with up to 35% of hybridization in regime P) (c) maintains a higher diversity of male types in the population 65 . Likewise, many theoretical studies have found cases where sexual selection favors partial choosiness 16,[24][25][26] . Our study adds to this quest the possibility of preference cycling in situations where choosiness evolves as a quantitative trait (see Supplementary Note 2). In particular, if choosiness only evolves to a partial degree, preference cycling may not occur; interestingly, this may favor speciation.
Our predictions are not limited to the context of emerging reproductive isolation among diverging populations in sympatry. We can expect similar coevolutionary dynamics of female and male preferences in more advanced stages of reproductive isolation, e.g., after secondary contact. Indeed, disruptive viability selection may be caused by genetic incompatibilities among more distantly related taxa. In this context, preference cycling could likewise temporally increase hybridization rate, and conceivably, explain the formation of "hybrid swarms" and subsequent genetic introgression 5 or hybrid speciation 66 .
Overall, our theoretical model adds support to the idea that premating isolation may often be readily reversible 27,56,67 . Intriguingly, we show that premating isolation should be particularly unstable under conditions that favor mutual mate choice. We highlighted some factors that could inhibit preference cycling (strong selection against hybrids, high carrying capacity, imperfect choices, and extensive reallocation of courtship effort). The geographical context of speciation and a more detailed look into alternative genetic architectures (e.g., "two-allele mechanisms" 68 , physical linkage among choosiness loci) could conceivably change the modalities of preference cycling and should therefore be investigated in future theoretical studies. It will also be fruitful to consider preference cycling in a system with multiple potential barriers to gene flow, for if the time between episodes of hybridization during preference cycling is very long, other barriers might have time to evolve despite preference coevolution taking place. On the empirical side, the occurrence of preference cycling and its impact on reproductive isolation remains to be tested. More generally, our study should stimulate further research on the stability of barriers to gene flow, and on the link between mating systems and speciation.

Methods
Genotypes. Our population genetic model is based on three autosomal diploid loci. Alternative alleles at each locus are represented by small and capital letters. An ecological locus, A, is subject to disruptive selection and can be used as a basis for mate choice (so called "magic trait"). Additional loci F and M alter female and male choosiness (i.e., strengths of homotypic preference) before mating. We assume that choosiness alleles code for either no choosiness or strong choosiness, i.e., preferences vary from indiscriminate to almost fully assortative. We assume no physical linkage (i.e., loci are on different chromosomes or very far apart on the same chromosome) and alleles assort independently of one another in gametes following Mendel's second law. This assumption allows us to reduce the number of dynamic variables needed to describe our genetic system. There are three genotypes per locus (e.g., AA, Aa, and aa for the A locus), and, therefore, we track the frequencies of 3 3 ¼ 27 genotypes in the population.
Assuming discrete generations, we follow the evolution of genotype frequencies pðtÞ within an infinite population (deterministic simulations) or within a finite population (stochastic simulations). p ¼ fp i g is a vector consisting of 27 elements {p 1 ; p 2 ; :::; p 27 } referring to the frequencies of the 27 genotypes present in newborn offspring. The life cycle is as follow: census, viability selection, courtship/mating, zygote formation (and random sampling in stochastic simulations).
Disruptive viability selection on the ecological locus. Environmental/ecological pressures act on an adaptive ecological trait and favour sympatric divergence into two distinct ecotypes occupying niches of equal size. To ensure the maintenance of polymorphism, a parameter s 0 > 0 confers an advantage to the rarer of the two homozygotes AA and aa. We also assume that heterozygotes Aa suffer viability costs described by a parameter s. Following these assumptions, the genotype frequencies after viability selection are: The normalization factorN ensures that the genotype frequencies p S i sum up to 1. If s 0 ¼ 0, the ecological allele that is more frequent initially may outcompete the other allele ("gene swamping" 69 ), hampering divergence 23,70,71 . To prevent fixation of a single ecological genotype, we always include some negative frequency dependence by implementing s 0 > 0 (leading to p S AA % p S aa ).
Male choice and courtship. denotes the courtship effort of a male with genotype m towards females with genotype f (m and f 2 f1; 2; :::; 27g). Males with genotype mm or Mm at locus M are nonchoosy and court all females with the same intensity ( toward all females). Homozygous MM males are choosy (i.e., they express homotypic preferences), and their courtship depends on the match between the ecological trait (locus A) of the female and their own. In case of a mismatch (e.g., between a male with genotype AA and a female with genotype Aa or aa), choosy males reduce their courtship effort to a small fraction (small ϵ m thus reflects strong choosiness). In other words, choosy males reduce resources (e.g., time or energy) spent on courting unpreferred females. Saved courtship effort can be reallocated (totally, partially, or not at all) toward courtship of preferred (matching) females. The extent of this reallocation is described by parameter α. Overall, of all possible courtship events that could happen in the population, a fraction C m;f will occur between males of genotype m and females of genotype f : where p S m and p S f are the frequencies of males of genotype m and females of genotype f after viability selection. If α ¼ 1, choosy males reallocate all saved courtship effort toward preferred females and therefore enjoy a strong mating advantage over their competitors with these particular females. Contrary to previous models 31,34,39,72 , however, male preferences can induce lost courtship opportunities. If α < 1, only part of the saved courtship effort is reallocated ( P m P f C m;f < 1), and total courtship effort may differ between individual males ( P f C m;f ≠ p S m ). Equation (2) therefore differs from those previous models and is instead analogous to a model of female mating preferences with opportunity costs 25 .
Female choice and mating. We assume that males are "visible" to females (i.e., available as potential mates) proportionally to their courtship effort, defining a baseline mating rate which can then be adjusted downwards or upwards by female denotes the willingness of a female with genotype f to mate with males with genotype m. Females with genotype ff or Ff at the locus F mate indiscriminately ( with all males), leading to mating rates that are directly proportional to courtship efforts. Homozygous FF females are choosy (i.e., they express homotypic preferences). Their decision to mate depends on the match between the ecological trait of the male and their own. In case of a mismatch, choosy females reduce the probability of mating to a small fraction of the baseline (small ϵ f reflects strong choosiness). Thus, the overall proportion of matings M m;f that occur between males of genotype m and females of genotype f is given by ð3Þ This equation is analogous to previous population genetic models of mating with female preferences 16,31,33,34,39,72 . It ensures that all females, even the ones that are less preferred by males (or that prefer rare males), have the same mating success (no cost of choosiness and no sexual selection from male choice). Likewise, the mating success of females with and without preference is identical ( P m M m;f ¼ p S f ). These assumptions are realistic for many polygynous mating systems; relaxing them by implementing a weak cost of female choosiness ( P m M m;f < p S f for f 2 f FF g) does not change our conclusions (whereas female choosiness does not evolve if it associates with a strong cost; see Supplementary Note 1).
Zygote formation. Expected genotype frequencies pðt þ 1Þ of zygotes in the next generation are calculated by summing the appropriate mating frequencies M m;f , assuming Mendelian segregation and free recombination between all loci.
Random sampling for stochastic simulations. Based on the above deterministic model, we also perform stochastic simulations in finite populations to account for drift at loci under weak selection. To do so, for each generation, we first apply Eqs (1) to (3) to the vector pðtÞ, yielding the expected frequency distribution of genotypes in generation t þ 1. The new vector pðt þ 1Þ is then obtained by randomly sampling K offspring individuals from this distribution. In addition, we assume that mutation can occur in each offspring individual with a probability μ per diploid locus.
Parameters and initialization. Unless stated otherwise, we perform simulations with strong choosiness in genotypes FF and MM (ϵ f ¼ 0:01, ϵ m ¼ 0:01). We also set s 0 ¼ 0:5 to ensure that polymorphism at the ecological locus A is maintained. Simulations start with only alleles f and m present in the population. Once the population has reached ecological equilibrium, 1% of choosy males and females are introduced at Hardy-Weinberg equilibrium, such that choosiness alleles are in linkage equilibrium with each other and with alleles at the ecological locus.
Deterministic equilibrium is typically reached in <1000 generations. In stochastic simulations, we model populations of appreciable size K ¼ 500 with a probability of mutation μ ¼ 10 À3 per individual and per diploid locus. For each parameter combination, 40 stochastic simulations were run for 100,000 generations and statistics were calculated by averaging over time within run and over these runs.
Note that in previous models [23][24][25][26] , sexual selection created by female choosiness has been shown to impede speciation if mating is initially random (such that intermediate phenotypes have the highest mating success) and if choosiness evolves in small steps. This is not the case in our model, where alleles conferring strong choosiness are allowed to directly compete with alleles for random mating (Supplementary Fig. 4). This allows a strong linkage disequilibrium to develop between FF genotypes and AA or aa genotypes, such that choosy females are mostly homozygous at the ecological locus, increasing the mating success of homozygous males and indirectly favouring ecological divergence.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All relevant data are available from the authors.