Peter Czuppon

George W A Constable

Invasion and extinction dynamics of mating types under facultative sexual reproduction

Keywords: Mating types, Isogamy, Facultative sex, Self-incompatibility, Balancing selection, Frequency-dependent selection, Finite populations

In sexually reproducing isogamous species, syngamy between gametes is generally not indiscriminate, but rather restricted to occurring between complementary self-incompatible mating types. A longstanding question regards the evolutionary pressures that control the number of mating types observed in natural populations, which ranges from two to many thousands. Here, we describe a population genetic null model of this reproductive system and derive expressions for the stationary probability distribution of the number of mating types, the establishment probability of a newly arising mating type and the mean time to extinction of a resident type. Our results yield that the average rate of sexual reproduction in a population correlates positively with the expected number of mating types observed. We further show that the low number of mating types predicted in the rare-sex regime is primarily driven by low invasion probabilities of new mating type alleles, with established resident alleles being very stable over long evolutionary periods. Moreover, our model naturally exhibits varying selection strength dependent on the number of resident mating types. This results in higher extinction and lower invasion rates for an increasing number of residents.

Introduction

In isogamous species, the gamete size differentiation that defines the sexes in anisogamous species is absent [START_REF] Lehtonen | What do isogamous organisms teach us about sex and the two sexes? Phil[END_REF]. Despite their morphological similarity, the gametes 1 of isogamous species are typically not interchangeable, but rather fall into one of a number of genetically determined self-incompatible gamete classes, termed mating types. Syngamy can only occur between gametes of distinct and complementary mating types.

The evolutionary explanation for this self-incompatibility, which limits the number of potential mates available to an isogamous organism, has been the subject of debate, with a number of competing hypotheses proposed [START_REF] Billiard | Having sex, yes, but with whom? Inferences from fungi on the evolution of anisogamy and mating types[END_REF]. These include theories that such self-incompatibility alleles limit inbreeding depression [START_REF] Charlesworth | The evolution and breakdown of S-allele systems[END_REF], increase encounter rates between gamete pairs [START_REF] Hoekstra | On the asymmetry of sex: Evolution of mating types in isogamous populations[END_REF][START_REF] Hadjivasiliou | Cell-cell signalling in sexual chemotaxis: a basis for gametic differentiation, mating types and sexes[END_REF][START_REF] Hadjivasiliou | Evolution of the mating type locus with suppressed recombination[END_REF], allow for ploidy level detection and the instigation of the zygote developmental program [START_REF] Haag | Why two sexes? Sex determination in multicellular organisms and protistan mating types[END_REF][START_REF] Perrin | What uses are mating types? The "developmental switch" model[END_REF], or manage cytoplasmic conflict by promoting uniparental inheritance of organelles (UPI) [START_REF] Hurst | Cytoplasmic fusion and the nature of sexes[END_REF][START_REF] Hadjivasiliou | Selection for mitonuclear co-adaptation could favour the evolution of two sexes[END_REF].

Empirical observations of the number of mating types vary between species, ranging from 2 to many thousands [START_REF] Kothe | Tetrapolar fungal mating types: sexes by the thousands[END_REF]. For example, intermediate numbers of mating types are reported in ciliates of the Tetrahymena-species (3 -9 different mating types [START_REF] Doerder | High frequency of sex and equal frequencies of mating types in natural populations of the ciliate Tetrahymena thermophila[END_REF][START_REF] Phadke | Rapid diversification of mating systems in ciliates[END_REF]) and slime molds (2 -13 mating types [START_REF] Bloomfield | Sex determination in the social amoeba Dictyostelium discoideum[END_REF][START_REF] Clark | Reproductive systems in the myxomycetes: a review[END_REF]). Larger numbers can be found in fungal populations; in Coprinellus disseminatus the global population is estimated to contain 123 different mating types [START_REF] James | Evolution of the bipolar mating system of the mushroom Coprinellus disseminatus from its tetrapolar ancestors involves loss of mating-type-specific pheromone receptor function[END_REF], while Schizophyllum commune has a staggering 23, 328 distinct types [START_REF] Kothe | Mating types and pheromone recognition in the homobasidiomycete Schizophyllum commune[END_REF]. This naturally leads to the question which type of evolutionary pressures govern this diversity in mating type number? This has been the subject of much debate, motivated in part by a discrepancy between these empirical observations and simple evolutionary reasoning.

From a theoretical standpoint, one might naïvely expect to see a very large number of mating types within any given species, due to the "rare sex advantage" of novel types (Iwasa and Sasaki, 1987). Since mating types are self-incompatible, rare types have more opportunities for mating and thus each type experiences negative-frequency dependent selection. Therefore, each novel mating type produced by mutation should establish in the population and the number of mating types should consistently grow. However, this prediction stands in stark contrast to what is observed in the natural world. Although isogamous species with hundreds, or even thousands, of mating types are possible [START_REF] Kothe | Tetrapolar fungal mating types: sexes by the thousands[END_REF], examples of such species are very rare; the vast majority have very few (typically two) mating types [START_REF] Hadjivasiliou | Theoretical studies on the role and evolution of mating types and two sexes[END_REF]. Explaining this discrepancy between theory and empirical observation has been the focus of much work, and multiple theories have been proposed [START_REF] Billiard | Having sex, yes, but with whom? Inferences from fungi on the evolution of anisogamy and mating types[END_REF].

One prominent hypothesis is that UPI drives the evolution of two mating types [START_REF] Hadjivasiliou | Selection for mitonuclear co-adaptation could favour the evolution of two sexes[END_REF], with larger mating type numbers becoming less stable with the increased complexity of coordinating an organelle donor-receiver program [START_REF] Hurst | Cytoplasmic fusion and the nature of sexes[END_REF]. While recent modelling work has shown that when the frequency-dependent effects of UPI are accounted for, an invading donor-receiver program does not reduce the expected number of mating types [START_REF] Hadjivasiliou | Dynamics of mitochondrial inheritance in the evolution of binary mating types and two sexes[END_REF], perceived empirical support for the theory comes from the fact that many species with more than two mating types have developed mechanisms to ensure homoplasmy without such a program. For instance, in Paramecium bursaria, with up to eight mating types [START_REF] Phadke | Rapid diversification of mating systems in ciliates[END_REF], and S. commune, sexual reproduction is achieved following the exchange of nuclei between cells without cytoplasmic mixing [START_REF] Birky | Uniparental inheritance of mitochondrial and chloroplast genes: mechanisms and evolution[END_REF]. However this in turn leads to an opposing question; if it really is only a donor-receiver program that limits the number of mating types, why do ciliates and Agaricomycetes, with these alternative methods for ensuring homoplasmy, still feature species with very few types (as low as two and three respectively, see [START_REF] James | Why mushrooms have evolved to be so promiscuous: Insights from evolutionary and ecological patterns[END_REF])?

A further hypothesis [START_REF] Hoekstra | On the asymmetry of sex: Evolution of mating types in isogamous populations[END_REF][START_REF] Hoekstra | The origin of isogamous sexual differentiation[END_REF], with renewed attention (Hadjivasiliou and Pomiankowski, 2016), suggests that it is cell-cell signalling between gametes that limits the number of mating types. Here co-evolution between two resident types has been shown to potentially limit the evolutionary success of a third mutant mating type. If derived from a resident type, this new mutant must essentially cross a fitness valley before it can develop encounter rates with the residents comparable to their current pairwise encounter rate, limiting its invasion potential (Hadjivasiliou and Pomiankowski, 2016).

In each of these hypotheses, a biologically plausible emergent mechanism is sought that generates a selective advantage to a pair of mating types, thus limiting their number to two. A notable exception can be found in Iwasa and Sasaki (1987). Here it was demonstrated that under certain dynamics for the mating type encounter rate, the advantage to rare mating types could be suppressed. In the limit of infinitely long-lived gametes (that can always survive until a suitable partner becomes available) selection for mating type numbers greater than two could be eliminated entirely. It was verbally suggested that under such a scenario, genetic drift would purge new mating types, limiting their number. This led to a bimodal prediction for the number of types; populations would either have two (given a particular set of encounter rate dynamics and immortal gametes), or have infinitely many otherwise (with each at infinitely low frequency).

