TEMPERATURE‐DEPENDENT TURNOVERS IN SEX‐DETERMINATION MECHANISMS: A QUANTITATIVE MODEL

Sex determination is often seen as a dichotomous process: individual sex is assumed to be determined either by genetic (genotypic sex determination, GSD) or by environmental factors (environmental sex determination, ESD), most often temperature (temperature sex determination, TSD). We endorse an alternative view, which sees GSD and TSD as the ends of a continuum. Both effects interact a priori, because temperature can affect gene expression at any step along the sex‐determination cascade. We propose to define sex‐determination systems at the population‐ (rather than individual) level, via the proportion of variance in phenotypic sex stemming from genetic versus environmental factors, and we formalize this concept in a quantitative‐genetics framework. Sex is seen as a threshold trait underlain by a liability factor, and reaction norms allow modeling interactions between genotypic and temperature effects (seen as the necessary consequences of thermodynamic constraints on the underlying physiological processes). As this formalization shows, temperature changes (due to e.g., climatic changes or range expansions) are expected to provoke turnovers in sex‐ determination mechanisms, by inducing large‐scale sex reversal and thereby sex‐ratio selection for alternative sex‐determining genes. The frequency of turnovers and prevalence of homomorphic sex chromosomes in cold‐blooded vertebrates might thus directly relate to the temperature dependence in sex‐determination mechanisms.

