Martijn A Schenkel,Jean-Christophe Billeter,Leo W Beukeboom,Ido Pen
Divergent evolution of genetic sex determination mechanisms along environmental gradients
clinal distribution,environmental sex determination,female determiner,genetic sex determination,polymorphism,housefly,male determiner,sex chromosomes,temperature,transformer
Abstract Sex determination (SD) is a crucial developmental process, but its molecular underpinnings are very diverse, both between and within species. SD mechanisms have traditionally been categorized as either genetic (GSD) or environmental (ESD), depending on the type of cue that triggers sexual differentiation. However, mixed systems, with both genetic and environmental components, are more prevalent than previously thought. Here, we show theoretically that environmental effects on expression levels of genes within SD regulatory mechanisms can easily trigger within-species evolutionary divergence of SD mechanisms. This may lead to the stable coexistence of multiple SD mechanisms and to spatial variation in the occurrence of different SD mechanisms along environmental gradients. We applied the model to the SD system of the housefly, a global species with world-wide latitudinal clines in the frequencies of different SD systems, and found that it correctly predicted these clines if specific genes in the housefly SD system were assumed to have temperature-dependent expression levels. We conclude that environmental sensitivity of gene regulatory networks may play an important role in diversification of SD mechanisms.
Divergent evolution of genetic sex determination mechanisms along environmental gradients

Sex determination (SD) is a crucial developmental process, but its molecular underpinnings are very diverse, both between and within species. SD mechanisms have traditionally been categorized as either genetic (GSD) or environmental (ESD), depending on the type of cue that triggers sexual differentiation. However, mixed systems, with both genetic and environmental components, are more prevalent than previously thought. Here, we show theoretically that environmental effects on expression levels of genes within SD regulatory mechanisms can easily trigger within-species evolutionary divergence of SD mechanisms. This may lead to the stable coexistence of multiple SD mechanisms and to spatial variation in the occurrence of different SD mechanisms along environmental gradients. We applied the model to the SD system of the housefly, a global species with world-wide latitudinal clines in the frequencies of different SD systems, and found that it correctly predicted these clines if specific genes in the housefly SD system were assumed to have temperature-dependent expression levels. We conclude that environmental sensitivity of gene regulatory networks may play an important role in diversification of SD mechanisms.

