High elevation increases the risk of Y chromosome loss in Alpine skink populations with sex reversal

The view that has genotypic sex determination and environmental sex determination as mutually exclusive states in fishes and reptiles has been contradicted by the discovery that chromosomal sex and environmental influences can co-exist within the same species, hinting at a continuum of intermediate states. Systems where genes and the environment interact to determine sex present the opportunity for sex reversal to occur, where the phenotypic sex is the opposite of that predicted by their sex chromosome complement. The skink Bassiana duperreyi has XX/XY sex chromosomes with sex reversal of the XX genotype to a male phenotype, in laboratory experiments, and in field nests, in response to exposure to cold incubation temperatures. Here we studied the frequency of sex reversal in adult populations of B. duperreyi in response to climatic variation, using elevation as a surrogate for environmental temperatures. We demonstrate sex reversal in the wild for the first time in adults of a reptile species with XX/XY sex determination. The highest frequency of sex reversal occurred at the highest coolest elevation location, Mt Ginini (18.46%) and decreased in frequency to zero with decreasing elevation. We model the impact of this under Fisher’s frequency-dependent selection to show that, at the highest elevations, populations risk the loss of the Y chromosome and a transition to temperature-dependent sex determination. This study contributes to our understanding of the risks of extinction from climate change in species subject to sex reversal by temperature, and will provide focus for future research to test on-the-ground management strategies to mitigate the effects of climate in local populations.


Introduction
The primary signal directing embryos to develop phenotypically as either male or female is known as sex determination. The diversity of sex determination systems within vertebrates is remarkable, and has classically been considered to occur via one of two mechanisms-genotypic sex determination (sex determined at the time of fertilization by genetic factors independent of environmental influence, GSD) or environmental sex determination (sex determined by environmental factors that act after fertilization, ESD). In reptiles, these alternatives were thought to be incompatible and to form a mutually exclusive dichotomy. However, reversal of sexual genotype by environmental factors is common in fish (Wang et al. 2019;Hattori et al. 2020;Miyoshi et al. 2020) and amphibians (Flament 2016;Lambert et al. 2019;Nemesházi et al. 2020). Sex reversal was suspected in turtles based on screening for H-Y antigens (Zaborski et al. 1982;Servan et al. 1989), but it was not conclusively demonstrated in reptiles until recently (Quinn et al. 2007;Holleley et al. 2015;2016;Wiggins et al. 2020). Sex reversal in reptiles, as with temperaturedependent sex determination or TSD (Mitchell and Janzen 2010), has profound implications for our understanding of how this diverse group of species can respond to climate change.
Theoretical studies show that as the frequency of sex reversal increases in a population, a likely response is a reduction and possible elimination of the Y or W chromosome under Fisher's frequency-dependent selection (Fisher 1930;Bull 1981;Grossen et al. 2011;Holleley et al. 2015;Bókony et al. 2017;Schwanz et al. 2020;Geffroy and Wedekind 2020). Indeed, at least one population of central bearded dragon (Pogona vitticeps, (Ahl 1926)), in which the sex-reversed ZZ female phenotype is fertile, appears to be on the brink of such a transition (Fig. S2 of Holleley et al. 2015). The implications of this work are that species with GSD may not be entirely immune from the demographic destabilization that climate change potentially brings to species with TSD. For some GSD species at least, global warming may well drive a transition to TSD (already demonstrated in the lab in one generation- Holleley et al. 2015); further warming increases the risk of extinction because of insufficient time for them to evolve and optimize the threshold for sex reversal and avoid the adverse demographic consequences of extreme sex ratio bias.
The temporal scale of these combined demographic and evolutionary changes make these processes challenging to study. However, it is possible to use altitudinal or latitudinal climatic trends as a surrogate for how species might respond to temporal trends (Doody et al. 2006;Pen et al. 2010;Castelli et al. 2020). The widespread dragon lizard, P. vitticeps, has a ZZ/ZW system of sex determination (Ezaz et al. 2005) as in birds (Smith et al. 2007) but the ZZ genotype is reversed to present a viable female phenotype by high incubation temperatures (Quinn et al. 2007;Holleley et al. 2015). Furthermore, modelling the consequences of sex reversal has shown that the W chromosome (in this case) can persist indefinitely at low frequency in the population, particularly if there is immigration from adjacent populations that remain predominantly GSD (Schwanz et al. 2020). So the transition from GSD to TSD is not a unidirectional process, with the residual and cryptic presence of the W chromosome enabling a rapid halt to the transition and reversion to GSD if and when conditions change.
The scincid lizard Bassiana duperreyi (Gray 1838) also has genotypic sex determination, but with an XX/XY system (Donnellan 1985;Matsubara et al. 2016) and reversal of the XX female genotype to a male phenotype by low incubation temperatures (Shine et al. 2002;Radder et al. 2008;Quinn et al. 2009). Despite having sex chromosomes, hatchling sex is influenced by the temperature during incubation, both in the lab and the field (Shine et al. 2002;Radder et al. 2008;Holleley et al. 2016); it is not influenced by hydric variation over the range of soil water potentials recorded in natural nests (Flatt et al. 2001). High rates of sex reversal occurred in nests of an Alpine population (36°5 ′ 8.15″ S 148°13′ 1.93″ E, Jagumba in Kosciuszko National Park, Australia; 28% XX male- Holleley et al. 2016), but whether this translates to equally high rates of sex-reversed adults, and whether those sex-reversed adults are reproductively viable is not known. Here we use a new and reliable PCR test  to show that, unlike P. vitticeps (Castelli et al. 2020), the frequency of sex reversal in adults of B. duperreyi does vary in accordance with expectation along an ambient thermal gradient (elevational gradient) extending from the Australian Alps (1640 m a.s.l.) to the Victorian lowlands (20 m a.s.l.). Our modelling of the impact of sex reversal on the frequency of the XY genotype suggests that, at the upper altitudinal extremes of the distribution of the species, the Y chromosome may be lost, precipitating a transition to TSD.