Birds and mammals display a strictly genotypic sex determination (GSD), with highly differentiated sex chromosomes (Graves 2008). The XX/XY male-heterogametic system of mammals has remained stable since the master male-determining gene SRY first appeared close to 200 million years ago (Potrzebowski et al. 2008;Veyrunes et al. 2008). Loss of recombination has since induced strong differentiation and degeneration of the Y chromosome. A similar process occurred in birds, where females are the heterogametic sex and thus carry the decayed W chromosome (ZW/ZZ system).
Patterns are strikingly different in other vertebrates. First, sex determination is often extremely labile. Different sex-determination systems may be found in closely related taxa (Hillis and Green 1990;Ezaz et al. 2006; Baroiller et al. 2009) sometimes even in different populations from the same species (Miura et al. 1998). Transitions between systems seem frequent on an evolutionary time scale (Hillis and Green 1990;Janzen and Krenz 2004;Ezaz et al. 2006;Mank et al. 2006;Graves 2008) with the result that sex chromosomes (here defined as chromosomes bearing sexdetermining genes, independent of recombination patterns), are usually homomorphic or poorly differentiated (e.g., Schmid 1983;Hayes 1998;Devlin and Nagahama 2002;Eggert 2004).
Second, the environment also plays a role, temperature in particular. Sex may be determined mostly by temperature (such as found in turtles or alligators) or by the joint action of temperature and genetic factors (Bull 1980;Conover and Heins 1987;Conover et al. 1992;Janzen and Phillips 2006). There is mounting evidence for mixed sex determination systems (where GSD may be overridden by temperature within the natural range) in lizards (Shine et al. 2002;Quinn et al. 2007;Radder et al. 2007) and fish (reviewed in Ospina-Álvarez et al. 2008; Baroiller et al. 2009), including species with differentiated sex chromosomes. Sex reversal also occurs spontaneously in species considered to have purely GSD (e.g., Crew 1921;Witschi 1929b;Aida 1936;Kawamura 1977;Nagler et al. 2001;Matsuba et al. 2008). Even when sex determination is strictly genotypic under natural conditions, temperature has been shown to play a role in experimental conditions (Witschi 1929a;Wallace and Wallace 2000;Eggert 2004;Ospina-Álvarez et al. 2008). Temperature effects are probably not adaptive in such cases, but more likely the side effect of a temperature-dependence in underlying molecular processes.
Two highly contrasted views exist regarding interactions between, respectively, temperature sex determination (TSD) and GSD. The conventional view sees sex determination as a dichotomous process. As formulated by Valenzuela et al. (2003), the sex of organisms is determined by two distinct and mutually exclusive mechanisms. In GSD, sex is determined at conception by genes, whereas in environmental sex determination (ESD, as for instance TSD), sex is determined after fertilization by environmental factors. To account for intermediate situations (where both genes and environment interact), Valenzuela et al. (2003) propose that (1) some species can be categorized as GSD + EE: that is , sex is determined at conception by genes, but may be then reversed by environment during the embryonic development; (2) TSD species may harbor genetic variance in their sensitivity to environment, but are still to be considered as TSD, because at the individual-level sex is determined by environment; (3) even when TSD and GSD coexist within species or populations, at the individual level, sex is still determined either by genes or by the environment (Valenzuela et al. 2003).
The alternative view, by contrast, sees GSD and TSD as the two ends of a continuum (e.g., Sarre et al. 2004). The observed continuity in phenotypic patterns of sex determination reflects the extraordinary conservatism in the gonadal developmental pathways of vertebrates, as revealed by recent molecular studies. In both TSD and GSD species, the same genes are involved in the cascade of processes that result in sex determination (Sarre et al. 2004). At each stage along this path, temperature may affect molecular processes, and thereby the final outcome. Temperature is known, for instance, to affect the activity of enzymes (such as aromatase, which transforms testosterone into estradiol; e.g., Desvages et al. 1993;Crews and Bergeron 1994;Crews et al. 2001;D'Cotta et al. 2001;Lance 2009) or the expres-sion of genes (such as DMRT1) involved in the sex-determining cascade.
This may be illustrated by the case of medaka fish (Oryzias latipes). New sex chromosomes recently evolved in this lineage through the duplication of DMRT1 (Matsuda et al. 2002;Nanda et al. 2002), a gene encoding a transcription factor with a DMdomain playing a crucial role in the male determination cascade throughout all vertebrates (Hodgkin 2002;Haag and Doty 2005), including both GSD Ounap et al. 2004;Yoshimoto et al. 2008;Smith et al. 2009) and TSD species (e.g., Kettlewell et al. 2000;Shoemaker et al. 2007). At normal temperatures (25 • C), DMRT1 expression from the autosomal copy is too low to reach the threshold required to induce male development, so that XX individuals develop into females. By contrast, XY individuals develop into males due to the additional expression of the duplicated copy on the proto Y (DMRT1bY). Higher temperatures (32 • C), however, upregulate the autosomal DMRT1 copy, so that a significant fraction of XX individuals develop into males, which are perfectly functional and phenotypically indistinguishable from XY males (e.g., Sato et al. 2005;Hattori et al. 2007). Sex is thus determined by the amount of DMRT1 transcription factor, which itself depends on a combination of genetic and environmental effects. Interestingly, DMRT1bY is first expressed at the neurula stage (stages 17-18; Nanda et al. 2002), which coincides temporally with the thermally sensitive period of masculinization (Hattori et al. 2007). Hence, opposing the conventional definition (Valenzuela et al. 2003), the timing of sex determination does not necessarily differ between TSD and GSD.
In a quantitative-genetics perspective, sex can be seen as a threshold trait that depends on an underlying liability factor (e.g., DMRT1 expression in medaka), which is under both genetic and environmental influence (e.g., Bull 1981;Bulmer and Bull 1982;Roff 1996). Under this concept, an individual trait value cannot be assigned to either genes or environment. Why does an XX medaka fish become male at 32 • C? Sex determination is neither purely genetic (as temperature plays an obvious role), nor purely environmental (as with different autosomal DMRT1 alleles, the individual might have developed into a female) but results from an interaction between the two factors (temperature dependence of DMRT1 expression). The relevant question in quantitative genetics is that of the apportionment of phenotypic variance within populations, and this directly relates to the definition of sex-determination systems. In line with the alternative view, we propose to define a system as GSD if all (or most) of the variance in sex determination within the normal environmental range can be assigned to genetic factors, and ESD if all (or most) of this variance can be assigned to environmental factors. These are to be seen as special cases of a more general situation, where phenotypic variance has both genetic and environmental components. This definition obviously differs from the one proposed by the conventional view. At the XX pivotal temperature (T XX ), the system is purely TSD ( A XX,T = 0) with even sex ratio. At higher temperatures, sex-ratio selection will favor a feminizing mutation (W), leading to a femaleheterogametic system ( A XW,T < 0 and A XX,T > 0). Similar transitions are expected for temperature decreases.
An appealing way to formalize gene-environment interactions in this unifying view is to represent genotypes as reaction norms in the space defined by phenotype versus environment (Fig. 1). In any given environment, different genotypes may express large differences (major genes) or small differences (minor genes) in liability-trait values. Environmental effects define the shape of the norm (e.g., thermal up-regulation of DMRT1 expression) which may be linear or not, and may or not vary among genotypes (G × E interactions). Environmental effects are to be seen as constraints stemming from the laws of thermodynamics. Thermal upregulation of aromatase inhibitors in tilapia (D'Cotta et al. 2001), or DMRT1 expression in medaka (Hattori et al. 2007), for instance, are seen as the necessary consequences of increased kinetic energy in the underlying physico-chemical processes. In GSD species, this effect of temperature on the liability trait might be hidden by genes with major effects (i.e., the liability-trait values expressed by different genotypes are so far apart from the threshold that sex is under complete genetic control within the natural range). TSD species, by contrast, have capitalized on this constraint to make temperature dependence a functional system. This evolution clearly required a set of specific co-adaptations (regarding e.g., nest-choice behavior or the shape of reaction norms) to render this strategy adaptive. The evolution of TSD is certainly a fascinating question, but one outside the scope of the present article.
The first aim of our formalization effort is to show that the alternative paradigm, in which sex determination depends on both genes and environment, directly accounts for the frequent transitions occurring between sex-determination systems. If a population of medaka fish were to live in a consistently warm habitat, where a majority of XX individuals develop into males, then the ensuing sex-ratio bias will necessarily induce a strong selective pressure for a new sex-determining system. This might be achieved, for example, by the appearance of a new allele on the autosomal DMRT1 copy that downregulates its expression. Individuals with this copy would develop into females, so the system would have evolved toward female heterogamety through a sex-chromosome turnover.
The second point we want to make is that this process may help explaining why transitions seem more frequent in coldblooded than in warm-blooded vertebrates. In the formers, the developing embryos are under the direct influence of external temperature, which will necessarily interfere with sex determination owing to underlying thermodynamic constraints. Hence, temperature drifts (induced e.g., by climatic changes or range expansions) will necessarily cause turnovers. In the warm-blooded birds and mammals, by contrast, embryonic temperature is kept constant by parental control during the sensitive period of development, which will prevent sex-reversal and the ensuing sex-ratio selection for turnovers. These groups are thus expected to display evolutionary stable GSD.

PHENOTYPIC DISTRIBUTION OF THE LIABILITY TRAIT
Let us consider one master sex-determining locus and two alleles X and Y with additive effects, X having a feminizing effect and Y a masculinizing effect (i.e., Y produces more of the liability trait A than does X). Depending on the strengths of these effects relative to the threshold value, the model may account for either male heterogamety (XX females, XY males) or female heterogamety (XY females, YY males). In both cases the two genotypes involved constitute a recurrent pair (sensu Bull 1983), that is, two genotypes of opposite sex, which when mated produce only the same two parental genotypes, in a 1:1 ratio.
Masculinization or feminization effects also depend on temperature. Each genotype IJ (I, J = X, Y) is characterized by a norm of reaction, defining its liability-trait value, A IJ,T (Table 1), as a function of temperature T. Such a norm can be represented as a curve in the space phenotype-environment, assumed here to be linear with slope β (Fig. 1). The genetic component of the phenotypic variance in A within populations comes from the coexistence of different genotypes, producing different A values at any given temperature. The environmental component stems from the fact that different individuals with the same genotype

A IJ,T
Liability-trait value for genotype IJ (I, J=X, Y) at temperature T ζ threshold value for A, such that individuals develop into males for A>ζ and into females otherwise proportion of males produced by genotype IJ at temperature T R IJ,T =r IJ,T / number of males per female for genotype (1−r IJ,T ) IJ at temperature T p IJ,k frequency of individuals with genotype IJ and sex k (=m, f ) within the population β slope of reaction norms (increase in the production of the liability factor A with a unit increase in temperature) T IJ Pivotal temperature for genotype IJ (defined by R IJ,T = 1).
experience slightly different microenvironments, and thus have different values of the liability trait ( Fig. 1). Assuming normal distributions for environmental deviations from genotypic averages (with standard deviation σ E ), the phenotypic distribution of the liability trait A within populations, at any average temperature T, is the sum of several normal distributions (one per genotype) weighted by genotypic frequencies.

RATIOS
Sex is considered as a threshold trait, that is, there is a threshold value (ζ) for the liability trait A, such that for each genotype IJ, a proportion r IJ,T exceeds the threshold and develops into males, whereas the complementary proportion 1 − r IJ,T develops into females. Pure GSD will result if, at given T, liability-trait values (A IJ,T ) are far apart from the threshold, and σ E is low. In such a case, phenotypic sex perfectly correlates with genotypes. In contrast, pure TSD will result if a single genotype is fixed, and phenotypic variance in A only results from environmental variation. In between these two extreme cases, the correlation between phenotypic sex and genotypes will be imperfect, with a variable amount of sex-reversed individuals. The relevant quantity to determine the amount of sex reversal for genotype IJ at temperature T is the distance of its liability value to the threshold ζ, in units of σ E . Hence it is useful to standardize these values as α IJ,T = AIJ,T −ζ σE . Using this metrics, the slope of the reaction norm is expressed in units of standard errors (β/σ E ) and the proportion of males among individuals of genotype IJ becomes r IJ, x z=0 e −z 2 dz is the so-called error function. This sex ratio can thus be written where α represents the standardized amount of liability trait (α = A−ζ σE ). It can be easily checked from equation (1) that genotype IJ will produce both sexes in equal proportions at the temperature for which its liability-trait value exactly matches the threshold (i.e., A IJ,T = ζ and α IJ,T = 0), referred to as its pivotal temperature (T IJ ).

PARAMETER ESTIMATION
Equation (1) can be used to predict the sex ratio produced by a given genotype at a given temperature. Reciprocally, it can also be used to infer α IJ,T values from observed sex ratios. For illustration ( Fig. 2), we calculated α IJ,T values at different temperatures for four TSD lineages, namely the fish Odonthestes bonariensis (Strussmann et al. 1997) and Poeciliopsis lucida (TSD strain; Sullivan and Schultz 1986), the lizard Niveoscincus ocellatus (Wapstra et al. 2009) and the turtle Chrysemys picta (Ewert et al. 2004). The excellent fit supports the assumption of a single underlying genotype in all four cases (i.e., pure TSD). Reaction norms appear linear, with positive slopes in fish (β/σ E = 0.51 • C −1 and 0.28 • C −1 respectively) and negative slopes in reptiles (β/σ E = −0.40 • C −1 and −1.38 • C −1 ). The pivotal temperatures (i.e., temperatures providing equal sex ratios), obtained as the intercept of the reaction norms with the abscissa (α IJ,T = 0), amount to 24.3 • C, 25.3 • C, 17.7 • C and 26.9 • C, respectively.
From the normal-distribution assumption, sex ratio is expected to be a sigmoid function of temperature in pure TSD systems such as in O. bonariensis (Fig. 3A). In heterogametic species with mixed sex determination, by contrast, the two genotypes are expected to generate a double sigmoid, as illustrated for two male-heterogametic species in Fig. 3B, C, namely the newt Triturus cristatus (Wallace et al. 1999) and the fish Patagonina hatcheri (Strussmann et al. 1997). In both cases, data were fitted assuming linear and parallel norms of reaction. In P. hatcheri, the XX and XY genotypes are more differentiated on the α axis (21.2, as opposed to 6.4 in T.cristatus) but display closer pivotal temperatures (25.1 • C and 15.4 • C respectively, as opposed to 31.2 • C and 8.4 • C for T. cristatus), because the slope (in terms of β/σ E ) is much steeper (2.17 • C −1 vs. 0.28 • C −1 for T. cristatus).
Intersex individuals sometimes appear around the pivotal temperature, due to disorders of sexual development arising when the liability-trait value lies too close to the threshold. Their frequency allows estimating an intersex range along the α axis. In the female-heterogametic salamander Pleurodeles waltl, for instance, the ZW genotype produces at 30 • C 44% females, 46% males and 10% intersex individuals (Dournon et al. 2007, Table 3). Thus, assuming a normal distribution for the liability trait, intersex occurs in this species for α values ranging −0.13 to + 0.13.

Evolutionary Responses to Temperature Change
All of the systems illustrated here are expected to show evolutionary response to temperature change, stemming, for example, from climatic change or range expansion. The short-term response will be a sex-ratio bias: for example, rearing medaka fish at 32 • C will produce an excess of males in the first generation, because some XX individuals develop into males. Within a few generations, however, sex-ratio selection will restore a Fisherian sex ratio by lowering the frequency of the Y chromosome. In Appendix A, we derive the equilibrium frequency values for all three genotypes involved (XX, XY, and YY) and for the Y chromosome as a function of temperature and σ 2 E . Specific values are also derived at the several pivotal temperatures. In the following, we discuss in more intuitive terms the evolutionary outcomes predicted from consistent increases or decreases in temperature.

TEMPERATURE INCREASE
Let us start with an initial situation (T = 0) where the liability-trait values for the XX and XY genotypes define a male-heterogametic system. Due to sex-ratio selection, the two genotypes are in equal proportions and Y frequency approximates 0.25 (Fig. 4A). As temperature increases, a few XX individuals turn into males, inducing a mixed system with genetic and environmental compo-nents. Sex-ratio selection will automatically lower the frequencies of XY individuals (and thus of Y) so as to maintain an even sex ratio at the population level. Note, however, that sex ratios are biased within families, because mating between XY males and XX females produces an excess of males, whereas mating between XX males and XX females produces an excess of females.
At the XX pivotal temperature (e.g., 25.1 • C in P. hatcheri or 31.2 • C in T. cristatus, Fig. 3), half of the XX individuals develop into males, so that within a few generations Y is lost and sex determination becomes purely environmental. As temperature further increases, the sex ratio becomes progressively male biased (Fig. 4A) and, were the process to continue, extinction would occur as soon as females become too rare to sustain a viable population. This sex-ratio bias will favor the emergence of a new sex-determination mechanism, selecting for any feminizing mutation, either on the same locus (XW/XX) or on any other locus involved in the sex-determining cascade (XXfF/XXff ). The system will thus evolve toward female heterogamety.

TEMPERATURE DECREASE
If, from the initial conditions (T = 0), temperature is decreasing rather than increasing, then at some point a few XY genotypes will develop into females. Mating with XY males will produce 25% of YY progeny, which may be viable or not, depending on the load of deleterious mutations that have accumulated on the nonrecombining segment of the Y chromosome.
In the case YY are viable and develop into fertile males, this genotype will increase in frequency, and XX females decrease genotypes develop more likely into females. Sex-ratio selection will thus increase the frequency of YY genotypes, and decrease that of XX genotypes. At mid distance between the XY and YY pivotal temperatures, the system will have reached a new state of pure GSD with female heterogamety. The genotypes YY and XY will produce males and females, respectively, in equal quantities, and the equilibrium frequency of Y tends toward 0.75 (eq. A3, Fig. 4A).
If temperature keeps on decreasing, the system will behave symmetrically to the high-temperature conditions. A pure ESD situation will be reached at the YY pivotal temperature, where allele Y will become fixed in the population (Fig. 4A). Thereafter sex ratio will progressively depart from even (Fig. 4A) so that the population will eventually go extinct due to absence of males, unless a masculinizing mutation, either on the same locus (YY/YZ) or on a different one (YYmm/YYmM), rescues the system (which will then turn back to male heterogamety).
Unviable YY genotypes will induce an earlier female bias (i.e., before the XY pivotal temperature is reached), and a strong segregation load (because crosses between XY males and XY females keep on producing 25% of lethal YY genotypes). Relative to the YY-viable situation, both XY and XX genotypes will reach higher frequencies (because the only males are XY). Extinctions are thus expected to occur earlier, unless the system is rescued by a masculinizing mutation (either XX/XZ or XXmm/XXmM).

ENVIRONMENTAL VARIANCE EFFECTS
The equilibrium Y frequency is a stepwise function of temperature (Fig. 4A), decreasing from 1 at the YY pivotal temperatures to 0 at the XX pivotal temperature. Steps are strongly marked when the environmental component of phenotypic variance is low (i.e., the slope β/σ E is steep) because the standardized liability-trait values lie then far apart from the threshold for most of the temperature ranges, so that sex reversal is rare. As a result, GSD systems (either male or female heterogamety) are stable over large ranges of temperatures, and show abrupt changes around critical temperature values. In contrast, a high environmental variance (i.e., shallow β/σ E slope) makes changes in equilibrium Y frequency smooth and continuous, resulting in a mixed sex-determination system over large temperature changes.

RATIONALE AND SETTINGS
Whether temperature changes result in extinctions or turnovers will depend on drift and mutations, in addition to the selective pressures outlined here above. To cross-validate our analytical predictions and investigate the interplay between these evolutionary forces, we performed individual-based simulations, exploring a subregion of the domain of climatic change plotted in Figure 1.
Simulations (detailed description in Appendix B, see also Fig. S1) were run with a modified version of the program quantiNemo 1.0.3 (Neuenschwander et al. 2008). Starting with a male heterogametic system (XX/XY) at T = 0, temperature was increased by progressive steps after a burn-in of 400 generations. We varied the environmental variance σ 2 E (from 10 −7 to 40) and the population size N (from 50 to 10,000) to investigate their effects on transition probabilities. A first set was run without any mutations to cross-validate our analytical expectations. In the other sets, we allowed mutations, either to a strongly feminizing allele W (set 2), or to randomly sampled allelic values from a large range of masculinizing and feminizing mutations (sets 3 and 4; further details in Appendix B). Mutations affected only the intercept (genetic up-or down-regulation of the liability trait), not the slope (assumed to reflect environmental effects stemming from physiological constraints), but we also performed simulations with a flat slope (β = 0) as a control for endothermy (body temperature assumed constant, independent of environmental change). We also ran simulations assuming YY lethality and/or some intersex sterility. Intersex sterility was implemented as a fitness decrease for genotypes with a liability-trait value close to the threshold (ranging −1.0 to 1.0 Appendix B, Fig. S2), and YY lethality by assigning a fitness of zero for YY genotypes. For each combination of parameter values, we ran 500 replicates over 3,000 generations.

CROSS-VALIDATING ANALYTICAL RESULTS
The results obtained from the first set match remarkably well our analytical expectations (see e.g., Fig. 4B for N = 1,000), notwithstanding slight departures stemming from genetic drift at small population sizes (N ≤ 500) and large environmental variances (σ 2 E > 5; data not shown). At starting conditions (T = 0), the Y frequency slightly exceeds 0.25 for high σ E values (because the few sex-reversed XY females produce some YY males). As expected from the analytical calculations, Y was quickly lost (and ESD achieved) as soon as the XX pivotal temperature was reached. At higher temperatures, increases in the frequency of XX males could not be accommodated anymore by a decrease in Y frequency, so that male biases progressively accrued (Fig. 4B) and extinctions eventually occurred at low σ 2 E values. Both gene dynamics and extinction probabilities were affected by environmental variance σ 2 E in interaction with population size. As expected from analytical results, a large environmental variance (shallow β/σ E slope) increased the proportion of sex-reversed XX males below the pivotal temperature (α XX,T < 0), and thus lowered equilibrium Y frequency (Fig. 4). When genetic drift was large (i.e., small population size N), this selective pressure accelerated the loss of Y, and thus the transition to ESD. Above the pivotal temperature by contrast (α XX,T > 0), a large environmental variance σ 2 E increased the proportion of females among XX individuals, and thereby moderated the biases in sex ratio (Fig. 4). Extinction risk (stemming from such biases) was thus large at small σ 2 E (steep β/σ E slope), but declined rapidly with increasing σ 2 E , and the more so in large populations.

MUTATIONS AND TRANSITIONS TO NEW GSD
The remaining sets of simulations allowed for environmental change to be combined with mutations to new alleles. These sets were used to investigate how sex-ratio selection interacts with drift and environmental variance in determining the timings and probabilities of transitions (or extinctions). Under the assumption of endothermy (β/σ E = 0), the initial XX/XY system was usually maintained throughout (e.g., 63% to 75% of simulations at N = 500, depending on the mutation model). Multiallelic polymorphisms (XYW, akin to the system found e.g., in Xiphophorus Orzack et al. 1980) sometimes occurred (20-35% of simulations at N = 500), while transitions to another recurrent pair were rare (2% to 5% at N = 500, depending on the mutation model).
Under the assumption of ectothermy (β/σ E = 0), by contrast, the initial XX/XY system was always overturned. As all mutation models provided the same qualitative results, we present only the first one in the main text (set 2: mutation to one strongly feminizing allele W; Fig. 5) and the others as Supporting information (Fig. S3). The results can be classified in four qualitatively distinct processes, depending on whether the final outcome was extinction or transition to a new GSD, and whether this end result was preceded or not by an ESD period. (1) Transitions to ESD (by loss of Y) followed by extinction occurred when no feminizing allele was available during the ESD phase. (2) Transitions to ESD, followed by a new GSD (most often female heterogametic, XW/XX) occurred when mutation to a feminizing allele occurred during the ESD phase. (3) Direct transition to a new GSD (most often female heterogametic, XW/XX), occurred when a feminizing allele appeared before the loss of Y, and was maintained throughout by chance to be then favored by sex-ratio selection. In outcomes (2) and (3), transition to a new male heterogametic GSD sometimes occurred (e.g., WW/WY), when Y was maintained at low frequency throughout, or a new masculinizing mutation arose.
(4) In a few cases, finally, a transition to a new GSD was followed by extinction when the feminizing effect of the new allele was insufficient to cope with climatic change.
Outcome frequencies depended on population size N and environmental variance σ 2 E ( Figs. 5A and S3A). Direct transitions (outcome 3) were very rare in small populations, because feminizing mutations were both less likely to appear (Nμ being up to 200 times smaller), and less likely to be maintained (due to strong drift). Environmental variance had no noticeable effect on direct transitions, but a strong effect on extinction rate (outcome 1), in interaction with population sizes: extinctions were very frequent at small N and low σ 2 E . Outcome 4 only occurred when both N and σ 2 E were low (Fig. S3A).

INTERSEX STERILITY AND YY LETHALITY
Intersex sterility had a drastic influence on the outcomes, independent of whether YY was viable (Fig. 5B) or not (data not shown). ESD was rare and only observed for low population size and high σ 2 E (shallow β/σ E slope), whereas direct transitions to a new GSD (always male heterogametic, WW/WY) were frequent outcomes. Outcome 1 (ESD then extinction) was rare and outcome 4 (new GSD followed by extinction) occurred at low population size and low σ 2 E . Surprisingly, intersex sterility generally lowered extinction risks in small populations.
By contrast, YY lethality only had a slight negative effect on the probability of direct transitions. This was due to selection against the WY/YY recurrent pair, which prevented the feminizing allele W to spread before attainment of ESD (data not shown). More serious effects are obviously expected under temperature decrease (see above).

SEX-DETERMINATION IN A QUANTITATIVE-GENETICS FRAMEWORK
The conventional definition of sex-determination systems (e.g., Valenzuela et al. 2003;Ospina-Álvarez et al. 2008) is undoubtedly applicable in the extreme cases of pure GSD or pure ESD, when all the phenotypic variance of the liability trait in a population is either genetic or environmental: If for example, all individuals in a population share the same genotype, the reason why a focal individual develops into a male or a female can easily be assigned to its environment. In intermediate cases, however, such a dichotomous assignment is not possible, because genetic and environmental effects cannot be partitioned at the individual level: the liability-trait value necessarily depends on the interplay between genes and environment.
What can be partitioned, however, is the amount of variance in a population. The relevant question thus becomes that of the apportionment of phenotypic-sex variance into genetic and environmental factors. We think the alternative definition of sexdetermination systems we propose here, and its formalization in a quantitative-genetics framework, offer better opportunities to account not only for intermediate situations, but also for transitions between sex-determination systems. Both seem common in cold-blooded vertebrates. The dynamics and diversity of sexdetermination systems observed in nature, including transitions (e.g., Volff et al. 2007), mixed systems (e.g., Quinn et al. 2007), or multiallelic polymorphisms (where X, Y, and W coexist, as observed in some fish and amphibians; Orzack et al. 1980;Miura et al. 1998), were actually regular outcomes in our simulations.

COLD-BLOODED VERTEBRATES?
Several mechanisms have already been proposed to explain turnovers. From these models, new sex determiners may spread and invade if (1) some of the new sex genotypes have a higher intrinsic adaptive value Orzack et al. 1980;Basolo 2001;Kraak and Pen 2002), (2) changes are driven by a genetic conflict arising from sex-chromosome meiotic drive or cytoplasmic sex-ratio distorters (e.g., Wolbachia, Hamilton 1967;Werren and Beukeboom 1998;Caubet et al. 2000), (3) sex-ratio selection is affected by changes in the production costs of male and female offspring (Kozielska et al. 2006), or (4) the new sex determiners are linked to a locus with strong sex-antagonistic effect (Van Doorn and Kirkpatrick 2007). None of these mechanisms, however, seem to specifically apply to cold-blooded vertebrates, and hence to account for the striking contrast offered with birds or mammals.
From our formalization, turnovers in cold-blooded vertebrates are the mere consequences of environmental changes. The underlying rationale is simple: Because temperature affects sex determination, environmental changes will induce sex-ratio biases, which in turn should favor the emergence of new sexdetermination genes or alleles. Sex-ratio selection is a powerful force, already invoked to explain turnovers between GSD systems (e.g., Caubet et al. 2000;Kozielska et al. 2006) or transitions from GSD to TSD (e.g., Bull 1981). As our results show, sex-ratio selection is triggered in ectotherms by environmental changes, due to the thermodynamic constraints shaping their reaction norms, and is thereby expected to induce turnovers in sex-determination systems.
This process does not affect warm-blooded vertebrates, because thermoregulation prevents the expression of temperature effects. Besides birds and mammals, some fish (e.g., tunas) do also show some levels of endothermy, but usually no parental thermoregulation during the thermally sensitive window of embryonic development. Carcharinid sharks apparently possess the relevant combination of viviparity and endothermy, but sexdetermination mechanisms are unfortunately poorly known in Elasmobranchs (although many seem to display heteromorphic sex-chromosomes; Maddock and Schwartz 1996).
Our model directly relates the prevalence of homomorphic sex chromosomes in cold-blooded vertebrates to the temperature sensitivity of sex determination imposed by ectothermy. The driving forces behind turnovers are to be found in climatic changes or range expansions. Note that the present argument differs from Perrin (2009), who also proposes a link between homomorphic sex chromosomes and the sex-reversal events induced by temperature changes, but without any turnover. From his model, X-Y recombination may occur in sex-reversed XY females (because recombination depends on phenotypic sex, not on genotypic sex), and this should prevent the decay of Y chromosomes even in absence of turnovers. As turnovers do occur in cold-blooded vertebrates (e.g., Hillis and Green 1990;Ezaz et al. 2006;Volff et al. 2007), our model provides a (nonexclusive) alternative to Perrin (2009) to explain the prevalence of homomorphic sex chromosomes in these groups.