Sex determination (SD) is a crucial aspect of the development of sexually reproducing organisms, yet the regulatory mechanisms underlying SD are very diverse and prone to evolutionary change (Bachtrog et al., 2014; Beukeboom & Perrin, 2014). SD mechanisms have traditionally been classified as either environmental (ESD) or genetic (GSD) depending on the type of signal that initiates the determination of an individual’s sex. Under ESD, such signals include temperature, salinity, and acidity during a sensitive period in embryonic development (reviewed in Devlin & Nagahama, 2002; Janzen & Paukstis, 1991). Under GSD, the signal is a specific gene (or set of genes) present in zygotes, leading to either male or female development, such as the male-determining Sex-determining Region Y (SRY) gene in mammals or transformer (tra) in many insects (Berta et al., 1990; Goodfellow & Lovell-Badge, 1993; Pane et al., 2005; Sinclair et al., 1990; Verhulst et al., 2010). Diversification of SD mechanisms occurs via the evolution of novel SD mechanisms that replace their predecessors, a process called SD transition. SD transitions are prevalent in some taxa but not in others (Bachtrog et al., 2014; Beukeboom & Perrin, 2014), indicating variable evolvability of SD systems. What enhances the evolvability of some SD systems but not others and what causes SD transitions is still poorly understood. One often-overlooked aspect of the evolution of SD is how environmental and genetic factors may simultaneously affect SD. Rather than being mutually exclusive, ESD and GSD could instead be considered as two extremes of a continuum (Beukeboom & Perrin, 2014; Pen et al., 2010; Uller & Helanterä, 2011), with mixed systems occurring in several organismal groups, such as amphibians, fish, and insects (Devlin & Nagahama, 2002; Doums et al., 1996; Hamm et al., 2015; Ma et al., 2016; Mezzasalma et al., 2021; Nigro et al., 2007; Pen et al., 2010). In such mixed systems, SD may reflect a delicate balance between genetic effects that bias the process of SD toward male or female development, counteracted by environmental effects that push the system in the opposite direction (Alho et al., 2010; Holleley et al., 2015). GSD has repeatedly evolved in species which previously had ESD (e.g., Ezaz et al., 2009), and several theoretical models have been developed to predict when such turnovers should occur (Muralidhar & Veller, 2018; Pen et al., 2010; Van Dooren & Leimar, 2003). However, such models typically do not explicitly consider the underlying molecular mechanisms of sex determination. Although sex is determined by genes in species with GSD, environmental conditions can affect SD by modifying the expression levels of these genes (Hodson et al., 2017; Schmidt et al., 1997). Despite clear evidence that such environmental effects may perturb the action of GSD mechanisms (e.g., Devlin & Nagahama, 2002), their effect on the evolution of GSD systems is still unknown. In many species, SD involves hierarchical gene regulation (Kopp, 2012; Pomiankowski et al., 2004; Wilkins, 1995), where an initial signal targets a small number of regulatory genes that in turn regulate downstream targets. Evolutionary transitions between GSD systems are thought to occur primarily by the displacement of the initial signal gene by another gene with a similar function, or by the recruitment of a new regulatory gene on top of the ancestral SD regulatory pathway (Bull & Charnov, 1977; Pomiankowski et al., 2004; Wilkins, 1995). Genes downstream of the top regulatory genes are considered to be more constrained in terms of evolutionary change, as mutations in such genes may interfere with their pre-existing function in regulating SD. Nonetheless, they are not fully prohibited from undergoing further evolution, and some changes may still occur (Herpin et al., 2013). A new evolutionary framework needs to integrate the views that (1) SD is not solely environmentally or genetically determined but is to varying degrees affected by both types of cues; and (2) changes in SD cascades do not only occur at the top but may also occur via changes in the underlying genetic network. Here, we formalize this framework by developing a theoretical model of the evolution of SD systems in the presence of spatial environmental variation that affects expression of SD genes. The model is inspired by the polymorphic SD system of the housefly Musca domestica, but can be applied more broadly to other systems as well (Table 1). In this system, temperature is likely to act as an environmental cue because (1) variation in SD systems is correlated with variation in temperature between natural populations, and (2) temperature affects SD in several M. domestica strains harboring mutant SD genes (see also Box 1). Like in M. domestica SD, the model features two types of SD genes: an environmentally-sensitive gene F induces femaleness when active, and one or more M-genes that induce maleness by inhibiting F (Figure 1, see details below). We investigate here how the (co-)evolution of F and M can yield novel SD systems under the influence of environmental sensitivity of F. We then use the model to explain how the multifactorial SD system of houseflies has evolved (Box 1).

