Epistatic interactions between sex chromosomes and autosomes can affect the stability of sex determination systems

Abstract Sex determination (SD) is an essential and ancient developmental process, but the genetic systems that regulate this process are surprisingly variable. Why SD mechanisms vary so much is a longstanding question in evolutionary biology. SD genes are generally located on sex chromosomes which also carry genes that interact epistatically with autosomes to affect fitness. How this affects the evolutionary stability of SD mechanisms is still unknown. Here, we explore how epistatic interactions between a sexually antagonistic (SA) non‐SD gene, located on either an ancestral or novel sex chromosome, and an autosomal gene affect the conditions under which an evolutionary transition to a new SD system occurs. We find that when the SD gene is linked to an ancestral sex‐chromosomal gene which engages in epistatic interactions, epistasis enhances the stability of the sex chromosomes so that they are retained under conditions where transitions would otherwise occur. This occurs both when weaker fitness effects are associated with the ancestral sex chromosome pair or stronger fitness effects associated with a newly evolved SD gene. However, the probability that novel SD genes spread is unaffected if they arise near genes involved in epistasis. This discrepancy occurs because, on autosomes, SA allele frequencies are typically lower than on sex chromosomes. In our model, increased frequencies of these alleles contribute to a higher frequency of epistasis which may therefore more readily occur on sex chromosomes. Because sex chromosome–autosome interactions are abundant and can take several forms, they may play a large role in maintaining sex chromosomes.