While the conditions required for two mating types in Iwasa and Sasaki (1987) were stringent (and indeed, the possibility of intermediate numbers of mating types impossible) it was notable for suggesting that genetic drift may have a key role to play in determining mating type number. In a similar vein, recent work has emphasized the relevance of finite population size null models (i.e. models in which all mating types are phenotypically similar) for addressing the distribution of mating type numbers observed in nature (Constable and Kokko, 2018;[START_REF] Czuppon | Evolution of mating types in finite populations: the precarious advantage of being rare[END_REF]. These studies stress that even in the absence of species-specific biological processes, the number of mating types in any real finite population cannot be infinite. Instead the expected number of types will arise from balance between mutations (which introduce new mating types) and extinctions (which decrease the number of types), leading to a number of mating types well below the studied population size.

As an explanation for the low number of mating types often observed in isogamous species, this mutation-extinction balance hypothesis may seem at first improbable, particularly in the light of a number of classic population genetic studies of self-incompatibility (SI) alleles in plants [START_REF] Wright | The distribution of self-sterility alleles in populations[END_REF][START_REF] Wright | On the number of self-incompatibility alleles maintained in equilibrium by a given mutation rate in a population of a given size: a re-examination[END_REF][START_REF] Wright | The distribution of self-incompatibility alleles in populations[END_REF][START_REF] Ewens | On the problem of self-sterility alleles[END_REF][START_REF] Nagylaki | The deterministic behavior of self-incompatibility alleles[END_REF][START_REF] Yokoyama | Population dynamics of sex-determining alleles in honey bees and self-incompatibility alleles in plants[END_REF][START_REF] Yokoyama | The expected number of self-incompatibility alleles in finite plant populations[END_REF], see also [START_REF] Clark | Self-incompatibility: theoretical concepts and evolution[END_REF] for a review. From a modelling perspective, the dynamics of SI alleles in gametophytic species, those where SI is determined at the haploid pollen stage [START_REF] Bod'ová | Evolutionary pathways for the generation of new self-incompatibility haplotypes in a nonself-recognition system[END_REF], closely resemble those of SI mating type alleles in isogamous species. (A comparison with the dynamics of SI alleles in sporophytic species, in which SI of the haploid pollen is determined by the diploid parent, can be more complicated due to the complex dominance relationships that regulate SI in some of these species [START_REF] Thompson | Non-linear dominance relationships between S alleles[END_REF][START_REF] Prigoda | Phenotypic and genotypic expression of selfincompatibility haplotypes in Arabidopsis lyrata suggests unique origin of alleles in different dominance classes[END_REF][START_REF] Billiard | A general model to explore complex dominance patterns in plant sporophytic self-incompatibility systems[END_REF].) In a seminal paper [START_REF] Wright | The distribution of self-sterility alleles in populations[END_REF], Wright predicted that for the model plant Oenothera organensis, a population of around 500 individuals could sustain approximately 13 SI alleles (see also [START_REF] Jack | Self-incompatibility alleles in the population of Oenothera organensis[END_REF]). As the number of predicted SI alleles would rise with increasing effective population size (which intuitively reduces genetic drift and hence extinction rates), many more SI mating type alleles might naïvely be expected in isogamous species, for which effective population sizes can be of the order 10 6 [START_REF] Baranova | Extraordinary genetic diversity in a wood decay mushroom[END_REF].

While in the above context the plausibility of the mutation-extinction balance hypothesis might seem doubtful, an important biological feature of isogamous species with less prominence in plants is the potential to reproduce asexually. In many isogamous species, long periods of asexual reproduction are punctuated by rare bouts of facultative sex. For instance, in the single celled green algae, Chlamydomonas reinhardtii (with two mating types [START_REF] Goodenough | Sex determination in chlamydomonas[END_REF]), recent genomic estimates have placed the rate of sexual reproduction to be once in every 770 asexual generations [START_REF] Hasan | The genomic landscape of recombination rate variation in Chlamydomonas reinhardtii reveals effects of linked selection[END_REF]. Sex in yeast, which also typically have two mating types [START_REF] Butler | The evolution of MAT: the ascomycetes[END_REF], appears to be rarer still, with estimates of once in every thousand to three thousand asexual generations reported in some species [START_REF] Tsai | Population genomics of the wild yeast Saccharomyces paradoxus: Quantifying the life cycle[END_REF].

In Hadjivasiliou et al. (2016), a computational model was used to show that long periods of asexual reproduction would lead to substantial drift in mating types frequencies. Using a population genetic model, Constable and Kokko (2018) demonstrated further that as sex becomes increasingly rare, the relative strength of genetic drift over selection for novel mating types is amplified. This leads to a lower expected number of mating types in mutationextinction balance (a similar observation has been made of Solanaceae gametophytic SI alleles [START_REF] Vallejo-Marín | On the Evolutionary Modification of Self-Incompatibility: Implications of Partial Clonality for Allelic Diversity and Genealogical Structure[END_REF]). This was true even in the absence of any species-specific selective mechanisms. A survey of available empirical data further supported the view that the rate of facultative sex is a key predictor of the number of mating types in isogamous species.

Analytic bounds on the number of expected mating types under a mutation-extinction balance were calculated in Constable and Kokko (2018), using the stationary distribution of a simple Moran-type model of mating type dynamics. In [START_REF] Czuppon | Evolution of mating types in finite populations: the precarious advantage of being rare[END_REF] an approximation of the expected number of mating types was calculated under the assumption that sex is obligate. [START_REF] Muirhead | Modeling multiallelic selection using a Moran model[END_REF] used a similar model to calculate the stationary distribution of the frequency spectrum of mating type alleles, however estimates on the number of expected mating types were not explicitly calculated. In addition, facultative sex was not modelled (although the mathematical analysis employed can account for this factor).

While such estimates of the number of mating type alleles in the stationary distribution are informative, they obscure the dynamical processes that drive and maintain this equilibrium, the establishment and extinction of alleles. These quantities are important for two reasons.

Firstly, they are both independent of the arrival rate of new SI alleles. Null models of mating type dynamics often equate the arrival rate of new SI alleles with an effective mutation rate, and in doing so take a course grained approach to the genetic, molecular and physiological mechanisms underlying mating type determination. However this effective mutation rate is a confounding parameter, as it is difficult to estimate empirically. (A similar problem was encountered by Wright, who conceded that there could be argument about what a "reasonable" rate should be [START_REF] Wright | The distribution of self-sterility alleles in populations[END_REF]). In contrast, the effective mutation rate affects neither the establishment probability of a mutant mating type nor the expected extinction time of a resident type. These quantities thus provide mutation rate independent measures of the evolutionary dynamics of mating types, that can be used to test the feasibility of the mutationextinction hypothesis.

Secondly, whereas the stationary distribution of mating type number can only be observed over evolutionary time periods, the establishment probability of new mating types and the mean extinction time of mating types can be observed over significantly shorter time periods. These quantities are thus important for providing empirically testable insight into the evolution of mating type number.

Extinction rates have previously been studied in related systems featuring negative frequency dependent selection, which results in high levels of polymorphism. In [START_REF] Takahata | A simple genealogical structure of strongly balanced allelic lines and trans-species evolution of polymorphism[END_REF] a one-dimensional stochastic diffusion was used to identify the dynamics of the polymorphic major histocompatibility complex (MHC) system with a time re-scaled neutral coalescent. This allowed for estimating the extinction rates of MHC alleles. These results were verified numerically in [START_REF] Takahata | Allelic genealogy under overdominant and frequency-dependent selection and polymorphism of major histocompatibility complex loci[END_REF] and later in [START_REF] Slatkin | Overdominant alleles in a population of variable size[END_REF]. The same theoretical approach was adapted to gametophytic SI in plants obtaining results showing that the time to the most recent common ancestor of a sex allele might even exceed the speciation time [START_REF] Vekemans | Gene and allelic genealogies at a gametophytic self-incompatibility locus[END_REF]. In the specific context of mating type alleles, such results have thus far been absent. However recently in [START_REF] Czuppon | Evolution of mating types in finite populations: the precarious advantage of being rare[END_REF] the establishment probability of newly arising mating type was calculated for populations in which sex is obligate. We will rely on this quantity to obtain our estimate on the mean extinction time of resident mating types.

In this paper we complement this literature by exploring scenarios of self-incompatibility under facultative sex. We begin by calculating an analytic expression for the stationary distribution of the number of mating types, extending the results of Constable and Kokko (2018) (where only bounds of the mode of this distribution were calculated). We then provide an expression for the establishment probability of a novel mating type in facultatively sexual populations, generalizing the results of [START_REF] Czuppon | Evolution of mating types in finite populations: the precarious advantage of being rare[END_REF]. Finally, combining these expressions, we calculate the mean time to extinction of a resident mating type allele in a novel way. To be more precise, instead of relying solely on the one dimensional dynamics of a focal mating type (in mating type frequency space) we make use of the dynamics on the number of mating types (in mating type number space). We conclude by discussing the biological implications of these results.

Model

We consider a population in which the self-incompatible mating type of an individual is determined by one of an infinite set of potential alleles at a single locus. Such single locus determination systems are for example found in the social amoebae Dictyostelium discoideum (with three mating types [START_REF] Bloomfield | Sex determination in the social amoeba Dictyostelium discoideum[END_REF]) and Didymium squamulosum (with up to 12 mating types [START_REF] Clark | Reproductive systems in the myxomycetes: a review[END_REF]). The evolutionary dynamics are determined by a Moran process; generations are overlapping and the population is assumed to be at an ecological equilibrium at which the number of individuals, N , is constant over time. We denote the number of individuals carrying mating type i by n i , such that ∞ i =1 n i = N . Asexual reproduction happens with probability c ∈ [0, 1]. In this case a randomly chosen individual produces a clone which replaces one of the individuals in the population, also chosen uniformly at random.

During sexual reproduction, which occurs with probability 1c, two randomly chosen individuals mate if they express different mating types. In the case of a successful mating, an offspring individual is produced that replaces one of the other individuals in the population. The mating type of the offspring is chosen at random from among the two parental alleles.

Additionally, we consider the emergence of novel mating types through mutations which occur at rate m. We implement mutations according to the infinitely many alleles model [START_REF] Kimura | The number of alleles that can be maintained in a finite population[END_REF], where the mutated individual expresses a completely new mating type not previously present in the population. For simplicity we consider the case in which mutation events are decoupled from reproduction. Mutants possess the same characteristics as any other mating type in the population; that is they are self-incompatible and mate with non-self types at the same per-capita rate as the resident types.

Mathematically implementing the model described above, the probability per unit time for a type i to increase by one and a type j to decrease by one through a birth-death event is given by

T i j = c n i N + (1 -c) 2 n i N k =i n k N reproduction n j N death , for n i > 0. ( 1
)
The reproduction term is split into an asexual component (the first term) and a sexual component (the second term). We note, that the sum in the sexual reproduction term goes over all non-i mating types present in the population, generating a reproductive advantage to rare mating types. The probability per unit time that a novel mating type i is generated from an ancestral mating type j is given by

T i j ∝ m n j N , for n i = 0, (2)
which can only occur when type i is not already present in the population. As described in Eq. (S6) in the Supplementary Information, the full expression for the probability transition rate given in Eq. ( 2) features an additional normalization constant, included to ensure that the mutation rate is independent of the population composition (i.e. that ∞ i =1 T i j | n i =0 = m). The deterministic limit (N → ∞) of this model is described by a system of ordinary differential equations where the dynamics (neglecting mutations) of a the i t h mating type, with frequency

x i = lim N →∞ n i N , is given by d x i d t = (1 -c) 2 x i j =i x j (x j -x i ) .
(3)

This system has been studied in Iwasa and Sasaki (1987) (see their Mating kinetics I). The dynamical system possesses an internal stable fixed point where all mating types are present at equal frequencies. Analogous finite population size models lead to stochastic differential equations and have been analyzed in [START_REF] Czuppon | Evolution of mating types in finite populations: the precarious advantage of being rare[END_REF]. Note that in that paper mutation was implemented in a distinct manner, with the assumption that mutation occurs with reproductive events. While that implementation is biologically more reasonable, the choice leads to only minor quantitative differences in the dynamics.

Data availability:

The code and simulation outputs can be found on github: https://github.com/gwaconstable/InvExtDynMatTypes.

Results

In this section we mathematically analyze the null model just presented. We begin by characterizing the long-time equilibrium behavior of the model, which is described by its stationary distribution. This allows us to answer the key question of how many mating types, M , are predicted by the null model as a function of the population size N , mutation rate m and, importantly, relative rate of asexual reproduction, c. However, while this results in an important benchmark for the expected number of mating types in real populations, it provides little insight into the dynamics of mating type number. To this end we further calculate the probability that a novel mutant mating type establishes in the population and the expected extinction time of a resident mating type allele. These quantities provide a deeper insight into the ongoing dynamics of mating type number in real populations.

The stationary distribution of individuals carrying each mating type allele

We denote by P st n the stationary distribution of the number of individuals of each mating type. In Constable and Kokko (2018) it was shown that an analytic solution for P st n is accessible as the probability transition rates T i j in Eqs. (1-2) can be decomposed into the product of a birth function b(n i ) and death function d (n j ) that each depend only on the number of each mating type reproducing and dying respectively;

T i j = b(n i )d (n j ) , (4)
where

b(n i ) = c n i N + (1 -c) 2 n i N N -n i N , if n i ≥ 1 , b(n i ) ∝ m , if n i = 0 , (5) d (n j ) = n j N , for all n j .
Under this decomposition, the stationary distribution of mating type alleles takes the exact form

P st n = 1 N M -1 i =1 n ↓ i -1 k=0 b(k)d (Φ i -k) b (Φ i -(k + 1)) d (k + 1) (6)
where n ↓ is the vector n reordered with its largest entries first, N is a normalization constant that enforces n P st n = 1 and Φ i is defined as

Φ i = N - i -1 j =1 n ↓ j . (7)
The full derivation of Eq. ( 6) is given in Constable and Kokko (2018). In terms of birth-death processes, it is interesting to note that Eq. ( 6) is the product of standard stationary distributions obtained for a two-allele (single variable) population [START_REF] Karlin | A First Course In Stochastic Processes[END_REF], moderated by a series of effective population sizes. A similar observation has been made in a related model of multi-allelic selection [START_REF] Muirhead | Modeling multiallelic selection using a Moran model[END_REF].

The stationary distribution of the number of mating type alleles

The stationary distribution of populations with M mating types is defined as P st M . This gives the probability of finding a population with a given set of parameters (N , c and m) in a state with M mating types at long times. P st M is related to the distribution P st n through the following summation;

P st M = n∈S (M ) P st n , (8)
where S (M ) is the set of all vectors n that represent a population with M present mating types (that is, all vectors n with M non-zero elements).

While Eq. ( 8) can be expressed neatly, evaluating this quantity is problematic as it involves summations of Eq. ( 6) over sets of the infinite vector n. To make analytic progress we note that when the population size N is large and the per-generation mutation rate m g = mN is small, at intermediate times the population resides in a quasi-stationary distribution in the region of a fixed point of the dynamics in the deterministic limit. More precisely, when the population is comprised of M mating types, this fixed point is given by n i ≈ N /M for each of the i present mating types, and n i = 0 otherwise. Following an extinction or mutation event (which changes the number of mating types in the population) the population quickly relaxes to a new quasi-stationary distribution in the region of an alternate fixed point with M ± 1 mating types, i.e. n i ≈ N /(M ± 1). Within the limit of large N and small m g then, P st n can be approximated by a superposition of these quasi-stationary distributions.

Briefly, the quasi-stationary distribution of the population in the region of a fixed point can be calculated by conducting a diffusion approximation on the underlying Moran model and linearizing the resultant advection-diffusion equation about a deterministic fixed point (the van Kampen approximation). The full calculation is conducted in the Supplementary Information, where we show that these quasi-stationary distributions are Gaussian (see Eqns. (S18)-(S22)).

Each of the quasi-stationary distributions can now be renormalized, using Eq. ( 6) to 'pin' the height of each quasi-stationary distribution to the height of Eq. ( 6) in the region of the deterministic fixed point (see [START_REF] Peter | Intrinsic noise in systems with switching environments[END_REF][START_REF] Vítor | Stochastic dynamics through hierarchically embedded markov chains[END_REF] for similar approaches). The full calculation is detailed in Section S3 in the Supplementary Information. M . For any given rate of asexual reproduction, c, larger populations (larger N ) contain more mating types on average. Meanwhile increasing the rate of asexual to sexual reproduction (increasing c) decreases the expected number of mating types. Analytic results (solid lines) are obtained by evaluating Eq. ( 9). Simulation results (corresponding shaded regions) are obtained using the Gillespie algorithm [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF], averaged over 5 × 10 6 generations, sampled every 10 2 generations and following a relaxation period of 10 3 generations from initial conditions with one more mating type than the mode number predicted analytically by the model (see Eq. ( 9)). The mutation rate is m = 10 -8 in all plots. The analytic description can be seen to capture the simulated behavior of the model with very high accuracy.

Substituting the resulting approximation for P st n into Eq. (8) and taking the limit of large N we find

P st M ≈ (2π) M 2 -1 M 2m 1 + c M -1 M M -1 2 M ! θ M -1 θ -M -1 1/2 × N (M -1)/2 θ -M -1 θ θ -M -1 M θ N , ( 9
)
where M is a normalization constant such that ∞ M =1 P st M = 1 and

θ = 1 + c 1 -c . ( 10
)
Note that when sexual reproduction is obligate (c = 0), θ = 1, while when sex is facultative and rare θ becomes large.

Comparing our approximate expression for the distribution of the number of mating types, Eq. ( 9), against the results of stochastic simulation of the population, we find excellent agreement (see Figure 1). We note that in a facultatively sexual population, where the role of genetic drift is amplified by high clonal reproduction rates, c, extinctions can even lower the number of mating types to one (see Supplementary Information, Section S8). In this scenario no further sexual reproduction events can take place, and the population reproduces purely asexually until a new mating type is generated by mutation.

The mode number of mating types

Since the distribution P st M is unimodal, obtaining an estimate for its mode is straightforward numerically. We define by r M the ratio of probabilities of having M and M -1 mating types;

r M = P st M P st M -1 . ( 11
)
For any given set of parameters, this ratio is independent of the normalization constant M . Substituting Eq. ( 9) into Eq. ( 11), we obtain a simplified analytic expression for r M (see Supplementary Information, Eq. ( S42)). Since the distribution P st M is unimodal, finding its mode is equivalent to finding the value of M for which P st M starts to decrease. The approximate mode of P st M , which we denote M o , can then be obtained as the solution to the equation

r M o = 1 . ( 12
)
The function r M -1 has a single root, and thus solving r M o = 1 for M o is numerically straightforward.

In Figure 2 we plot the mode number of mating types as a function of the per-generation mutation rate, m g and the population size, N , for four different rates of asexual to sexual reproduction. Here we see that facultative sex has a strong influence over the number of mating types, with even a ratio of 9 : 1 asexual to sexual divisions (panel (b)) reducing the expected number of mating types by an order of magnitude compared to the obligately sexual scenario (panel (a)).

The establishment probability of a new mating type allele

Having approximated the mode of the number of mating types we now proceed to study the evolutionary dynamics of the number of mating types. We begin by computing the probability of a successful establishment, Q Est M , of a newly arising mutant in a population of M resident mating types. We define Q Est M as the probability that this novel mating type allele (initially at a frequency 1/N ) reaches the new stationary frequency of mating types (1/(M + 1)) before any of the resident mating types goes extinct.

For c < 1 the deterministic equilibrium is an internal fixed point. This is crucial for our approximation of the establishment probability since we will identify it by the survival probability of a corresponding branching process (see Section S5 in the Supplementary Information). Furthermore, we assume that N M , i.e. the population size is sufficiently large that the resident mating types are not quickly lost through genetic drift. Under these assumptions we can extend the computation from [START_REF] Czuppon | Evolution of mating types in finite populations: the precarious advantage of being rare[END_REF] for obligate sexual reproduction and find In the case of obligate sex (c = 0), Q Est M reduces to 1/M , recovering the result found in [START_REF] Czuppon | Evolution of mating types in finite populations: the precarious advantage of being rare[END_REF]. For facultative sex (0 < c < 1) we see that the establishment probability is reduced. Since asexual reproduction increases the time to reach the stationary state in the deterministic system (i.e. the selection strength for even mating type ratios is reduced, see also Eq. ( 3)), newly arising mutants spend more time at low frequencies where they are susceptible to purging by genetic drift. This is reflected by Q Est M decreasing with c (see Figure 3). Our definition of the establishment probability includes that none of the resident mating types become extinct during the invasion process. This impedes a straightforward comparison with the case of obligate asexual reproduction, c = 1 (i.e. a neutral Moran model). The corresponding establishment probability might naïvely be assumed to be 1/(M + 1); however this value ignores the survival of all resident mating types and therefore overestimates the actual establishment probability. 13).

