## Abstract

Evolutionary transitions between male and female heterogamety are common in both vertebrates and invertebrates. Theoretical studies of these transitions have found that, when all genotypes are equally fit, continuous paths of intermediate equilibria link the two sex chromosome systems. This observation has led to a belief that neutral evolution along these paths can drive transitions, and that arbitrarily small fitness differences among sex chromosome genotypes can determine the system to which evolution leads. Here, we study stochastic evolutionary dynamics along these equilibrium paths. We find non-neutrality, both in transitions retaining the ancestral pair of sex chromosomes, and in those creating a new pair. In fact, substitution rates are biased in favor of dominant sex determining chromosomes, which fix with higher probabilities than mutations of no effect. Using diffusion approximations, we show that this non-neutrality is a result of “drift-induced selection” operating at every point along the equilibrium paths: stochastic jumps off the paths return with, on average, a directional bias in favor of the dominant segregating sex chromosome. Our results offer a novel explanation for the observed preponderance of dominant sex determining genes, and hint that drift-induced selection may be a common force in standard population genetic systems.

IN most animals, sex is determined genetically (Bull 1983). Among those animals with genetic sex determination, the majority exhibit heterogametic sex determination: the presence or absence of a sex-specific chromosome triggers sexual differentiation (Bull 1983; Beukeboom and Perrin 2014). Depending on whether the sex-specific chromosome is in males or in females, the system is, respectively, male heterogamety (XX females and XY males) or female heterogamety (ZW females and ZZ males).

The system of heterogamety is a fundamental genetic property of a species. It is therefore surprising that it is evolutionarily very labile, with transitions between male and female heterogamety having occurred frequently in amphibians (Hillis and Green 1990), reptiles (Ezaz *et al.* 2009; Pokorna and Kratochvíl 2009), and fishes (Ezaz *et al.* 2006; Mank *et al.* 2006; Mank and Avise 2009), as well as in invertebrates (Kaiser and Bachtrog 2010; Vicoso and Bachtrog 2015; Becking *et al.* 2017). A striking example of a recent transition is found in the frog *Rana rugosa*; populations in northern Japan exhibit female heterogamety, while populations in southern Japan exhibit male heterogamety (Nishioka *et al.* 1993; Miura *et al.* 1998).

In a classic theoretical study, Bull and Charnov (1977) showed that continuous paths of population genetic equilibria exist between male and female heterogamety. States along these paths are equilibria in the sense that the evolutionary dynamics of an infinite, randomly mating, population are stationary at them when all genotypes are equally fit. Intermediate states along the paths involve the presence of multiple genotypes for each sex, a situation observed in several species, *e.g.*, the platyfish *Xiphophorus maculatus* (Kallman 1965, 1968), the blue tilapia *Oreochromis aureus* (Lee *et al.* 2004), a Lake Malawi cichlid *Metriaclima pyrsonotus* (Ser *et al.* 2010), the western clawed frog *Xenopus tropicalis* (Roco *et al.* 2015), and the housefly *Musca domestica* (Franco *et al.* 1982; Feldmeyer *et al.* 2008; Meisel *et al.* 2016).

Two of these equilibrium paths are of particular interest. The first, which we shall refer to as “model 1,” governs those transitions between male and female heterogamety that involve the same pair of sex chromosomes (Bull and Charnov 1977) (see Figure 1 and Figure 3A). The second, which we shall refer to as “model 2,” governs those transitions between male and female heterogamety that involve different pairs of sex chromosomes, *i.e.*, where the sex chromosome pair in one system is autosomal in the other, and vice versa (Scudo 1964, 1967; Bull and Charnov 1977) (see Figure 2 and Figure 3B).

The existence of these deterministic equilibrium paths has led to a belief that neutral drift along them in finite populations could be responsible for transitions between male and female heterogamety (Bull 1983, 1987; van Doorn 2014)—an important baseline model for such transitions (van Doorn 2014). Moreover, arbitrarily small fitness differences between the sexual genotypes can eliminate the equilibrium paths under deterministic evolutionary dynamics, and render one of the heterogametic systems stable and the other unstable (Bull and Charnov 1977). This has led to a belief that small fitness differences alone can determine which transitions are possible (Bull 1983, 1987; van Doorn 2014).

Here, we examine these claims in the context of finite-population evolutionary dynamics. Using both stochastic simulations and analytical approximations based on the removal of fast variables, we estimate the fixation probabilities of the various sex determining mutations along the two equilibrium paths, starting from a state of simple heterogamety. We find that evolution along these “neutral” equilibrium paths is in fact not neutral, instead showing a clear bias in favor of dominant sex-determining mutations. Selection for otherwise deterministically neutral genotypes has previously been recognized in other settings (Verner 1965; Gillespie 1977; Taylor and Sauer 1980; Parsons *et al.* 2010; Lin *et al.* 2012; Kogan *et al.* 2014; Constable *et al.* 2016; Chotibut and Nelson 2017). Perhaps the most prominent example is Gillespie’s criterion (Gillespie 1977), which states that if the reproductive rates of two genotypes have equal arithmetic mean but different variance, then the genotype with lower variance will be selected for, owing to higher geometric mean fitness. A key difference is that, in the neutral models we study, we assume no *a priori* differences in the reproductive rates of the various genotypes; these instead emerge naturally in our analysis of the dynamics along, and around, the equilibrium paths.

When all genotypes are equally fit, we find in both models that the substitution rates of the dominant sex-determining mutations (in directions *a* and *c* in Figure 3) are higher than the substitution rates of recessive mutations (in directions *b* and *d* in Figure 3). Thus, in finite populations, stochastic evolutionary dynamics have a clear directionality along the equilibrium paths, owing to a drift-induced selective force. However, this selective force is a weak one—its effect on fixation probabilities in the neutral case is of order and therefore can be dominated by direct selective forces such as viability, fecundity, or fertility differences between genotypes.

We also study the case where the sex-specific chromosome (the Y or the W) has accumulated deleterious recessive mutations, as it is expected to do over time in its nonrecombining region (Charlesworth 1978; Charlesworth and Charlesworth 2000). When the sex-specific chromosome is so degraded that homozygosity for it is lethal (as in mammals and birds), most heterogametic transitions are impossible (Bull 1983). However, in many other taxa, the sex-specific chromosome has accumulated only few recessive deleterious mutations—this state is common, for example, in amphibians, reptiles, and fish (Bachtrog *et al.* 2014). In such cases, homozygosity for the sex-specific chromosome is not expected to be lethal, and in fact can be associated with only a small fitness cost, if any at all. When taking into account these natural selective forces, we find, in model 1, that the drift-induced selective bias in favor of dominant sex-determining mutations is amplified, though the overall substitution rates in both directions are reduced. In model 2, the drift-induced bias in favor of dominant sex determining mutations is reduced when these selective forces are taken into account. Thus, in model 2, when selection is weak and the population is small, dominant sex-determining mutations substitute at a higher rate; when selection is strong and the population is large, recessive sex-determining mutations substitute at a higher rate.