Here, we briefly describe the model; a more detailed description of the model and simulation techniques is provided in the Supplementary Methods. We developed an individual-based simulation model, where individuals occupy a linear array of subpopulations (demes) arranged along a temperature gradient. The life cycle is as follows: adults reproduce sexually and then die; their offspring undergo sexual development and viability selection; a fraction of the surviving offspring migrates from their natal subpopulation to a neighboring subpopulation; finally, individuals mature and the next round of reproduction begins. Motivated by the SD mechanisms of the housefly (Box 1), individual sexual development was modeled to result from interactions between a single feminizing gene F, one or more masculinizing genes M and the local temperature of an individual’s environment; note that labels such as “female” and “feminizing” are interchangeable for “male” and “masculinizing” to suit other SD mechanisms, e.g., ZW systems instead of XY systems. The F gene produces a temperature-dependent amount of product which is partially inhibited or degraded by the products from the M genes; the remaining or net amount of F product, denoted by , determines sex: if exceeds a threshold value , the individual develops as female, whereas if is below a second threshold value it develops into a male. If is between the two thresholds, , then the individual develops into an infertile intersex. The value of is obtained by summing up the net expression levels of both F alleles in an individual. The quantitative value of can vary between F alleles and depends on (1) the temperature-dependent expression level of the allele, (2) the allele’s sensitivity to M product, and (3) the amount of M product. Specifically, normalized temperature varies between at one end (the “north”) of the array of subpopulations and at the other end (the “south”) and it affects the net amount of F product according to The first factor on the right, , represents the initial amount of F product, before partial degradation by M product, where is the F allele’s baseline expression level at , quantifies the linear rate of increase in F expression with temperature, and represents random variation in expression levels due to environmental and/or developmental noise. The second factor, , represents the proportion of F product that is not degraded by M product, where is the F allele’s sensitivity to M product and is the cumulative amount of M product produced by the individual’s M alleles. The baseline expression level and sensitivity of F alleles are evolvable trait values, as are the expression levels of M alleles. Whenever a gamete is produced, each allelic trait incurs a mutation with a certain probability and its trait value is modified. Most mutations are “regular” mutations that modify traits by adding a small Gaussian amount with mean zero, but a small fraction of mutations are “null mutations” such that the resulting trait values (allelic expression levels or sensitivity to M product) are zero and cannot evolve any further. See Supplementary Table 2 for all model parameters and their default values used in the simulations. The initial populations all have an XY male heterogametic system: all individuals carry two F alleles on an autosomal F locus and males additionally carry a single M allele (designated MY) on their Y-chromosome. Initially there are no M alleles on autosomes, but we assume that during meiosis sometimes a new M allele (designated MA) is created on an autosomal locus; this is assumed to occur via de novo evolution of a novel M allele, but can also occur via transposition from a Y-chromosomal M allele. Individuals carry at most four M alleles: two on the Y-chromosome (if they have two Y chromosomes) and two on an autosome. In natural systems, many loci may be capable of evolving a male-determining function (Bopp, 2010), but we consider here only a single autosomal locus for simplicity. Thus, the initial XY population has the potential to evolve into a population with a new system of male heterogamety where an autosomal chromosome with an M allele has replaced the original Y-chromosome. Evolution of populations with female heterogamety is also a potential outcome, if females become heterozygous for an insensitive F allele (i.e., with ) with a sufficiently high expression level. We also allowed for Y-chromosomal fitness effects: (1) individuals homozygous for the Y-chromosome will have reduced viability and (2) Y-chromosomes carry sexually antagonistic alleles that are beneficial to males and detrimental in females, with additive effect males and in females (where is the antagonistic effect). The combined effects of Y-homozygosity and sexual antagonism are assumed to be multiplicative, i.e., a male with two Y-chromosomes will have his expected fitness modified by a factor , whereas a female with two Y-chromosomes gets .