Q Est M = 1 -c 1 + c 1 M . ( 13
Simulation results are obtained from Gillespie simulations; in each plot N = 10 5 and results are averaged over 10 3 runs (c = 0), 10 4 runs (c = 0.9 and c = 0.99), and 10 5 runs (c = 0.999). Error bars are ±2 times the standard-deviation of the sample binomial distribution.

The mean time until the extinction of a mating type allele

Our final goal is to calculate the mean extinction time of a mating type allele. We assume that initially the population is close to its equilibrium (i.e. all mating types are approximately at equal frequencies). We then model the arrival of mating types (through mutation events) and their extinction (due to genetic drift) as a birth-death process on the number of mating types. Let β M be the probability per unit time that a population with M resident mating types gains a new mating type, and δ M be the probability per unit time that a mating type goes extinct.

Although we are interested in calculating the time between extinction events, we first turn our attention to the stationary probability distribution of this effective 'birth-death' process in the number of mating types; P st M . Solving the equations for the stationary distribution of the effective birth-death process, we obtain the well-known result (see for instance (Allen, 2011, Chapter 6))

P st M = β M -1 δ M P st M -1 . ( 14
)
Rearranging for δ M , we find

δ M = β M -1 P st M -1 P st M = β M -1 r M , ( 15
)
where we have used the definition of r M given in Eq. ( 11). The mean time between extinction events is the inverse of this effective death rate, i.e.

T Ext M = 1 N 1 δ M = r M N β M -1 , ( 16
)
where the factor N accounts for the fact that we are measuring the time until an extinction event in time units of generations. Further, we note that since T Ext M is determined using P st M (which itself accounts for covariances in allelic frequencies in the stationary state), our expression for T Ext M captures the effect of allele frequency fluctuations in all M types. As we have already analytically calculated r M (see Eq. ( 11)), all we now need to evaluate the mean time until the extinction of a mating type allele is an expression for β M (see Eq. ( 16)). We assume that the effective 'probability birth rate', β M (the probability per unit time that the number of mating types in the system increases) is given by the product of the mutation rate and the establishment probability of a mutant allele (see Eq. ( 13)):

β M ≈ mQ Est M = m 1 -c 1 + c 1 M . ( 17
)
This approximation relies on the frequencies of the resident types being close to their deterministic equilibrium (an assumption made in the derivation of Q Est M ). In general this will be true if the state containing M types is sufficiently stable and if mutations are occurring sufficiently infrequently that the population has time to relax to this quasi-stationary state. If these conditions are not met, it is possible that a new mutant mating type may arise in a population with a highly uneven mating type distribution, where Q Est M will no longer provide an accurate estimate of the establishment probability. We discuss these points further below.

Inserting this effective birth rate into Eq. ( 16) and substituting r M from Eq. ( 11) (see Eq. ( S42)), we can express the mean extinction time as

T Ext M = 2π 2(M -1) 1 + c θ 3 2 M M -1 M -3 2 θ -1 M θ -1 M -1 -1 2 × 1 N 3 2   θ θ θ -1 M θ -1 M -1 1-M θ θ - 1 M -1 -θ   N . (18)
We first note that as we would expect, Eq. ( 18) is independent of the mutation rate of novel mating type alleles (the mean time until the extinction of a mating type allele should not depend on the time until a new mating type allele arrives). Further, we can see that when the population size, N , is large the mean time to extinction is dominated by the term raised to the power of N in Eq. ( 18) (explaining the linear growth of log(T Ext M ) for large N in Figure 4). Although Eq. ( 18) is lengthy, it is useful as it allows us to quickly evaluate the mean extinction time for arbitrary parameter values. We can then rapidly explore parameter regimes that would be prohibitively time consuming to simulate (see Figure 5).

The expression derived in Eq. ( 18) becomes increasingly accurate with increasing population size, N , decreasing rates of asexual reproduction, c, and smaller numbers of initial mating types, M . The approximation can break down however when either N is small, or c or M are large. In this latter range, as the population dynamics become increasingly dominated by genetic drift, various assumptions involved in the derivation of Eq. ( 18) can become invalid. Most importantly the Gaussian approximation for the quasi-stationary distribution of a focal mating type breaks down. The true distribution becomes increasingly flat, while the Gaussian approximation does not respect the boundary conditions requiring non-negative number of 18) (right of dashed line) and the neutral limit of Eq. ( 19) (left of dashed line), defined by Eq. ( 20).