The potential for drift-induced selection in transitions between sex-determining systems has previously been noted by Vuilleumier *et al.* (2007), who use simulations to investigate the stochastic dynamics of model 1, focusing predominantly on the effects of population structure. Our work overlaps with theirs in one particular case, that of a single deme, with no viability differences between the various genotypes. In that case, they too find non-neutral fixation probabilities for new dominant and recessive sex determining mutations, but find fixation probabilities substantially above the neutral expectation for both classes of mutation. In contrast, we find in this case that dominant sex determining mutations fix with probability above the neutral expectation, but recessive mutations fix with probability lower than the neutral expectation. The fixation probabilities that they report are also orders of magnitude different from those we report, especially for large populations [compare, for example, their figure 1(a) with our Figure 4A]. They also find mean conditional fixation times that are invariant, and, in some cases, even decreased, as population size increases; we find mean conditional fixation times that increase linearly with population size, consistent with drift-like dynamics. We have independently simulated their population model for the case that overlaps with ours, and obtain results consistent with ours, but different from those they report. The validity of our results is supported by mathematical analysis, which also allows us to explain them in analytical detail. In addition, we also consider transitions between sex-determining systems in model 2, where the sex chromosome locus is changed in the course of a transition; empirically, this scenario is possibly even more common than model 1 (Ezaz *et al.* 2006). Our results thus suggest that biases favoring dominant sex-determining mutations may be general to transitions between male and female heterogamety.

### Characterization of the equilibrium paths

#### Model 1:

Consider an initial male-heterogametic system, XX/XY. Suppose now that a mutation occurs on an X chromosome that renders the feminizing tendency of the resulting chromosome—call the new chromosome X′—dominant to the masculinizing tendency of the Y (so that X′Y individuals are female). Allowing all possible male–female matings between the genotypes results in a closed system of five sexual genotypes: females can be XX, X′X, or X′Y, while males can be XY or YY (Figure 1A, system arises in direction *a*). Clearly, this system could also arise in the reverse direction: starting from a female-heterogametic system, X′Y/YY, an X′ chromosome can mutate to an X chromosome with recessive feminizing tendency (Figure 1A, direction *b*). (It should of course be noted that the usual sex chromosome labels—X, Y, Z, and W—are arbitrary, so that we could just as validly label those in a female-heterogametic system X′ and Y. We also note that “sex chromosomes” are defined simply by the presence of a locus at which genes of major sex-determining effect segregate—it is possible, for example, that X and Y chromosomes are identical along their entire length, and recombine along their entire length, except at the major sex-determining locus.)

Assuming the genotypes all to have equal fitness (*i.e.*, that the system is neutral), enumerating them in the above order (XX, X′X, X′Y, XY, and YY), and letting be the population frequency of genotype *i*, Bull and Charnov (1977) showed that, for any value the population state(1)is an equilibrium in an infinite, randomly mating, population (Figure 1B). When all males are XY and all females XX, so that a system of male heterogamety operates; when males are YY and females X′Y, and the system is female heterogamety. For intermediate values all five genotypes are present at positive frequency.

A symmetric, though distinct, path exists, where, from an initial female-heterogametic system ZW/ZZ, a dominant masculinizing Z′ arises from a mutated Z, and, if it reaches high enough frequency, establishes a male-heterogametic WW/WZ′ system. The reverse transition along the same path involves fixation of the recessive Z chromosome from an initial WW/WZ′ system. Intermediate states along this path involve two female and three male genotypes.

These symmetric paths, one with three female and two male genotypes (as illustrated in Figure 1), and the other with two female and three male genotypes (described in the previous paragraph), are illustrated in a general format in Figure 3A. There, among the transitions involving fixation of dominant sex-determining mutations, those from male to female heterogamety are labeled *a* (as in Figure 1A), while those from female to male heterogamety are labeled *c*. Among the transitions involving fixation of recessive sex-determining mutations, those from female to male heterogamety are labeled *b* (as in Figure 1A), while those from male to female heterogamety are labeled *d*. We shall use this labeling throughout for model 1 transitions.

Notice that, if there are no demographic differences between males and females, then the dynamics in directions *c* and *d* are identical to those in directions *a* and *b*, respectively, up to a relabeling of males and females.

For the non-neutral case, the equilibrium path connecting systems of male and female heterogamety no longer exists (Bull and Charnov 1977). The resultant dynamics in this scenario depend on the selective forces acting on each of the genotypes. Transitions in either direction involve the production of individuals homozygous for a previously sex-specific chromosome (the Y for transitions in direction *a*, and the X′ and X for those in direction *b*). Sex-specific chromosomes are expected to accumulate recessive deleterious mutations in a heterogametic system (Charlesworth 1978; Charlesworth and Charlesworth 2000), causing individuals homozygous for them to be selected against weakly in young heterogametic systems (as in many reptiles, amphibians, and fishes), and to be inviable in old systems (as in mammals and birds). In the latter case, model 1 transitions are impossible. In the former case, we expect YY genotypes to be selected against in transitions in direction *a*, and X′X, and XX genotypes to be selected against in transitions in direction *b* (in direction *b*, the X chromosome is created simply by mutation at a major sex determining locus on the X′ chromosome; the X is therefore expected to carry the same deleterious mutations that the X′ does).

#### Model 2:

Begin with a male-heterogametic system XX,AA/XY,AA, where A is initially autosomal. Now, suppose that a mutation occurs on an A chromosome that confers on the resultant chromosome, A′, a feminizing tendency dominant to the masculinizing tendency of the Y (so that XY,AA′ and YY,AA′ individuals are female). All possible matings then yield a closed system of six sexual genotypes, females being XX,AA, XX,AA′, XY,AA′, or YY,AA′, and males being XY,AA or YY,AA (Figure 2, direction *a*). Again, this system could arise in the reverse direction as well, starting from the female-heterogametic system YY,AA′/YY,AA and introducing, as a mutated Y chromosome, a recessive feminizing X (Figure 2, direction *b*).

Scudo (1964, 1967) and Bull and Charnov (1977) showed that, when all genotypes are equally fit (*i.e.*, when the system is neutral), enumerating them in the above order (XX,AA, XX,AA′, XY,AA′, YY,AA′, XY,AA, YY,AA), and writing for the frequency of genotype *i*, a continuous path of equilibria,(2)connects male heterogamety at one end (XX,AA/XY,AA; ) with female heterogamety at the other (YY,AA′/YY,AA; ), with intermediate equilibria () involving all six genotypes.

Notice that, if the former system transitions to the latter (Figure 2, direction *a*), the Y chromosome becomes autosomal, and the previous autosome A becomes a sex chromosome. Notice too that this transition involves the production of YY individuals. On the other hand, the reverse transition (Figure 2, direction *b*) does not involve the production of individuals homozygous for the previously sex-specific chromosome (here, the A′ chromosome)—an important difference.

Again, a symmetric, though distinct, equilibrium path exists that connects a female-heterogametic system AA,ZW/AA,ZZ with a male-heterogametic system AA′,WW/AA,WW via an intermediate system with two female and four male genotypes. Transitions from female to male heterogamety along this path involve production of WW individuals, the W having been female-specific in the original system of female heterogamety. But transitions from male to female heterogamety along this path do not involve the production of A′A′ individuals, the A′ having been male-specific under male heterogamety. That heterogametic transitions along standard equilibrium paths are possible without the production of WW or YY individuals (using the standard labeling) is an important fact often forgotten in the literature on evolutionary transitions between sex-determining mechanisms.