As an initial exploration of our model, we performed a set of 10,000 simulations without temperature-dependent overexpression, but for different values for the SD thresholds and , as well as different MA activation rates (sampled from a uniform distribution with range [0, 0.05]). Here, de novo activation of MA causes a male-biased sex ratio, thereby promoting the invasion of female-determining alleles (Bull & Charnov, 1977; Wilkins, 1995), similar to, e.g., sex chromosome meiotic drive (Jaenike, 2001; Kozielska et al., 2010). In addition, for these simulations, we did not incorporate fitness effects associated with MY. For each simulation, we determined the most prevalent genotype among females and males at equilibrium, based on the expression and sensitivity levels of their F alleles as well as the number of expressed MY and MA alleles. We categorized simulations according to the activity of a single F allele () relative to and to determine which SD systems can evolve under different levels of F activity. We find that under different relationships between the maximum potential activity of a single F allele () and the SD thresholds, different SD systems can evolve (Supplementary Table 3). When , nearly all simulations yield an SD system where both sexes have two sensitive and expressed F alleles, and males are heterozygous for a single M (either MY or MA; females F/F; +/+, males F/F; +/M, with + indicating the absence of an expressed M allele). In contrast, when or , we additionally encounter systems in which F alleles evolved to become insensitive and/or unexpressed, but the distribution of these alleles over the sexes differs between these two scenarios: when , females carry two insensitive F alleles, whereas males carry a single insensitive F and one unexpressed F; here, insensitive F alleles can be regarded as recessive female-determining alleles, whereas the sensitive and expressed F allele (in presence of M) or unexpressed F plays the role of a dominant male-determining allele. When , females may carry a single insensitive F allele, whereas males carry none, suggesting that insensitive F alleles are dominant female-determining alleles. This is corroborated by the presence of expressed M alleles in both sexes in the simulations, which would otherwise induce maleness in their carriers. In simulations where insensitive F alleles have evolved, we find that the remaining F alleles have become unexpressed, and in some cases that expressed M alleles are lost as well. Here, insensitive F alleles spread first, along with an increased frequency of expressed M alleles in both sexes (Supplementary Figure 1), so that an equal sex ratio is maintained. Next, sensitive F alleles become unexpressed (and subsequently insensitive due to the constant mutation pressure affecting F sensitivity); the insensitive expressed F allele is retained as it now performs the SD function. Along with the increase in unexpressed F alleles, M alleles similarly become unexpressed. Based on these dynamics, we infer that evolution of F insensitivity indirectly renders the loss of expression selectively neutral for sensitive F alleles, which in turn renders the loss of M expression neutral; both genes may therefore decay via mutation accumulation in a stepwise order (Supplementary Figure 2). In addition to the systems identified in our simulations, we speculate some other systems may also be evolutionary stable (Figure 2). Their absence from our simulations may be because they only rarely arise through evolution, or because they represent intermediate states between some of the systems that are observed. The latter for example applies to systems where females have one insensitive and one sensitive F allele, males have two sensitive F alleles, and M is fixed in both sexes (females FI/F; M/M, males F/F; M/M in Figure 2). This system is prevalent in some simulations prior to the accumulation of unexpressed F alleles and loss of M (e.g., Supplementary Figure 1). In absence of fitness effects related to MY, recurrent de novo mutation of MA results in MA replacing MY as the male-determining factor; without evolution of F, this represents a transition from one male heterogamety system to another, as MY and MA perform equivalent functions. Male heterogamety systems with MY or MA as the male-determining gene are observed regardless of the activity of F relative to the SD thresholds (Supplementary Table 3). Moreover, male heterogamety with MY or MA as a dominant male-determiner is observed in all simulations when , suggesting this is the only stable SD system under these conditions, so that insensitive F alleles cannot be part of any SD mechanism under these conditions.

The above results show how F, M, and the SD system as a whole evolve in relation to the SD thresholds and . Most importantly, we find that F can evolve into a dominant female-determining gene by becoming insensitive to M, provided that a single F allele generates enough F product to induce feminization, i.e., the threshold for feminization is sufficiently low relative to F expression . We refer to such insensitive dominant feminizing F alleles as FI. Invasion of FI alleles in our model is primarily driven by the de novo evolution of MA, which causes a slight male-biased sex ratio, thereby causing feminizing mutations to be favored (but see Discussion) (Bull & Charnov, 1977). Provided that the feminization threshold is constant under all conditions, local variation in temperature could then lead to divergent evolution of SD systems as temperature effects on F expression could then allow for local evolution of a female heterogamety system. To test this reasoning, we used simulations in which we varied two parameters: the rate at which temperature increases F expression (see Equation 1) and the threshold value for feminization. Here, we exclude Y-chromosomal fitness effects, but explore their impact on SD evolution in the following section. We found that FI can spread in the entire metapopulation when the feminization threshold is sufficiently small (Figure 3A), regardless of how strongly temperature affects F expression. For sufficiently high values of , FI was unable to spread in colder demes because it would result in intersexual development (Supplementary Figure 3), but could still spread in warmer demes. Under these conditions, a geographical cline in the frequency of FI evolves (Figure 3B). These results underline that the distribution of FI is constrained by the expression level of F, and show that temperature-dependent effects on gene expression can establish gradients when temperature varies. Due to its feminizing effect, an FI allele is transmitted as if it were a female-limited W-chromosome in a ZW system, wherein males have a ZZ genotype, whereas females have a ZW genotype (in contrast to XY systems where females are XX and males are XY). As a result, it occurs only in FI/F females (or non-reproducing intersexes). In the presence of M product, the FI product is not broken down but the product of regular sensitive F alleles is degraded, so that regular F alleles do not contribute to the total F product. This scenario becomes increasingly probable as either MY and/or MA frequency increases, as more FI-bearing individuals will also bear MY and/or MA. Therefore, feminization of developing embryos under these conditions is achieved solely by the activity of the FI allele. Because FI is insensitive to M, the total F activity is determined by its expression (see Equation 1). When no temperature-dependent expression occurs, feminization is only achieved when the baseline genetic expression level exceeds the feminization threshold, but in the presence of temperature effects is less constrained. Therefore, when is sufficiently low, FI can spread everywhere independent of temperature (Figure 3A, left panel), whereas otherwise FI frequencies depend on the temperature-dependent expression rate (Figure 3A, right panel) and the local temperature (Figure 3B).