As the rate of asexual reproduction is increased (from c = 0.1 in panel (a) to c = 0.9 in panel (b)), the mean extinction time drops rapidly (note differing scales on x-axis).

mating types n i . Since this results in erroneous approximations of P st M , (the crucial ingredient calculating Eq. ( 14)) when the number of mating types is highly unstable (i.e. when M M o ), we do not expect our analytically derived mean extinction time to provide reasonable results in this regime. This reasoning is assessed in more detail and validated in the Supplementary Information (see Section S6).

We now seek a more quantitative measure of when we expect Eq. ( 18) to remain valid. We first note that naturally T Ext M should increase monotonically with the population size, and decrease monotonically with the rate of asexual reproduction and the number of resident mating types. However, on evaluating Eq. ( 18), we find that these expectations are violated as the relative strength of genetic drift increases (i.e. as the initial number of mating types, N /M , becomes small) or, conversely, as the strength of selection for equal numbers of mating types decreases (i.e. c becomes large). In these regimes, the dynamics of the mating type alleles approach neutrality, and the extinction time is better approximated by standard results on the 20)). The blue dashed line indicates the parameter regime in which extinction rates become approximately neutral (Eq. ( 19)). Extinction times in the gray shaded region exceed 10 11 generations. As these time are approximately longer than the evolutionary history of fungi, they have been omitted for clarity.

neutral multi-allelic Moran model Baxter et al. (2007) (see also Eq. (S51) in the Supplementary Information);

T Neutral M = -N M -1 s=1 (-1) s-1 M s s M log s M . ( 19
)
Combining T Ext M and T Neutral M , the mean time to extinction of a mating type allele can then be approximated by

T M = T Ext M while ∂T Ext M ∂N > 0, ∂T Ext M ∂c < 0, ∂T Ext M ∂M < 0, T Neutral M otherwise . ( 20
)
In Figure 4, we can see that this captures the results of simulations very well.

In Figure 5 we plot predicted extinction times for regions of parameter space that are prohibitively time consuming to investigate numerically. The gray regions which dominate the parameter space for low asexual reproduction rates depict extinction rates larger than 10 11 generations. To set this into context, we compare this to the evolutionary history of fungi. Fungi first evolved around 1.5 billion years ago [START_REF] Daniel | Divergence time estimates for the early history of animal phyla and the origin of plants, animals and fungi[END_REF]. Assuming approximately 100 generations a year, a conservative estimate considering the doubling time of yeast which is around 90 minutes, this leads to an order of magnitude guess of 10 11 generations for the entire evolutionary history of fungi. Hence, in the context of obligate sex (see Figure 5 top-left panel), mating types alleles are expected to remain fixed in a population while their number remains less than 6. This condition is severely relaxed in populations in which sex is facultative (see for example Figure 5 bottom-right panel). Even though extinction times remain high for realistic effective population sizes, some turnover of mating types might be expected over long time periods, e.g. the evolutionary history of fungi. Moreover, this observation indicates that mating types loss is most likely be driven by extreme population bottlenecks decreasing the effective population size N .

Discussion

The evolutionary mechanisms that drive the number of mating type alleles observed across species has been the topic of numerous theoretical studies. Our analysis in the context of haploid self-incompatibility adds new results and perspectives on the evolutionary dynamics of this type of balancing selection.

A specificity of mating type SI, as studied here, when compared to gametophytic SI, as prevalent in plants, is the possibility to reproduce asexually. Our findings support and extend the previous result that switching between an asexual and a sexual life cycle significantly reduces the number of mating types in a population (Constable and Kokko, 2018). Empirically, available data appears to support the view that more frequent sex is correlated with more mating types (see Table 1 and again Constable and Kokko (2018)). However we are hampered from making any strong empirical claims in this area as a result of a paucity of empirical estimates of rates of sex in natural populations. In most species such estimates are absent, however with new methodologies and understanding for estimating the rate of sex arising increasingly frequently [START_REF] Nieuwenhuis | Repeated evolution of self-compatibility for reproductive assurance[END_REF][START_REF] Ennos | Estimating the number of sexual events per generation in a facultatively sexual haploid population[END_REF][START_REF] Hartfield | Coalescence and linkage disequilibrium in facultatively sexual diploids[END_REF], it is our hope that this gap in the literature will soon be filled.

In Table 1 we test our model quantitatively against four species where estimates for the rate of sex are available. We compare the number of mating types observed empirically with the mode number predicted theoretically. Since estimates for the mutation rate of new mating type alleles and the effective population size parameters are difficult to obtain, we considered a range of values. We find that while facultative sex explains much of the variation in mating type number, there are quantitative disagreements. In particular, for a range of parameters (particularly large effective population size) the number of mating types is overestimated in Saccharomyces cerevisiae, C. reinhardtii and Tetrahymena, while it is underestimated in S. commune.

Given the simplicity of the model that we have proposed, these quantitative disagreements are not unexpected. An interesting open problem is to explore other null models that are capable of explaining these discrepancies. Essentially this task translates to identifying ad- 20)).

Species Estimated c Empirical M Null model M o Q Est M T M (Generations) S. cerevisiae * 1 -1/2000 [1] 2 [5] 9 ≥ M o ≥ 1 Q Est 2 = 1.25 × 10 -4 10 274 ≥ T 2 ≥ 10 6 C. reinhardtii 1 -1/770 [2] 2 [6] 14 ≥ M o ≥ 1 Q Est 2 = 3.25 × 10 -4 10 707 ≥ T 2 ≥ 10 9 Tetrahymena † 1 -1/100 [3] 3 -9 [7] 38 ≥ M o ≥ 4 Q Est 3 = 1.68 × 10 -3 10 1822 ≥ T 3 ≥ 10 20 Q Est 9 = 5.
Table 1: A comparison of empirical values for the number of mating types in various isogamous species with the predictions of the null model. We assume an effective population size of between 10 7 ≥ N ≥ 10 5 in our calculations of M and T M (note that Q Est M is independent of N ). In our calculation for M , we further assume a per-generation mutation rate, m g = mN , of 10 -6 ≥ m g ≥ 10 -8 . For references to empirical data, see the list at the end of the manuscript. For a visualization of the full theoretical distributions of M , P st M , for each parameter set, see Section S8 in the Supplementary Information.