These symmetric equilibrium paths, one with four female and two male genotypes (as in Figure 2), and the other with two female and four male genotypes (as described in the previous paragraph), are illustrated in a general way in Figure 3B. There, among model 2 transitions involving fixation of dominant sex determining mutations, those from male to female heterogamety are labeled *a*, while those from female to male heterogamety are labeled *c*. Among the reverse transitions involving fixation of recessive sex determining mutations, those from female to male heterogamety are labeled *b*, while those from male to female heterogamety are labeled *d*. We shall use this labeling throughout for model 2 transitions. Transitions *b* and *d* do not involve the production of individuals homozygous for a previously sex-specific chromosome, while transitions *a* and *c* do.

Again, notice that, if there are no demographic differences between males and females, then the dynamics in directions *c* and *d* are, respectively, identical to those in directions *a* and *b* up to a relabeling of males and females.

## Methods

The equilibrium paths described in the previous section arise under deterministic evolutionary dynamics. Our aim is to study stochastic evolution along and around these paths. To do so, we employ Monte Carlo simulations to estimate the substitution rates of the various mutations along the paths, and approximation techniques to analytically investigate the results of these simulations.

For computational efficiency, our simulations are of a nonoverlapping generations Wright-Fisher process (Fisher 1930; Wright 1931; Hartl and Clark 2007). For mathematical tractability, on the other hand, our analytical treatment considers models of the overlapping generations Moran type (Moran 1958); mapping between these processes in the diffusion limit simply requires rescaling the rate of genetic drift (Ewens 2004). The agreement of the results under the two processes will demonstrate their robustness to the consideration of overlapping or nonoverlapping generations.

### Monte Carlo simulations

For both model 1 and model 2, we simulate a population of constant size *N*, which comprises males and females, and evolves according to a sexual Wright-Fisher process (Fisher 1930; Wright 1931; Hartl and Clark 2007). Each generation, males and females form mating pairs, *N* in total. An individual can be in more than one pair, and the probability that an individual is in a given pair is proportional to that individual’s fitness relative to other members of its sex (and is independent across pairs). Each mating pair produces a single offspring, whose sexual genotype (and therefore sex) is determined by randomly choosing a gamete from each of its parents. The sex chromosomes of heterogametic individuals are assumed to segregate in a Mendelian fashion. After offspring production, all individuals in the parental generation die. (This model is equivalent to one in which each male contributes a large number of sperm, proportional to his fitness, to a common sperm pool, each female contributes a large number of eggs, proportional to her fitness, to a common egg pool, and each of *N* offspring is then formed by drawing a random sperm and a random egg from the respective pools.)

As a baseline for both models, we consider the case where each genotype is equally fit (the “neutral” case), *i.e.*, where each individual within a sex is equally likely to be chosen to be in a given mating pair. Thereafter, we focus on cases where individuals that are homozygous for a previously sex-specific chromosome suffer a selective disadvantage: each such individual is a fraction as likely as other members of its sex to be chosen to be in a given mating pair.

In all cases, we assume no population structure (mating is random), and no demographic differences between males and females in our simulations. Relaxing these assumptions is an important direction of future work, and would be aided by parallel theoretical developments in the general study of drift-induced selection. In Supplemental Material, Section S1.3 in File S1, we give some numerical picture of how demographic differences between males and female (*viz*. greater variance in male reproductive success) affect our results (also see *Discussion* section).

### Diffusion approximations

In the Moran formulation, we likewise consider a discrete population of *N* individuals. Males and females of each genotype encounter each other with a probability per unit time proportional to their frequency in the population. On encountering each other, a pair produces a single offspring, which inherits its sexual genotype from its parents in a Mendelian fashion. In order to keep the population size constant at *N*, another individual is chosen to die. In the neutral case, each individual in the population is equally likely to be chosen to die. In the case with selection, the probability that an individual is chosen to die is weighted by a normalized death probability, the inverse of that genotype’s fitness. This implementation of selection will give similar dynamics to the Wright-Fisher model to leading order in the selection strength (see Section S2.1 in File S1).

In order to simplify the probabilistic model, we make use of a diffusion approximation (Crow and Kimura 1970). Denoting the frequency of each genotype by and assuming *N* to be large but finite, the evolution of the pdf for the system is approximately governed by(3)(Gardiner 2009; McKane *et al.* 2014), where the forms of the vector and the matrix can be directly calculated for both models 1 and 2 from their respective probability transition rates (see Section S2.1 in File S1). The term primarily controls the time-evolution of the mean of and can thus be interpreted as a deterministic selective term. Indeed, in the deterministic limit (), the dynamics of the system are described by the ordinary differential equations Meanwhile, the term controls the diffusion of thus capturing the effect of genetic drift.

In comparing this diffusion equation with one for the analogous Wright-Fisher process, we must take account of two scalings. First, genetic drift is larger by a factor of two in the Moran formulation than in the Wright-Fisher formulation (Felsenstein 1971; Ewens 2004). This leads to an additional prefactor of appearing before the in the Wright-Fisher formulation of Equation (3). Second, reproduction occurs times faster in the Wright-Fisher model. This is because individual reproduction events in the sexual Moran model occur at an average rate proportional to the product of the frequency of males and females in the population [which is at equilibrium sex ratios], whereas, in the Wright-Fisher model, *N* reproduction events occur every time-step. Since we have already taken into account a factor *N* in timescale in obtaining Equation (3) (), the Wright-Fisher formulation of Equation (3) contains an additional factor of 4 preceding all terms on the right-hand side of Equation (3).

We wish to determine the probabilities of transitions from male to female heterogamety, and vice versa. The calculation of these quantities is not straightforward; the diffusion Equation (3) governing the dynamics is nonlinear in either four or five variables (for models 1 and 2, respectively). However, analytical progress can be made by exploiting a separation of timescales that is present in both models. Such methods have long been successfully employed in population genetics to solve a wide range of problems (Ethier and Nagylaki 1980; Stephan *et al.* 1999; Newberry *et al.* 2016). The key to progress is in noting that, in the deterministic neutral systems, a trajectory starting from any initial point will quickly collapse to a point on the equilibrium path of the system, Equation (1) for model 1 (Figure 1B) and Equation (2) for model 2, where it will then stay. In the current notation, this line is the set of solutions to the equation When genetic drift and selection are taken into account, the system will no longer reach and subsequently remain at a position on this line. However, if selection is weak and *N* is large [such that the rate of genetic drift, is small], then the system will quickly collapse to a subspace in the vicinity of this line; it will then slowly move along this “slow subspace” until the system fixes in a state of either male or female heterogamety. We exploit this separation of timescales by removing the fast transient dynamics to obtain an approximation for the system dynamics in the slow subspace. Since this approximate description in the slow subspace is in terms of a single variable, *q* [see Equation (1) and Equation (2)], fixation probabilities are then straightforward to calculate.