CONSERVATION ISSUES
Whether the processes underlying sex-determination are dichotomous or continuous is not only a semantic question. It will affect the dynamics of turnovers, but also the risks of extinction. It was recently argued that TSD species are put at significant extinction risk by climatic changes, due to the progressive biases in sex ratios induced by temperature rises (e.g., Janzen 1994;Janzen and Morjan 2001;Mitchell et al. 2008). GSD species are often implicitly considered as protected against such changes, because genotypic systems (e.g., XY or ZW) are expected to necessarily produce even sex ratios. From our formalization, however, environmental change will also induce sex-ratio biases in such systems. Mass events of sex reversal and sex-ratio biases have already been documented in natural populations of amphibians and fish supposed to display pure GSD (e.g., Nagler et al. 2001;Matsuba et al. 2008).
From our simulations, extinctions are possible outcomes from such events. Small populations, in particular, were often locked in one recurrent pair or in ESD, which prevented adaptive transitions during climatic changes, and induced strong genetic loads associated with YY lethality, intersex sterility, or biased sex ratios. This often led to extinction when combined with low environmental variance in the liability trait A (steep β/σ E slope). Further negative effects of small population size stemmed from the lower probability of appearance of a new mutation, and a higher demographic stochasticity in sex ratios, inducing a risk of loosing by chance all members of one sex. These processes are also likely to hinder range expansions, because colonizing populations reaching new (and climatically different) areas often stem from rare long-distance propagules and have small effective population sizes.
The environmental component of variance (σ 2 E ) also mattered, with overall positive effects. High σ 2 E values (shallow β/σ E slope) made any genotype more likely to produce individuals of both sexes, which provided insurance against extinction, particularly under ESD. Interactions with drift were important, as in many instances large σ 2 E values could counteract the negative effects of small population sizes. As a consequence, extinctions are expected to occur mostly when both N and σ 2 E are low.