For the simulations discussed above, we assumed that MY is not associated with any fitness effects. Under these conditions, MY was always lost due to de novo evolution of MA, which causes a male-biased sex ratio so that both MY and MA are selected against. Recurrent mutation of MA ensures that it nonetheless remains present in the population, but this does not hold for MY; in absence of any associated selective effect, MY is therefore purged. However, as MY is the ancestral SD gene, it may have induced the chromosome on which it is located to undergo Y-chromosome evolution (reviewed in Bachtrog, 2013; Schenkel & Beukeboom, 2016). If so, this chromosome is expected to become enriched for sexually antagonistic (SA) genetic variants as well as recessive deleterious mutations. In effect, this would cause MY to be associated with higher fitness in heterozygous +/MY males, but to induce a fitness cost to MY-bearing females (who also carry FI) as well as all homozygous MY/MY individuals. Both sexually antagonistic genetic variation and a cost of homozygosity can in theory strongly affect evolutionary transitions in SD (reviewed in van Doorn, 2014). We therefore performed additional simulations where we considered a sexually antagonistic fitness effect of the Y-chromosome (see Supplementary Table 1), causing the Y-chromosome to affect individual fitness during the mating phase or during development (positively in males, negatively in females). In addition, we include costs of MY/MY homozygosity as a form of viability selection during embryogenesis. When MY is under sexually antagonistic selection, we find that it is maintained over MA (Figure 4). Sexually antagonistic selection on MY can also inhibit invasion of FI if selection is sufficiently strong, so that negative effects of MY in females reduces the fitness of FI/F females relative to non-MY-carrying F/F females. However, when the rate of introduction of novel MA alleles is sufficiently high, FI can always invade even if the selective effects associated with MY are strong. When MA originates more frequently and therefore reaches higher frequencies, a more male-biased sex ratio results (similar to Y-chromosomal meiotic drive; Jaenike, 2001; Kozielska et al., 2010) and hence the selective benefit for FI as a female-determining factor increases. Y-chromosomal fitness effects can also enable maintenance of both MY and MA in the population, albeit in different locations. This occurs when MY is disfavored through reduced survival of YY homozygotes in subpopulations with FI, in contrast to subpopulations without FI, where MY is favored over MA via sexually antagonistic selection in +/MY heterozygotes. Such MY versus MA polymorphism is highly similar to the distribution of Y-chromosomal versus autosomal M-factors in the housefly M. domestica. In this species, autosomal M-factors are more prevalent at lower latitudes and coincide with a dominant feminizing allele traD (which is equivalent to FI as described above), resulting in three latitudinal gradients in Y-chromosomal M-factors, autosomal M-factors, and the traD allele. We find that Y-chromosomal fitness effects enable the evolution of this complex system in our model (see Box 1). This provides an adaptive explanation for the evolution of this system in contrast to existing models of SD evolution, which have been unable to predict when stable polymorphisms for tra versus traD along with Y-chromosomal M-factors versus autosomal M-factors may occur in general, and in particular when these coincide along environmental clines as observed under natural conditions (e.g., Kozielska et al., 2006, 2010, but see Meisel et al., 2016, for an alternative explanation that invokes linkage between different SD genes and sexually antagonistic loci).