There is no unique way to mathematically implement the approach outlined conceptually above. While methods similar to the projection operator formalism described in Gardiner (2009) have historically found most favor in the population genetics literature (*e.g.*, Ethier and Nagylaki 1980; Stephan *et al.* 1999), in this paper we implement the approach described in Parsons and Rogers (2015), which is more intuitive in the present case. The full calculation is provided in Section S2.2 in File S1. Here, we simply state the key results. On removing the fast variables, the dynamics of Equation (3) can be approximated by(4)Here, is the probability density function for *q* along the slow subspace. In a similar fashion to Equation (3), the terms and control the time-evolution of the mean of and can thus be interpreted as selective terms for *q* along the slow subspace. Likewise, the term controls the diffusion of and can thus be interpreted as capturing the effect of genetic drift along the slow subspace.

The equation for the slow subspace itself can be approximated by the equation for the equilibrium line, since they lie close to each other in the weak selection limit, and coincide when selection is absent. The equations for and can be determined from and along with the equation for the slow subspace. The terms and are simply the components of and along the slow subspace (*i.e.*, respectively, the components of deterministic selection and genetic drift in the subspace). More interesting is the term which is a selective term induced by genetic drift. Its origin can be interpreted in various ways. First, it can be graphically understood as resulting from a bias in how fluctuations taking the system off the slow subspace return to the slow subspace (*i.e.*, fluctuations, on average, do not return to the point from which they originated on the slow subspace) (Parsons and Rogers 2015; Constable *et al.* 2016). Second, it can be mathematically understood as the result of making a nonlinear change of stochastic variables (Risken 1989) into the system’s slow variable. Finally, it can be understood as the result of a selective pressure favoring genotypes with a lower variance in their reproductive output—Gillespie’s Criterion (Gillespie 1974; Parsons *et al.* 2010; Hansen 2017). In Equation (4), differences in the variance of the reproductive output of the genotypes emerge naturally from the confinement of the system to the slow subspace.

### Data availability

The authors state that all data necessary for confirming the conclusions presented in the article are represented fully within the article.

## Results

### Model 1: transitions using the same chromosomes

In this section, we study transitions between male and female heterogamety where the sex chromosome locus is the same under both systems (Figure 1 and Figure 3A).

#### Monte Carlo simulations:

We begin with a male-heterogametic system, XX/XY, in a population of constant size *N*, initially with females (all XX) and males (all XY). We consider a mutation to one of the X chromosomes, rendering it an X′ chromosome, X′Y and X′X bearers of which are female. YY individuals are male. If the X′ chromosome “fixes” in the population, a female-heterogametic system, X′Y/YY, is established (direction *a* in Figure 1A and Figure 3A).

Influenced by whether the original X′ mutation occurs in oogenesis or spermatogenesis, the X′ could initially find itself in an X′X or an X′Y female. We consider both cases in our simulations, and in both begin with a population that is females (one of which carries an X′ chromosome) and males.

Initially assuming all five sexual genotypes to be equally fit, what is a reasonable null expectation for the fixation probability of this X′ chromosome? The fixation probability of a mutation on the X chromosome of no effect is simply the inverse of the initial census count of the X chromosome, *i.e.*, This we take to be the neutral expectation for the fixation probability of the X′ chromosome, so that

Instead, we find in our simulations that the fixation probability of the X′ chromosome is higher than this neutral expectation, with (Figure 4A). This is irrespective of the background on which the initial X′ chromosome finds itself (Figure S1A in File S1). The average conditional fixation time of the X′ chromosome is close to for each *N* considered in our simulations, and, again, this is irrespective of the mutation’s initial background (Figure S1B in File S1). A fixation time that scales linearly with *N* is strongly suggestive of drift-like evolutionary dynamics.

We now consider evolution in the other direction along the path (direction *b* in Figure 1A and Figure 3A). We begin with an established X′Y/YY female-heterogametic system, and consider a mutation to one of the X′ chromosome that renders it an X. If this X chromosome subsequently fixes in the population of X and X′ chromosomes, an XX/XY male-heterogametic system would be established. Here, the X mutation must occur in oogenesis, since it derives from an X′ chromosome that must have been borne by a female. Therefore, the first individual to carry the new X chromosome must be an XY male. For consistency, we begin our simulations with females and males (one of which carries the X chromosome).

Similar to before, our neutral expectation for the fixation probability of the X is just the inverse of the number of X′ chromosomes initially in the population, *i.e.*, so that Instead, in our simulations we find the fixation probability of the X chromosome to be much lower than this neutral expectation, with (Figure 4B). The average conditional fixation time of the X chromosome is, like that of the X′, close to for each *N* considered (Figure S2B in File S1).

Thus, the fixation probability of the X′, starting from an initial XX/XY system, is almost twice as great as the fixation probability of the X, starting from an initial X′Y/YY system. However, to properly determine whether evolution along the equilibrium path is biased in favor of the X′ chromosome, we must consider the *substitution rates* of the two chromosomes. In doing so, we assume symmetric mutation rates in the two directions, *i.e.*, that the probability of an X chromosome mutating to an X′ is the same as the probability of an X′ chromosome mutating to an X. We also assume this mutation rate to be sufficiently small that at most one mutation segregates in the population at any given time [a common assumption (McCandlish and Stoltzfus 2014)].

Suppose this mutation rate to be *u* per chromosome per generation. An XX/XY population with even sex ratio produces X′ mutations at a rate of per generation, while an X′Y/YY population produces X mutations at a third of this rate, per generation. Under our naive neutral expectations, the substitution rates of the X and X′ chromosomes should both be *u* (). Instead, the substitution rate from an XX/XY to an X′Y/YY system (in direction *a* of Figure 1A and Figure 3A) is about while the substitution rate from an X′Y/YY to an XX/XY system (direction *b* of Figure 1A and Figure 3A) is about That is, the substitution rate of the X′ is more than *five* *times higher* than that of the X.

Since, in the population model we have simulated, there are no demographic differences between males and females, the transitions symmetric to those above (*i.e.*, in directions *c* and *d* in Figure 3) occur with fixation probabilities and substitution rates equivalent to those we have estimated above (*i.e.*, with reference to Figure 3, the fixation probabilities in directions *a* and *b* match those in directions *c* and *d*, respectively).

Therefore, the most likely trajectory that the neutral dynamical system will follow in the long term is the recurrent invasion and fixation of successive *dominant* sex-determining mutations, flipping the system repeatedly between male and female heterogamety. This corresponds to a bias in favor of rightward transitions in Figure 3.

We now consider the possibility that some genotypes are fitter than others. In particular, we study the case where individuals homozygous for a previously sex-specific chromosome (including mutated versions of it) are of lower fitness.

The transition from a male-heterogametic XX/XY system to a female-heterogametic X′Y/YY system, via fixation of the dominant X′ chromosome and displacement of the X, involves the production of YY males. Assume that YY males have relative fitness while all other sexual genotypes are of equal fitness 1. Figure 4C gives the fixation probability of the X′ chromosome for various values of *s* and *N*. As expected, the fixation probability decreases as selection against the YY genotype increases. Nonetheless, it is clear that the YY genotype can suffer appreciable fitness reductions with the X′ chromosome still fixing with probability higher than the neutral expectation. This is especially true in small populations, in which selection acts less efficiently (Lanfear *et al.* 2014); we would expect this to carry over to structured populations of larger size, *e.g.*, those that are subdivided into many small demes in which drift is an important force (Laporte and Charlesworth 2002; Whitlock 2003).

