Sexual reproduction with variable mating systems can resist asexuality in a rock–paper–scissors dynamics

While sex can be advantageous for a lineage in the long term, we still lack an explanation for its maintenance with the twofold cost per generation. Here we model an infinite diploid population where two autosomal loci determine, respectively, the reproductive mode, sexual versus asexual and the mating system, polygynous (costly sex) versus monogamous (assuming equal contribution of parents to offspring, i.e. non-costly sex). We show that alleles for costly sex can spread when non-costly sexual modes buffer the interaction between asexual and costly sexual strategies, even without twofold benefit of recombination with respect to asexuality. The three interacting strategies have intransitive fitness relationships leading to a rock–paper–scissors dynamics, so that alleles for costly sex cannot be eliminated by asexuals in most situations throughout the parameter space. Our results indicate that sexual lineages with variable mating systems can resist the invasion of asexuals and allow for long-term effects to accumulate, thus providing a solution to the persisting theoretical question of why sex was not displaced by asexuality along evolution.


Introduction
The evolutionary maintenance of sexual reproduction has been considered the main unsolved problem in evolutionary biology [1][2][3][4] and it remains challenging after more than 40  genetic associations ( [5,6] see also [7,8]), decoupling male and female fitness [9], expending time and energy to find suitable mates or incurring risks of predation or parasitism associated with mating [10]. However, the major issue in understanding the evolutionary maintenance of sex is how to balance the disadvantages of sexuality in terms of the cost of producing males who provide no resources to offspring, the so-called 'twofold cost of sex', which may reduce the reproductive potential of a lineage to one-half per generation ( [1][2][3][4][11][12][13][14][15]).
The twofold cost of sex entails disadvantages at two levels. One is the competition between species and lineages. Mean fecundity of an asexual lineage may be twofold per generation compared with a sexual one so that, all others being equal, the former may invade the niche of the latter ( [2,16] but see [17]). A number of theories have explored how the mean fitness of a sexual population could be higher than that of an asexual one in the long run to compensate the twofold cost of sex. Considerable research effort has been put on the potential benefits of recombination and segregation in hastening the accumulation of beneficial mutations and in the elimination of deleterious mutations [12,14,15,[18][19][20][21][22][23][24][25][26], as well as on the lower ability of asexual lineages to adapt to heterogeneous or rapidly changing environments [4,[27][28][29], or to maintain evolutionary arm races with parasites (reviewed in [30], see also [28]). As a consequence, asexual lineages may have disadvantages in the long term compared with sexual ones (e.g. [18,[31][32][33], but see: [34][35][36][37][38]).
The other, main type of sex disadvantages is at the individual level. Even if asexuals are short-lived, this does not explain why parthenogenetic individuals do not outcompete sexual breeders in the short term as group selection is normally regarded as much weaker than individual selection (e.g. [39,40]). Indeed, much of the literature on the evolution of sex attempts to weigh the short-term advantage of asexual reproduction with the long-term benefits of sex, as well as asking whether sex has significant benefits that already act over shorter time-scales (see e.g. [12,13,17,19,21,24,[41][42][43][44][45]). Models show that for the short-term benefits of a sexual female to be twofold they require rather extreme conditions of environmental variation, which are still too extreme even under strong selection from parasites [46][47][48], when parasites are considered together with the accumulation of deleterious mutations [49][50][51], or when females are assumed to mate with genetically diverse males [52]. Modifier alleles that produce changes from asexuality to small increases in the frequency of sex may be shown to spread within a small, finite population if genetic drift produces negative linkage disequilibrium, although even the range of population sizes that favours sex are rather small when the twofold cost of sex is included [14,53,54]. Therefore, it seems that in many cases the twofold cost may actually not being compensated in the short term. Hence, the question is still open on how alleles promoting sex are not displaced in the short term by those promoting asexuality [15,55].
The twofold cost of sex originated through competition over mates [4,56,57]. The full twofold cost only occurs when sexual females invest 50% of their resources [11] to produce males that do not contribute with resources to the production of offspring, normally because they invest their entire reproductive budget in competition with other males to secure mates. Therefore, it is polygyny, and the associated male competition for mates, which causes the twofold cost of sex.
Polygynous, costly sex must have been favoured by selection because of their advantages over earlier sexual strategists [56][57][58][59]. Alleles that promoted in males the allocation of reproductive effort into mating behaviour rather than into contributing to the production of offspring should have benefited from the higher mating success of their carriers compared with those of less polygynous males. Thus, one evident short-term benefit that might have contributed to compensate for the twofold cost of sex for a sexual female is the likely increment in the ability of her sons to obtain multiple mates. However, as the polygynous system spreads, mean population fitness may have decreased, as well as that of the individual females compared with similar monogamous or parthenogenetic ones, which would have made the population more vulnerable to invasion by asexuals. Thus, our understanding of the maintenance of costly sex should benefit by considering its interactions not only with asexuals but also with low-cost sexuals (when males contribute with variable amounts of resources to offspring production). To our knowledge, no model has so far explored: (i) whether the interactions between several reproductive strategies operating simultaneously in the population may have a role in explaining the maintenance of costly sex, and (ii) especially, whether costly sex may persist under such circumstances without the need of recombination benefits that outweigh the twofold cost.
Our goal in this article is to study the dynamics of alleles that promote costly sex when they coexist with alleles promoting low-cost sexual modes and asexuality, when recombination or genetic segregation cannot offset the full twofold cost of sex. We use game theory in a sexual selection genetic model that investigates the evolutionary dynamics of asexual and sexual individuals that use different reproductive strategies in an infinite population.

The model
We considered a hypothetical population composed of infinite diploid individuals, without overlapping generations. Individuals may reproduce asexually (parthenogenetically) or sexually. When sexual, they produced even proportions of male and female offspring. The reproductive strategy of individuals was determined by two autosomal two-allele loci with independent segregation. The first locus S, influenced whether the reproduction was asexual (allele a) or sexual (allele s). Genotype aa determined an asexual (or parthenogenetic) reproductive strategy S 1 ; genotype ss produced sexual reproduction; and heterozygote as individuals may choose between sexual and asexual reproduction with a probability α and 1 − α, respectively. A second locus M, influenced the mating system and its associated cost-of-sex, it being non-costly (allele n) or costly sex (allele c) and only expresses in the presence of allele s in locus S. Genotype nn provided a non-costly sexual strategy S 2 , typically due to males contributing as much as females to the resources to offspring (e.g. symmetrical parental care), although it may also encompass negligible expenditure in males or male function in hermaphrodites; and genotype cc determined a costly sexual strategy (polygyny) S 3 , where males are polygynous and do not provide parental care at all, so that sexual females invest on average half of their resources to produce males that do not contribute to resources to offspring. The behaviour of heterozygotic individuals nc depended on the dominance relationship between n and c alleles indicated by parameter k (k = 1 when allele c is dominant, k = 0 when c is recessive, and 0 < k < 1 when n and c have some degree of codominant relationship).
All females, irrespective of their genotype and the resources provided by males, allocated equal parental resources to offspring. Recombination was assumed to influence fitness, because of differences in any characteristics of sexually produced offspring with respect to asexually produced ones, although we did not deal with the nature of possible recombination benefits. We introduced parameter R as the relative benefit of recombination with respect to asexuality. R does not include the twofold cost, as the latter is only considered as a reduction in number of offspring.
We used parameter b m for the male-relative-to-female potential contribution to offspring production (0 ≤ b m ≤ 1 for sexuals and b m = 0 for asexuals). Female breeding strategy might not allow symmetrical male contribution (consider, for instance, egg production in many invertebrates or mammalian pregnancy and lactation). Thus, we introduced parameter b f as the proportion of potential male care that finally contributes to the production of offspring, let us call it the proportion accepted by the female (0 ≤ b f ≤ 1), so that the product b m · b f is the actual male contribution relative to female (0 ≤ b m · b f ≤ 1 for sexuals and b m · b f = 0 for asexuals). Table 1 summarizes the parameters used in the model.
We assumed fitness (W) of each individual displaying the sexual or asexual strategy being determined by five parameters as follows: r being the degree of relatedness from parents to offspring (r = 0.5 for all sexual strategies and r = 1 for asexuals); N, the number of descendants produced (we assumed N = 1 for asexuals and for every sexual individual allocating the whole parental budget to offspring); R, the relative benefits of recombination with respect to asexuality; and the product b m b f was the contribution of male-relative-tofemale-to-offspring production. In summary, equation (2.1)) yields for sexual (W s , considering the added contributions of male and female breeders) and asexual (W a ) strategists: Thus, male contribution to offspring production (in short, parental care), b m · b f , depended on male and female genotypes (i.e. number of c alleles), and on the dominance relationship between n and c alleles in the M locus, as follows. Males with genotype cc did not provide parental care at all (b m = 0) and females of the same genotype only accepted a proportion b (0 ≤ b ≤ 1) of the male care to offspring. Males and females of the genotype nn invested equally to offspring and thus b m = b f = 1. The behaviour of heterozygotic individuals nc depended on the dominance relationship between n and c alleles indicated by parameter k: b m = 1 − k for nc males, and b f = 1 − k(1 − b) for nc females.
We used parameter α (0 ≤ α ≤ 1) to simulate a mixed sexual/asexual strategy for heterozygotic genotypes, as, in the S locus. Parameter α represents the proportion of individuals of these genotypes that, irrespective of being generated by a sexual or an asexual strategy, reproduce using a sexual strategy,    the rest 1 − α being the proportion of asexual ones. Therefore, we assumed that 1 − α proportion of individuals of genotypes (as nn), (as nc) and (as cc) used asexual reproductive strategies (α = 1, 0 or 0.5, when allele a is recessive, dominant or codominant, respectively, with respect to allele s).
For sexual individuals, differences in male mating success between genotypes were addressed by two positive functions that depended on allelic frequencies at the M locus, η and π (η ≥ 0 and 0 ≤ π ≤ 1, see below). Functions 1 + η and 1 − π represented the mean number of mates for homozigotic cc (polygynous) and nn (monogamous) males, respectively. Therefore, the mean numbers of mates were represented as deviations from unity. The mean mating success of heterozigotic nc males was θ = 1 − π (1 − k) + ηk, and thus 1 + η ≤ θ ≤ 1 − π depending on the type of dominance relationship between n and c alleles determined by parameter k.
The ability in mating competition was assumed to be higher for polygynous cc males compared with monogamous nn males, and hence the presence of polygynous males reduced the mating success of monogamous ones below unity. Therefore, depending on phenotype frequencies of male types in the population, 1 ≤ 1 + η ≤ m and 0 ≤ 1 − π ≤ 1 for polygynous and monogamous males respectively, where m (m ≥ 1) represents the maximum potential number of mates for cc males (only achieved when polygynous males are scarce in a population composed mainly of monogamous males). Thus, mean mating success for polygynous males ranged from m (when polygynous males are rare) to 1 (when monogamous males are rare), and for monogamous males it ranged from 1 (when polygynous males are rare) to 0 (when monogamous males are rare).
We computed mating functions 1 + η, θ and 1 − π , assuming that all females in the population mated, as follows.
Let P i (when i varies from 1 to 9) be the frequencies of genotypes (G 1 = aa, nn), (G 2 = aa, nc), (G 3 = aa, cc), (G 4 = as, nn), (G 5 = as, nc), (G 6 = as, cc), (G 7 = ss, nn), (G 8 = ss, nc) and (G 9 = ss, cc), respectively. P 1 , P 2 and P 3 and 1 − α proportion of P 4 , P 5 and P 6 represent the frequencies of asexual individuals in the population. Thus, the proportion of sexual individuals in the population is  Table 2. Genetic contribution of asexual individuals to the next generation. (Only a proportion 1 − α of heterozigotic individuals in the R locus (as -) choose a parthenogenetic strategy S 1 . The genotype of the filial generation is an exact copy of the parental generation. The asexual contribution to filial genotypes are L Ai = P i when i = 1, 2 or 3 and L Ai = (1 − α)P i when i = 4, 5 or 6.) half of them being females. Thus, the accumulated proportion of mates for all males must equal the total proportion of females in the population, τ /2: In the limit, when the maximum polygynous potential is m → ∞, only polygynous males can reproduce and monogamous males have no access to females. Then and lim π m→∞ = 1. (2.8) A solution that satisfies equations (2.7) and (2.8) is .
(2.10) Therefore, the mean number of mates per male in relation to his genotype in the M locus and the dominance relationship between n and c alleles is (2.12) , (2.13) where 1 + η, θ and 1 − π being decreasing functions with the frequency of the allele c in the population. The composition of the filial generation in the t + 1 generation arises through the reproductive effects of all sexual and asexual individuals in the t generation. Asexual individuals are the homozygotic ones aa at the R locus and a proportion 1 − α of heterozygotic individuals as. The asexual contribution to filial genotypes, L Ai (i = 1-6), is obtained by replicating the asexual genotype with fitness equal to unity (table 2):              The proportion of descendants generated by a strategy S 1 is 1 − τ . However, this value must be corrected by changes in proportions of genotypes owing to the action of sexual selection (competition for mates) operating in sexual strategies (table 3; see below). The contribution of sexual individuals to the filial generation includes a chain of effects: (i) proportion of male-female pairs per genotype and mean number of mates by males of this genotype, (ii) differences in fitness of male-female pairs, and finally, (iii) filial segregation of these pairs. The sexual contribution to the filial genotype i, L Si (i = 1-9), is the sum of the effects of pairs that segregate into this genotype. However, these values need to be corrected after the sexual selection and asexual reproduction processes have operated (table 3).
For sexual strategies, 36 male × female crossing genotypes are possible (combinations of as nn = G 4 , as nc = G 5 , as cc = G 6 , ss nn = G 7 , ss nc = G 8 and ss cc = G 9 genotypes; table 3). We do not consider departures from random mating (i.e. panmixia) other than differences in the degree of polygyny between males in relation to their genotypes. Two offspring are produced by each sexual pair, with recombined (Mendelian) genetic material and even sex ratio. The proportion of each of the 36 combinations of male × female pairs can be generated by the terms of a squared polynomial: [α(P 4 + P 5 + P 6 ) + P 7 + P 8 + P 9 ] 2 = τ 2 . (2.16) The relative contribution for each male × female pair to the next generation can be determined by the proportion of this combination in the population (term in the squared polynomial, equation (2.16)), its fitness (equation (2.2)), the mean number of mates of the male in relation to his genotype (equations (2.11), (2.12) and (2.13)) and the Mendelian segregation of genotypes (table 3). For example, the proportion of pairs with as nc male (θ being its mean number of mates, equation (2.12)) and as cc female (i.e. G 5 × G 6 mating in table 3) is 1 2τ Thus, the fitness of their descendants (equation (2.2)) is Finally, after Mendelian segregation of genotypes, the relative contribution of G 5 × G 6 pairs for each of their six types of offspring genotypes (G 2 , G 3 , G 5 , G 6 , G 8 and G 9 , see table 3) is 1/4 or 1/8 of the expression The frequencies for the nine possible genotypes in the t + 1 generation can be obtained after a twostep process. First, we add all relative contributions coming from sexual and asexual strategies, L Ai (t) and L Si (t). Let us now consider L i (t) = L Ai (t) + L Si (t), i = 1, 2, . . . 9, the sum of all partial contributions (tables 2 and 3) to filial genotypes. However, owing to the action of the sexual selection pressures, the sum of all contributions to the nine filial genotypes does not equal 1. Second, we correct L i (t) proportions to obtain the new genotype frequencies: where 9 i=1 P i (t + 1) = 1. The new allelic frequencies are

Results
We derived analytical solutions when feasible and carried out computer simulations to study the fate of allele mutants or migrants when introduced in a population where other alleles predominated. Allele frequencies increased/decreased according to fitness functions, with parameter values corresponding to the behaviour of the interacting strategists. We first explored the interaction for pairs of strategies; that is, when one of the alleles was not present in the population (i.e. strategies S 1 versus S 2 ; S 1 versus S 3 ; and S 2 versus S 3 ). Thereafter, we investigated the interaction between all the strategies when the four alleles at both R and M loci were present. Simulations were run for the relevant range of R (1-2) and m (≥2) parameters. Variations in initial conditions for α, k and b as well as in initial allele frequencies were explored but they are not fully presented here for simplicity. Interested readers can explore themselves all desired variations in the provided R script (see Data accessibility section).

Asex versus monogamous, non-costly sex
Here we consider the competition between non-costly sex (strategy S 2 ) and asexual reproduction (strategy S 1 ). For this particular case, allele frequency of the costly sex in the M locus is c = 0 and fitness return (equation (2.2)) becomes W s = 2R. The frequency of the allele a between two successive generations, t and t + 1, has an analytical solution:

Asex versus polygynous, costly sex
Let us now consider the interaction between costly sex (S 3 ) and asexual reproduction (S 1 ). For this case, allele frequency of the costly sex in the M locus is c = 1 and fitness return (equation (2.2)) becomes W s = R. The problem is similar to the previous S 1 versus S 2 and also has an analytical solution: where τ (t) = αP 6 + P 9 (see equation (2.4)) is the proportion of sexual individuals in the parental generation t; P 6 and P 9 being the proportions of genotypes as cc and ss cc, respectively, in the same generation t. Equation (3.2) has the following possible solutions: (i) pure-asexual strategy S 1 (i.e. allele frequency a = 1); (ii) pure costly sex strategy S 3 (i.e. allele frequency s = 1); and (iii) equilibrium between costly sex and asexual strategy (S 1 /S 3 ), irrespective of all possible values of allele frequencies a and s, when R = 2. Parameters m and b have no effect here as all sexual individuals are equally competitive, hence mean number of mates (equation (2.11)) becomes 1 + η = 1.
ESSs are asexuality (S 1 ) when R < 2, polygynous sex (S 3 ) when R > 2, and the equilibrium between asexuality and polygynous sex (S 1 /S 3 ), for all possible values of alleles a and s, when R = 2. This result is in agreement with the common view that costly sex needs twofold benefits of recombination to spread in a population of asexuals. As for the previous S 1 versus S 2 problem, the dynamics of the invasion of a pure-asexual population S 1 by polygynous c mutants is dependent on the value of parameter R; the higher the value of R above 2 the more speedy is the change from S 1 to the ESS S 3 .

Monogamous versus polygynous sex
The interaction between costly (S 3 ) and non-costly (S 2 ) sex is represented by c and n alleles, when asexuals are not present in the population (i.e. s = 1 in S locus and τ = 1). The mean number of mates of homozygote polygynous cc males and those of heterozygote polygynous cn males compared to monogamous homozygote nn males are m and 1 + k(m − 1), respectively (from equation (2.11-2.13).
The evolutionary game between these two strategies is analysed elsewhere (Polo & Carranza [60]). Analytical solutions are possible at both extremes of the dominance relationship between alleles n and c, and only numeric solutions for intermediate values of dominance between alleles a and c, (i.e. when 0 < k < 1). The game drive model shows that the mating potential for polygynous males, m, is the most important factor explaining the dynamics of the competition between n and c alleles, and that k, and the proportion of male parental care accepted by females, b, have secondary, although important effects (see below the interaction among the three strategies).
In summary, the model for the paired interaction between both sexual strategies shows: (i) when allele n is dominant or recessive (k = 0 or 1, respectively) non-costly sex (S 2 ) is the only ESS if m < 1 + b, costly sex (S 3

Asex versus monogamous versus polygynous sex
When the four alleles at both loci were present, our results confirmed that asexuality, S 1 , was the only ESS when the relative benefit of recombination was R < 1, and that for high values of the relative benefit of recombination (R > 2), parthenogenetic individuals were replaced by either costly or non-costly sexual ones, the interaction between both sexual strategies depending again on parameters m, k, b and initial allele frequencies.
Indeed, the most interesting results are those for intermediate values of the relative benefit of recombination, 1 < R < 2. For this situation we found that the dynamics of the three strategies (S 1 , S 2 , S 3 ) that resulted from the combined effects of alleles a, s, n and c, represented a rock-paper-scissors' type model based on the following, intransitive fitness inequalities between strategies: W(S 2 ) > W(S 1 ) (as R > 1; see §3.1 above); W(S 1 ) > W(S 3 ) (as R < 2; see §3.2 above) and W(S 3 ) > W(S 2 ) when m was above a minimum value, also depending on b and k parameters (see §3.3 above).
To exemplify these interactions, we may observe the dynamics of a population where the parthenogenetic S 1 strategy was dominant but alleles s, n and c were also present (figure 2). The frequency of the allele a in locus S tended to decrease as non-costly sex S 2 (alleles s and n in loci S and population. This result arises under these conditions because higher m makes alleles from polygynous males (c) spread very quickly at the expense of low-cost alleles n. Then, asexuality is not so effectively displaced because: (i) sexual reproduction predominates in its costly form, and (ii) the low R value provides very little advantage to sexuals. Other parameters affected more the duration of the transitory phases than the equilibrium balance between strategies. Thus, a decrease in b highly increased the transitory phase, and also had a secondorder effect of decreasing the final proportion of asexuals at the attractor point ( figure 5). On the other hand, deviations from perfect codominace between alleles (i.e. k > 0. 5   -non-costly sex, S 2 , was the only ESS when 1 < R < 2 and m ≤ 2; -the costly sex allele c was never completely eliminated for 1 < R < 2 except when nallele fixates at m ≤ 2; -as R approaches 2 (approx. above 1.4, table 4) and for m values only moderately higher than 2 (approx. below 2.2, table 4) the rock-paper-scissors dynamics showed transitory cycles between regions of dominant asexuality S 1 (a approx. 1) followed by a short or non-existent region of dominant non-costly sex S 2 (s and n approx. 1) and finally by a region of dominant costly sex S 3 (s and c approx. 1) and then a new S 1 region. In many cases, the dominant region S 2 was not present and the dynamics became mainly a transition between asexuality and costly sex strategies with the amplitude of transitory cycles increasing so that the frequencies of either a or c  Table 4. Parameter space data (i.e. allele frequencies for alleles s and c) at generation 10 000 as a function of R and m values, when α = b = k = 0.5 and initial frequencies were s = c = 0.01. (Table shows allele frequencies at t = 10 000. Cells represent cases that reached stable equilibria or instances in oscillations (*), some of which (**) reached values very close to fixation or extinction (1 or 0) with only traces of the competing alleles that can resume the cycles thanks to the infinite size of the modelled population. Note that transitions between these areas within the parameter space are not discrete, so that for some cells near the borders of the equilibrium area allele frequencies can still experience slight variations after generation 10 000, also depending on variations in initial conditions, although in all explored cases within the area without asterisks they finally reached fixation. Likewise, cases with oscillations can reach steady states after enough number of generations, the sooner the closer to the equilibrium area. See also figure 6 and     alleles may be very close to 1. This means that they would probably fixate in any finite population leading to either pure-asexual or costly sexual populations for R > 1 and m > 2, the odds for sexual ones being higher the closer they are to the R = 2 and m = 2 corner in the parameter space; and -most of the parameter space for 1 < R < 2 showed three coexisting strategies S 1 , S 2 and S 3 (either in stable equilibrium S 1 /S 2 /S 3 or in cycles), the main exception being the fixation of non-costly sex allele n against costly sex c when polygyny potential was low (m ≤ 2), rather than the displacement of sexuals by asexuals ( figure 6 and table 4; see also figure 2 for the dynamics of some particular cases within this parameter space).

Discussion
Our model can provide an answer to the old, persistent question of why sex is not displaced by asexuality in the short term: the diversity of sexual breeding systems that can occur during the evolution of a lineage, with variable degrees in the twofold cost-of-sex and in the opportunity for sexual selection, prevent the invasion and fixation of asexuality. In particular, our results show that with very small benefits of recombination (well below twofold), alleles promoting costly sex cannot be displaced by those promoting asexuality if sexual reproduction is represented by two modes that correspond to breeding systems commonly occurring in nature, such as monogamy and polygyny, which can frequently occur as flexible options even within a population (e.g. [62][63][64][65][66][67][68][69]; see below). Our model does not deal with the nature of the benefits of recombination. Rather, it simply explores a scenario in which recombination may be assumed to provide some benefit in the short term, although not necessarily twofold relative to asexuality. This is important in formulating the problem of sex [1][2][3][4]: some benefit of recombination may be necessary at the origin of sexual reproduction [70], but a twofold advantage (either from single or multiple causes: [12]) may not be required for costly sex to be maintained.
The three strategies included in our model have intransitive fitness relationships that follow the conditions of the rock-paper-scissors game [71,72], so that all the alleles involved can be maintained either in equilibrium or cycling depending on parameter conditions. This situation can solve the question of why sex is not displaced by asexuality before long-term effects have time enough to accumulate.
Besides the presence of the three strategies, the remaining conditions in our approach are rather unfavourable for the maintenance of costly sex. For example, populations in nature are not infinite, and finite populations have been proved to favour the persistence of sex [14,15,53,54]. In our model there are no density-dependent effects, and density may favour sex against clones [17]. We assumed global dispersion of individuals throughout the population while local dispersion and geographical structure can make sexuals to resist asexual invasions [73].
Our results for pairs of strategies simply confirm previous findings of models that investigated sex versus non-sex [4,14,21,24,46]. Non-costly sex and asexuality can maintain in equilibrium when there is no net benefit of recombination (R = 1), irrespective of all possible values of allele frequencies a and n; but non-costly sex is the only ESS when 2 > R > 1. Likewise, a pure costly sex strategy (S 3 ) needs twofold benefits of recombination to spread when only asexual individuals are present in the population. However, when alleles promoting the three basic strategies are allowed to interact, our model does not support the classical assumption that asexuality can replace sexual reproduction when recombination does not compensate for the twofold cost of sex [1][2][3]11]. Rather, the nature of the problem of the maintenance of sex appears to be cyclical, simply by considering the intransitive interactions between reproductive modes that may vary with respect to the twofold cost instead of the contest between only two extremes as in the original proposal.
With some value of R not so high as to compensate the twofold cost, variations in the polygyny potential m can largely influence the dynamics of the interaction between alleles promoting different reproductive strategies. Polygyny potential provides advantages to c alleles with respect to n alleles, causing a cycling dynamic between the twofold cost, non-costly sex and the asexuality tending to a continuous rock-paper-scissors game or to equilibrium between alleles. Therefore, an important result in our model is that such prevalence of asexuality caused by the twofold cost of polygyny does not completely eliminate the sexual strategies whenever any non-costly sexual mode remains as a possible alternative to proliferate in the long term. Although reproductive modes such as fair hermaphroditism, monoecy, nuptial gifts or biparental care are acknowledged to reduce the twofold cost of sex ( [2,5], but see [74]), their role in the maintenance of sexual reproduction has been much overlooked. Our model shows that these strategies that do not include the full twofold cost of sex may be crucial for the evolutionary maintenance of polygynous sex because they can buffer the evolutionary interactions between asexual and costly sexuals.
Data for current polygynous systems reflect that successful males can inseminate many females without contributing to parental care (see e.g. [10,75]), but a gradualistic variation occurs in many systems with reducing paternal care while increasing mating effort. For example, a proportion of multiple inseminations by unidirectional sperm transfer are present in most hermaphrodites [70,77], and males of most monogamous species engage in extrapair copulations [78] at the expense of reducing their contribution to parental care [10,79]. For a number of organisms the closest alternative to costly sex is not asexuality but biparental care. Female preference for good genes versus paternal investment represents the flexible balance at the individual level between fecundity (by paternal contribution) and costly sex (with stronger selection on male mating ability) [10,58,59,80]. Displacements along the gradient between low-cost and more-costly sex may have repeatedly occurred during the evolution of any lineage, for example by adjusting the mating system. Temporal and spatial variations in ecological conditions in the distribution range of a species promote the variability of mating systems (e.g. [69,81]), and phylogenies reflect that changes in mating systems should have occurred during the evolution of lineages, affecting the gains and losses of sexual traits [82].
Changes between breeding/mating strategies may not be equally probable. Our model has treated all changes equally, although in nature the maintenance of a reproductive strategy against invasion by alternative strategies may not only depend on parameter conditions but also on the ease of shifting between strategies, which is also likely to be different according to the direction of change. Thus, the step from asexuality to sexual reproduction can still occur in nature [66], and most predominately asexual organisms keep some sexual functioning [34,67,68,83,84]. Displacements between low-cost and costly sex can be achieved in both directions simply by changing the mating system and the relative contribution of the sexes to parental care (e.g. [10,85]). By contrast, mutations to parthenogenesis may be common in simpler organisms but less likely as they become more complex ( [21,43,46,64,65,83,86,87], but see [88] and citations herein). Some organisms have blocks to asexuality such as genomic imprinting in mammals [86] or paternally derived organelles in gymnosperms [21], which may represent actual 'traps' that impede asexual reproduction [89]. Livnat et al. [87] have proposed that selection for mixability and modularity promoted by sex from its origins should have had enormous influence in the architecture of sexual genomes, which may also have implications in the likelihood and success of an eventual shifting to asexual reproduction.
On the other hand, blocks to sexuality of asexual organisms may also be possible, and despite the fact that they may exist in obligate asexuals, they are rare in nature [35,63,83] perhaps because of the short longevity of pure-asexual lineages [25,37,83]. Evolutionary equilibria and cycles that result in our model show how sexual and asexual lineages may have successive opportunities to evolve blocking traits as side products of adaptive complexity [21,87], as well as the accumulation of other populational benefits and effects [12,13,[19][20][21][22][23]25,30,47,[90][91][92], which may help to explain how an evolving sexual lineage may persist up to accumulating complex differences without being replaced by asexual clones. Data accessibility. Interested readers can explore the model and introduce all desired variations by using the R-script [61] in electronic supplementary material and at: https://github.com/JRubalcaba/Carranza-Polo.
Authors' contributions. J.C. elaborated the idea and the theoretical basis for the model. V.P. made the mathematical model and the computer simulations. Both authors wrote the paper and approved submission.