ditional mechanisms that decrease selection for more mating types (where their number is overestimated) or increase this selection strength (where their number is underestimated). There exist a number of biologically reasonable potential candidates.

In Eq. ( 1) we have assumed mass action encounter rate dynamics between gametes, leading to a linear relationship between a mating type's frequency and its probability of finding a sexual partner (i.e (1x i )). This implementation neglects active mate search. In reality, in species such as C. reinhardtii or the diatom Ditylum brightwellii, that have developed active methods to increase the encounter rate between complementary mating types [START_REF] Snell | Flagellar adhesion, flagellar-generated signaling, and gamete fusion during mating[END_REF][START_REF] Waite | Role of sinking and ascent during sexual reproduction in the marine diatom Ditylum brightwellii[END_REF], this term might be more accurately described by a decreasing concave up function (see for example [START_REF] Ashby | Parasitic castration promotes coevolutionary cycling but also imposes a cost on sex[END_REF]). This would qualitatively recapitulate the results of Iwasa and Sasaki (1987) (Mating Kinetics 3 and 4) and lead to a decrease in the predicted number of mating types.

Further, our model assumes that a mating type is determined by a single locus. Thus, our model is restricted to bipolar mating type systems, as opposed to tetrapolar systems where the mating type is determined by alleles at two loci [START_REF] Nieuwenhuis | Evolution of uni-and bifactorial sexual compatibility systems in fungi[END_REF]. We expect that additional loci may lead to larger numbers of mating types maintained in the population (as mating type allele combinations can be regenerated through recombination, the extinction rate of mating types will be reduced). Indeed, this would be consistent with the larger number of mating types empirically observed in the tetrapolar S. commune, as well as the general observation that tetrapolarity is associated with an increase in mating type diversity [START_REF] Nieuwenhuis | Evolution of uni-and bifactorial sexual compatibility systems in fungi[END_REF]. While it is interesting to note that the number of alleles predicted by our model capture the number of A and B self-incompatibility complexes in S. commune separately (see Table 1), theoretically addressing the mutation-extinction hypothesis for tetrapolar systems in a systematic way, along with the effects of encounter rate dynamics, will be interesting areas for future investigations.

Thus far we have focused on modelling considerations that may well strengthen the mutationextinction balance hypothesis. Of course, non-neutral differences between mating types must be considered as well, and a pluralistic view that incorporates the remaining hypotheses for observed mating type numbers may be necessary. For instance, we have not accounted for the fact that newly arising mutants are unlikely to be fully compatible with (or equally as fit as) resident mating types, a biologically important consideration [START_REF] Power | On forces of selection in the evolution of mating types[END_REF]Hadjivasiliou and Pomiankowski, 2016;[START_REF] Krumbeck | Fitness differences suppress the number of mating types in evolving isogamous species[END_REF] that would lower the number of types relative to our idealized estimate. In addition, we have not considered how forms of homothalism (i.e. the evolution of self-compatibility) may affect our results. For instance most yeasts have evolved the ability to switch mating types between sexual generations [START_REF] Nieuwenhuis | The evolution of mating-type switching for reproductive assurance[END_REF][START_REF] Nieuwenhuis | Repeated evolution of self-compatibility for reproductive assurance[END_REF], while among ciliates some species exhibit probabilistic mating type expression [START_REF] Paixao | Sex ratio evolution under probabilistic sex determination[END_REF]. Indeed, as demonstrated in Hadjivasiliou et al. (2016), the rapid loss of mating type diversity when sex is rare and local population sizes are low (see Figure 5) can drive selection for mating type switching in order to maintain the capacity for sexual reproduction. Such factors make a comparison of our model with these species problematic. Further mechanisms that may limit the number of mating types that are not included in our model, including UPI and costly mate search are reviewed in [START_REF] Billiard | Having sex, yes, but with whom? Inferences from fungi on the evolution of anisogamy and mating types[END_REF].

Our model for haploid SI is conceptually similar to extensively studied gametophytic SI systems observed in plants. Yet we take a novel mathematical approach in estimating the number of mating types supported by a finite population. Instead of using the extinction boundaries of the one-dimensional diffusion approximation of a focal mating type as in most previous studies [START_REF] Wright | The distribution of self-sterility alleles in populations[END_REF][START_REF] Wright | On the number of self-incompatibility alleles maintained in equilibrium by a given mutation rate in a population of a given size: a re-examination[END_REF][START_REF] Wright | The distribution of self-incompatibility alleles in populations[END_REF][START_REF] Ewens | On the problem of self-sterility alleles[END_REF][START_REF] Yokoyama | Population dynamics of sex-determining alleles in honey bees and self-incompatibility alleles in plants[END_REF][START_REF] Yokoyama | The expected number of self-incompatibility alleles in finite plant populations[END_REF][START_REF] Vallejo-Marín | On the Evolutionary Modification of Self-Incompatibility: Implications of Partial Clonality for Allelic Diversity and Genealogical Structure[END_REF], we study a birth-death process on the total number of mating types. Utilizing just the local description of the stationary distribution around an interior stable fixed point, we circumvent the problem of the diffusion approximation for systems with a stable deterministic behavior being inaccurate at the boundaries [START_REF] Assaf | WKB theory of large deviations in stochastic populations[END_REF]. This leads to robust predictions for the stationary probability distribution P st M (Figure 1), and, importantly, the mean extinction time of a mating type allele, T M (Figure 4).

In calculating this mean extinction time, we deviate from the approaches taken before which rely on the one-dimensional dynamics of a focal mating type [START_REF] Takahata | A simple genealogical structure of strongly balanced allelic lines and trans-species evolution of polymorphism[END_REF][START_REF] Vekemans | Gene and allelic genealogies at a gametophytic self-incompatibility locus[END_REF]. This previous method compares the one-dimensional diffusion to a timerescaled neutral coalescent, thus obtaining estimates for the diversification rate (i.e. the establishment rate of a novel mating type). This comparison is conducted by assuming constant selective strength. In our model, this would be equivalent to assuming that the term j x i x j is constant. However, with a varying number of present mating types (and indeed, also with fluctuations in the mating type frequencies themselves) this value is not constant in time. Hence, a comparison with a time-rescaled neutral coalescent does not seem appropriate since the selection strength would need to be re-assessed after each coalescence event which is associated with a loss of a mating type. This variation in selection strength has been empirically observed by analyzing the genome of various Coprinus cinereus (mushroom fungus) populations worldwide that show that selective forces in the mating type dynamics are dependent on the number of distinct lineages [START_REF] May | The signature of balancing selection: Fungal mating compatibility gene evolution[END_REF]; see also [START_REF] Richman | Evolution of balanced genetic polymorphism[END_REF] for other examples in the context of balancing selection. Our technique, explicitly computing the mean extinction time by the effective birth-death process on the number of mating types, avoids this problem. Since the varying selective strength enters in both the establishment probability and the stationary distribution, our computed extinction time accounts for the fluctuating selective pressure.

These explicit estimates on establishment probabilities and typical extinction times yield new insight into the dynamics that drive the low number of mating types predicted by our model in the rare sex regime. One might initially suspect that these numbers are the result of a high turnover of mating type alleles, i.e. frequent extinction and invasion events. However we find that in fact they are a result of very low invasion probabilities for novel mating types, combined with rapidly decreasing extinction times as a function of resident mating type diversity. It is worth mentioning that although our approximations for the extinction time break down when the resident state becomes highly unstable (i.e. when M M o ), these represent states that would be very short lived in nature, and thus are not relevant from a biological perspective.

Finally, in Table 1 we link our theoretical observation of small numbers of mating type alleles yielding very large extinction times (see Figure 4) to some species examples. Indeed, the large extinction times found by the quantitative study are in line with previous empirical observations of long terminal branches in allelic genealogies under negative-frequency dependent selection, e.g. in fungi [START_REF] May | The signature of balancing selection: Fungal mating compatibility gene evolution[END_REF] and Solanaceae [START_REF] Uyenoyama | Genealogical structure among alleles regulating self-incompatibility in natural populations of flowering plants[END_REF]. Similar to previous simulation studies [START_REF] Slatkin | Overdominant alleles in a population of variable size[END_REF][START_REF] Gervais | Origin and diversification dynamics of self-incompatibility haplotypes[END_REF] we find that the number of resident mating types strongly influences the diversification rate; the larger the resident number, the lower the diversification rate. This slowdown is ubiquitously observed in natural systems under balancing selection, such as self-incompatibility in fungi [START_REF] May | The signature of balancing selection: Fungal mating compatibility gene evolution[END_REF] or Solanaceae [START_REF] Uyenoyama | Genealogical structure among alleles regulating self-incompatibility in natural populations of flowering plants[END_REF] and MHC-systems [START_REF] Solberg | Balancing selection and heterogeneity across the classical human leukocyte antigen loci: A meta-analytic review of 497 population studies[END_REF]. We find that long terminal branches in self-incompatibility systems, corresponding to old mating type allele ages (with some being older than the corresponding species age [START_REF] Richman | Evolution of balanced genetic polymorphism[END_REF]) emerge naturally under balancing selection for two reasons: (i) the stability of the internal equilibrium can lead to extremely large extinction times of alleles thus enabling such long branches and (ii) the lower establishment rate due to a large number of resident SI alleles slows down the diversification rate and hence decreases the number of newly established mating types.

In conclusion, we analyzed the evolutionary dynamics of self-incompatible mating types in facultatively sexual isogamous species. Our results refine previous estimates on the number of mating types maintained in a finite population as well as on the establishment probability of a newly arising mating type allele. Furthermore, we have used these results to compute the mean extinction time of a focal mating type via an effective birth-death process describing the number of mating types. This estimate naturally incorporates variation in selection strength due to varying numbers of resident mating types, a fact that previous studies have failed to incorporate. We are therefore able to qualitatively explain the empirically and numerically observed slowdown of allelic diversification rates in populations under balancing selection. The here presented methodology is theoretically extendable to other systems exhibiting negativefrequency dependent selection.

Acknowledgements