The reverse transition, from the X′Y/YY to the XX/XY system, involves the production of X′X and XX females. The substitution rate of the X, in direction *b* of Figure 1A, was very low even when no genotypes were of reduced fitness; selection against the XX and X′X genotypes (recall that the X′ chromosome is just an X mutated at the sex-determining locus) severely exacerbates this disadvantage (Figure 4D and Figure S4 in File S1). Again, given the symmetry between transitions in directions *a* and *c* in Figure 3, and between transitions in directions *b* and *d*, the above results imply a bias toward *c* transitions over *d* transitions, with this bias exacerbated by selection. That is, selection generally exacerbates the bias in favor of rightward (dominant) transitions in Figure 3.

### Analytical results

To better understand these simulation results, we now study the system analytically in the diffusion limit. We first consider the neutral case. Applying the fast-variable elimination described in *Methods* and detailed in full in the Section S2 in File S1, we find, for the model 1 system described in Equation (1), that the system dynamics can be approximated by Equation (4), with and(5)(calculations in Section S2.2.1 in File S1). As we expect in this neutral scenario, there are no deterministic contributions to selection along the equilibrium path []. However, there is a drift-induced selection term []. The strength of the drift-induced selection is of order generated as it is by demographic fluctuations. Recalling that corresponds to male heterogamety [XX/XY; in Equation (1)], while corresponds to female heterogamety [X′Y/YY; in Equation (1)], we find that the drift-induced selection selects for the fixation of the dominant sex determining mutation *at every point* on the equilibrium path [ for all ].

How does this drift-induced selection emerge? Essentially, in this “neutral” stochastic system, demographic fluctuations continually perturb the system away from the equilibrium path (to which the deterministic system is constrained). There is then a selective pressure for the system to return to the equilibrium path. However, the nonlinear trajectories along which these fluctuations return, combined with the curvature of the equilibrium path, give rise to a bias in the average position to which fluctuations return. In other words, fluctuations arising at a point *q* return on average to a point on the equilibrium path.

Mathematically quantifying this bias requires taking account of the probability distribution of fluctuations in each genotype, the form of trajectories back to the equilibrium path, and the curvature of the equilibrium path itself, each of which varies as *q* is varied. However, further intuition can be gained by decomposing into two components: The first term, quantifies contributions to arising from the nonlinearity of trajectories that take the system back to the equilibrium path. The second term, quantifies contributions to arising from the curvature of the equilibrium path itself. Plotting these terms together in Figure 5A, we see that it is the curvature of the equilibrium path that contributes most to the observed drift-induced selection.

Standard methods can be used to numerically calculate the fixation probability of either male or female heterogamety for any initial condition, on the slow subspace [see Gardiner (2009); Risken (1989); and Section S2.3 in File S1]. Our final task is to calculate how the initial conditions described in the previous section, *i.e.*, those of single mutants invading a resident population, map onto initial conditions on the equilibrium path. That is, for a given we wish to determine This calculation is given in Section S2.3.1 in File S1; with it, the neutral fixation probabilities can be calculated.

The results for the fixation probabilities of the mutants are given in Figure 4A, in which we see excellent agreement between theory and simulations. In particular, we find that the fixation probability of the X′ chromosome under the diffusion approximation is (Figure 4A). This is irrespective of the background on which the initial X′ chromosome finds itself, as both of these initial conditions lead to the same initial condition on the equilibrium path (Section S2.3.1 in File S1). We may also compute the mean conditional fixation time of the X, which agrees well with our simulation estimates when the differences in variance between the Wright-Fisher and Moran processes are taken into account (Figure S1B in File S1). Meanwhile, the fixation probability of an X chromosome is again in close agreement with our simulation results (Figure 4B), as is the computed mean conditional fixation time (Figure S2B in File S1).

We now consider the possibility that some genotypes are fitter than others. For transition direction *a*, assume that YY males (with frequency ) have relative death rate while all other sexual genotypes have death rate 1. When *s* is small, we can still utilize fast-variable elimination to arrive at the approximate description of the system dynamics given by Equation (4). The functions and retain the forms given in Equation (5), but now where(6)and the superscript “” denotes that this is evaluated for transitions in direction *a* (Figure 5B; calculations in Section S2.2.1 in File S1). This term is of order *s* as this is the deterministic contribution to the dynamics along the slow subspace.

For the reverse transition *b*, assume that XX and X′X females (with frequencies and ) have relative death rates and respectively, while all other sexual genotypes have death rate 1. Again utilizing fast variable elimination, we find where(7)and the superscript “” denotes that this is evaluated for transitions in direction *b* (Figure 5B; calculations in Section S2.2.1 in File S1).

We are now in a position to calculate the respective fixation probabilities of mutations arising in directions *a* and *b*. We find that, while the fixation probability of X′ in direction *a* is of course lower than in the neutral case (since YY is selected against), the relative reduction of the fixation probability of X in direction *b* is even higher (Figure 4D).

### Model 2: transitions that change the sex chromosome pair

In this section, we study transitions between male and female heterogamety where, in the course of the transitions, a pair of chromosomes that are initially autosomal are co-opted as new sex chromosomes, while one of the old sex chromosomes becomes autosomal (Figure 2 and Figure 3B).

#### Monte Carlo simulations:

Beginning with an XX,AA/XY,AA male-heterogametic system (where X and Y are sex chromosomes and A is an autosome), assume that a mutation appears on an A chromosome, rendering it an A′ such that XX,AA, XX,AA′, XY,AA′, and YY,AA′ individuals are female, while XY,AA and YY,AA individuals are male (Figure 2). If the A′ chromosome reaches sufficiently high frequency in the population, the X chromosome is eliminated, and a YY,AA′/YY,AA female-heterogametic system establishes (direction *a* in Figure 2 and Figure 3B).

We initially assume all six sexual genotypes to be equally fit. Unlike for the case of transitions involving the same chromosome pair, we do not propose a null *expectation* for the fixation probability of the A′ mutation. This is because, if it fixes, it displaces the unlinked X chromosome from the population: this is not the “population” in which the A′ arises, being a mutated A chromosome. (In contrast, in model 1, the X′ chromosome for example is a mutated X chromosome, and, if it fixes, it displaces the X chromosome.) Therefore, we shall focus predominantly on comparing the substitution rates of the A′ and X chromosomes (*i.e.*, the substitution rates, respectively, in directions *a* and *b* of Figure 2). We may, however, take as a *reference* neutral fixation probability for both mutations, that of a mutation of no effect occurring on an autosome:

In our simulations, we find the fixation probability of the A′ chromosome to be (Figure 6A), substantially higher than our reference value of This value is insensitive to the genetic background of the initial mutation (Figure S5A in File S1). The average conditional fixation time of the A′ in our simulations is ∼ for each *N* considered, again regardless of the initial background of the mutation (Figure S5B in File S1), and again suggestive of drift-like dynamics.