Study species
The eastern three-lined skink, B. duperreyi, is a mediumsized (60-80 mm snout-vent length, SVL) oviparous lizard widely distributed at cool high-elevation sites close to the upper elevational limits for oviparous reproduction by Australian lizards (Shine and Harlow 1996); its distribution extends through an altitudinal gradient to coastal regions of New South Wales (NSW) and Victoria in south-eastern Australia (Cogger 2014). B. duperreyi appeared in the Upper Miocene and/or Upper Pliocene period (4.7-6.6 Mya), during which its distribution would have been strongly influenced by paleoclimatic conditions (which fluctuated between warm-wet and cool-dry) and other biogeographical processes which have prevailed since (Dubey and Shine 2010). Female B. duperreyi lay a single clutch of 3-11 eggs from early Australian summer (late November) to December (Shine and Harlow 1996).

Study area
Ten sites along an elevational gradient were selected in mainland south-eastern Australia spanning a large portion of the distribution of B. duperreyi (Fig. 1A The climate in south-eastern Australia is temperate. We compiled climatic data for each study location using the Scientific Information for Land Owners database maintained by the Queensland Department of Natural Resources and Water (SILO, https://www.longpaddock.qld.gov.au/ silo/, last accessed 8 Apr 2020) (Jeffrey et al. 2001). The SILO data drill provide daily temperatures, rainfall and other climate variables for 0.050°grids across Australia, interpolated from point measurements made by the Australian Bureau of Meteorology. Historical monthly climate data from January 1895 to January 2019 were downloaded for analysis. The mean annual maximum temperature (Tmax) was 15.4 ± 6.02°C and mean annual minimum (Tmin) of 4.5 ± 4.25°C at Mt Ginini (n = 1561). Anglesea Tmax was 18.0 ± 3.95°C and Tmin was 9.6 ± 2.66°C (n = 1561). The highest mean annual rainfall was recorded in the Mt Kosciuszko National Park field location (113.9 ± 70.32 mm), while the highest mean annual solar radiation was recorded in Piccadilly Circus field location (17.3 ± 6.39 MJ/m 2 ). Detailed climatic data for each location are provided in supplementary materials (Table S1), which we used to analyze the frequency of sex reversal with an elevational gradient.

Sample collection
We conducted our fieldwork in areas where the lizards were most abundant, i.e., in natural open areas but also often in areas artificially cleared for overhead hydroelectric power lines, fire banding trails, roads and or tracks inside the forest areas in above-selected field locations. Adults tend to aggregate in these open areas during the nesting season (Radder and Shine 2007). Adult specimens were captured by hand either when active or when inactive under rocks or logs. We measured each SVL (tip of snout to the anterior end of cloaca) and tail length (the anterior end of cloaca to tip of the tail) for each lizard. All measurements were taken using digital Vernier callipers on living specimens (nearest 1 mm). The phenotypic sex of males was identified by hemipene extrusion (Harlow 1996); those for which prominent hemipenes could not be everted were identified as female. Female lizards were checked by abdominal palpation to see if they were gravid (Holmes and Cree 2006). Breeding coloration was used to corroborate the determination of sex by hemipenile extrusion-males have a prominent reddish-orange throat when in breeding condition, lacked by females. All animals with reddish-orange coloration on the throat extruded hemipenes when probed, so there were no cases of conflict between the two criteria when present. Some animals without breeding coloration also extruded hemipenes, but this was typically outside the breeding season, so it was possible that some individuals who did not extrude hemipenes and did not have breeding coloration were nevertheless males. However, in no cases did individuals presenting as phenotypic females deliver a male outcome (XY) in the sex testing. Indeed, our approach to phenotypic sexing was considered highly reliable because in no case was there a mismatch in phenotypic sex and gonadal sex (testes or ovary) in dissections for other studies, and no mismatches occurred in the panel of 20 males and 20 females used as the validation panel for our sex-linked marker . In seven instances, lizards escaped capture before their phenotypic sex could be determined.
Tail tips (4-5 mm) were removed with a sterile blade and the free-flowing blood drop collected onto a labelled Whatman FTA™ Elute Card (WHAWB12-0401, GE Healthcare UK Limited, UK); tail tips were collected into labelled 1.5 ml tubes containing 95% ethanol. All animals were immediately released at their point of capture. All experimental protocols were conducted with the permission of Animal Ethics Committees at the University of Canberra and the CSIRO. All experiments were conducted in accordance with guidelines and regulations established by these committees.

SNP genotyping
Tail tissue samples of a total of 100 B. duperreyi from the 10 populations (10 samples per location) were passed to Diversity Arrays Technology Pty Ltd (Canberra) (DArT) for reduced representational sequencing. Genomic DNA was extracted using a NucleoMag ® 96 Tissue kit (Macherey -Nagel, Düren, Germany) coupled with NucleoMag SEP (Ref. 744900) to allow automated separation of high-quality DNA on a Freedom Evo robotic liquid handler (TECAN, Männedorf, Switzerland). Four combinations of restriction enzymes were evaluated for B. duperreyi (PstI enzyme combined with either HpaII, SphI, NspI, and MseI) and the restriction enzyme combination of PstI (recognition sequence 5′-CTGCA|G-3′) and SphI (5′-GCATG|C-3′) was selected for complexity reduction by double digestion (Kilian et al. 2012). A full account of the process of generating SNP genotypes for these samples is given by Georges et al. (2018). The data comprised a matrix of SNP loci by individuals, with the contents stored as 0, homozygote, reference state; 1, heterozygote; and 2, homozygote, alternate state.

Additional SNP filtering and visualization
The SNP data and associated locus metadata were read into a genlight object (R Package adegenet-Jombart 2008) to facilitate processing with package dartR v.1.5.5 . Only loci with >99% repeatability (repAvg) were chosen for subsequent analysis. Further filtering was undertaken based on the call rate (>90%) and where multiple SNPs occurred within a single sequence tag, only one was retained at random. The population sample sizes were small (n = 10), so we could not filter loci for departures from Hardy-Weinberg equilibrium or linkage disequilibrium; the sparse sampling of loci across the genome allows the reasonable assumption of little or no linkage between loci. We regard the data remaining after this additional filtering as highly reliable. Genetic similarities for individuals and populations were visualized using principal coordinates analysis (PCoA) as implemented in the gl.pcoa and gl.pcoa.plot functions of dartR. A scree plot of eigenvalues (Cattell 1966) guided the number of informative axes to examine, taken in the context of the average percentage variation explained by the original variables (using the diagnostics provided by gl.pcoa function in dartR).

Fixed-difference analysis and genetic diversity
To examine the possibility that more than one taxon (Operational Taxonomic Unit, OTU) might be contributing to the altitudinal cline, potentially confounding the comparisons, a fixed-difference analysis was done using the scripts gl.fixed.diff and gl.collapse in dartR. A fixed difference between two populations at a biallelic SNP locus occurs when all individuals in one population are fixed for the reference allele and all individuals in the other population are fixed for the alternate allele. Accumulation of fixed differences between two populations is a clear indication of lack of gene flow . Expected heterozygosity, a measure of genetic diversity, was obtained for each population from allele frequencies using the gl.report. heterozygosity function of dartR. We used gl.ibd function in dartR to calculate isolation by distance (F ST /1 − F ST versus Euclidean distance), tested with a Mantel test.

Genotypic sex and sex reversal frequency
For sex testing, DNA was extracted from tail tips using a Gentra Puregene commercial kit (Qiagen Science, Maryland, USA) following the manufacturer's protocols; DNA was extracted from blood samples (FTA TM Elute Micro Card) following the manufacturer's protocols. DNA purity was determined using a NanoDrop 1000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA) and quantified using the Qubit 2.0 Fluorometric Quantitation (Invitrogen, Life technologies, Sydney, NSW, Australia). Genotypic sex was identified using a PCR test comprising seven Y-specific markers validated against samples of known sex individuals from across the study sites (i.e., Anglesea to Piccadilly Circus) of B. duperreyi . Briefly, in applying the test we used 1× MyTaq TM HS Red mix (Bioline USA Inc. USA), 4 µM of each primer and 25 ng of genomic DNA. The PCR cycling conditions used an initial touchdown phase to increase the specificity of amplification: denaturing at 95°C, annealing temperature stepping down from 70°C by 0.5°C per cycle for 10 cycles, extension at 72°C. This was followed by 30 cycles at 65°C annealing and 72°C extension. PCR products were visualized on a 1.5% agarose gel using SYBR Safe (Life Technologies, Carlsbad, USA). The samples that showed an amplified band for each of the seven markers are recognized as XY individuals, whereas as the samples for which a band was not amplified in all seven markers were recognized as XX individuals. The seven markers always concurred in their determination in the original study  and in this study, which renders false negatives highly unlikely (they would present as some but not all markers failing). There were also no feminized individuals. Phenotypic male lizards showing genotype-phenotype discordance were classified as sex reversed ). All molecular sex tests were conducted blind.
Coree Flat West population was assigned to the Coree Flat East population owing to the low number of samples that we acquired from these two locations. Simple linear regression was used to determine the relationship between the proportion of sex-reversed males (arcsin transformed) in each population with respect to elevation, and to characterize the trends in yearly mean maximum (Tmax) and yearly mean minimum temperature (Tmin) against elevation in each field location. The Pearson correlation was used to evaluate the magnitude and direction of the association between the frequency of sex reversal and climatic variables [Tmax, Tmin, total rain (mm), evaporation (mm), radiation (Mj/m 2 ) and vapour pressure (hPa)] including elevation recorded at each field site. For these correlations, we used only data from the Alpine populations.