We thank the Department of Evolutionary Theory at the Max Planck Institute for Evolutionary Biology in Plön for inviting GWAC, resulting in this collaboration. PC received funding from the Agence Nationale de la Recherche, grant number ANR-14-ACHN-0003-01 provided to Florence Débarre. GWAC thanks the Leverhulme Early Career Fellowship provided by the Leverhulme Trust for funding. We also thank Deborah and Brian Charlesworth for their insightful feedback when presenting this work. We are also grateful to Sylvain Billiard for directing our attention towards diversification rates of mating types and for comments on the manuscript which increased its readability. Lastly, the suggestions of two anonymous reviewers helped to clarify certain aspects of the manuscript. Most notably this improved the discussion around equation ( 18).

6 Citations for Table 1.

[1] As addressed within the discussion, our model is only strictly appropriate for populations in which the mating type does not switch between generations. Mating type switching is prevalent among many Saccharomyces, however most populations feature non-switching strains [START_REF] Nieuwenhuis | Repeated evolution of self-compatibility for reproductive assurance[END_REF]. We assume S. cerevisiae has a similar ecology to S. paradoxus, where sex has been estimated to occur between once every thousand and three thousand generations [START_REF] Tsai | Population genomics of the wild yeast Saccharomyces paradoxus: Quantifying the life cycle[END_REF] (indeed, estimates for the rates of outcrossing in both species are comparable).

[2] The rate of sex in C. reinhardtii has been estimated to once in every 770 generations [START_REF] Hasan | The genomic landscape of recombination rate variation in Chlamydomonas reinhardtii reveals effects of linked selection[END_REF].

[3] We assume that synclonal Tetrahymena species (such as T. americanis, see [7]) feature similar rates of sex to T. thermophila. Our estimate on the rate of sex is based on a minimal sexually immature period of 60 -100 fissions in T. thermophila [START_REF] Doerder | High frequency of sex and equal frequencies of mating types in natural populations of the ciliate Tetrahymena thermophila[END_REF]. We note however that the frequency of cellular conjugation events observed in other ciliates indicates that lower rates of sex may be appropriate [START_REF] Lucchesi | How often does conjugation in ciliates occur? Clues from a seven-year study on marine sandy shores[END_REF].

[4] Molecular analysis suggests that S. commune features some of the highest rates of sex within the fungal kingdom [START_REF] Nieuwenhuis | The frequency of sex in fungi[END_REF]. In addition, most of the Agaricomycotina (the class to which S. commune belongs) are known to be obligately sexual [START_REF] Nieuwenhuis | Sexual selection in fungi[END_REF]. We therefore assume a rate of asexual reproduction of c = 0 for S. commune. We note however that the life-cycle of S. commune features more complicated dynamics than accounted for by our model; it is multicellular and can exhibit vegetative growth as a haploid mycelium.

[5] The mating type of Saccharomyces species is determined by one of two alleles at the MAT locus, termed α and a [START_REF] Butler | The evolution of MAT: the ascomycetes[END_REF]. We note that often these variants are referred to as idiomorphs (rather than alleles) because of large differences in their size and genetic composition [START_REF] Butler | Evolution of the mat locus and its ho endonuclease in yeast species[END_REF].

[6] C. reinhardtii has two mating types, denoted mt + and mt -, each determined by one of two alleles at the MT locus. Similar mating type determination is found among other isogamous algae [START_REF] Hamaji | Distribution of the sex-determining gene mid and molecular correspondence of mating types within the isogamous genus gonium (volvocales, chlorophyta)[END_REF][START_REF] Sekimoto | Sexual reproduction and sex determination in green algae[END_REF].

[7] Tetrahymena species T. americanis, T. hegewischi, T. hyperangularis and T. pigmentosa have mating type numbers nine, eight, four and three respectively [START_REF] Phadke | Rapid diversification of mating systems in ciliates[END_REF]. We take care to focus only on these synclonal species, defined as those in which mating type is inherited deterministically based on parental genotype. In caryonidal species, such as T. thermophila, the mating type of progeny is determined stochastically [START_REF] Phadke | Rapid diversification of mating systems in ciliates[END_REF], and thus a comparison of our model with these species is not appropriate.

[8] The mating type of haploid S. commune is determined by two complementary pathways, controlled by unlinked genetic complexes A and B . Each of these regions consists of two weakly recombining loci, leading to a total of four mating type loci. These loci, denoted Aα, Aβ, B α, and B β, have respectively 9, 32, 9 and 9 alleles [START_REF] Stankis | The Aα mating locus of Schizophyllum commune encodes two dissimilar multiallelic homeodomain proteins[END_REF]. Full compatibility between two mating types is achieved when both the A and B complexes of the mates are of different specificities; semi-compatibility occurs when the only one of these complexes is of different specificity [START_REF] Raudaskoski | Basidiomycete mating type genes and pheromone signaling[END_REF].

Supplementary Information: Invasion and extinction dynamics of mating types under facultative sexual reproduction In this section we will calculate an approximate expression for the quasistationary distribution of mating type alleles around one of the fixed points (deterministic equilibria). We expect this approximation to be valid on timescales shorter than those of mutation or extinction events.

The approximation that we will use is equivalent to the Linear Noise Approximation (LNA) or a central limit theorem for Markov processes, see [START_REF] Stewart | Markov processes: characterization and convergence[END_REF][START_REF] Ng Van Kampen | Stochastic processes in physics and chemistry[END_REF]. The procedure is as follows: First we apply a diffusion approximation to Eq. ( S3), to obtain a non-linear advection-diffusion equation for the frequency of mating type alleles. Second we will linearize this resultant equation about its fixed point value, to obtain a form of the equation that is amenable to analytical simplifications. Since we are interested in the distribution of mating type alleles close to the fixed point of a focal mating type allele, we will assume both that m = 0 and that there are no extinctions (so that M is fixed to its initial value) for the remainder of this section.

We begin by applying the diffusion approximation to Eq. ( S3); we introduce variables x i = n i /N that measure the frequency of a given mating type allele in the population. The variables x i are approximately continuous when the population size N is large. We transform into these variables and conduct a Taylor expansion of Eq. (S3) in the small parameter 1/N . Truncating at next to leading order, we obtain the following non-linear advection-diffusion equation, the Fokker-Planck equation (FPE),

∂p(x, t) ∂τ = - M i=1 ∂ ∂x i [A i (x)p(x, t)] + 1 2N M i,j=1 ∂ 2 ∂x i ∂x j [B ij (x)p(x, t)] , (S9)
describing the time-evolution of the continuous probability distribution of the variables x, p(x, t). The expressions for the advection vector A(x) and diffusion matrix B(x) can be uniquely defined from the underlying probability transition rates Eq. (S1) using well-practiced standard methods [START_REF] Ng Van Kampen | Stochastic processes in physics and chemistry[END_REF].

In the current notation we find that the advection vector in Eq. (S9) is given by

A i (x) = M j =i T ij -T ji , = M j =i x i c + (1 -c) 2 (1 -x i ) x j -x j c + (1 -c) 2 (1 -x j ) x i , = 1 -c 2 M j=1 x i x j (x j -x i ) . ( S10
)
As ẋ = A(x) in the N → ∞ limit (see [START_REF] Mckane | Stochastic pattern formation and spontaneous polarisation: the linear noise approximation and beyond[END_REF]), this vector essentially gives the deterministic dynamics of the population. Therefore solving A(x * ) = 0 allows us to calculate the deterministic fixed point of the dynamics, which we find lies at x * i = 1/M for i = 1, . . . , M . Meanwhile the diagonal elements of the diffusion matrix are given by

B ii (x) = M j =i T ij + T ji , = M j =i x i c + (1 -c) 2 (1 -x i ) x j + x j c + (1 -c) 2 (1 -x j ) x i , = M j =i x i x j 2c + 1 -c 2 (2 -x i -x j ) , (S11)
while the off-diagonal entries are given by

B ij (x) = -(T ij + T ji ) , ∀i = j , = -x i c + (1 -c) 2 (1 -x i ) x j + x j c + (1 -c) 2 (1 -x j ) x i , = -x i x j 2c + 1 -c 2 (2 -x i -x j ) . (S12)
With the non-linear FPE for p(x, t) now defined, we proceed to linearize the system about its deterministic fixed point.

As addressed, the deterministic fixed point of the system with M mating types is given by x * i = 1/M for i = 1, . . . , M . We assume that the population relaxes to the vicinity of this fixed point before any mutation or extinction events have had time to occur. If N is large, the population will then approach a quasi-stationary distribution around this fixed point on an intermediate timescale. Fluctuations about the fixed point will be of order 1/ √ N due to the central limit theorem for density dependent Markov processes [START_REF] Stewart | Markov processes: characterization and convergence[END_REF].

The LNA [START_REF] Ng Van Kampen | Stochastic processes in physics and chemistry[END_REF] utilizes this fact by making the change of variables

x i = x * i + 1 √ N ξ i (S13)
to linearize the FPE around the fixed point. Neglecting terms of order 1/N or lower, we obtain the following FPE for φ(ξ, t), the probability distribution of ξ (see [START_REF] Ng Van Kampen | Stochastic processes in physics and chemistry[END_REF]; Eq. ( 6.4) and surrounding discussion):

∂φ(ξ, t) ∂τ = - M i,j=1 J ij ∂ ∂ξ i [ξ j φ(ξ, t)] + 1 2 M i,j=1 B ij (x * ) ∂ 2 ∂ξ i ∂ξ j [φ(ξ, t)] , (S14)
where J is the Jacobian matrix of A(x),

J ij = ∂A i ∂x j x=x * (S15)
As addressed, we are interested in obtaining the stationary distribution of fluctuations about the fixed point (the quasi-stationary distribution of x around x * ). Our first step is to obtain the stationary distribution φ st (ξ) that is the solution to Eq. (S14) given a constant number of mating types M (the quasi-stationarity assumption);