Turning our attention to transitions in the other direction along the equilibrium path (direction *b* in Figure 2 and Figure 3B), we begin with a female-heterogametic YY,AA′/YY,AA system, and assume that a mutation on a Y chromosome occurs, rendering the chromosome a recessive feminizing X. If this mutation reaches sufficient frequency, the male-heterogametic system XX,AA/XY,AA establishes.

We estimate in our simulations that the fixation probability of the X is (Figure 6A), which is fractionally higher than the reference value of The background on which the mutation arises has no effect on its fixation probability (Figure S6A in File S1). The average conditional fixation time of the X chromosome is also close to for all values of *N* considered (Figure S6B in File S1).

To see if there is a directional bias in one direction or the other along the equilibrium path, we assume that these mutations occur at the same rate, *u* per chromosome per generation, and calculate the substitution rates of the two transitions.

The male-heterogametic system XX,AA/XY,AA generates A′ mutations at rate so that the substitution rate from an XX,AA/XY,AA system to a YY,AA′/YY,AA system (in direction *a* of Figure 2 and Figure 3B) is about Similarly, the female-heterogametic system YY,AA′/YY,AA generates X mutations at rate and so the substitution rate from a YY,AA′/YY,AA system to an XX,AA/XY,AA system (direction *b* of Figure 2 and Figure 3B) is about or about half that of the reverse transition.

Again, with no demographic differences between males and females, the transitions symmetric to those above (*i.e.*, in directions *c* and *d* in Figure 3B) occur with fixation probabilities and substitution rates equivalent to those we have estimated above (*i.e.*, the fixation probabilities in directions *a* and *b* match those in directions *c* and *d*, respectively).

We now consider the role of selective differences between the genotypes. This is an especially important question here because, unlike in model 1 heterogametic transitions, model 2 transitions are possible without the production of individuals homozygous for a previously sex-specific chromosome. In particular, the transitions in directions *b* and *d* of Figure 3B change the heterogametic system, but do not involve the production of individuals homozygous for the previously sex-specific chromosome. In contrast, the reverse transitions (in directions *a* and *c*) do involve the production of individuals homozygous for previously sex-specific chromosomes. Since we have found these latter transitions to have higher substitution rates than the reverse transitions in the neutral case, we should expect selection to reduce, and, when strong enough, to overturn, this bias.

We focus on the transition from the male-heterogametic XX,AA/XY,AA system to the female-heterogametic YY,AA′/YY,AA system (direction *a* in Figure 2 and Figure 3B), involving the substitution of a dominant female-determining A′ chromosome as a mutated A. The Y chromosome is sex-specific in the original XX,AA/XY,AA system, and we assume that it has accumulated deleterious recessive mutations such that the two YY genotypes are of fitness relative to all other genotypes’ fitness of 1. Figure 6B gives the fixation probability of the A′ chromosome for various values of *s* and *N*. Naturally, the fixation probability decreases as selection against the YY genotypes increases, and this effect is stronger in larger populations [in which selection acts more efficiently (Lanfear *et al.* 2014)]. Indeed, in large populations, even very small degrees of selection against the YY genotypes are enough to overturn the substitution rate bias in favor of dominant sex-determining mutations.

#### Analytical results:

We begin by considering the dynamics of the neutral model. Once again, fast-variable elimination can be used to calculate the effective dynamics of the system along the equilibrium path [see Equation (4)]. For model 2 (as with model 1), (that is, there is no deterministic contribution to the dynamics along the equilibrium path), but there is a drift-induced selection term, which now takes the form(8)which, along with the expression for diffusion along the equilibrium path,(9)approximates the stochastic dynamics (calculations in Section S2.2.2 in File S1). Recalling that corresponds to male heterogamety [XX,AA/XY,AA; in Equation (2)], while corresponds to female heterogamety [YY,AA′/YY,AA; in Equation (2)], we find that the drift-induced selection selects for the fixation of the dominant sex determining mutation *at every point* on the equilibrium path [ for all ].

As in model 1, we find that demographic fluctuations away from the equilibrium path return to the equilibrium path on average with a bias, described by the drift-induced selection term Once again, can be split into two components, and that respectively capture the contribution to arising from the nonlinearity of trajectories taking the system back to the equilibrium path, and the curvature of the equilibrium path itself. These terms are plotted in Figure 7A, in which we see that, as with model 1, it is the curvature of the equilibrium path that contributes most to drift-induced selection in model 2.

Using this single-variable description of the dynamics, the fixation probability for any initial condition on the equilibrium path can be calculated. In order to determine the fixation probability of mutants in the system starting from a single copy, we need to calculate the mapping, for each initial mutation scenario, from to This calculation is given in Section S2.3.2 in File S1. We can now evaluate our numerical expression for the neutral fixation probability at these initial conditions. Recall that simulations of the Wright-Fisher process showed the fixation probability of the A′ chromosome to be higher than our reference value of Our analytical prediction slightly underestimates this fixation probability ( Figure 6A). The background on which the mutation arises has no effect on its fixation probability, because the two scenarios initially have an identical component along the slow manifold (Section S2.3.2 in File S1). The computed mean conditional fixation time of the A′ agrees well with our simulation estimates (Figure S5B in File S1).

We next consider the fixation probability of an X mutation occurring on an initially autosomal Y. Our Wright-Fisher simulations returned an estimated fixation probability of Once again, there is a small discrepancy with our analytical results, which overestimate this value ( Figure 6A). The fixation probability is again the same irrespective of the background on which the X mutation occurs, because the initial conditions on the equilibrium path are the same (Section S2.3.2 in File S1). Again, we may compute the mean conditional fixation time of the X, and find good agreement with our simulation estimates when the variance difference between the Wright-Fisher and Moran processes is taken into account (Figure S6B in File S1).

We now consider the role of selective differences between the genotypes. As we have noted, comparing transitions in directions *a* and *b*, only those in direction *a* involve the production of individuals homozygous for a previously sex-specific chromosome, and so we focus here on transitions in direction *a*. We assume that the two YY genotypes (with frequencies and ) have elevated death rates and relative to all other genotypes’ death rates of 1. The term in Equation (4) now becomes where(10) for all and so this term acts in the opposite direction to (Figure 7B; calculations in Section S2.2.2 in File S1). There is thus an antagonism between the deterministic contribution to selection (favoring transitions in direction *b*) and the drift-induced selection (favoring transitions in direction *a*). Which of these dominates depends on the strength of selection and the population size: increases with *s*, while decreases with *N*.

In contrast, since transitions in direction *b* do not produce individuals homozygous for a previously sex-specific chromosome, none of the genotypes is selected against. Therefore, in direction *b* is

Whereas in model 1, selection exacerbated the directionality of switching between sex determining systems (while decreasing the overall switching rate), in model 2 we see that deterministic selection and drift-induced selection work in opposite directions. Thus, for small populations with weak selection, transitions in directions *a* and *c* in Figure 3B occur more often than transitions in directions *b* and *d*, while for large populations with strong selection, transitions in *b* and *d* occur more often than transitions in directions *a* and *c*.