| INTRODUC TI ON
In sexually reproducing species, the process of sex determination (SD) is an essential part of an individual's development, but the manner in which the sexual phenotype is set is far from conserved. An astounding variety of SD mechanisms has been described (Bachtrog et al., 2014;Beukeboom & Perrin, 2014); among organisms with genetic sex determination systems (GSD), there exists large variation in the genes and mechanisms that control the sexual phenotype.
In most GSD systems, the primary SD gene lies on a sex chromosome, resulting in either male heterogamety (males XY, females XX) or female heterogamety (females ZW, males ZZ). In some organismal groups, the SD gene (and by extension, the sex chromosome pair) that determines sex is strongly conserved, such as the SRY gene and the XY system of therian mammals (Graves, 2006).
However, other organismal groups exhibit substantially more variation, with different sex chromosome systems present in different species (Vicoso, 2019), such as in lizards (Ezaz et al., 2009;Pokorná & Kratochvíl, 2016), teleost fishes (Mank, 2009) and flies (Vicoso & Bachtrog, 2015). In addition to interspecific variation in SD mechanisms, intraspecific SD variation exists in several species, such as the southern platyfish Xiphophorus maculatus, in which X-, Y-and W-chromosomes are found (Orzack et al., 1980), and the housefly Musca domestica, in which some populations have an XY system and others a ZW system (Feldmeyer et al., 2008;Hamm et al., 2015).
The variability of SD mechanisms between and within organismal groups suggests that evolutionary turnovers between SD systems occur readily (Meisel, 2020;Vicoso, 2019).
Various population genetic models have been developed for evolutionary turnovers in SD systems (reviewed in van Doorn, 2014), of which two are of most interest here. First, sex ratio selection can favour a new SD gene when it induces development into the sex with the higher fitness, typically the minority sex (Fisher, 1930;Wilkins, 1995, but see Pen, 2006. Sex ratios can be biased due to, for example, sex chromosome meiotic drive (Jaenike, 2001;Kozielska et al., 2010), and selection can then favour a new SD gene that brings the sex ratio closer to 50:50. However, sex ratio selection can also favour rather than counteract deviations from equal sex ratios (Uller et al., 2007), and SD genes may also evolve when they actually cause such deviations (Kuijper & Pen, 2014). Second, linkage with sexually antagonistic (SA) loci has been proposed as a selective force in SD turnovers. As the regions flanking an SD locus are transmitted through males and females at different rates, SA loci can become genetically differentiated between the sexes. For example, a male-determining allele might become linked to a malebeneficial allele (on a primordial Y-chromosome) whereas chromosomes lacking the male-determining allele can become enriched for female-beneficial alleles (X-chromosome) (Charlesworth et al., 2014;Jordan & Charlesworth, 2012;Rice, 1984). Effectively, the SA locus and the SD locus evolve to form a co-adapted gene complex, and depending on the fitness effects and degree of linkage of SA and SD loci, the new gene complex may spread (van Doorn & Kirkpatrick, 2007. The acquisition of an SD gene on a chromosome initiates a process of sex chromosome differentiation (reviewed in Bachtrog et al., 2011;Charlesworth et al., 2005;Schenkel & Beukeboom, 2016). SA genes are expected to accumulate on the sex chromosomes along with the evolution of suppressed recombination on the Y-chromosome (or the W-chromosome in ZW systems) (Rice, 1987(Rice, , 1996. Subsequent degradation and masculinization of the Y-chromosome can help stabilize the SD system, by preventing it from becoming either fixed or lost (Marin & Baker, 1998). Overall, the stability of an SD mechanism can be affected by the association between the SD gene and nearby linked genes, and depending on the function of these linked genes different selective pressures may act on the SD gene.
Models on the evolution of SD mechanisms often focus on direct selection on the SD gene or the sex chromosome on which it is located. However, sex chromosomes represent only a fraction of the genome and the autosomes typically make up the majority.
Besides direct effects on the individual (e.g., by determining its sex), sex chromosomes may also have indirect effects through interactions with other (autosomal) genes, such as in humans (Bellott et al., 2014) and Drosophila melanogaster (Jiang et al., 2010;Lemos et al., 2008); in both species, the Y-chromosome harbours multiple genes that extensively regulate X-chromosomal and/or autosomal gene expression, and thereby eventually affect fitness. The evolution of gene expression differences and dosage compensation in recently formed sex chromosome systems suggests that even from an early point on sex chromosomes may interact with autosomes to affect fitness (Archer et al., 2017;Lachance et al., 2011;Zhou & Bachtrog, 2012). This is not surprising as sex chromosomes are thought to originate from autosomes (Ohno, 1967) and may prior to becoming sex chromosomes have been involved in autosome-autosome epistatic interactions. Although SA genes may accumulate on the sex chromosomes, they could also remain on the autosomes but become regulated by sex-chromosomal genes that control their expression (Parsch & Ellegren, 2013). Thus, although the sex chromosomes represent a specialized part of the genome, they can have crucial effects on autosomal gene expression and individual fitness by interacting with other components of the genome. Whether and how these interactions can influence the stability of SD mechanisms have however not been investigated yet.
Our aim is to investigate whether epistasis between autosomes and sex chromosomes can affect the stability of SD systems. We build on previous work by Van Doorn andKirkpatrick (2007, 2010) who investigated the influence of SA loci on transitions between SD mechanisms. Their models focus on two unlinked SD genes, each of which is linked to an SA locus. This mimics a situation in which the ancestral sex chromosome pair has begun differentiating into a full-fledged sex chromosome as described above, but has not yet undergone extensive genetic differentiation; the novel SD gene then arises near an autosomal SA locus. Depending on the selective pressures acting on the SA loci, the new SD gene may then invade or not.
Such transitions can be between identical sex chromosome systems (e.g., between different male heterogamety systems; van Doorn & Kirkpatrick, 2007) or between different types of sex chromosome systems (e.g., male heterogamety to female heterogamety or vice versa; van Doorn & Kirkpatrick, 2010). We focus here specifically on how epistasis alters the scope for turnover as predicted by these previous models. Thus, we investigate how epistatic interactions can affect the occurrence of SD transitions.

| Model overview
We provide here a conceptual description of our model; a more technical treatment is presented in the Appendix. We work with discrete, non-overlapping generations and random mating in a population with an infinite size. Offspring genotypes are determined based on Mendelian segregation whilst accounting for recombination, followed by viability selection based on their relative fitness. Our model features a diploid genome consisting of four different linkage groups ( Figure 1a). The first three linkage groups (XY, I A and II W ) each carry one SD locus and one SA locus that recombine at a rate r that can vary per linkage group. The fourth linkage group carries a single locus, called EPI, that interacts epistatically with the SA locus on XY, I A or II W to affect male fitness. Each locus has two possible alleles (referred to as the non-focal and focal alleles). The non-focal allele corresponds to a recessive allele without phenotypic effects (generally denoted +), whereas the focal allele affects the sex (for SD loci) or fitness (for SA loci) of an individual. We refer to the focal alleles by the name of their respective loci; all allele frequencies reported represent the frequencies of these focal alleles.
For the SD loci on XY and I A , the focal allele constitutes a maledetermining factor (Y and A, respectively), whereas the SD locus on II W corresponds to a dominant female determiner (W) which overrides the action of Y. The terms male and female are interchangeable, and hence, the model also applies to, for example, competing female-determining genes or transitions from female to male heterogamety. The SD genotypes that can be formed in either Y→A and Y→W transitions and their corresponding sex are listed in Table 1.
Genotypic fitness is defined as the relative viability of individuals carrying a particular genotype. An individual's fitness is determined by the genotype at the SA loci, whose effects depend on the individual's sex. Fitness effects of the focal allele at a single SA locus are determined by its fitness effect sizes in males and females s M and s F in homozygotes and additionally the sex-specific dominances for these effects in heterozygotes h M and h F (for details see Table 2). SA Y and SA A both have positive effects in males (s M > 0) and inversely negative effects in females (s F < 0). Conversely SA W has positive effects in females but negative effects in males (See Table S1). A female's total fitness is given by the product of the fitness scores of all SA loci, that is: F I G U R E 1 Model overview. (a) Genetic components of the model. All loci are labelled with their focal allele. Recombination rates between the SD and SA loci are given by r XY , r A and r W for linkage groups XY, I A and II W , respectively. (b) SD transitions scenarios considered. Epistatic interactions between loci are indicated in green. Only linkage groups which harbour SD genes involved in the transition and the linkage group carrying the EPI locus are depicted. All scenarios start out with a population where Y is the ancestral SD locus into which we introduce a new SD allele (either A or W). (c) Male and female karyotypes before (left) and after (right) transitions. Coloured chromosomes indicate the presence of an SD gene (Y, A or W) whereas white chromosomes indicate absence of an SD gene Here, w SA Y, w SA A and w SA W refer to the locus-specific fitness scores at the SA Y , SA A and SA W loci. In males, epistasis can further affect fitness, and therefore their fitness is given by: Here w EPI represents the fitness effect of epistasis. This is the outcome of interactions between an SA locus and the EPI locus.

| Epistasis scenarios and epistatic fitness effects
We let EPI interact with different SA genes to reflect situations in which either the established SD gene or the novel invading SD gene is linked to an SA locus that interacts epistatically with an autosomal gene. Although we are not aware of specific examples in which SA loci are indeed involved in epistasis, the existence of such loci is highly likely given that in several species, sex chromosomes are both enriched for sexually antagonistic genetic variation and for genes that play important roles in regulating autosomal gene expression (e.g., Innocenti & Morrow, 2010;Lemos et al., 2008). Alternatively, if such functions are performed by independent but linked genes, these may segregate as a single supergene to the same effect.
This scenario may be particularly relevant for sex chromosomes in which recombination suppression has recently begun to evolve (Rice, 1987). Epistatic interactions between EPI and an SA locus only occur in males, and their effects depend both on the genotype at the SA locus and the genotype at the EPI locus. We do not incorporate epistatic effects in females to limit the complexity of the model, whereas including it would likely only affect the dynamics of the model marginally. This is because the frequency of the focal allele at the SA Y locus is reduced in females and hence epistasis involving this locus would already occur at very low rates in females. For SA A and SA W , the focal alleles may be present at higher frequencies in females than SA Y . However, these loci are still autosomal prior to the SD transition, and therefore, the frequency of their focal alleles may only be slightly higher than the focal SA Y allele, so that epistasis involving these loci is similarly rare. Nonetheless, when the assumption that epistasis is sex-specific is not met, the scope for turnover from Y to A or W may differ slightly from that predicted by our model. Epistasis in general represents a situation in which the effect of one gene is modulated by another gene, and the manner in which such gene-gene interactions influence fitness may be modelled via numerous different approaches (reviewed in Wade et al., 2001). To explore all possibilities is therefore infeasible, and instead, we consider three standardized scenarios which we refer to as dominance, overdominance and coadaptation (see also The W/W genotype at II W cannot be obtained in our model as the W allele cannot be transmitted through males. b A low frequency of A alleles is introduced by mutation across all genotypes present in the population at that time; this results in small numbers of Y/+; A/+ individuals that decrease in frequency over time due to producing a 75% sex ratio (compared with favoured 50% sex ratios for males with a single Y or a single A allele).

TA B L E 2
alternative epistasis types are conceivable, they ultimately conform to minor variations to those considered here in that they share an underlying selective scenario. Interactions between the SA locus and EPI affect male fitness multiplicatively according to the factor w EPI = 1 + , where denotes the epistasis effect size, and the binary factor determines whether or not epistasis occurs or not. Table 3 lists the values of for every genotype combination in the different epistasis scenarios.

| Model initialization and sex determination transition types
In each scenario, we start with a standard XY system with a single male-determining allele Y which is fixed on the paternally inherited copy in males. New SD genes are not present in the ancestral popu-  Table S1. The selective effect parameters for the SA loci linked to SD loci involved (SA Y and SA A for the Y→A transition; SA Y and SA W for the Y→W transition) and the epistasis effect size were randomly sampled from uniform distributions with range (0, 0.05) for each independent simulation. For each combination of the four SD transitions and the three epistasis types, we ran 1000 independent simulations.

| RE SULTS
In our analysis, we focussed on the fitness effects of the SA loci linked to the ancestral and novel SD gene, as well as the effect size of epistasis on whether or not the new SD gene could invade or not.
We additionally varied the type of SD transition (male heterogamety to male heterogamety (Y→A) or male heterogamety to female heterogamety (Y→W)), the type of epistatic interactions (coadaptation, dominance and overdominance), and which SA locus engaged in epistatic interactions. We find that the SA effects of the linked loci remain a key determinant of whether or not SD transitions may take place as described by Van Doorn andKirkpatrick (2007, 2010).
However, epistatic interactions of different types affect the range of parameter values for which transitions take place. We find that the parameter range resulting in an SD transition is differently affected depending on (1) the type of SD transition, (2) the type of epistasis and (3) the gene which interacts with EPI.
For the Y→A scenarios, we find that interactions between SA Y and EPI tend to have a stabilizing effect on the Y allele as the male determiner for dominance and overdominance epistasis (Figure 2).
More specifically, the minimal sexually antagonistic fitness effect of SA A that results in an SD transition from Y to A is higher when SA Y interacts with EPI. The stabilizing effect is more pronounced when the effect of epistasis is higher; that is, epistasis has a stronger effect.
This effect however does not apply for overdominance epistasis involving SA A , where we instead observe that the scope for turnover is virtually unaffected (Figure 2, lower right panel). In contrast to the stabilizing effect of epistasis for dominance and overdominance epistasis, we find that for coadaptation epistasis, the effect of epistasis tends to be destabilizing except for when SA A is involved in epistasis ( Figure 2, Figure S1). When SA Y interacts with EPI, we find that A can invade for a large range of parameter values except for when epistasis is weak. Similarly, when SA A interacts with EPI, we find that A fails to invade and instead Y is maintained.
In the Y→W transitions, we find that the effects of epistasis on the scope for turnover are comparable with those for Y→A transitions ( Figure 3). Some differences do however exist; first, the effects of epistasis are much weaker for overdominance and dominance epistasis when SA Y is involved. We find that overdominance epistasis involving SA W has virtually no effect on the invasive capacity of W, which is similar to the case in Y→A transitions where overdominance epistasis involving SA A does not affect the scope for turnover. Taken together, overdominance epistasis involving the SA locus linked to the novel SD gene appears to have no effect on the conditions which permit this new SD gene to invade. For coadaptation epistasis, we again find that when SA Y is involved, this tends to promote turnover to W (Figure 3, Figure S2). When SA W is involved, the dynamics are slightly more complicated; when the sexually antagonistic fitness effect of SA Y is relatively weak, the effect of epistasis tends to favour its maintenance as the SD gene for higher values of W. However, as the selective effects associated with SA Y are higher, the scope for turnover becomes larger. Interestingly, the strength of epistasis appears not to have a major effect on the scope for turnover; rather, it is only the presence/absence of epistasis that affects the outcome.

| DISCUSS ION
We investigated whether epistatic interactions can affect the stability of and explain turnovers in SD mechanisms. Our model builds on previous work by Van Doorn andKirkpatrick, 2007, 2010), who showed that SA selection can drive evolutionary transitions between SD mechanisms. Our model is an extension in that male fitness can be affected by an epistatic interaction between an SA locus on the ancestral or the novel pair of sex chromosomes and a neutral autosomal locus. We considered transitions between different male heterogamety systems and from male heterogamety to female heterogamety in combination with three different types of epistatic interactions. We furthermore varied the strength of epistasis and the SA loci involved, and whether the ancestral sex chromosome or the invading sex chromosome is involved in epistasis.
We found that epistasis can affect the scope for SD transitions, but the manner in which it does so depends on a variety of factors.
For dominance and overdominance epistasis, epistasis tends to have either very little effect on the outcome of SD transitions (e.g., overdominance involving SA A or SA W ) or tends to promote maintenance of the ancestral sex chromosome pair. A possible explanation is that the allele frequencies of the SA locus on the ancestral sex chromosome pair have already diverged between the X-and Y-chromosome.
As the frequency of SA Y increases on the Y-chromosome, interactions between EPI and SA Y occur more frequently than interactions between EPI and SA A or SA W , who start out as autosomal SA genes and hence have a lower frequency in males. Effectively, under these conditions epistatic interactions are capable of enhancing stability of an ancestral SD system, but fail to enhance the invasive capacity of a new SD system; therefore, the effects of epistasis do not equally affect all SD genes. Instead, differentiation of the established sex chromosome pair leads to enrichment for alleles that engage in epistasis, thereby promoting its stability. Autosomal loci cannot become differentiated, so that they are not enriched for alleles involved in epistasis, and therefore, novel autosomal SD alleles do not experience the same benefit from epistatic interactions.
For coadaptation epistasis, the effects of epistasis tended to be destabilizing so that it facilitates turnover. In this scenario, doubly homozygous males (e.g., SA Y /SA Y ; EPI/EPI or +/+; +/+ genotypes) experience a fitness benefit from epistasis. When the epistasis ef-

F I G U R E 2
Maintenance of Y male heterogamety versus transition to A male heterogamety in Y→A transitions. Y may be maintained as the sex-determining locus depending on the strength of SA effects associated with SA Y (horizontal axis) and SA A (vertical axis) as well as the effect of epistasis (differently coloured lines). Lines indicate boundaries for the maintenance of Y, with Y generally being maintained when parameter values are below the boundary line and A invading when they are above the line (see indications in the plots). An exception applies for the coadaptation epistasis scenario involving SA A ; when = 0.05, A invades below the boundary line rather than above it (see also Figure S1). Horizontal bars indicate different epistasis types, whereas vertical bars indicate the SA locus involved in epistasis with EPI F I G U R E 3 Maintenance of Y male heterogamety versus transition to W female heterogamety in Y→W transitions. Y may be maintained as the sex-determining locus depending on the strength of SA effects associated with SA Y (horizontal axis) and SA W (vertical axis) as well as the effect of epistasis (differently coloured lines). Lines indicate lower boundaries for the invasion of W, with W being unable to invade and therefore Y being maintained as the sex-determining gene when parameter values are below the boundary line and W invading and Y being fixed when they are above the line (see indications in the plots). An exception applies for the coadaptation epistasis scenario involving SA Y ; when = 0.01, W invades above the top-left boundary line and below the bottom-right boundary line (see also Figure S2). Horizontal bars indicate different epistasis types, whereas vertical bars indicate the SA locus involved in epistasis with EPI coadaptation epistasis in our model favours fixation of the SA locus involved for either allele, it nullifies the ability for SA selection to favour the spread or maintenance of the linked SD locus. This results in the destabilization of the existing SD system or the inability of novel SD genes to invade.
The effects of epistasis provide another explanation for the apparent stability of some sex chromosome systems such as those of most mammals. Here, the sex chromosome system may not be stable solely due to the characteristics of a given SD gene (e.g., being insensitive to becoming regulated by a newly evolved upstream SD gene), but rather because of genetic differentiation of the region linked to the SD gene. This includes for example male-essential genes on the Y-chromosome that prevent its loss (e.g., as in the case of the Y→A transitions) or the decayed nature of older Y-chromosomes preventing fixation of Y, as homozygous YY individuals experience severe fitness costs (e.g., as in the case of Y→W transitions) (Bull & Charnov, 1977;Graves, 2006;van Doorn, 2014). Both of these effects however mostly apply to older sex chromosome pairs that have already persisted for extended periods of time. In contrast, the stabilizing effect of epistasis as reported here can apply from the very onset of sex chromosome evolution, as we show that a single locus involved in epistasis may already affect the stability of the sex chromosome system. This means that the stabilizing effect of epistatic interactions may occur on relatively undifferentiated sex chromosomes.
Although these effects are less substantial, we find they can be sufficient to prevent early displacement of an SD gene once established (although a sufficiently strong selective pressure on the new SD gene may still enable a transition). Over time, other factors such as acquisition of male-essential genes or recessive deleterious mutations may then further enhance the stability of ancestral sex chromosomes so that these can persist over extended periods of time.
As a caveat to the above, it must be noted that the effects of epistasis appear to depend on the type of SD transition considered, with the effects of epistasis being more pronounced in Y→A transitions as compared to Y→W transitions. A possible explanation is that in the latter, Y is fixed rather than lost. If the Y-bearing chromosome has become enriched for SA Y alleles (as described in Jordan & Charlesworth, 2012;Rice, 1987), the frequency of male-beneficial epistatic interactions does not decrease directly as W spreads in the population. This instead only drops later as the Y-bearing chromosome is no longer male-restricted, and therefore, the frequency of SA Y on this chromosome decreases as well. Even then, the frequency of SA Y in this new 'quasi-autosomal' state may still be higher than the frequency of SA Y on the ancestral X-chromosome (i.e., the non-Y-bearing chromosome that existed prior to the spread of W and fixation of Y), which instead had been enriched for the female-beneficial non-focal allele at the SA Y locus. In Y→A transitions, A-bearing males must also bear two such 'Xchromosomes' which severely reduces their odds of experiencing the benefits of epistasis. This poses an additional burden to the invasion of A that does not apply to invasion of W.
We focussed here specifically on a model involving SA loci, but other mechanisms capable of driving SD transitions may likewise be modulated by the effects of epistatic interactions (e.g., meiotic drive (Kozielska et al., 2010)). The benefit of Y-chromosomal differentiation with regard to SA loci, resulting in an increased frequency of epistasis, may apply more broadly to other genes as well, with the only requirement being that the Y-chromosome becomes enriched for an allele that engages in epistatic interactions. Examples of this include the evolution of Y-chromosomal regulating genes such as those regulating the expression of autosomal SA genes (Ågren et al., 2019).
Y-chromosomes (or W-chromosome in ZW systems) of several species have essential regulatory functions (Lahn & Page, 1997;Wright et al., 2014), as evident from their gene content and the impact of Y-chromosomal genetic variation in a variety of species (e.g., Bellott et al., 2014;Lemos et al., 2008). Therefore, given that Y-autosome interactions are prevalent in many species, the effects of such interactions on SD transitions may likewise apply in many organisms.
In this study, we have explored the effect of epistatic interactions between a sex chromosome (either ancestral or novel) and an autosome on the scope of turnover from an ancestral to a novel sex chromosome system. Our results demonstrate that such interactions can confer additional stability to an ancestral sex chromosome system for some types of epistatic interactions, whereas other interactions can reduce the stability of the ancestral sex chromosomes.
The capacity for sex chromosomes to become genetically differentiated relative to autosomes here enables epistatic effects to become more prevalent and/or pronounced, thereby resulting in increased stability of established systems. When a novel SD gene evolves on an autosome, no such genetic differentiation has occurred and therefore epistasis benefits the spread of novel SD genes to a lesser extent. The effect of epistasis on transitions in SD is further largely dependent on the type of sex chromosome transition considered and the strength of epistasis. In conclusion, the stability of a sex chromosome pair does not depend solely on its own characteristics, but instead should be considered as part of an interactive network with the remainder of the genome.

ACK N OWLED G EM ENTS
We would like to thank the Center for Information Technology of the University of Groningen for providing access to the Peregrine highperformance computing cluster. MAS is supported by an Adaptive Life grant awarded to LWB, IP and Jean-Christophe Billeter by the University of Groningen.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

AUTH O R CO NTR I B UTI O N S
MAS, LWB and IP conceived and designed the study; MAS collected and analysed the data; MAS drafted the initial version of the manuscript; MAS, LWB and IP contributed to later versions of the manuscript.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.13939.

O PE N R E S E A RCH BA D G E S
This article has been awarded Open Data, Open Materials Badges.
All materials and data are publicly accessible via the Open Science

Model initialization
To initialize the model, we first determine the frequency of different haplotypes at each linkage group based on a series of starting frequencies for the focal allele 1 (as described in the main text) on the maternally inherited (first) and the paternally inherited (second) copy, which are used to calculate the frequencies of the possible haplotypes on each linkage group at the first and the second copy.
For example, a 11 haplotype (which indicates a focal SD allele linked to a focal SA allele) on the maternal copy of linkage group XY is given by P(Y M 1 ) × P(SA Y 1 ). For linkage group III, which only carries the EPI locus, the haplotype frequencies are given simply by the allele frequencies of alleles 0 and 1. We define an array F i,j,k which contains the frequencies of the different haplotypes linkage groups XY, I A and II W , where i denotes linkage group (1 = XY; 2 = I A ; 3 = II W ), j denotes haplotype (1 = 00; 2 = 01; 3 = 10; 4 = 11), and k denotes the allele copy (1 = maternally inherited allele; 2 = paternally inherited allele).
For linkage group III, we similarly define an array F EPI i,j where i and j indicate which allele is present on respectively the maternal and the paternal copy (for both i and j, 1 = allele 0; 2 = allele 1).
Based on F i,j,k and F EPI i,j , we can define our initial population, which is given by an array P with dimensions H × H, where H represents a 4 × 4 × 4 × 2 array where each element represents a particular combination of haplotypes. Each element p in P represents a particular genotype, and the value of element p gives its frequency.
The frequency of each genotype p in P i,j,k,l,m,n,o,p (hereafter we use the subscript i to indicate an array with dimensions i through p) is given by F 1,1,i × F 2,1,j × F 3,1,k × F 1,2,m × F 2,2,n × F 3,2,o × F EPI l,p , where i and m indicate the haplotype on XY at the maternally and paternally inherited copies respectively, and similarly j and n, and k and o, for linkage groups I A and II W ; l and p indicate the genotype at the maternal and paternal copies of the EPI locus. Similarly, we define an array N h,i which counts the number of focal alleles 1 at locus h (1 = Y, 2 = SA Y , 3 = A, 4 = SA A , 5 = W, 6 = SA W , 7 = EPI) for genotype p in P ; values in N can be 0 (homozygous 0/0), 1 (heterozygous 1/0 or 0/1) or 2 (homozygous 1/1).
Because we randomly select the parameter values for the SA loci and the epistasis effect, the population upon initiation does not conform to a population at equilibrium. We therefore incorporate 10,000 generations of burn in during which allele frequencies can reach equilibrium prior to introducing new SD genes via mutation. This mutation procedure is described under 'Gametogenesis and reproduction'.

Sex determination
Sex determination takes place based on the number of focal alleles 1 at loci Y, A and W in each genotype (see Table S1). To do so, we define two binary arrays S M andS F , which describes for each genotype in whether that genotype is male (S M = 1 andS F = 0) or female (S M = 0 andS F = 1); we use superscripts F and M to identify arrays containing female-and male-specific components of the model throughout. A genotype in P i is male when (1) N 5,i = 0 and (2) N 1,i > 0 and/or N 3,i > 0. The frequencies of different genotypes in males (P M ) are given by the entrywise product of P • S M , and the frequencies in females (P F ) are given byP • S F . To keep track of the frequencies of the focal alleles per sex, the frequency of the focal allele at locus h in males is given by the inner productP M ⋅ N hi , whereas those in females are given byP F ⋅ N hi . Let a indicate the number of focal alleles at the SA Y locus of genotype p in P, which is given by N 2 , and similarly b (given by N 4 ) and c (given by N 6 ) be the number of focal alleles at loci SA A and SA W , respectively. The male fitness score of genotype p in P (W M ) is then given by

Sex-specific fitness effects and epistasis
and similarly its female fitness score (W F ) is given by w F XY a+1 × w F A b+1 × w F W c+1 (note that the +1 here is necessary to ensure that the proper index is used, i.e., a male with 0/0 genotype at SA Y has 0 focal alleles at this locus, but its fitness score for this locus corresponds to w F XY 1 = 1). The genotype frequencies in females after selection (G F ) are given by P F • W F . Here, k indicates which SA locus is involved in epistasis (k = 2 for SA Y , k = 4 for SA A and k = 6 for SA F ).

Gametogenesis and reproduction
Reproduction takes place by gametogenesis in males and females to yield pools of sperm and oocytes, respectively. To calculate the frequencies of the haplotypes amongst the sperm and oocytes, we first define an array U i,j,k,l , in which i denotes linkage group (1 = XY; 2 = I A ; 3 = II W ), jand kdenote haplotype on the maternally inherited and paternally inherited chromosome copies, respectively, and l denotes the haplotype to be sampled (for j, kand l : 1 = 00; 2 = 01; 3 = 10; 4 = 11), such that element u i,j,k,l gives the probability of sampling a haplotype l from a diploid genotype consisting of haplotypes jand kon linkage group i . In defining U, we also account for recombination which may yield novel haplotypes; for example, a genotype 00/11 can yield 01 and 10 haplotypes if recombination occurs (see Table S2). Similarly, we define an array V i,j,k,l , where j denotes the haplotype at the paternally inherited allele at the EPI locus, kdenotes that of the maternally inherited allele, and i denotes the type of allele to be sampled (for i , jand kalike, 1 = the non-focal allele 0, 2 = the focal allele 1), such that element v i,j,k,l gives the likelihood of sampling an allele of type i from a diploid genotype of jand k. Uand Vcan be used to construct an array T , which is defined as T i,j,k,l,m,n,o,p,q,r,s,t = U 1,i,m,q × U 2,j,n,r × U 3,k,o,s × V l,p,t (note that the order of subscripts in Uand Vhere deviates from that used above). The matrix product of G F Tyields an array H F containing the frequencies of each haplotype among the oocytes; similarly, the matrix product of G M Tgives an array H M which describes the haplotype frequencies among sperm. Genotypes are formed by fusion of sperm and oocytes, with genotype frequency being given by the product of their respective frequencies, such that the frequencies of genotypes in the next generation are given by the Kronecker product H F ⊗ H M , yielding an array Q. Qhas identical dimensions to Pand effectively represents the offspring produced by this ancestral population P. Therefore, we can update P = Qto represent moving forward one generation in our model. All simulations are carried out for 200,000 generations, during which we track the frequency of the focal allele 1 at each locus.
We introduce either A and W by mutations in gametes by manipulation of the H F and H M arrays. To introduce A into the population, in both the H F and H M arrays, we redefine H i3kl = H i1kl × A and H i4kl = H i2kl × A to convert a proportion A of 00 and 01 gametes into 10 and 11 gametes; we also redefine H i1kl = H i1kl × (1 − A ) and H i2kl = H i2kl × (1 − A ). Similarly for introducing W, we redefine H ij3l = H ij1l × W and H ij4l = H ij2l × W to the same effect, and redefine H i1kl = H i1kl × (1 − W ) and H i2kl = H i2kl × (1 − W ). In Y→A transitions, we introduce A at a frequency A = 10 −4 whereas in Y→W transitions, we introduce W at an identical frequency W = 10 −4 . A and W are introduced after 10,000 generations of burn in during which the initial population has been allowed to evolve to reach an equilibrium as dependent on the selective effect parameter values for the various SA genes and the epistatic effect.

Parameter value selection
In this manuscript, we focus on the fitness effects of the SA genes and the effect of epistasis primarily, and therefore do not vary some other aspects described here such as the recombination rates between the SD and SA genes or the dominances of the SA genes in males and females. The parameter range for the SA effects is taken to reflect a fitness effect between 0% and 5% in SA/SA homozygotes (where SA reflects the focal allele for a given SA locus) relative to +/+ homozygotes. To do so, we sampled the fitness effects for each SA locus from a uniform distribution with range (0, 0.5); the epistasis effect size is sampled from an identical uniform distribution. Fitness effects are identical in males and females but with inverse sign, that is, a fitness effect s in males is associated with a fitness effect − s in females and vice versa. In line with this, we assume a low rate of recombination (1%) which generally results in a case of strong linkage as defined by van Doorn and Kirkpatrick (2007).
Dominances for the fitness effects of SA alleles are set to 0.6 in the sex in which they have a beneficial effect (males for SA Y and SA A , females for SA W ) and 0.4 in the sex in which they have a deleterious effect. This configuration for recombination rate and dominance parameters should promote the maintenance of SA polymorphism for a large range of fitness effects (Jordan & Charlesworth, 2012). For both Y→A and Y→W transitions, we only consider fitness effects of SA alleles for those loci that are linked to one of the potential SD genes, for example, in Y→A transitions we include a fitness effect for SA Y and SA A , but not for SA W whose fitness effect is instead taken to be zero.