- M i,j=1 J ij ∂ ∂ξ i ξ j φ st (ξ) + 1 2 M i,j=1 B ij (x * ) ∂ 2 ∂ξ i ∂ξ j φ st (ξ) = 0 . (S16)
Since Eq. ( S16) is linear, φ st (ξ) is normally distributed [START_REF] Ng Van Kampen | Stochastic processes in physics and chemistry[END_REF], with mean 0 and a covariance matrix, Σ, that is the solution to the following Lyapunov equation (see [START_REF] Horn | Topics in Matrix Analysis[END_REF]); JΣ + ΣJ + B(x * ) = 0 .

(S17)

We now must determine the form of the matrices J and B(x * ), which can be calculated from Eqs. (S10-S12). Substituting x M = 1 -M -1 j=1 x j we find that the Jacobian matrix is diagonal with

J ii = - 1 -c 2 1 M , ∀i , J ij = 0 , ∀i = j , (S18)
where we have again absorbed any terms independent of M into the normalization factor M. We can see that Eq. (S37) matches the full expression in Eq. (S35) very well (see also Figure 1 in the Main Text).

S6 Breakdown of the mean extinction time calculation

In the main text, we calculate an approximation for the mean extinction time as a function of the population size, N , rate of asexual reproduction c and number of resident mating types, M (see Eq. ( 18)). We also show that in certain regions of parameter space, this approximation breaks down (see Eq. ( 20)). In these parameter regions, in which genetic drift dominates the dynamics, the mean time to extinction of a resident allele is better approximated by the neutral multi-allelic Moran model (see below).

In the following we show that it is the inaccuracies of our quasi-stationary linear approximation for the distribution of mating type allele frequencies about deterministic fixed points (see Section S2) that drives this break down. This approximation assumes that a given fixed point of the deterministic system is sufficiently stable that the distribution of allele frequencies around this fixed point can be well captured by a Gaussian (normal) distribution (see Eq. ( S25)). Clearly this local description of the distribution does not account for the behaviour of the system at the extinction boundaries. When population sizes are high, rates of asexual reproduction low, and the number of resident mating types small, the probability mass predicted by the Gaussian approximation at these extinction boundaries is negligible. In this regime our approximation continues to be accurate, as illustrated in Figure S1. However, outside this range (i.e when N is low, c high, or M large), the variance of the Gaussian approximation becomes sufficiently large that nonnegligible probability mass is predicted at the boundary (see Figure S3). These parameter regions correspond to areas where drift dominates the dynamics.

In regions of parameter space where the Gaussian approximation becomes inaccurate, we expect our expression for the distribution of mating types, P st M , (which relies on the Gaussian approximation, see Eq. (S28) and the subsequent calculations) to also become inaccurate. However, as these inaccuracies occur in biologically less interesting regions of parameter space where mating types are frequently becoming extinct, they do not affect the dominant modes of P st M (see Figure 1 in the Main Text). In contrast, in our investigation of extinction times we are sometimes probing very unstable configurations of the resident mating types (see Figure 5 in the Main Text), and thus the approximation breaks down in these parameter regions.

We now show that it is indeed the Gaussian approximation for the quasi- S27)). We consider two mating types in a fairly small population, M = 2, N = 200. For c = 0, the approximation remains physically reasonable. As c is increased however, a non-negligible amount of probability mass builds up on the boundary relative to that at the fixed point around which the distribution is centered (that is, the ratio of P qst(2) (η (1) ) to P qst(2) (η (2) ) tends to one). A similar pattern can be observed as N decreases and M increases.

stationary distributions of mating type allele frequencies about deterministic fixed points that drives the breakdown of the extinction time calculation. In Figure S4 we plot the probability mass predicted by the Gaussian distribution at an extinction boundary relative to the distribution's value at the corresponding fixed point. This quantity should be very low for the approximation to remain physically reasonable (i.e. there should be very little mass at the extinction boundary), while it approaches one as the predicted distribution becomes increasingly flat and inaccurate. We observe that high values of predicted probability mass at the extinction boundary coincide with regions where our prediction for the mean time to extinction, Eq. ( 18), break down, and the neutral theory becomes more appropriate (see Eq. ( 20)). S27)). The colorbar indicates the probability mass predicted by the quasi-stationary distribution at the nearest extinction boundary, relative to the probability mass at the center of the quasi-stationary distribution (i.e. at the deterministic fixed point). Low values (in white) are associated with regions of parameter space where the quasi-stationary distribution remains reasonable. High values (in deep blue) are associated with regions where the quasi-stationary distribution becomes a poor approximation (see Figure S3). Areas to the top left of the red dashed lines are those where the conditions in Eq. ( 18) are violated (that is, our approximation for the mean time to fixation breaks down).

Figure 1 :

 1 Figure 1: Stationary distribution of mating types, P stM . For any given rate of asexual reproduction, c, larger populations (larger N ) contain more mating types on average. Meanwhile increasing the rate of asexual to sexual reproduction (increasing c) decreases the expected number of mating types. Analytic results (solid lines) are obtained by evaluating Eq. (9). Simulation results (corresponding shaded regions) are obtained using the Gillespie algorithm[START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF], averaged over 5 × 10 6 generations, sampled every 10 2 generations and following a relaxation period of 10 3 generations from initial conditions with one more mating type than the mode number predicted analytically by the model (see Eq. (9)). The mutation rate is m = 10 -8 in all plots. The analytic description can be seen to capture the simulated behavior of the model with very high accuracy.

Figure 2 :

 2 Figure 2: Analytic results on the mode number of mating types, M o , as a function of the population size, N and the per-generation mutation rate, m g = N m, under differing rates of asexual reproduction, c. As the rate of sexual to asexual reproduction is decreased [panels (a-d)] so too does the expected number of mating types. When sex is obligate [panel (a)], of the order of hundreds of mating types are expected. When sex is facultative and very rare, occurring approximately once to every thousand asexual reproduction events [panel (d)] far fewer mating types are observed. Results are obtained by numerically solving Eq. (12).

Figure 3 :

 3 Figure 3: Establishment probability of a novel mating type as a function of M for varying rates of asexual reproduction, c. Analytic results are obtained by evaluating Eq. (13).Simulation results are obtained from Gillespie simulations; in each plot N = 10 5 and results are averaged over 10 3 runs (c = 0), 10 4 runs (c = 0.9 and c = 0.99), and 10 5 runs (c = 0.999). Error bars are ±2 times the standard-deviation of the sample binomial distribution.

Figure 4 :

 4 Figure 4: Mean extinction time of the first mating type allele, T M , as a function of N for populations with initially M = 2, M = 3 and M = 4 mating types. Analytic results (see Eq. (20)) are plotted as solid lines. Simulation results are averaged over 10 3 Gillespie simulations, with error bars indicating the standard deviation of results. Vertical dashed lines indicate the transition between the approximation of Eq. (18) (right of dashed line) and the neutral limit of Eq. (19) (left of dashed line), defined by Eq. (20).As the rate of asexual reproduction is increased (from c = 0.1 in panel (a) to c = 0.9 in panel (b)), the mean extinction time drops rapidly (note differing scales on x-axis).

Figure 5 :

 5 Figure 5: Analytic predictions for the mean extinction time of a resident mating type allele as a function of N and M (see Eq. (20)). The blue dashed line indicates the parameter regime in which extinction rates become approximately neutral (Eq. (19)). Extinction times in the gray shaded region exceed 10 11 generations. As these time are approximately longer than the evolutionary history of fungi, they have been omitted for clarity.

Figure S3 :

 S3 Figure S3: Figure illustrating the breakdown in the quasi-stationary approximation for the distribution of alleles around a fixed point, P qst (n) (see Eq. (S27)). We consider two mating types in a fairly small population, M = 2, N = 200. For c = 0, the approximation remains physically reasonable. As c is increased however, a non-negligible amount of probability mass builds up on the boundary relative to that at the fixed point around which the distribution is centered (that is, the ratio of P qst(2) (η (1) ) to P qst(2) (η (2) ) tends to one). A similar pattern can be observed as N decreases and M increases.