Our model has been inspired by the common housefly Musca domestica, in which wild populations harbor different SD systems (reviewed in Hamm et al., 2015). Here, we discuss how our model can explain the evolution of this system. In houseflies, sex is determined by a linear cascade of genes. First, transformer (tra) induces female development when active in developing embryos (Hediger et al., 2010). Its activity depends on an autocatalytic feedback loop where female-specific TRAF protein directs the splicing of Tra pre-mRNA into the female-specific TraF mRNA, thereby ensuring the production of novel TRAF protein to maintain the loop. Possibly, temperature effects on tra work by affecting pre-mRNA splicing, as splicing is sensitive to temperature as well as other stressors (Palusa et al., 2007). Second, masculinizing factors (M-factors) such as Mdmd (Sharma et al., 2017) trigger male development by inhibition of the tra loop through splicing of tra pre-mRNA to a male-specific variant. Intriguingly, the M-factors in M. domestica are found on different chromosomes in different populations; most of these correspond to orthologs of Mdmd, and the genomic regions harboring these M-factors consist of multiple (partial) repeats of the Mdmd gene, suggesting that these have originated through gene duplication and/or translocation (Sharma et al., 2017, Li et al., in prep). There also exists an insensitive variant of tra, traD, that induces female development in all carriers regardless of whether they carry any M-factors. In our model, F represents tra and likewise traD corresponds to the dominant FI as discussed in the main text; MY and MA in our model represent M-factors. High-latitude M. domestica populations have a male heterogamety (XY) system in which the Y-chromosomal Mdmd gene induces maleness, and all individuals carry two regularly sensitive tra alleles. At lower latitudes, however, females usually carry the insensitive traD allele, and both sexes can be homozygous for an autosomal copy of Mdmd; hence, these populations have a female heterogamety (ZW) system (Supplementary Figure 4; Hamm et al., 2015; Kozielska et al., 2008). The geographical transition between these two SD systems is gradual, so that clines exist in the frequencies of Y-chromosomal Mdmd (decreasing toward lower latitudes), autosomal Mdmd and traD (both increasing toward lower latitudes). Given that a single traD allele is sufficient for feminization (even in the presence of multiple M-factors; Hamm et al., 2015) suggests that at least in natural populations it would mimic an FI allele with a net expression level that exceeds . Therefore, in our model this system would likely mimic a transition from to . Temperature likely plays a causal role in maintaining these gradients by affecting the SD process (Adhikari et al., 2021; Delclos et al., 2021; Feldmeyer, 2009; Feldmeyer et al., 2008). Temperature effects on housefly SD have been reported in the form of biased sex ratios produced in wildtype crosses (Feldmeyer, 2009) as well as in females carrying the masculinizer (man) mutation, another variant of tra (Hediger et al., 2010; Schmidt et al., 1997). The man mutation represents a maternal-effect male-determining gene, where man-carrying females can produce all-male progeny even if the progeny do not carry an M-factor. However, this effect is incomplete and temperature-sensitive (Schmidt et al., 1997), with offspring sex ratios more male-biased at higher temperatures. Altogether, temperature seems to have an important influence on SD in M. domestica, but the underlying mechanisms are not yet fully understood. The housefly with its different SD systems is represented in our model by gradients in the frequencies of MY (decreasing with temperature), MA, and FI (both increasing with temperature). Presumably, traD is limited to warmer localities for similar reasons as FI in our model, i.e., because a single traD allele may not be sufficient to induce feminization at low temperatures. Y-chromosomal Mdmd and autosomal Mdmd may follow similar dynamics as MY and MA. Genomic analyses of the Mdmd loci on the Y-chromosome and various other chromosomes suggests that Mdmd may exhibit (or have exhibited) transposon-like activity, leading to tandem duplications of transposition of these duplications as a complex to other chromosomes (Sharma et al., 2017, Li et al., in prep). Y-chromosomal fitness effects can yield gradients in MY and MA, but may also prevent the spread of FI, particularly when novel MA alleles enter the population at a low rate. The evolution of a housefly-like polymorphic SD system therefore depends on a balance between the Y-chromosomal fitness effects and the rate at which new autosomal Mdmd arises. In our model, we find that a housefly-like system can evolve under various rates of MA de novo evolution (Figure 5). Higher rates of MA evolution require stronger SA effects for MY to be maintained in low-temperature demes. Costs of YY homozygosity do not appear essential for MY to be lost in the presence of FI, although they may increase the likelihood of MY being lost in favor of MA in the presence of FI by reducing the fitness of YY homozygotes. Possibly, costs of MY in females due to sexually antagonistic selection may suffice to promote the loss of MY. Our model therefore can explain the evolution of complex SD systems such as those found in houseflies. Thus, we propose a novel hypothesis for the evolution of the housefly system in which the sex determination cascade is shaped by a combination of environmental influences on tra, recurrent evolution of autosomal Mdmd, and fitness effects associated with the Y-chromosomal copy of Mdmd (Figure 6). This hypothesis, along with the results in this manuscript, highlights several promising aspects of housefly SD and its evolution that require further investigation. First and foremost, the effect of temperature on tra functionality needs to be verified, as current data are inconclusive (Adhikari et al., 2021; Feldmeyer, 2009). Our model predicts that tra function is enhanced at high temperatures, allowing for the evolution of a dominant feminizing allele traD. Given that tra regulates its own activity via an autoregulatory feedback loop (Dübendorfer et al., 2002; Hediger et al., 2010), this enhanced function may be effectuated in a variety of manners, including pre-mRNA expression, mRNA splicing, and protein catalytic activity, so that a purely molecular analysis of tra function may be inconclusive. These issues are aggravated by the fact that whether or not the loop is established is typically determined in the zygotic stage, so that analyses of tra function in the adult stage may not be informative (in contrast to Adhikari et al., 2021). Given that tra is expressed even in presence of Mdmd, it seems likely that the hypothesized gain-of-function of traD must act post-transcriptionally. Instead, an inverse approach may be more fruitful; if traD can only evolve at high temperatures, then traD-bearing flies reared at low temperatures should be more prone to sex reversal or intersexuality. Second, the evolution of Mdmd must be further explored to evaluate the transpositional activity of Mdmd itself or the supergene complex containing this gene. The complex arrangement of Mdmd copies in M-loci strongly suggests Mdmd may at one point have transposed frequently (Li et al., in prep), but it is undetermined whether this is still the case. Third and final, fitness assays of houseflies with Y-chromosomal versus autosomal Mdmd must be made to determine if there are any fitness effects linked to these genes. In vitro studies on houseflies strains that carry Mdmd on different chromosomes found that the Mdmd-bearing chromosome may exhibit some genetic differentiation (Meisel et al., 2020), and that these may be associated with fitness effects (Son et al., 2019). In vivo characterization is however required to validate these effects, and to determine their impact on the maintenance of MY over MA as predicted in our model.

