Establishment of a new sex-determining allele driven by sexually antagonistic selection

Abstract The turnover of sex-determining loci has repeatedly occurred in a number of species, rather than having a diverged pair of sex chromosomes. We model the turnover process by considering a linked locus under sexually antagonistic selection. The entire process of a turnover may be divided into two phases, which are referred to as the stochastic and deterministic phases. The stochastic phase is when a new sex-determining allele just arises and is still rare and random genetic drift plays an important role. In the deterministic phase, the new allele further increases in frequency by positive selection. The theoretical results currently available are for the deterministic phase, which demonstrated that a turnover of a newly arisen sex-determining locus could benefit from selection at a linked locus under sexually antagonistic selection, by assuming that sexually antagonistic selection works in a form of balancing selection. In this work, we provide a comprehensive theoretical description of the entire process from the stochastic phase to the deterministic phase. In addition to balancing selection, we explore several other modes of selection on the linked locus. Our theory allows us make a quantitative argument on the rate of turnover and the effect of the mode of selection at the linked locus. We also performed simulations to explore the pattern of polymorphism around the new sex-determining locus. We find that the pattern of polymorphism is informative to infer how selection worked through the turnover process.


Introduction
Recent genome analyses have demonstrated that genetic systems that determine sex are more labile than previously thought and that the turnover of sex-determining loci has repeatedly occurred. In some clades such as teleost fish and amphibians, sex-determining loci differ among closely related species or even within a species (Bachtrog et al. 2014;Beukeboom and Perrin 2014). Frequent turnover should occur in species especially when sex is determined by a single sex-determining locus, rather than a pair of highly diverged sex chromosomes. Turnover should be initiated by mutation at another potentially sex-determining locus, which could become a new sex-determining locus while the polymorphism at the original sex-determining locus disappears.
Many theoretical studies have investigated the evolutionary process of such turnover (Bull and Charnov 1977;van Doorn andKirkpatrick 2007, 2010;Kozielska et al. 2010;Blaser et al. 2013Blaser et al. , 2014Veller et al. 2017;Scott et al. 2018;Saunders et al. 2018Saunders et al. , 2019. A consensus has been established that, if a new sexdetermination system has higher fitness than the old one, the new system could potentially override the old one. This explains why turnover hardly occurs in species with a diverged pair of sex chromosomes, such as the X/Y system in mammals and the W/Z system in birds. Theoretical examinations of turnover usually involves a two-locus system, under which turnover has to pass through a phase in which dimorphic sex-determining alleles segregate at both of the two sex-determining loci. A deterministic theory assuming a population with an infinite size (Bull and Charnov 1977) demonstrated that the system with higher fitness can be stably maintained by selection. It is indicated that the fitness advantage of the new locus over the existing one would be an important factor for the turnover of sex-determining loci.
A possible scenario that a new sex-determining locus confers a fitness advantage is that the new locus arises in close linkage to a locus under sexually antagonistic selection (van Doorn and Kirkpatrick 2007Kirkpatrick , 2010. This process is explained in Figure 1A. Initially, sex is determined by the original X/Y locus, while the A/a locus is another potential sex-determining locus on a different chromosome (autosome). That is, the A/a locus is monomorphic (fixed for allele a) so that it does not play a role in sex determination, but when allele A arises from allele a by mutation, it creates a new sex-determination system. We refer to allele A as a sex-determining allele. There is another polymorphic locus (B/b) that is located close to the A/a locus. The B/b locus is assumed to be subject to sexually antagonistic selection; for example, allele B is beneficial in males and allele b is beneficial in females. At this point, the B/b locus does not confer any advantage or disadvantage because there is no physical linkage to the sex-determining locus X/Y. If a sex-determining allele (A) arises by mutation at the A/a locus, then one of the possible outcomes is that this new dimorphic locus takes over the role of sex determination and the original X/Y locus becomes monomorphic.
We are here interested in how often such turnover of sex-determining loci occurs. To understand this, it is crucial to theoretically describe the entire process, from the birth of a new sex-determining allele to its stable establishment. However, previous studies on this topic mainly by van Doorn and Kirkpatrick (2007Kirkpatrick ( , 2010 focused on the second half of the process (as explained below). The purpose of this work is to provide an analytical description for the first half, which largely determines how often turnover occurs under what conditions.
The entire process may be divided into two phases, which are referred to as the stochastic and deterministic phases. The stochastic phase starts when a new sex-determining allele A arises and continues until the frequency becomes high enough to escape from extinction due to random genetic drift. Then, the deterministic phase follows, in which allele A further increases in frequency by positive selection and genotype Aa becomes fixed in the heterogametic sex; in this way, the new A/a locus takes over the role in sex determination. The theoretical results currently available are for the deterministic phase, in which analytical treatment is quite straightforward because random genetic drift may be ignored. van Doorn and Kirkpatrick (2007Kirkpatrick ( , 2010 used a deterministic approach to describe the rate of increase in the frequency of allele A during the early stages of the deterministic phase. In contrast, the behavior of allele A in the initial stochastic phase is much more complicated and many factors are involved in it, which is the scope of this work. One such factor is the effect of random genetic drift, which is the major evolutionary force acting when the allele frequency is low. The linkage to B/b is also a very important factor. If allele A arises in linkage with allele B, allele A immediately benefits from positive selection because the haplotype A-B confers a selective advantage. On the other hand, if allele A arises in linkage with allele b, selection works against allele A because the haplotype A-b is deleterious. Recombination plays a crucial role in determining the rate at which the advantageous haplotype A-B is created and broken up. Another important factor is how selection works on the B/b locus. Depending on the parameter setting for the effect of the B/b locus on fitness, selection operates in various modes (see Figure 2): selection works for or against allele B, or even in some parameter space, balancing selection works to maintain the two alleles at intermediate frequencies.
It is easy to imagine that this parameter setting largely affects the fate of the sex-determining mutation.
The purpose of this work is to theoretically understand how often turnover of a new sex-determining A/a locus occurs. We here mathematically describe the probability that a newly arisen sex-determining locus turns over the old one, which is referred to as the establishment probability. We found that the establishment probability is markedly high in the case of balancing selection, while it is very low in other modes of selection. It is indicated that the mode of selection at the B/b locus is critical in determining the fate of a newly arisen sex-determining allele. Therefore, to understand the rate of turnover, it is necessary to know how many sexually antagonistic loci are under balancing selection. Our simulations provide insight into how to distinguish the mode of selection from the pattern of polymorphism around the sex-determining locus.

Model
We use a discrete-generation model of a diploid species with population size N. Three loci are considered in the model ( Figure 1A). One is the ancestral sex-determining locus, which is located on the sex chromosome where males have genotype XY and females have genotype XX. Although male heterogamety (XY system) is assumed here, our model can also handle female heterogamety (ZW system) by swapping males and females. The model includes another sex-determining locus A/a on the autosome, at which the initial state is that allele a is fixed, so that the locus does not play a role in sex determination. Then, we consider that masculinizing or feminizing allele A just arises by mutation at locus A/a. We are interested in how this new sex-determining allele A behaves in a finite population and how often it spreads and eventually enables the A/a system to override the old X/Y system. In our model, if allele A has a masculinizing effect, its turnover does not change the heterogamety sex and involves four genotypes (case 1, Figure 1B). If allele A has a strongly feminizing effect (i.e. strong enough to make XYAa a female) (van Doorn and Kirkpatrick 2010), the turnover changes the heterogamety system and involves six genotypes (case 2, Figure 1C). In either case, to evaluate the effect of the B/b locus, we assume that the levels of fitness of the new and old systems are identical. Therefore, the fate of allele A is mainly determined by the selection effect of another linked locus B/b. It is assumed that the recombination rate is r between the A/a and the B/b loci. We ignore mutation at the A/a locus to trace the fate of a single mutation, whereas at the B/b locus, recurrent mutation is allowed between alleles B and b such that this locus should be under selection-mutation balance. The mutation rate from allele B to allele b is denoted by u and v denotes the reverse mutation rate from allele b to allele B. The Relationship between sexes and genotypes when allele A has a feminizing effect. The red and blue stars are given to the genotypes that determine sex under the ancestral and new systems, respectively. The genotypes with no stars arise in the phase of transition from the ancestral to the new system. frequency of allele B is denoted by p. The fitness of genotypes BB, Bb, and bb is given by 1 þ s m ; 1 þ h m s m , and 1 in males and 1 þ s f ; 1 þ h f s f , and 1 in females, respectively. As we assume allele B is beneficial for allele A, we set s m > 0 and s f < 0 when allele A has a masculinizing effect (case 1), and s m < 0 and s f > 0 is set when allele A has a feminizing effect (case 2).

Data availability
The authors state that all data necessary for confirming the conclusions presented in the manuscript are represented fully within the manuscript. Codes used for numerical analyses and simulation experiments are available at https://github.com/TSakamotoevo/sex_codes_2020. Supplementary material is available at figshare DOI: https:// doi.org/10.25387/g3.13300967.

Results
We derive the probability that a single sex-determining mutation at the A/a locus spreads in the population so that the A/a locus becomes the new sex-determining locus and the old X/Y locus becomes monomorphic under the three-locus model illustrated in Figure 1A. This probability is essentially identical to the probability that the new sex-determining mutation successfully increases in frequency by avoiding extinction immediately after its introduction, as pointed out by van Doorn and Kirkpatrick (2007Kirkpatrick ( , 2010. Once the frequency of the new sex-determining mutation increases, it is very likely that it further increases to an intermediate frequency so that the sex-determining locus transitions from the X/Y locus to the A/a locus. This is because our model assumes that the presence of the B/b sex antagonistic locus provides a benefit only for the new A/a system. In the following, we derive this probability of the successful spread of allele A, which is referred to as the "establishment probability". The "establishment" means that the new system is used by all individuals in the population (i.e. the new system is fixed but alleles A and a coexist stably), rather than the old and new systems coexisting in the population. The establishment probability is derived separately for the two cases (cases 1 and 2).
Case 1: turnover without changing the heterogametic sex We first consider case 1, where allele A has a masculinizing effect ( Figure 1B). For allele B to be beneficial for allele A, we assume s m > 0 and s f < 0 in this section. We derive the probability of allele A escaping immediate extinction by using the branching process. Let u i ðp 0 Þ ði 2 fB; bgÞ be the establishment probability of allele A that arises in linkage with allele i at the B/b locus when the frequency of allele B is p 0 . If allele A links with malebeneficial allele B, allele A is favored by linked selection. On the other hand, if allele A links with female-beneficial allele b, allele A is disfavored by linked selection. By denoting the frequencies of haplotypes A-B and A-b by x B and x b , respectively, and ignoring second-order terms in s m ; r; u and v, the expected changes of the frequencies of A-B and A-b in one generation are given by where aðpÞ and bðpÞ are functions defined such as aðpÞ ¼ h m s m þ ðs m À 3h m s m Þp þ ð2h m s m À s m Þp 2 and bðpÞ ¼ Àh m s m p þ ð2h m s m À s m Þp 2 (for details, see Appendix B). It is interesting to point out that Equation (1) involves only selection parameters in males, not those in females. This means that E½Dx i directly Balancing selection works such that the average fitness is in the following order: depends on selection among males but does not necessarily mean that selection does not work in females. Selection in females is involved because it affects p, the frequency of allele B.
In the following, we first derive the establishment probability for a special case where selection at the B/b locus does not provide any systematic pressure on p. In this case, p does not change rapidly, so we can treat it as a constant, at least in the timescale of a newly arisen allele escaping from initial extinction (see Appendix C). With this assumption, we can obtain the establishment probability as a solution of a cubic equation, which is given by a function of p 0 , the frequency of allele B when allele A arises. Next, we consider a more general case where p changes, from the initial value p 0 to the equilibrium value p Ã . In this case, by contrast, we show that the establishment probability is given by a function of both p 0 and p Ã , of which the effect of p Ã is quite large. In either case, we first obtain the establishment probability conditional on p 0 , and then we derive the unconditional establishment probability by incorporating the stationary distribution of p 0 .

Establishment probability when a constant p can be assumed
We consider a case where we can assume p does not change significantly in the timescale in which we are interested in. In other words, the expected change of p in one generation, M p E½Dp, is assumed to be as small as $ 1 N , where M p is given by: if the second-order terms of s m ; s f ; u, and v are ignored. This assumption holds when h m % h f and s m % Às f , so that the levels of selection in the two sexes are of the same strength, but work in opposite directions. This mode of selection is referred to as "equalizing selection". We assume that s m ; s f ; u; v, and r are so small that their second-order terms can be ignored. Then, following the branching process approximation (Barton 1987), it is straightforward to show that the two establishment probabilities satisfy: (for details, see Appendix C). Notably, similar models were used to analyze the effect of a linked allele on the establishment of locally beneficial alleles (Aeschbacher and Bü rger 2014) or the effect of population structure on the establishment of a beneficial mutation (Pollak 1966;Barton 1987; Tomasini and Peischl 2018; Sakamoto and Innan 2019). Equation (3) can be reduced to a cubic equation and be solved analytically (for details, see Appendix D).
It is important to note that establishment is promoted by selection if u B ðpÞ is positive; otherwise, allele A is likely selected against and goes extinct. Whether u B ðpÞ > 0 depends on the leading eigenvalue of the matrix in Equation (1). When the leading eigenvalue is positive, u B ðpÞ and u b ðpÞ are given by: and where When r is small, they can be approximated in quite simple equations: We performed simple forward simulations in a Wright-Fisher population to check the accuracy of our derivation (for details of the simulations, see Appendix A). We confirmed that Equation (3) is in excellent agreement with the simulations for a wide range of the parameters. Some of the results are shown in Figure 3, where the two establishment probabilities, u B ðp 0 Þ and u b ðp 0 Þ, are plotted by assuming s m ¼ Às f ¼ 0:02; N ¼ 10; 000. Through this work, the mutation rates at locus B/b are fixed to be quite small values, u ¼ v ¼ 1:0 Â 10 À6 , unless otherwise mentioned. Note that the effect of mutation rate is small in the establishment process unless the mutation rate is very large (but the mutation matters when the stationary distribution of allele B is considered, as demonstrated below).
We first focus on the case of r ¼ 0, that is the A/a and B/b loci are completely linked, to investigate the effect of dominance (h). Suppose that allele B is recessive (h ¼ 0, left panel in Figure 3A). Then, a newly arisen allele A can benefit from the B/b locus only when it arises in a BB homozygote. In such a case, u B increases as p 0 increases to p 0 ¼ 0:5 (plotted in red), where the effect of sexually antagonistic selection is maximized. u B decreases as p 0 decreases from 0.5 to 1, making a symmetric function. With the assumption of no recombination, it is obvious that u b % 0 for any p 0 (plotted in blue). The weighted average of u B and u b (i.e. In contrast, in the dominant case (h ¼ 1, left panel in Figure 3C), u B for a small p 0 is quite high because a newly arisen allele A in linkage with allele B is immediately selected for, regardless of the genotype at the B/b locus. This selection works particularly efficiently when B is so rare that the selective advantage of A-B haplotypes is large in comparison with the population fitness. Therefore, u B is given by a monotonically increasing function with decreasing p 0 , but as p 0 decreases the probability that allele A arises in linkage with allele B decreases; therefore, the weighted average has a peak in the middle. An intermediate pattern is observed in a case of partial dominance (h ¼ 0.5, left panel in Figure 3B). Similar results were obtained for other values of selection intensity as long as s m ¼ Às f (not shown).
We next consider the effect of recombination. The approximations given by Equation (6) agree overall with the simulation results as long as r is small. As the recombination rate increases, u B decreases because the association with allele B becomes weaker. When r is relatively small, u b increases as the recombination rate increases because recombination gives a chance to link with allele B and the A-B association may be preserved if further frequent recombination does not break the association. However, this benefit does not hold for a large r, and u b decreases as the recombination rate increases because further recombination prevents the stable linkage of allele A and allele B and reduces the benefit of linkage.

Establishment probability when p changes
We next consider the case where p can change during the establishment process due to selection and mutation (i.e. M p ) 1 N ). To derive u B ðp 0 Þ and u b ðp 0 Þ, we use a continuous time approximation of the branching process (Barton 1995). With a deterministic approximation for M p (M p 6 ¼ 0), the two establishment probabilities satisfy: (for details, see Appendix C). This differential equation can be solved numerically by setting an initial condition, for which we use the establishment probability of allele A that arises when p is at stable equilibrium p Ã . u B ðp Ã Þ and u b ðp Ã Þ can be numerically computed by Equation (3) because we can assume p does not change significantly around p Ã . Technically, we cannot use the exact values of u B ðp Ã Þ and u b ðp Ã Þ as an initial condition because they violate the assumption of M p 6 ¼ 0. To avoid this problem, assuming a very small e, we use u B ðp Ã 6eÞ % u B ðp Ã Þ and u b ðp Ã 6eÞ % u b ðp Ã Þ as an initial condition, which does not markedly affect the numerical solutions as long as e is small.
When p is not constant, both p 0 and p Ã play an important role in determining u B ðp 0 Þ and u b ðp 0 Þ. The relative contribution of p Ã can be large when selection is strong and p 0 is not far from p Ã , so that p quickly approaches p Ã . In such a case, u B ðp 0 Þ and u b ðp 0 Þ cannot be large when u B ðp Ã Þ and u b ðp Ã Þ are very small. This argument can be explained as follows. Let us consider the fate of the descendants of an allele A that arises. After the mutation arises, p changes and finally reaches around p Ã . Denote the numbers of haplotypes A-B and A-b when p reaches p Ã as X B and X b , respectively. Then, the establishment probability is approximately given by: unless the population is very small. Because p Ã largely depends on the mode of selection on the B/ b locus (see Figure 2), we here consider u B ðp 0 Þ and u b ðp 0 Þ under three modes of selection separately (i.e. balancing, negative, and positive selection on allele B). For each mode of selection, we performed extensive forward simulations to check the performance of Equation (7). It was found that u B ðp 0 Þ and u b ðp 0 Þ computed by Equation (7) well agreed with the simulation results for a wide range of parameter space, and representative cases are shown in Figures 4-6 (see also Supplementary Figures S1 and S2).
First, we consider the case of balancing selection. To make a realization of balancing selection, it is necessary to set h m large enough to secure the average fitness of B/b heterozygote higher than 1. When h m is large, a newly arisen allele A would be immediately selected for together with allele B. Therefore, the overall behavior of u B ðp 0 Þ and u b ðp 0 Þ could be quite similar to that in the case of a high h when p does not change (e.g. Figure 3C where h m ¼ 1 is assumed). This is demonstrated in Figure 4, where we assume s m ¼ 0:02; s f ¼ À0:02; h m ¼ 1:0; h f ¼ 0:0 (as a consequence, p Ã ¼ 0:5), while all other parameters are the same as those used in Figure 3. To emphasize the difference from the case of a constant p, this figure also shows u B ðp 0 Þ and u b ðp 0 Þ computed by Equation (3) for comparison. The major difference from Equation (3) is that, when p 0 < p Ã ; u B ðp 0 Þ is lower than that in the case of constant p (the red broken lines in Figure 4) because the selective Error bars on the red and blue circles represent the 95% confidence interval, but they are too small to be seen.
advantage of haplotype A-B would be reduced by the rapid increase in the frequency of allele B by selection (i.e. both haplotypes A-B and a-B increase). u b ðp 0 Þ is overall very low, similar to the case of a constant p ( Figure 3C). On the other hand, when p 0 > p Ã ; u B ðp 0 Þ is larger than that in the case of constant p because the number of males having allele B decreases, with which haplotype A-B has to compete. Supplementary Figure S1, B and C, is for weaker selection coefficients (s m ¼ Às f ¼ 0:008 and 0.002), where u B ðp 0 Þ and u b ðp 0 Þ are overall reduced and the difference from Equation (3) is small perhaps because p does not change rapidly with weak selection.
We next consider the case of negative selection, where p would move from p 0 to p Ã , which is usually very low. If a very low p Ã is assumed, from Equation (3), we can approximate u B ðp Ã Þ and u b ðp Ã Þ as: These equations mean that the establishment of allele A is very unlikely when r > h m s m . Thus, the major difference from the case of a constant p is that, as the recombination rate increases, the establishment rate decreases to $0 around the threshold r > h m s m (see Appendix E for the behavior for large r). This is demonstrated in Figure 5, where s m ¼ 0:02; s f ¼ À0:04; h m ¼ h f ¼ 0:5 are assumed such that the comparable result for the case of a constant p is Figure 3B with the same selection parameters for males (all other parameters are identical). In Figure 5, r ¼ h m s m holds at r ¼ 0.01, which works as the threshold. When r is smaller than this threshold (r ¼ 0.0 and 0.005 in Figure 5), u B ðp 0 Þ and u b ðp 0 Þ are roughly in agreement with those in Figure 3B, although u B ðp 0 Þ and u b ðp 0 Þ are slightly higher than the case of a constant p, especially when p 0 is not small. The situation changes dramatically as the recombination rate exceeds the threshold (i.e. r ¼ 0.02 in Figure 5): u B ðp 0 Þ and u b ðp 0 Þ decrease to as low as $2=N (i.e. the neutral expectation for a sex-determining allele) where there is no benefit of linked selection and only drift-driven establishment occurs in a nearly neutral fashion. In contrast, in Figure 3B, u B ðp 0 Þ; u b ðp 0 Þ ) 2=N unless p 0 is close to 0 or 1. Such strong reduction of u B ðp 0 Þ when r exceeds h m s m ( Figure 5) can be explained as follows. A newly arisen allele A benefits from linked selection when it arises in association with allele B. Once the linkage is broken by recombination, allele A has almost no chance to recombine back to link to allele B because the frequency of allele B is very low. Therefore, after allele A loses linkage with allele B, there would be no selection for allele A so that the establishment of allele A has to rely on random genetic drift. Supplementary Figure S2, B and C, shows the results for weaker selection (s m ¼ 0:015 and 0.008), where the general pattern is similar to that in Figure 5, while u B ðp 0 Þ and u b ðp 0 Þ are overall reduced.
We finally consider the case of positive selection, where p Ã is generally very large (i.e. %1). Unless p 0 is small, p increases very quickly to p Ã % 1, that is the B/b locus is almost fixed for allele B. Therefore, even when allele A arises in association with allele B, allele A does not benefit from linked selection on locus B, resulting in a very low u B ðp 0 Þ. Particularly when allele A at p % p Ã , random genetic drift plays a role in the establishment process. The theoretical treatment for this situation is shown in Appendix E. Figure 6 shows u B ðp 0 Þ and u b ðp 0 Þ for the case of positive selection assuming s m ¼ 0:02; s f ¼ À0:01; h m ¼ h f ¼ 0:5. It is demonstrated that allele A can spread efficiently with the help of linked selection for allele B only when p 0 is small and r is so small that the initial association between A and B can be maintained for a while. The performance of Equation (7) is not as good as that in the cases of balancing selection and negative selection. It appears that Equation (7) underestimates the establishment probability because our derivation based on the branching process ignores establishment events occurring in a nearly neutral fashion.

Unconditional establishment probability
In the above, we consider the establishment probability as a function of the initial frequency of allele B, p 0 . We are here interested in the unconditional establishment probability, which is the weighted average over the stationary distribution of p 0 . Following Wright's formula (Wright 1931), the stationary distribution of p 0 , g(p), is given by: where V p ¼ pð1ÀpÞ

2N
and C is determined such that Ð 1 0 gðpÞ ¼ 1 Figure 4 Establishment probability for the case of balancing selection on allele B. Parameters are assumed to be s m ¼ Às Error bars on the red and blue circles represent the 95% confidence interval, but they are too small to be seen. For full version, see Supplementary Figure S1. Figure 5 Establishment probability for the case of negative selection against allele B. Parameters are assumed to be s m ¼ 0:02; s f ¼ À0:04; h m ¼ h f ¼ 0:5, N ¼ 10,000, and u ¼ v ¼ 1:0 Â 10 À6 . Error bars on the red and blue circles represent the 95% confidence interval, but they are too small to be seen. For full version, see Supplementary Figure S2. (Connallon and Clark 2012). It is well recognized that this formula works very well when the selection intensities, s f and s m , are relatively small so that their second-order terms are negligible (Crow and Kimura 1970;Ewens 2004). An intriguing exceptional case is when selection is weak but the absolute values of the selection intensities are not small. In the previous section, we considered a case where selection in males and that in females are well balanced when h f % h m and s f % Às m , irrespective of the absolute values of s f and s m . In such a special case of equalizing selection, the stationary distribution may be better given by: Then, the unconditional establishment probability of allele A, u, is given by: where u B ðpÞ and u b ðpÞ are given by Equation (3) for h f % h m and s f % Às m , and otherwise by Equation (7).
We can obtain an approximation of the establishment probability in a simple form for several special cases of r ¼ 0. When sexually antagonistic selection works as balancing selection, the establishment probability for r ¼ 0 is approximated by: u %2p Ã aðp Ã Þ %2p Ã ½h m s m þ ðs m À 3h m s m Þp Ã þ ð2h m s m À s m Þp Ã2 : It is implied that u is on the same order of magnitude as s m .
In the case of negative selection, assuming that selection is much stronger than mutation (u and v) and h m is not very small, we have p Ã % 0 and u $ u b ð0Þ. Furthermore, because u B ðpÞ ) u b ðpÞ and M p % 0, Equation (7) for r ¼ 0 is roughly simplified: Under these approximations, u is given by indicating that u is on the order of ffiffiffiffiffiffiffiffiffiffiffiffiffiffi vh m s m p . For selection to be dominant over random genetic drift, N 2 vh m s m ) 1 is required.
If sexually antagonistic selection works as positive selection and the mutation rate is low, linked selection no longer works and establishment is driven by random genetic drift, as discussed above. Then, the establishment probability is roughly given by 2 N .
In such cases, sexually antagonistic selection does not markedly increase the establishment probability. The unconditional establishment probability computed by Equation (12) is plotted for the four modes of sexually antagonistic selection in Figure 7. The approximations for r ¼ 0 [Equations (13) and (14)] are also presented by triangles on the y-axis, which show excellent agreement with the exact formula and simulation results. Three different values of the mutation rate are considered (u ¼ v ¼ 10 À5 ; 10 À6 ; 10 À7 ).
The establishment probability is in general highest when balancing selection works ( Figure 7A). This is because an intermediate p provides both the benefit of linkage and a higher chance to link with allele B to allele A. Because the mutation rate at locus B/b does not markedly influence the stationary distribution of p, the establishment probability does not depend on the mutation rate either.
When negative selection works at locus B, the establishment probability is quite low ( Figure 7B) because negative selection keeps allele B very rare. Mutation from allele b to B contributes to the establishment of allele A in two ways. First, it increases the frequency of beneficial allele B, resulting in a higher chance that allele A acquires linkage with allele B by recombination. Second, it increases the probability that haplotype A-b mutates to A-B. As a consequence, as the mutation rate v increases, the establishment probability becomes larger.
When positive selection works, the establishment probability is very low for all three mutation rates ( Figure 7C). This is because allele B is already fixed and linkage with allele B is no longer beneficial. Allele A may benefit from linked selection only when u is extremely high.
Thus, if selection is directional, p Ã is very low or high (under negative or positive selection, respectively), which does not markedly help the establishment of allele A. On the other hand, if balancing selection maintains alleles B and b in an intermediate frequency, allele A can take advantage of it. An intermediate situation is when equalizing selection works ( Figure 7D, where the process is nearly neutral because s f % Às m and h f % h m ). Because p fluctuates between 0 and 1 by random genetic drift, allele A can likely benefit from locus B/b if it arises when p is intermediate. This is why the establishment probability is largely affected by the mutation rate at locus B/b, which determines how often mutation is produced.
It is interesting to note that very tight linkage does not necessarily increase the unconditional establishment probability (e.g. Figure 7, B and D). With an increasing the recombination rate, in general, u b increases while u B decreases because recombination enhances exchange between haplotype A-B and haplotype A-b. Because the unconditional probability is given by the weighted average of u B and u b , there appears to be an optimal recombination rate to maximize u.

The pattern of neutral polymorphism after turnover
To investigate the effect of the mode of selection at locus B/b on the pattern of neutral polymorphism in a surrounding region, we Figure 6 Establishment probability for the case of positive selection for allele B. Parameters are assumed to be s m ¼ 0:02; s f ¼ À0:01; h m ¼ h f ¼ 0:5, N ¼ 10,000, and u ¼ v ¼ 1:0 Â 10 À6 . Error bars on the red and blue circles represent the 95% confidence interval, but they are too small to be seen.

Figure 8
The temporal dynamics of nucleotide diversity after turnover of the sex-determining locus. Results for a single simulation run with N ¼ 10,000 are shown for each mode of selection. The spacial distributions of the normalized p AÀA ; p aÀa , and p AÀa are shown in red, blue, and black, respectively, along the relative position in the (0,1) interval (also rescaled in units of 4Nr in parentheses). In the simulation, locus A/a is located at relative position 0.25 (vertical orange lines) in the (0,1) simulated region, while locus B/b is located at relative position 0.75 (vertical green lines). t is denoted as the number of generations since the new sex-determination system was fixed (i.e. since the old sex-determining locus X/Y became monomorphic). The recombination rate of the entire region was assumed to be 0.0002 such that the recombination rate between loci A/a and B/b ¼ 0.0001. The selection parameters used for each mode of selection were (A) s m ¼ Às performed forward simulations under the infinite-sites model. The spatial distributions of the nucleotide diversity after turnover in a typical run are plotted for each of the four modes of selection in Figure 8. The amount of nucleotide variation is measured by the average pairwise differences (p) standardized by the neutral expectation (h ¼ 4Nl where l is the neutral mutation rate). Haplotype i (i 2fA, ag) is denoted as the haplotype with allele i, and Figure 8 plots three p values: p AÀA for p within haplotype A (red line); p aÀa for p within haplotype a (blue line); and p AÀa for p between haplotypes A and a (black line). Let us first focus on the case of t ¼ 0 (i.e. just after the turnover has completed). When balancing selection works, p AÀa has a striking peak around the locus B/b ( Figure 8A) and a weaker peak is observed for p AÀA and p aÀa due to recombination between the two loci. This is because the B/b polymorphism had been maintained for a very long time by balancing selection before the turnover occurred, which is not observed in the other modes of selection (Figure 8, B-D).
As time passes, in all four cases, neutral mutations start to accumulate around the A/a locus, making a novel peak of divergence between haplotypes A and a (i.e. p AÀa ). The growth of the peak at the A/a locus proceeds while maintaining the peak at the B/b locus in the case of balancing selection ( Figure 8A). The pattern is as predicted by Kirkpatrick and Guerrero (2014). In contrast, in the other three cases (Figure 8, B-D), a peak newly arises at the B/b locus and grows along with the peak at the A/a locus. Thus, the patterns are similar in the cases of negative, positive and equalizing selection, whereas balancing selection is unique in that the peak at the B/b locus already exists, which is much higher than that at the A/a locus shortly after the turnover.
Case 2: turnover with changing heterogametic sex We next consider case 2, where allele A has a feminizing effect so that the heterogametic sex changes from male to female. That is, case 2 assumes s f > 0 for allele B to be beneficial for allele A. As in case 1, allele A can spread if allele A successfully avoids extinction shortly after it arises (van Doorn and Kirkpatrick 2010). Therefore, we again use the approximation of the branching process, following case 1. By ignoring the second-order terms in s f ; r; u, and v, the expected change of the frequency in one generation is given by: where By comparing Equation (15) with Equation (1), we notice that the two equations are identical if we replace h f by h m and s f by s m . It is indicated that the above arguments for case 1 are also applicable to case 2 by changing these variables. When equalizing selection works and p does not change rapidly (i.e. h m % h f and s m % Às f ), the establishment probability is given by Equation (3). On the other hand, when p changes when allele A is rare, the establishment probability depends on both p 0 and stable equilibrium p Ã . If balancing selection works, the establishment probability is given by Equation (7). When allele B is subject to negative selection, the establishment probability depends on whether r is larger or smaller than the threshold that is given by h f s f . When h f s f > r, the establishment probability is given by Equation (7). Whereas, when h f s f < r; u B ðp Ã Þ % 0 so that the establishment of allele A would be drift-driven once p reaches p Ã : therefore, u B ðp 0 Þ could be as small as $1=N. When allele B is subject to positive selection, we have u B ðp Ã Þ % 0. Consequently, the establishment process is driven by random genetic drift after p reaches p Ã . Therefore, as long as p 0 is near p Ã and positive selection is strong, the establishment probability is as low as $1=N (for more details, see Appendix E). Thus, the process in case 2 can be well described by the equations developed for case 1, as is demonstrated in Supplementary Figures S3-S6. The unconditional establishment probability is also given by Equation (12) (see Supplementary Figure S7).

DISCUSSION
In some clades such as teleost fish and amphibians, sex is often determined by a single locus rather than heteromorphic sex chromosomes. In such species, turnover of the sex-determining locus occurs so frequently that genetic divergence around the locus does not proceed. There are many factors that potentially promote the turnover of sex-determining loci, such as random genetic drift (Bull and Charnov 1977;Veller et al. 2017;Saunders et al. 2018), deleterious mutation load (Blaser et al. 2013(Blaser et al. , 2014, sexually antagonistic selection (van Doorn and Kirkpatrick 2007Kirkpatrick , 2010, sex ratio bias (Kozielska et al. 2010), and haploid selection (Scott et al. 2018). Recently, van Doorn and Kirkpatrick (2007Kirkpatrick ( , 2010 pointed out that turnover of sex-determining loci could be enhanced by linked selection. That is, a new sex-determining allele (allele A in our model) can be beneficial when it arises near a locus under sexually antagonistic selection (locus B/b). This article is aimed at understanding how often turnover of sexdetermining loci occurs with the help of linked selection at a nearby locus. Previous studies mainly by van Doorn and Kirkpatrick (2007Kirkpatrick ( , 2010 focused on the deterministic phase to obtain the rate of increase in the frequency of allele A, which is not sufficient to address the question on the rate of turnover, as mentioned in the Introduction. We here provide a full theoretical description of the behavior of allele A from when it newly arises to its establishment, where both the initial stochastic and following deterministic phases are taken into account. We provide some technical comments in Appendix F to provide intuitive insights on how to understand the results of van Doorn and Kirkpatrick (2007Kirkpatrick ( , 2010 in our framework. Our theory shows that the establishment probability is given by a function of the initial frequency of allele B, p 0 , and the equilibrium frequency, p Ã . It is demonstrated that the establishment probability largely depends on the mode of selection on allele B, which determines p Ã (Figure 7). The establishment probability is relatively high when balancing or equalizing selection works because polymorphism at locus B/b increases the benefit of linkage. In contrast, when directional selection works (either positive or negative), linked selection does not significantly help establishment. When negative selection works, the establishment probability is low because the frequency of allele B is so low that allele A has difficulty linking with it. When positive selection works, the establishment probability is also low because beneficial allele B should be almost fixed. The effect of p 0 appears to be smaller unless p 0 is very different from p Ã .
Our results demonstrate that the fate of newly arisen sexdetermining mutation is mainly determined in the early phases (i.e. stochastic phase), where random genetic drift is the major evolutionary force. It is indicated that the stochastic phase plays an important role in understanding the evolutionary process of turnover of sex-determining loci. In the stochastic phase, p 0 is a very important factor affecting the initial behavior of the new mutation, and the density distribution of p 0 largely depends on the mode of selection at the B/b locus. The establishment probability of allele A is very low when directional selection works at the B/b locus, indicating that the presence of loci under sexually antagonistic selection could significantly enhance turnover of sex-determining loci only when sexually antagonistic selection works in the form of balancing selection.
Our theory has great implications for the rate of turnover of sex-determining loci in natural populations. An important biological question is how sexually antagonistic loci help the turnover of sex-determining loci on the genomic scale. One might think that most turnover occurs with the help of sexually antagonistic loci under balancing selection, because the establishment probability is markedly high when balancing selection works at the B/b locus (see also van Doorn andKirkpatrick 2007, 2010). On one hand, one might consider that there may not be a large number of sexually antagonistic loci under balancing selection in a genome, so that their relative contribution might be small. Alternatively, there might be a negligible contribution of sexually antagonistic loci under other modes of selection (i.e. equalizing and directional selection). Even if the establishment probability is not high for each, on a genomic scale, their cumulative effect may not be small. To answer the question on the relative contribution of linked selection, we need to know how many sexually antagonistic loci exist in the genome, and what mode of selection is working. While empirical studies may be powerful to address this, it is also interesting to look at polymorphism data surrounding the sex-determining locus. Our simulations (Figure 8) provide insight into how to distinguish the mode of selection at a linked locus. Shortly after turnover, if divergence between male and female haplotypes is restricted in a very narrow region around the sex-determining locus, it may be likely that the sex-determining allele has become established with no help from linked selection (i.e. Myosho et al. 2012;Kamiya et al. 2012;Koyama et al. 2019). On the other hand, if a highly diverged region spreads surrounding the sex-determining locus, then either the sex-determining allele has become established together with a linked allele at a sexually antagonistic locus, or sexually antagonistic loci arose after the establishment of the sex-determining locus.

Funding
This work was partly supported by grants from the Graduate University for Advanced Studies, SOKENDAI, and Japan Society for the Promotion of Science (JSPS) to HI and TS (19H03207, 20J21760). Communicating editor: J. Ross-Ibarra

Appendices
Appendix A: Details of simulation model The details of the simulation model are presented. The life cycle is assumed to be in the order of mating (random genetic drift), selection, gamete production (recombination), and mutation. First, the details of simulations for establishment probabilities are provided. Let s ijk and e ijk be the frequency of genotype ijk among sperm and egg population at the beginning of the generation, respectively, where i 2 fA, ag, j 2 fB, bg and k 2 fX, Yg. In a mating step, N individuals are produced by random mating. Denote by f 0 ðijk=lmnÞ the frequency of zygotes between a sperm of genotype ijk and an egg of genotype lmn. E½f 0 ðijk=lmnÞ is given by s ijk e lmn . According to this probability, N zygotes are produced by multinomial sampling. Then, sex is determined by a combination of loci A/a and X/Y (see Figure 1, B and C).
During the selection step, fitness of each genotype is determined by locus B/b. Let F be a set of genotypes that grow into a female. The fitness of genotype X ¼ ijk=lmn, W(X), is given by In gamete production, genotype ijk/lmn produce eight kinds of gametes, ijk; ijn; imk; imn; ljk; ljn; lmk, and lmn, with proportion 1Àr 4 ; 1Àr 4 ; r 4 ; r 4 ; r 4 ; r 4 ; 1Àr 4 , and 1Àr 4 , respectively. After that, proportion u of allele B mutates to allele b while proportion v of allele b mutates to allele B. Then, the frequency s ijk and e ijk at the beginning of the next generation is defined.
Initial condition of simulations is determined as follows. In simulations for conditional probability, initial frequency of each gamete is set as After the mating step of first generation, an allele A is introduced. To investigate the fate of allele A arising in linkage with allele i (i 2 fB, bg), an allele i is randomly chosen and the linked allele a is changed into allele A. The simulation is run until either locus X/Y or locus A/a becomes monomorphic. In simulations for unconditional probability, we first run simulations without introducing allele A to obtain stationary distribution of s ajk and e ajk . After a burn-in period of 10 N generations, s ajk and e ajk are sampled. For each set of s ajk and e ajk , an allele a is randomly chosen and turned into allele A after the mating step of first generation. The simulation is run until either locus X/Y or locus A/a becomes monomorphic.
Next, the details of simulations for the pattern of neutral polymorphism are provided. A genomic region of relative length 1 is simulated where the locus A/a is located at relative position 0.25 and the locus B/b is located at relative position 0.75. For neutral sites, the infinite-site model is assumed. Although the basic assumptions are the same as the three loci model, implements are slightly different to improve efficiency of simulations such that all events are incorporated into the mating step.
In a mating step, a father and a mother are chosen from individuals in previous generation to form an offspring. To incorporate selection, the probability that an individual is chosen as a parent is proportional to its relative fitness among individuals of the same sex. For each parent, a haplotype is made through recombination where a constant rate of recombination across the simulated region is assumed. By combining two haplotypes, an offspring is formed. Mutation is finally introduced in both the locus B/b and neutral sites. These procedures are repeated N times to generate individuals of the present generation.
The dynamics of polymorphism are simulated as follows. First, to simulate the polymorphism just before turnover, we run simulations for 200N generations without introducing allele A. Then, allele A is introduced at the locus A/a. If turnover occurs, simulation continues until 50N generations pass since the turnover. the leading eigenvalue of the matrix in Equation (1). We also assume that kðpÞ ) 1 N such that allele A increases deterministically once its frequency becomes large. Note that e ) 1 N is required for kðpÞ ) 1 N . To derive the establishment probability, we focus on descendants of an allele A after one generation. Let a random variable n j i be the number of offsprings of haplotype A-j that is left by a haplotype A-i (i; j 2 fB, bg). Because we assume that allele A is rare, each n j i can be treated as independent and the probability distribution of n j i may be approximated by Poission distribution. Then, the probability generating functions of the number of offspring are given by where terms of Oðe 2 Þ are ignored in the exponent and x and y are indeterminates. Let Dp be the change of p in one generation. Then, establishment probabilities, u B and u b , satisfy following equation: Assume that u B ðpÞ; u b ðpÞ are the order of e. Then, Equation (C2) is expanded as ½ð1 À pÞr þ uE Dp ½u b ðp þ DpÞ þ ½1 þ aðpÞ À ð1 À pÞr À u E Dp ½u B ðp þ DpÞ À u B ðpÞ À 1 2 E Dp ½u B ðp þ DpÞ 2 ¼ Oðe 3 Þ ½1 þ bðpÞ À pr À vE Dp ½u b ðp þ DpÞ þ ½pr þ vE Dp ½u B ðp þ DpÞ Noting that evolutionary forces are so weak that the frequency of allele B changes slowly, E Dp ½u i ðp þ DpÞ and E Dp ½u i ðp þ DpÞ 2 can be approximated by E Dp ½u i ðp þ DpÞ%E Dp u i ðpÞ þ Dp du i ðpÞ dp þ ðDpÞ 2 2 Recall that e ) 1 N is assumed such that kðpÞ ) 1 N is satisfied. By substituting Equation (C4) into Equation (C3) and taking terms up to e 2 , Equation (7) is derived. By assuming M p is negligible (i.e. $ 1 N ), Equation (3) is also derived. The scope and the limitation of the approximation are discussed. Equations (3) and (7) are accurate if terms of order 1 N are ignored, as long as kðpÞ ) 1 N . Noting that u i ðpÞ $ 1 N when kðpÞ $ 1 N , they are generally accurate if the order of 1 N can be ignored. Since such small terms can be ignored in many cases, our approximations may be valid in broad situations.
However, there are exceptional situations, for which terms of order 1 N become too large to be ignored. Such situations occur when kðpÞ $ 1 N when p % p Ã while kðpÞ ) 1 N otherwise. Typical situations are that negative selection works on allele B and r is large, or positive selection works. In such cases, it is probable that the number of allele A may increase $N before p reaches p Ã , especially in small populations. Although u i ðpÞ are still proportional to 1 N , the absolute value of establishment probability may become large. For the details of these cases, see Appendix E. linked selection is very weak at p % p Ã We describe the behavior of the establishment probability when linked selection is very weak at p % p Ã . This is a typical situation when negative selection works against allele B and r is large, or positive selection works for allele B, so that the establishment probability is very small (unless p 0 is very different from p Ã ). We here explore such a case in detail because Equation (7) could underestimate the establishment probability (see Figures 5 and 6). This is because our derivation based on the branching process approximation ignores establishment events occurring in a nearly neutral fashion. In this appendix, we derive another approximation focusing on the establishment of allele A through a nearly neutral fashion. It is assumed that selection is strong and p 0 is close to p Ã .
First, we consider case 1 (i.e. the turnover without changing heterogametic sex) and negative selection on allele B is assumed. For linked selection to be very weak at p % p Ã ; r > h m s m is also assumed [see Equation (9)]. Under this condition, the establishment process does not occur simply by increasing the frequency of haplotype A-B: rather, because frequent recombination keeps breaking the linkage between alleles A and B, allele A has to increase without increasing allele B. In practice, negative selection works to reduce the frequency of allele B and allele b is almost fixed. Therefore, establishment occurs such that haplotype A-b increases. In this case, haplotype A-b is selectively neutral, so that we can derive an approximation of the establishment probability by following Equation (8). By using u b ðp Ã Þ % 2 N , where 2 N is the establishment probability of a neutral masculinizing allele, the establishment probability is approximated by where X b is the number of haplotype A-b when p reaches p Ã . We can calculate EðX b Þ numerically under the branching process approximation as described below. Next, positive selection on allele B is assumed in case 1. At equilibrium with p Ã % 1, haplotype A-B is selectively neutral because the frequency of allele b is very low. Similar to the previous case, we can derive an approximate formula for the establishmemt probability as 2 N Â EðX B Þ; where X B is the number of haplotype A-B when p reaches p Ã . We briefly explain how EðX B Þ and EðX b Þ can be computed. We here assume that mutation rate is low enough to be ignored. Let X i j ðp 0 Þði; j 2 fB; bgÞ denote the number of descendant A-j haplotypes at p ¼ p Ã originated from the focal single haplotype A-i that arose when p ¼ p 0 . EðX B Þ and EðX b Þ depend on p 0 and the allele at locus B/b with which the allele A initially links. By considering the change of the haplotype frequency in one generation, we can obtain the recursion equations: E½X B j ðpÞ ¼ ð1 þ aðpÞ À ð1 À pÞrÞE½X B j ðp þ DpÞ þ ð1 À pÞrE½X b j ðp þ DpÞ E½X b j ðpÞ ¼ prE½X B j ðp þ DpÞ þ ð1 þ bðpÞ À prÞE½X b j ðp þ DpÞ (E3) [for each coefficient, see Equation (1)] where Dp is the expected change of p in one generation. By ignoring the order of 1 N and using a continuous time approximation, Equation (E3) is rearranged as dE½X B j ðpÞ dt ¼ ÀðaðpÞ À ð1 À pÞrÞE½X B j ðpÞ À ð1 À pÞrE½X b j ðpÞ dE½X b j ðpÞ dt ¼ ÀprE½X B j ðpÞ À ðbðpÞ À prÞE½X b j ðpÞ dp dt ¼ h f s f þ h m s m 2 þ ð1 À 2h f Þs f þ ð1 À 2h m Þs m 2 ppð1 À pÞ; where t is an auxiliary variable. Then, by setting the initial condition at t ¼ 0, we can obtain E½X i j as a function of p. When negative selection is assumed, the initial condition is set by the values at equilibrium p Ã ¼ 0 as EðX B b Þ ¼ r rÀhmsm and EðX b b Þ ¼ 1. However, for a technical reason, we cannot use p ¼ 0 at initial state because we cannot calculate Equation (E4) under dp dt ¼ 0 (see the main text). Instead, we set p ¼ e ðe ( 1Þ as an initial state. When positive selection is assumed, EðX B B Þ ¼ 1; EðX b B Þ ¼ r rþð1ÀhmÞsm and p ¼ 1 À e is used as values at t ¼ 0. As long as e is small, its effect on the numerical value appears to be very subtle. For the case of negative selection, the establishment probability is plotted for different population sizes N ¼ 10,000 and 100,000 (Supplementary Figure S8). All other parameters are the same as those used for Supplementary Figure S2B, so that Supplementary Figure S8A is identical to Supplementary Figure S2B. In Supplementary Figure S8, simulation results are plotted together with Equations (3), (7), and (E1). Note that Equation (E1) can be applied only when r > h m s m . It is demonstrated that, when r < h m s m , as stated in the main text, establishment is driven by selection such that its probability does not depend on the population size (r ¼ 0:0; 0:005 in Supplementary Figure S8). On the other hand, when r > h m s m , the establishment probability is no longer determined by the selection intensity, and proportional to 1=N because random genetic drift is the dominant force involved in the establishment. This is shown in the inner panels for r ¼ 0:01; 0:02 in Supplementary Figure S8, where Equation (E1) agrees with the simulation results better than Equations (3) and (7) especially for a large r (i.e. r ¼ 0.02).
Supplementary Figure S9 shows the establishment probability for the case of positive selection. Simulation results are plotted together with Equations (3), (7), and (E2). Two population sizes (N ¼ 10; 000; 100; 000) are considered, and all other parameters are the same as those used for Figure 6 so that Supplementary Figure S9A is identical to Figure 6 except that the establishment probability is log-scaled in Supplementary Figure S9. For all recombination values, as the establishment process depends on random genetic drift, the establishment probability is on the order of 1=N (see Supplementary Figure S9). Equation (E2) agrees very well with simulations unless p 0 is very small.
Finally, we consider case 2 that involves a turnover of the heterogametic sex. Approximation formulae can be derived in a similar way to case 1. By noting that the establishment probability of a neutral feminizing allele is 1: 07=N (Veller et al. 2017), the establishment probability can be approximated by,