Figure S4 :

 S4 Figure S4: Figure illustrating the coincidence of the breakdown in the quasistationary approximation for the distribution of alleles around a fixed point, P qst (n) (see Eq. (S27)). The colorbar indicates the probability mass predicted by the quasi-stationary distribution at the nearest extinction boundary, relative to the probability mass at the center of the quasi-stationary distribution (i.e. at the deterministic fixed point). Low values (in white) are associated with regions of parameter space where the quasi-stationary distribution remains reasonable. High values (in deep blue) are associated with regions where the quasi-stationary distribution becomes a poor approximation (see FigureS3). Areas to the top left of the red dashed lines are those where the conditions in Eq. (18) are violated (that is, our approximation for the mean time to fixation breaks down).

  T 81 ≥ 10 4 * Considering only strains of S. cerevisiae that do not feature mating type switching, for which our model is appropriate. † Considering only Tetrahymena species with synclonal inheritance, for which our model is appropriate. ‡ Tetrapolar with mating type compatibility determined by genetic complexes A (with 288 variants) and B (with 81 variants). Lower-bound extinction time is effectively neutral (see Eq. (

				58 × 10 -4	10 153 ≥ T 9 ≥ 10 4
	S. commune	0 [4]	23, 328 ‡ [8]	510 ≥ M o ≥ 48 Q Est 23,328 = 4.29 × 10 -5 41 ≥ T 23,328 ≥ 0.5 Q Est 10 26 ≥ T 288 ≥ 58 288 = 3.47 × 10 -3 Q Est 81 = 1.23 × 10 -2 10 337 ≥

S1 Master equation of the null model

In this section we present the full mathematical description of the model. As discussed in the main text, the function T ij (see Eq. (1) in the Main Text) gives the probability per unit time that the number of individuals carrying the mating type i allele in the population increases by one and that of type j decreases by one from an initial state n = (n 1 , n 2 , . . . , n i , . . . , n j , . . .), e.g.

For ease of notation, we now introduce the probability transition rate T (n |n), which gives the probability per unit time that a population transitions to a state n from a state n. For our model, the two quantities are related as

The probability of being in a state n at time t, P n (t), evolves according to the master equation [START_REF] Ng Van Kampen | Stochastic processes in physics and chemistry[END_REF], which can be compactly expressed as This equation can be intuitively understood as follows: The probability of being in a state n increases with the probability that the population transitions into state n from a state n but decreases with the probability that the population was already in state n and transitioned out of it. For arbitrary initial conditions, it is difficult to solve Eq. ( S3), for P n (t). A simpler quantity is the stationary probability distribution P st n , to which the population relaxes on very long timescales. For a time-homogeneous process P st n is given by the solution to the set of difference equations

In [START_REF] Constable | The rate of facultative sex governs the number of expected mating types in isogamous species[END_REF], it was shown that an analytic solution for P st n was obtainable for the model defined by Eq. (S2) (see also, [START_REF] Constable | The rate of facultative sex governs the number of expected mating types in isogamous species[END_REF]: Supplementary Information). This is given in the Main Text, Eq. ( 6); denoting by n ↓ the vector n reordered with its largest elements first, the stationary distribution of the population composition can be expressed as

where M is the number of non-zero elements of n (i.e. the number of mating types in state n), and b(k) and d(k) are given by

for all n j , (see also, Eq. ( 4) in the Main Text) and

The parameter M max marks the length of the vector n (i.e. M max is the number of distinct possible mating types in the model, which may exceed N ).

The factor 1/(M max -M ) in the mutation term (which is suppressed in the main text) is an accounting term that ensures that the total population level mutation rate is m, i.e.

where δ(n i -0) is the Dirac delta function. As we take the limit M max → ∞ (see Section S3), the appearance of this factor in our transition rates does not alter our results. The validity of Eq. (S5) as a solution can be demonstrated by its direct substitution into Eq. (S4). while the diffusion matrix evaluated at the fixed point is

Since the Jacobian matrix can be expressed by

where I is the identity matrix, Eq. (S17) can be simplified to

We also note that since the original model is a Moran type model in which the number of individuals N is constant, J, B(x * ) and Σ are all (M -1) × (M -1) matrices. Therefore the determinant of Σ can be calculated to be

This can be verified by solving the matrix of interest, i.e. for a given M with a computer algebra program, e.g. Mathematica.

In summary, we find that the stationary distribution of ξ is given by

where 0 is the zero vector and Σ is given by Eq. (S22). Re-expressing this distribution in terms of our x variables we find

at intermediate times in the region of the fixed point, as illustrated in Figure S1.

Here we see that as the number of mating types in the population, M , increases, the mean of the quasi-stationary PDF, x * = 1/M , approaches the extinction boundaries. Simultaneously, increased M leads to increased variance in the quasi-stationary PDF (see Eqs. ( S19) and ( S22)), though the co-variance of the distribution decreased (fluctuations in each mating type frequency become less correlated, see Eqs. (S20) and ( S22)). The variance of the PDF also increases as the rate of asexual to sexual reproduction increases, and fluctuations around the fixed point frequency increase in magnitude (see Eqs. (S19) and (S22)). Our final step is to re-express Eq. (S25) in terms of the number of alleles of each type in the population. Let η (M ) be a vector giving an approximate, potentially non-discrete, value of n in the region of a deterministic fixed point with M mating types, e.g.

For clarity, we choose to write this quasi-stationary distribution as

where the superscript (M ) emphasizes that the form of this quasi-stationary distribution and the covariance matrix change with the number of mating types present at the fixed point.

S3 Approximation of the stationary distribution of the number of mating types, P st

M

In the previous two sections, we provided a solution for P st n , the stationary distribution of the population composition, and for P qst(M ) (n), the quasistationary distribution of a focal mating type. We will need both of these solutions to obtain an expression for the stationary distribution of the number of mating types, P st M -our actual object of interest. These stationary distributions of the population composition and the stationary distribution of the number of mating types are related by

where S (M ) is the set of all vectors n with M non-zero elements. In other words, to obtain P st M we need to sum our expression for P st n over all the states that contain M mating types. As this calculation is cumbersome, we seek an analytic approximation.

We consider the biologically reasonable parameter regime in which the population size, N , is large and the per-generation mutation rate, m g = mN , is small. Under these conditions the population quickly relaxes to a distribution in the region of a deterministic fixed point following an extinction or invasion event. We will use the expression in Eq. (S27) to approximate this quasistationary distribution in the region of the fixed point. By doing so, we can replace the sum in Eq. (S28) by a sum of Gaussian distributions (see Figure S2).

The calculation used to obtain Eq. ( S27) is based on a linearization of the population dynamics around a deterministic fixed point. As such the expression contains no information about how much more likely the population is to be in the region of a given fixed point with M mating types as opposed to in the region of a fixed point M + 1 mating types. In order to renormalize each of the normal distributions by the appropriate amount, we "pin" the peak of each normal distribution (see Eq. (S27)) to the height of the full distribution, P st n from Eq. (S5), evaluated in the region of the relevant fixed point. Therefore, we renormalize the height of the Gaussian distribution Eq. (S27) to the height of P st n at n = η (M ) (see Eq. ( S26)). Thus, the probability distribution in the region of a deterministic fixed point with M

+ . . . mating types can be described by the following function

(see Eq. (S27)). Now, the probability of being in the region of a specific fixed point with M mating types is simply the above function integrated over all n:

Eq. (S30) approximates the probability of being in the region of a particular fixed point with M mating types. The probability of the population having M mating types is therefore given by Eq. (S30) multiplied by the number of fixed points with M mating types. Recall that we introduced M max as a temporary parameter capturing the length of the vector n (we will shortly take the limit M max → ∞). The number of ways of choosing a fixed point containing M distinct present mating types from M max potential mating types is then given by the Binomial coefficient, such that

We now proceed to calculate P st η (M ) explicitly. By substitution of the functions b(k) and d(k) into Eq. (S5), we find after some algebraic simplification [see [START_REF] Constable | The rate of facultative sex governs the number of expected mating types in isogamous species[END_REF]: Eq. (S35)] that

It is therefore straightforward to show that

where

.

Upon substitution of Eq. (S32) into Eq. (S31), we find

As P st M is proportional to the function on the right hand side of Eq. (S33), we can divide through this function by a constant. We choose the binomial coefficient Mmax 1 . The terms involving M max in Eq. (S33) can then be considered separately, and the limit M max → ∞ taken;

Thus, for M max → ∞, our final expression for the approximate stationary distribution of the number of mating types is

where we have also absorbed any constant terms that do not involve M (e.g. [2cN/(1c)]!) into the normalization factor M, which is defined so that

Eq. ( S35) can be further simplified by noting that if N is large the terms in the final two factorials are large when M ≥ 2. We may then express these factorials using the Stirling approximation [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF];

Substituting these expressions into Eq. (S35), we find after some algebra that for M ≥ 2,

S4 The mode of the stationary distribution of the number of mating types

Since the distribution P st M is unimodal (see Figure 1, Main Text), determining the mode of P st M , amounts to obtaining the first value of M for which

This can be obtained with a simple numerical algorithm. Alternatively, we can assume that M is approximately continuous and solve

We introduce r M as

Solving Eq. ( S39) is equivalent to finding the root of the following expression

To this end we seek to simplify our expression for r M . Note that as r M involves the ratio of terms in P st M , the normalization factor M cancels. We obtain

Substituting Eq. (S42) into Eq. (S41), we find that a single root exists for M that gives the mode number of mating types that is straightforward to compute numerically with a standard root finding algorithm.

S5 Establishment probabilities

In the following we compute the establishment probability of a novel mating type in the population. We consider the case of facultative sex, i.e. the probability for asexual reproduction is given by c ∈ (0, 1]. Since the internal fixed point is stable in this case, see for example [IS87, Mating kinetics I], the calculation of the survival probability of the invading type gives a reasonable approximation of the establishment probability. Initially, when the mutant is still rare, the mutant individuals evolve independently. Hence, the mutant dynamics can be described by a branching process. For a general introduction on branching processes we refer to [START_REF] Haccou | Branching Processes: Variation, Growth, and Extinction of Populations. Cambridge Studies in Adaptive Dynamics[END_REF][START_REF] Linda | An Introduction to Stochastic Processes with Applications to Biology[END_REF].

The birth and death probability for a rare invading mating type in a population of M resident mating types can be written in terms of the transition probabilities given in Equation (S1). Assuming that the novel mating type has the index M + 1 and is present in k copies we find that it increases by one with rate

and decreases by one with rate

Assuming that the resident mating types are in the stationary state, i.e. setting n j = N M for all resident mating types and therefore having M j=1 n j /N = 1 we find

and

Using the formula for the probability of extinction in a birth-death process (see e.g. [START_REF] Linda | An Introduction to Stochastic Processes with Applications to Biology[END_REF]Theorem 6.2])

we see that it depends on the ratio of the rates given in equations ( S45) and (S46). We obtain

(S48) Then equation (S47), using the geometric series, simplifies to ) This results in a survival probability of

Since Q Est M equals the survival probability of this process we have found the expression in Eq. ( 13) from the Main Text.

We note that for c = 1, i.e. exclusively clonal reproduction, our model reduces to a multi-allelic Moran model. In this setting the establishment probabilities can not be calculated via the here implemented approach.

S7 Neutral mean extinction time

We derive an approximation for the neutral mean extinction time as given in Eq. ( 19) in the Main Text. Therefore, we identify our process with the multi-allelic neutral Moran model. This model has been analysed in detail in [START_REF] Gareth | Exact solution of the multi-allelic diffusion model[END_REF]. Using their Eq. ( 44) and plugging in our parameter values, i.e. x(0) = x * = 1/M and r = 1 we find

which is the result in Eq. ( 19) from the Main Text. Note, that there is a time-scale difference of 1/2 between their model (Wright-Fisher diffusion) and our implementation (Moran model) explaining. Furthermore, we consider the dynamics on the original time-scale resulting in the factor N in front of the sum (see also their comment preceding their Eq. ( 2)).

S8 The stationary distribution P st M for the parameters given in Table 1.