We have presented a model for the evolution of SD systems in a context where sex is determined by genetic factors in combination with environmental effects. In our model, female development is induced when the activity of a feminizing gene F exceeds a certain threshold, whereas male development is achieved when F activity is below a different and lower threshold. This can be caused by inhibition of F activity by the maleness-promoting gene(s) MY and/or MA, or by loss of F expression. We incorporated a positive effect of temperature on the expression of a feminizing locus F. We find that several different SD systems can be realized depending on the activity of an F allele relative to the threshold values for masculinization and feminization. Temperature-dependent effects on F expression may alter the relationships between F activity and SD thresholds, thereby enabling the evolution of different genetic SD systems. A particular prediction is the transition from male heterogamety to female heterogamety, which occurs in our model via the evolution of an insensitive dominant feminizing variant (FI) that induces femaleness even in the presence of MY and/or MA. FI can spread when activity of a single F or FI allele is sufficient to induce feminization; when activity is modulated by temperature, this can lead to local invasions of FI and subsequently differentiation between populations along temperature gradients. Differentiation can also occur for MY and MA, with MY being favored in absence of FI and MA in presence of FI. This occurs when MY is associated with certain fitness effects such due to linkage with sexually antagonistic variants or recessive deleterious mutations. Altogether, we show that this can lead to the coexistence of multiple gradients in SD genes as found in, e.g., the housefly M. domestica. Our model can be amended to other SD systems than the M. domestica system on which it is based, provided that they have a basic GSD framework influenced by an environmental effect (Table 1). Environmental effects on genetic sex determination systems are being reported in an increasing number of species (Edmands, 2021; Holleley et al., 2015). Although temperature-dependent effects are well-documented, other environmental effects may also influence SD in certain systems such as hormonal imbalance in fish due to pollution (Devlin & Nagahama, 2002). Other components and assumptions of the model that are based on the housefly system may be represented differently in other species but with similar effects. For example, the impact of MA evolution is not due to a specific mechanism of mutation, but more generally by causing a male-biased sex ratio, thereby promoting the invasion of FI. Male-biased sex ratios occur due to MA overrepresentation in the gene pool via its de novo evolution, but the same effect can be achieved by translocation of a Y-chromosomal male-determining gene or via association with meiotic drivers (Green, 1980; Kozielska et al., 2010). Similarly, although FI is favored due to sex ratio selection in our model, other selective processes such as sexually antagonistic selection (van Doorn & Kirkpatrick, 2010) as well as neutral processes such as genetic drift (Saunders et al., 2018; Veller et al., 2017) can also drive transitions from male to female heterogamety. Inversely, the spread of FI in our model may be impeded by selection against intersexes, but other factors may also limit its spread. For example, when FI is linked to a gene experiencing intralocus sexual conflict in one environment but not another (Plesnar-Bielak & Łukasiewicz, 2021; Schenkel et al., 2018), it may be favored in the first environment due to sex-specific benefits to females; if such benefits are absent in the second type, and instead it is similarly costly to males and females alike, this would restrict the spread of FI to the first environment. Additionally, we see that in absence of an association between MY and fitness effects, MA replaces MY altogether due to recurrent mutation at MA but not MY, yet still drives the invasion of FI, showing that our model does not strictly require a third locus. Inversely, it is likely that a more complex genetic basis generates similar evolutionary patterns, such as when various genes can evolve into a male-determining gene (Bopp, 2010). Bopp (2010) describes how a mutation that interferes with the key switch transformer, which regulates female development in houseflies and several other species (see also Box 1), is more likely to occur than a mutation that enhances the function of transformer, simply because there are many ways in which a gene’s function can be disrupted, but only a single way for it to function as intended. In our model, this mechanism may apply when F function can be disrupted in a variety of ways (though a similar reasoning could be applied to a gene with male determining function). In this scenario, many genes having a small chance to evolve into a male-determining gene may have the same net effect on sex ratio as a gene that is prone to evolving a masculinizing function. In this light, it will be interesting to test whether genes that have acquired sex-determining function in one species are prone to evolve a similar function in a related species, where it has no sex-determining function. Previous work has shed light on the evolution of ESD and GSD systems, and when transitions between these two may arise (Muralidhar & Veller, 2018; Pen et al., 2010; Schwanz et al., 2020). However, our understanding of the evolution of polymorphic SD systems and the potential for environmental heterogeneity to influence this process remains limited. Our results highlight environmental effects on GSD systems, and under which conditions this can lead to within-species polymorphism in GSD systems. Alternatively, these mechanisms may also explain why in some groups GSD mechanisms may be highly divergent between closely related species, such as in the mosquitofish Gambusia affinis and the closely related G. holbrooki (Kottler et al., 2020). In such cases, different species experience different ecological conditions in their respective niches, among which environmental factors that might impinge on the SD process as described here. Between-species divergence in SD mechanisms may then occur in much the same way as occurs in our model between populations. Our results furthermore suggest that hybridization need not be an impediment to such differentiation; potentially, such diversification may even enable the evolution of reproductive isolation via Haldane’s rule (Haldane, 1922). In conclusion, even in systems that appear to be fully GSD, the role of environmental influences on the SD processes must not be ignored as these may have played an important role in their evolution. In extension of this, changes in environment, e.g., due to global warming, may impose yet unforeseen selective pressures on species with GSD systems.

Click here for additional data file.