Modelling the consequences of sex reversal
We modelled the decline of the XY genotype resulting from frequency-dependent selection because of the overproduction of males through sex reversal of the XX genotype with decreasing temperature (Shine et al. 2002;Radder and Shine 2007) following the logic of Bull (1981). Sex reversed individuals are assumed to be viable and fitness within a sex is considered to be the same for all genotypes. Let the starting frequency of XY among zygotes be y, the starting frequency of XX among zygotes be 1-y, and let a fraction P[T] of XX become reversed to a male phenotype. In any one generation n, we have the proportion of male phenotypes (m n ) equal to the sum of the frequencies of normal XY males (y n ) and sex-reversed XX males, The frequency of the XY genotype in the next generation is given by We used previously published experimental data (Radder et al. 2008;Shine et al. 2002) and their estimate of 20°C for the threshold below which offspring sex is temperaturedependent and above which offspring sex is solely genetically determined (Telemeco et al. 2009) (not to be confused with the threshold at which a 50:50 sex ratio is observed), to generate the equation for P[T] (Fig. S1 and Table S2) given by, where T is in°C. We iterated for an equilibrium solution for y for various values of temperature T. Overlapping generations will delay the rate of convergence to equilibrium but will not affect the equilibrium value for a particular temperature.

SNP genotyping
A total of 97,417 polymorphic SNP loci were scored for 100 individuals from 10 populations of B. duperreyi in southeastern Australia (Fig. 1A). After stringent filtering on repeatability (10,161 loci <0.99) and call rate (61,524 loci <0.90), and after removing the secondary SNPs (7881), the number of SNP loci in the data set dropped to 17,851. Considerable genetic structure existed in B. duperreyi in southeastern Australia (Fig. 2). In particular, the lower elevation populations of Anglesea and Westernport Bay were distinctly different from the Australian Alps populations and probably represent a distinct taxon (OTU). One individual (AA080788) was intermediate between the coastal and Alpine populations, likely to represent contemporary hybridization, and was removed from the fixed difference analysis. Evidence for the coastal and Alpine populations being separate taxa was provided by the fixed difference analysis-the Anglesea population differed from the Alpine populations by an average of 280 (256-308) diagnostic fixed differences and the Westernport Bay population differed from the Alpine populations by an average of 336 (308-374) fixed differences (Table 1). All the Alpine populations collapsed into a single OTU based on corroborated fixed differences (option tpop = 1 in gl.collapse of dartR), after which the Anglesea population differed from the Alpine populations by 160 fixed differences (false positive expectation 76.1, p < 0.0001) and the Westernport Bay population differed from the Alpine population by 196 fixed differences (false positive expectation, 79.1, p < 0.0001). Anglesea and Westernport Bay differed by only 2 fixed differences (false positive expectation 15.1, p = 1.0000, ns) and so collapsed into a single OTU. Hence, two distinct taxa (OTUs) were identified, based on significant fixed differences, that of Anglesea/Westernport Bay (Coastal OTU) and the Alpine populations from Mt Ginini to Dartmouth (Alpine OTU), differing by 117 fixed allelic differences (false positive expectation 37.4) ( Table 1). The Alpine population showed evidence of isolation by distance (Mantel's test; r = 0.71, p < 0.002). The average expected heterozygosity varied from 0.10 at Westernport Bay to 0.16 at Piccadilly Circus (mean 0.13). Within the Alpine OTU, average expected heterozygosity varied from 0.19 at Cooma to 0.16 at Piccadilly Circus (mean 0.13).

Genotypic sex and sex reversal frequency
A total of 639 tail snips were collected during the field survey. Of them, 399 were phenotypic males and 233 were phenotypic females, and 7 escaped before they had their phenotypic sex identified, and those samples were removed from the PCR sex test. A total of 33 (5.22%) adult males were sex-reversed in the populations. Zero sex reversal was observed at the low elevation sites ranging from Dartmouth (340 m a.s.l) to Anglesea and Westernport Bay (40-20 m a.s.l), despite the markers having been validated for the putative distinct taxon from Anglesea/Westernport Bay . However, after identifying two OTUs in the B. duperreyi distribution within south-eastern Australia, we restricted our attention to the sex reversal frequency in the Alpine OTU (Dartmouth to Mt Ginini). The highest frequency of sex reversal occurred at the highest elevation, Mt Ginini (18.46%) (Fig. 1B) and a total of 6.58% adult males were sex-reversed across all the sampled Alpine populations (Table S3). Frequency of sex reversal was positively correlated with elevation (Pearson's correlation coefficient r = 0.97, p < 0.001) (Table S4); Tmax (R 2 = 0.66, p < 0.05) and Tmin (R 2 = 0.60, p < 0.05) were significantly negatively correlated with sex reversal frequency ( Fig. 3; Table S4).

Modelling the consequences of sex reversal
The frequency of the XY genotype is predicted to decline precipitously with decreasing incubation temperature as the system moves towards a 1:1 sex ratio equilibrium (Fig. 4). At incubation temperatures of 18°C and below, we would expect the complete loss of the Y chromosome from the population. Our wild populations, using the nest temperature data of Telemeco et al. (2009Telemeco et al. ( , 2010, indicate that some high-altitude sites (e.g., Mt Gingera, 1865 m a.s.l, not sampled in this study) are within the thermal range that would lead to loss of the Y chromosome (Fig. 4). Fig. 2 Graphical representation of genetic similarity between individuals using principal coordinates analysis. The two primary clusters (Alpine OTU and Coastal OTU) differed by 117 diagnostic fixed allelic differences. There were no significant fixed differences between populations within clusters. Axes not to scale.

Discussion
The eastern three-lined skink has heteromorphic sex chromosomes, with males as the heterogametic sex (XY) (Donnellan 1985). Individuals of the species with the XX genotype, normally destined to become female, can be sex reversed to a male phenotype at low incubation temperatures (<20°C) under laboratory conditions (Radder et al. 2008) and in natural nests ). Here we have presented a detailed account of naturally occurring sex reversal in adult B. duperreyi along an elevational gradient confirming that sex-reversed hatchlings are viable and survive to adulthood.
Our combination of seven sex-linked markers, concordant in their indication of genotypic sex in all determinations in this study and that of Dissanayake et al. (2020) eliminates the possibility that what we have observed is a gradient in the frequency of recombination. There is no reason to expect our seven independently derived markers to be tightly linked, so if recombination was varying with elevation, then the concordance of our seven markers would break down with elevation. There is no evidence of this. Seven novel Y-chromosome markers increases the confidence of chromosomal sex identification in B. duperreyi because it dramatically reduces the risk of a recombination event being misinterpreted as evidence of sex reversal.
No sex reversal was observed in the populations at the lowest elevation, on the coast at Anglesea and Westernport Bay. These populations appear to represent a different species from the Alpine populations, which opens the possibility that the coastal taxon does not have thermolabile sex for historical reasons, not as a direct contemporary response to climate. We therefore excluded from consideration the populations from Anglesea and Westernport Bay, and focused our analysis on the single taxon which we refer to as the Alpine OTU. Having done that, there was a significant positive relationship between sex reversal and elevation. The frequency of sex reversal in adults ranged from 18.46% at the highest, coolest location (Mt Ginini, 1640 m a.s.l., mean annual temperature 9.9 ± 7.58°C) to zero at the lowest, warmest location (Dartmouth, 340 m a.s. l., 13.4 ± 8.38°C). These results establish that the frequency of sex reversal varies as expected with prevailing climate.
Our modelling shows that the frequency of the Y chromosome in the population can be expected to decline as the frequency of sex reversal increases. The question arises as to whether temperatures exist within the species range, or are likely to exist under future climate change, that could lead to the complete loss of the Y chromosome driven by overproduction of males through sex reversal. Is a transition to TSD likely? Our modelling of the equilibrium state for the frequency of the XY genotype as nest temperature is decreased shows that under frequency-dependent selection Table 1 Matrix of euclidean genetic distances (below diagonal) and allelic genetic differences (above diagonal) between the ten populations sampled for Bassiana duperreyi.  the XY genotype begins to decline in frequency at 23°C and is lost at nest temperatures below 18°C (strictly constant temperature equivalents, CTE (Georges 1989;Georges et al. 1994)). Mean nest temperatures during the middle of incubation (weeks 5-8) were recorded as 17.1°C at Mt Gingera (1855 m a.s.l.), 19.3°C at Mt Ginini (1640 m a.s. l.), 20.6°C at Piccadilly Circus (1240 m a.s.l.) and 20.6°C at Coree Flat (1040 m a.s.l.) (Telemeco et al. 2009(Telemeco et al. , 2010Shine et al. 2003). Although these data have not been corrected for diel fluctuations in the nest, it is clear that at the highest elevations, there is the potential for loss of the Y chromosome and a transition to TSD provided there is no dispersal from adjacent regions where the Y is persistent (Schwanz et al. 2020).
Reduction of the Y chromosome under the influence of climate follows a neutral pathway, in that once the system achieves a 1:1 sex ratio, the frequency of the Y chromosome in the population is evolutionarily stable (Bull 1981).
Many theoretical models, such as genetic drift (Bull and Charnov 1977) and sex ratio selection (Bulmer and Bull 1982;Kozielska et al. 2010;Schwanz et al. 2020) emphasize the importance of natural selection for promoting transitions in sex determination systems (Natri et al. 2019). In P. vitticeps, sex reversal can lead to loss of the W chromosome under Fisher's frequency-dependent selection alone (Holleley et al. 2015), without the need to invoke an advantage under conventional natural selection (see also, Hurley et al. 2004;Grossen et al. 2011;Bókony et al. 2017). The same considerations apply to B. duperreyi. Nevertheless, the transition to a system free of the Y chromosome may be accelerated should a fitness advantage of sex reversal under conventional natural selection exist (Holleley et al. 2015;2016;Li et al. 2016) or retarded by a fitness cost (Cotton and Wedekind 2009), depending upon the species and the circumstances. This provides a great opportunity for future research on the evolutionary dynamics of sex reversal in B. duperreyi. P. vitticeps differs from B. duperreyi in that the direction of the reversal (at high temperatures rather than low) is aligned with the projections for global climate change. In P. vitticeps, sex reversal and the potential loss of the W chromosome occurs well within the range of temperatures experienced in nature (Holleley et al. 2015), and is likely to be exacerbated by climate change. Sex reversed individuals of P. vitticeps are fertile, and ZZ × ZZ crosses can be used to generate viable lines in which the W chromosome is absent. This raises a second unanswered question-the Fig. 3 The frequency of sex reversal in adult Bassiana duperreyi. A The rate of sex reversal in phenotypic males Bassiana duperreyi increases with elevation (F (1,5) = 71.39, p < 0.001; R 2 = 0.94) (green triangles). Numbers indicate the field locations, as described in Fig. 1. B Linear regression of yearly (January 1895 to January 2019) mean Tmax (°C) (red squares) (F (1,5) = 14.88, p < 0.05; R 2 = 0.74) and Tmin (°C) (blue circles) (F (1,5) = 7.82, p < 0.05; R 2 = 0.61) with elevation. Points represent mean ± SD (n = 1561). Temperature data obtained from SILO database (Jeffrey et al. 2001). Broken lines denote the 95% confidence interval. Fig. 4 Modelling of the decline of the XY genotype resulting from the frequency-dependent selection. Relative frequency of the XY genotype from the entire population declines precipitously with decreasing incubation temperature and requires only a small drop in environmental temperature to precipitate complete loss of the Y chromosome. Nest temperature of the highest altitude wild population sampled in this study is shown in green (Mt Ginini) (Telemeco et al. 2010) and resides on the precipice. One higher elevation population not sampled in this study (Mt Gingera, 1865 m a.s.l) (Telemeco et al. 2010) (red broken line) suggests conditions exist in alpine populations that are suitable for complete loss of the Y chromosome. Black broken line shows the actual recorded frequency of sex reversal in the adult population at the Mt Ginini. sex-reversed individuals of B. duperreyi survive to adulthood, but are they fertile? Establishing the fertility of sexreversed B. duperreyi is more challenging than for P. vitticeps because the reversal is to a male phenotype. The reproductive viability of sex-reversed male B. duperreyi is yet to be established. This places a caveat on any interpretations we may place on the evolutionary significance of sex reversal in this species.
A third question arises as to whether the demographic consequences of sex reversal-a predominance of males arising before the system can come to a 1:1 sex ratio equilibrium, and consequent reduction in effective population size-is a factor limiting the distribution of B. duperreyi at high altitudes, or whether the limitation arises through egg survivorship. If the climate is static, a 1:1 primary sex ratio is achieved locally despite a reduction in the frequency of the Y chromosome. The system comes to a series of equilibria across the landscape with no inherent limitations placed on the species distribution. Under climate change, however, the system is pushed to a disequilibrium state and there is the potential for a very substantial sex ratio skew until the populations come to evolutionary equilibrium. Such a skew could deliver a limit to the upper range of the species because of the demographic consequences of the gross overproduction of male phenotypes (Boyle et al. 2014). More detailed work needs to be done at the upper extremes of the distribution of this species (Mt Gingera, 1855 m a.s.l. and Mt Kosciuszko, 2020 m a.s.l.) to identify if these evolutionary and demographic processes are in play. The question of whether sex ratio skew arising from TSD or embryo physiology and viability establishes the distributional boundaries to species remains unresolved, complicated by the typical bias of sex ratio to overproduction of females at extreme temperatures in most species of reptile.
For B. duperreyi, global warming is likely to alleviate the demographic impact of sex ratio skew, because temperatures that produce a 50:50 sex ratio (23°C; Radder et al. 2008) will extend to higher altitudes leading to range expansion. However, projected changes in future climatic conditions (Benestad 2003;Tokarska et al. 2020) predict more extreme events, such as elevation of temperatures to greater extremes than have been experienced in the past, or dropping to exceptionally low temperatures during the nesting season with profound impacts on biological processes in many taxa (Easterling et al. 2000;Hari et al. 2006;Dunham et al. 2011;Scheffer et al. 2015;Gamelon et al. 2017;Valenzuela et al. 2019). Hence, at a finer temporal scale, extremely cold temperatures, even for a short time during the natural incubation period, might have a sex reversing effect in the Alpine populations and disrupt the contribution of a particular breeding year to the future population (see Schwanz et al. 2020-Fig. 3). Nest sites subjected to particularly cold temperatures in 1 year resulted in 28% of XX hatchlings to be sex-reversed .
Not many species have been examined for instances of sex reversal and the widespread occurrence of homomorphic sex chromosomes in reptiles means that instances of sex reversal would not come to notice incidentally. Sex reversal in the wild may also occur in the yellow-bellied water skink (Eulamprus heatwolei) and the spotted snow skink (Niveoscincus ocellatus), based on thermosensitive sex determination in both species, and on the discovery of XY chromosomes in E. heatwolei (Cornejo-Páramo et al. 2020) and sex-linked loci consistent with male heterogamety in N. ocellatus (Hill et al. 2018). Clearly, sex reversal under extremes of temperature experienced in natural nests or during gestation is reasonably common in lizards and likely a powerful evolutionary component in generating and maintaining lability and diversity in reptile sex-determining modes generally (Holleley et al. 2015;2016). This is true also of some amphibians (Rodrigues et al. 2014) and fish species (Baroiller et al. 2009;Honeycutt et al. 2019) in which climatic gradients and environmental temperature strongly correlate with sex chromosome differentiation. In the Atlantic silverside Menidia menidia, sex is determined by an interaction between genotype and temperature and sex ratio differs among populations from different latitudes in response to temperature (Duffy et al. 2015).
In conclusion, in species with TSD and GSD systems and systems where temperature and genotype interact, predicting evolutionary responses to climate change becomes complex. We have shown a significant correlation between environmental variables and sex reversal in the wild adult population of B. duperreyi for the first time in a reptile species with XX/XY sex determination. We provide evidence that sex reversal can be a biologically significant process in Alpine populations and provide an opportunity to reinterpret sex ratio trends observed in other Alpine reptiles (as has been the case with N. ocellatus and E. heatwolii). Theory suggests that frequency of the Y chromosome in wild populations decreases with higher elevation, and that under extreme cold conditions in Alpine Australia, the likelihood of complete loss of the Y chromosome and a transition of B. duperreyi from GSD to TSD is quite high under future climatic cooling. In contrast, lower elevation populations are likely to exist stably as a GSD population utilizing XY sex chromosomes in the absence of sex reversal. Studies of the frequency of sex reversal in the wild nests, reproductive viability of adults, and the fitness consequences of sex reversal in this remarkable species and others like it is a priority, if we are to fully establish the contribution of sex reversal to our understanding of sex ratio evolution and sex chromosome evolution under environmental change in heterothermic vertebrates.

Data availability
All data and materials are presented in the main paper and supplementary information. The SNPs data and custom scripts used in this study have been deposited are located on the Dryad Digital Repository: https://doi.org/10.5061/drya d.9p8cz8wdk. The Scientific Information for Land Owners data used in this paper can be accessed at https://www. longpaddock.qld.gov.au/silo/.