## Discussion

We have studied stochastic evolution along two “neutral” equilibrium paths connecting male and female heterogamety. We have shown that, even when all genotypes are equally fit, evolution along these paths is not neutral. Instead, it shows significant substitution rate biases in particular directions, specifically in favor of dominant sex-determining mutations. We have demonstrated these biases to be the result of drift-induced selection: random perturbations off the equilibrium path—inevitable in finite populations—return to the equilibrium path with an average directional bias in favor of the dominant segregating sex chromosome. The substitution rates of dominant sex determining mutations that switch the system of heterogamety are, in both of the cases we have studied, higher than those of truly neutral mutations occurring on the same chromosomes.

Evolutionary transitions between male and female heterogamety have been common in the evolutionary history of animals (Hillis and Green 1990; Ezaz *et al.* 2006, 2009; Mank *et al.* 2006; Mank and Avise 2009; Pokorna and Kratochvíl 2009; Kaiser and Bachtrog 2010; Bachtrog *et al.* 2014), despite the fact that, in any given heterogametic system, the sex-specific chromosome should degrade over time by the operation of Muller’s ratchet (Charlesworth 1978; Charlesworth and Charlesworth 2000). To counter this “evolutionary trap” (Bull and Charnov 1985; Rice 1998), mechanisms based on direct selective forces have been invoked to explain the frequency of heterogametic transitions [*e.g.*, the operation of sex-specific selection (van Doorn and Kirkpatrick 2010)]. Our demonstration of drift-induced selection in transitions between male and female heterogamety suggests that mechanisms based on direct selective differences between the sex chromosome genotypes may be unnecessary to explain empirical transitions. We have shown that small fitness reductions to individuals homozygous for previously sex-specific chromosomes are not enough to overturn the biases caused by drift-induced selection (Figure 4C and Figure 6B). Large fitness reductions to such individuals (and in the limit, their inviability) render most transitions impossible—the exception being the recessive transitions *b* and *d* in model 2.

The influence of drift-induced selection in heterogametic transitions is therefore best understood in terms of evolutionary timescale. From an initial heterogametic system with homomorphic sex chromosomes, the sex-specific chromosome gradually accumulates recessive deleterious mutations that would reduce the fitness of individuals homozygous for this chromosome. At some threshold period of time of accumulation of these mutations, the fitness reduction of the homozygote reduces the fixation probability of a dominant sex reversal mutation exactly enough to cancel this mutation’s drift-induced selective advantage (Figure 4C and Figure 6B). Prior to this time threshold, a transition via a dominant sex reversal mutation is likely—it becomes less likely as the threshold is neared. However, progress toward the threshold is reset every time a transition occurs, because each transition creates a new sex-specific chromosome from a chromosome that was previously not sex-specific.

A notable prediction of our findings is that heterogametic transitions should typically involve dominant sex-determining mutations. That is, we should expect transitions usually to occur in directions *a* and *c* in Figure 3, A and B, and not in directions *b* and *d*. This general prediction is not unique to our theory. For example, under the theory that heterogametic transitions are driven by linkage between novel sex-determining genes, and genes with sex-specific fitness effects (van Doorn and Kirkpatrick 2010), dominance of the sex-determining mutation is either required for a transition, or substantially increases the parameter range over which a transition may occur. However, our analysis does reveal that this prediction can be recapitulated with a minimal number of biological assumptions.

Evidence in favor of this prediction comes from intermediate “multi-factorial” systems (Bachtrog *et al.* 2011). In the platyfish *X. maculatus*, females are XX, WX, or WY, while males are XY or YY (Kallman 1965, 1968). In principle, such a system could have arisen either in directions *a* or *b* of Figure 3A, depending on the ancestral heterogametic system. Mapping known systems of heterogamety in the genus *Xiphophorus* (Tree of Sex Consortium 2014) onto a phylogeny of the clade (Cui *et al.* 2013) suggests this ancestral system to be male heterogamety, in which case the intermediate system of *X. maculatus* has arisen in direction *a* of Figure 3A, via a dominant sex-reversal mutation, consistent with our prediction. The western clawed frog, *X. tropicalis*, also has an intermediate system: females are ZW or WW, and males are ZZ, ZY, or WY (Roco *et al.* 2015). This system could have arisen in direction *c* or *d* in Figure 3A, depending on the ancestral system. The mechanism of sex determination has yet to be determined for most members of the genus *Xenopus* (Tree of Sex Consortium 2014; Roco *et al.* 2015), but the well-studied *X. laevis* is female-heterogametic (Chang and Witschi 1956), as is *X. borealis* (Furman and Evans 2016). If female heterogamety is ancestral to the intermediate system of *X. tropicalis*, then this intermediate system would have arisen in direction *c* of Figure 3A, again consistent with our prediction. Though it is possible that balancing selection operates to stabilize these observed instances of multi-factorial systems [*e.g.*, Orzack *et al.* (1980)]—with the important suggestion that, because observed instances are rare, most intermediate multi-factorial systems are transitional—the drift-induced selection that we have discovered operating at all points on the slow subspace near the line of equilibria will, even in this case, act so as to make invasion of dominant sex-determining mutations more likely.

The prediction that heterogametic transitions should usually involve dominant sex-determining mutations can also be tested by reciprocal crosses of species or populations on either side of a recent heterogametic transition, provided the ancestral system is known. In the frog *R. rugosa*, populations in northern Japan are female-heterogametic, while those in southern Japan are male-heterogametic (Nishioka *et al.* 1993; Miura *et al.* 1998). The sex chromosomes in these populations are all homologous (Uno *et al.* 2008), and so a model 1 transition appears to have occurred. Because the ancestral system is male heterogamety (Ogata *et al.* 2003), the candidate directions in Figure 3A are *a* and *d*. These two directions can be distinguished by crossing a homogametic male (from the north) with a homogametic female (from the south). If the transition occurred in direction *a*, all the offspring from this cross should be male, but if it occurred instead in direction *d*, all the offspring should be female. This test has been carried out using homogametic males from Hirosaki (in the north) and homogametic females from Kumano (in the south), the reciprocal cross of which yielded almost all male offspring (Nishioka *et al.* 1993). Again, this is consistent with our prediction. [Crossing heterogametic males and females yielded a sex ratio of (Nishioka *et al.* 1993), consistent with the model 1 transitions we have studied, though not informative of which direction the transition was in.]

The bias we have found in favor of substitution of dominant sex-determining mutations also carries predictions for how the genetic pathways underlying sex determination should look. We have referred throughout to “sex chromosomes,” but in reality we are talking about genes of major sex-determining effect, the presence or absence of which acts as a switch to direct development down separate molecular pathways, or sex-determining “cascades,” which then produce males and females. The downstream components of these cascades tend to be widely conserved (Beukeboom and Perrin 2014), but there is significant lability in the upstream components—through the addition of new sex-determining genes to the top of the cascade (Wilkins 1995), and the shuffling of genes already in a cascade (Schartl 2004)—suggesting that these cascades evolve “from the bottom up.” This is consistent with our findings: our model predicts successive transitions involving dominant sex-determining mutations, with comparatively few reversals involving fixation of recessive sex-determining mutations, and so we expect either the expansion of sex-determining cascades, or their shuffling, but seldom their contraction.