YY LETHALITY AND INTERSEX STERILITY
From our preliminary results, YY lethality limits the evolutionary potential of populations and enhances the risk of extinctions. This obviously prevents the evolution toward some specific systems (e.g., XY/YY, otherwise favored at low temperature), but also imposes segregation loads when YY genotypes are not directly involved in a recurrent pair (e.g., XX/XY or WW/WY), because sex reversal induced by high σ 2 E values produces unfit YY individuals. This may favor the transition to alternative recurrent pairs (e.g., XW/XX), but will also induce extinctions in small populations with reduced adaptive potential. YY lethality is likely to emerge during periods of evolutionary stasis, as functional genes become progressively involved into the expanding nonrecombining segment and accumulate deleterious mutations (Ohno 1967).
Frequent turnovers may thus allow maintaining the evolutionary potential and adaptability of populations, a point deserving further investigations.
Intersex sterility also had drastic influences on evolutionary paths, selecting for genotypes with liability-trait values far apart from the threshold. This might actually explain why sex reversal in cold-blooded vertebrates usually occurs outside the natural range of environmental conditions. In our simulations, this selective pressure often prevented the evolution of ESD (because XX genotypes were strongly counter-selected when producing sexfactor values within the intersex range) and favored direct transitions to a new male-heterogametic GSD system (e.g., WW/WY, with strongly masculinizing and feminizing effects of Y and W, respectively). It is worth noting that the increased selection for a feminizing allele W prevented its loss by drift in small populations, thereby lowering their extinction risk. This opposed our intuitive expectation that intersex sterility would increase extinction risks.