The drift-induced selective force that we have identified is a weak one: when all genotypes are equally fit, it shifts fixation probabilities away from the values expected under neutrality by amounts of order This raises two questions. First, are neutral transitions, even in the direction of the drift-induced bias, empirically relevant, given that they involve fixation probabilities of order ? Second, would direct selective forces, such as viability differences among genotypes, not overwhelm the drift-induced bias?

On the first question, as with the study of neutral substitutions elsewhere in the genome, this depends on how often the relevant mutations are produced. The extended sex-determining cascades discussed above present a large mutational target: mutations of major sex-determining effect can occur at many points along them. These mutations can also occur in many ways, in addition to standard sequence-mutation events: (i) translocation of a gene in a sex-determining pathway, and a resulting shift in expression or function (Traut and Willhoeft 1990; Charlesworth *et al.* 2005); (ii) duplication and subsequent neo-functionalization of such a gene (Schartl 2004; Bewick *et al.* 2011); and (iii) mutation of a major sex determining gene’s regulatory elements, such as transcription factors (Beukeboom and Perrin 2014; Takehana *et al.* 2014). The frequency of these mutations is evidenced by the heterogeneity of major sex-determining genes observed between and within clades (Beukeboom and Perrin 2014).

On the second question, it is true that the drift-induced bias we have found becomes very weak as population size increases, while direct selective forces do not. When such selective forces operate—for example, when there is selection against individuals homozygous for a previously sex-specific chromosome—we have argued that drift-induced biases are most likely to be relevant in populations with small effective sizes. One way in which effective population size can be reduced is through demographic differences between the sexes, for example if males exhibit higher variance in reproductive success. Notice that demographic differences between the sexes also eliminate the symmetry between dominant transitions from male to female heterogamety and from female to male heterogamety, and the same for recessive transitions, so that they also allow the possibility that male or female heterogamety may be favored. In Section S1.3 in File S1, we have introduced greater variance in male reproductive success for the case of model 1 transitions by altering our baseline population model so that, each generation, a reduced subset of males is randomly chosen to be candidate mates, with the other males denied the possibility of mating. In the neutral case, we find that dominant transitions from male to female heterogamety exhibit a substantially higher substitution rate than in the case with no demographic differences between the sexes (Figure S8A in File S1), while dominant transitions from female to male heterogamety have a slightly reduced substitution rate (Figure S8B in File S1). The substitution rate of recessive transitions from female to male heterogamety is increased marginally relative to the case of no sex differences (Figure S8D in File S1), while recessive transitions from male to female heterogamety have a significantly reduced substitution rate (Figure S8C in File S1). Therefore, the general effect of greater variance in male reproductive success in the neutral case is to exacerbate the bias in favor of dominant sex-determining mutations, and to bias evolution toward female heterogamety. When selection against individuals homozygous for a previously sex-specific chromosome is taken into account, we expect the results to differ from those under no sex differences for two reasons: because drift-induced selection operates differently in this regime, as for the neutral results just described, and because reducing the number of males eligible to mate reduces the effective population size, rendering the deterministic selection against individuals homozygous for a previously sex-specific chromosome less effective. To illustrate these effects, Figure S10 in File S1 displays the fixation probabilities of a dominant X′ chromosome in an initially XX/XY system, with varying strengths of selection against YY males. Compared with the case of no sex differences (Figure 4C and Figure S10C in File S1), the fixation probabilities of the X′ are greatly increased when males have greater variance in reproductive success (Figure S10, A and B in File S1), so that, even for substantial degrees of selection against YY individuals, the X′ fixes with non-negligible probability. Another effect of the reduced effective population size is that conditional fixation times are substantially smaller than in the case of no sex differences (Figure S9 in File S1).

In comparing the substitution rates of dominant and recessive sex-determining mutations, we have assumed that their respective mutation rates, *i.e.*, the rates at which they are generated, are equal. It is possible that this is not the case, and that one class of sex-determining mutations is generated more rapidly than the other (Hillis and Green 1990; Bachtrog *et al.* 2011). If this were the case, it would simply be a distinct mechanism by which we expect one class of sex-determining gene to be more prevalent than the other. We should note that this is not as simple as comparing the rates of generation of gain-of-function and loss-of-function mutations. Suppose, for example, that a system is initially XX/XY male heterogametic, and consider a transition in direction *a* of Figure 1A, involving fixation of a dominant feminizing X′ mutation. Depending on the molecular functioning of the initial XX/XY system, this dominant X′ could be gain-of-function or loss-of-function. If the Y is initially dominant male-determining [as, *e.g.*, in mammals (Koopman *et al.* 1991; Beukeboom and Perrin 2014)], then a mutation on the X chromosome that blocks the Y’s activity (a gain-of-function mutation) would be dominant feminizing, as we require. If, however, the initial system depends on the ratio of some gene (or genes) on the X chromosome with respect to autosomes [as in *Drosophila* (Bridges 1921; Beukeboom and Perrin 2014)], then a loss-of-function mutation to a relevant gene on the X chromosome could be dominant feminizing.

We have studied direct transitions between male and female heterogamety. These appear to have been common, at least in vertebrates (Ezaz *et al.* 2006). However, transitions are also possible between heterogametic and environmental sex determination (Bull 1981, 1983; Quinn *et al.* 2011; Holleley *et al.* 2015), so that transitions between heterogametic systems could also occur via an intermediate system of environmental sex determination.

We have identified the force driving our results to be drift-induced selection operating along the equilibrium path (in the neutral case) or in its near vicinity (in the case with selection). While drift-induced selection as a force in natural selection has analogs that have been known for some time [*e.g.*, Gillespie’s criterion (Gillespie 1974, 1977), and selection in favor of reduced variance in offspring sex ratios (Verner 1965; Taylor and Sauer 1980)], the demonstration that it acts endogenously in systems of biological interest is relatively recent (Parsons *et al.* 2010; Lin *et al.* 2012; Kogan *et al.* 2014; Constable *et al.* 2016; Chotibut and Nelson 2017). We suspect that drift-induced selection will come to be recognized as an important force in many dynamical systems in population biology.

## Acknowledgments

We are grateful to David Haig, Dan Hartl, Sally Otto, Naomi Pierce, and the referees for helpful comments. The computations in this paper were run on the Odyssey cluster supported by the Faculty of Arts and Sciences (FAS) Division of Science, Research Computing Group at Harvard University. P.M. is supported by an National Science Foundation (NSF) graduate research fellowship. G.W.A.C. thanks the Finnish Center for Excellence in Biological Interactions and the Leverhulme Early Career Fellowship provided by the Leverhulme Trust for funding.

## Footnotes

Supplemental material is available online at www.genetics.org/lookup/suppl/doi:10.1534/genetics.117.300151/-/DC1.

*Communicating editor: J. Hermisson*

- Received May 24, 2017.
- Accepted August 15, 2017.

- Copyright © 2017 by the Genetics Society of America