EMPIRICAL ISSUES AND MODEL EXTENSIONS
The quantitative predictions stemming from our simulations obviously depend on specific model assumptions and parameter values. As underlined above, empirical data can be used to calibrate key parameter values in specific cases. Values for α IJ,T (standardized liability-trait value for genotype IJ at temperature T) can be estimated from the proportion of sex reversals at this temperature. The shapes of reaction norms (and evidence for mixed systems) can be obtained by plotting α values as a function of temperature (Fig. 2). From our few empirical examples, parameter estimates fall well within the range used for simulations; empirical estimates of |β/σ E |, for instance, lie between 0.28 and 2.2 ( Figs. 2 and 3). Our assumptions of linear reaction norms and normal distribution of environmental variance are also well supported in the examples provided. Note however that a similar formalization might be provided for curvilinear reaction norms, such as found in lizards and turtles where males are only produced at intermediate temperatures (Ewert et al. 2004;Quinn et al. 2007). A larger scale literature survey to estimate the range of shapes and parameter values for reaction norms along the lines presented here (together with other relevant features such as intersex sterility or YY lethality) would constitute welcome empirical extensions of the present work.
Regarding theoretical extensions, the interaction between turnovers and Y decay will also constitute an important avenue for future research. Our preliminary results suggest that frequent turnovers might allow maintaining the evolutionary potential of populations, but a detailed formalization is required to precisely account for the dynamics of deleterious mutations. Sex-antagonistic genes are bound to interfere with this process, being the ultimate cause for the evolution of nonrecombination (e.g., Bergero and Charlesworth 2009). The fixation of male-beneficial sex-antagonistic alleles on the Y chromosome is expected to counter-select sex-reversed XY females, and thereby to affect the dynamics of turnovers.
Finally, a similar approach might be used to address the evolution of TSD as an adaptive strategy (e.g., Janzen and Phillips 2006). In our model, TSD only occurred as a side result when homozygous genotypes reached their pivotal temperature, and was counter-selected by sex-ratio selection at other temperatures. Evolving TSD as an adaptive strategy clearly requires behavioral adaptations allowing to fine-tune embryonic temperature so as to produce desired sex ratios. TSD might then outcompete GSD when optimal sex ratios differ from even (Hamilton 1967;Freedberg and Taylor 2007), or when temperature also affects fitness in a sex-specific way Bull 1981;Bulmer and Bull 1982;Conover 1984). Intersex sterility is also bound to play a crucial role in this context, because genotypes are mostly affected close to their pivotal temperature (Bull 1981). This should favor the evolution of reaction norms with extreme sensitivity to temperature (steep β/σ E slope) across the intersex sterility domain. Detailed investigations along the lines sketched out in the present article might help shedding light on this complex and fascinating issue.