Evolving thermal thresholds explain the distribution of temperature sex reversal in an Australian dragon lizard

Species with temperature‐dependent sex determination (TSD) are particularly vulnerable to climate change because a resultant skew in population sex ratio can have severe demographic consequences and increase vulnerability to local extinction. The Australian central bearded dragon (Pogona vitticeps) has a thermosensitive ZZ male/ZW female system of genetic sex determination (GSD). High incubation temperatures cause reversal of the ZZ genotype to a viable female phenotype. Nest temperatures in the wild are predicted to vary on a scale likely to produce heterogeneity in the occurrence of sex reversal, and so we predict that sex reversal will correlate positively with inferred incubation conditions.


| INTRODUC TI ON
In most vertebrates, the outcome of sexual differentiation is binary, with individuals developing phenotypically as male or female following distinct developmental trajectories. Although the outcomes of sexual differentiation are highly conserved among vertebrates, great diversity exists in the genetic and environmental factors that have acquired master or influential roles in the determination of sexual fate (Bachtrog et al., 2014). These signals can be genetic or environmental, but rather than a dichotomy, a continuum occurs from genetic sex determination (GSD) to environmental sex determination (ESD) in that the critical influence can be genetic, environmental or an interaction between the two (Sarre et al., 2004).
Transitions between genetic and environmental sex determination systems have occurred frequently among vertebrates on an evolutionary timescale (Bachtrog et al., 2014;Van Doorn, 2014;Sarre et al., 2011;Warner, 2011), and underlying these transitions is evolution in the threshold temperature for sex reversal (Quinn et al., 2011). The intrinsic and extrinsic factors which are capable of driving evolution towards one or the other extreme of the continuum have been widely modelled computationally and demonstrated experimentally (Sarre et al., 2004). Temperature sex determination (TSD) is the most common mode of ESD and is observed widely among fish, reptiles and amphibians (Bachtrog et al., 2014). Mechanistic models under different climatic scenarios have identified the importance of factors other than temperature alone in generating primary sex ratios, including maternal behaviour, nesting phenology, fecundity, and threshold temperature (Boyle et al., 2014;Schwanz et al., 2020). Species with TSD are particularly vulnerable to climate change because skews in the population sex ratio in response to extreme conditions can have severe demographic consequences, increasing vulnerability to local extinction (Mitchell & Janzen, 2010). Significantly skewed sex ratios caused by climate warming already occurs in contemporary populations of TSD reptiles (Jensen et al., 2018) and fish (Honeycutt et al., 2019). The tipping-point at which skewed sex ratios lead to demographic collapse is not yet known.
Species with GSD are typically assumed to be resilient to sex ratio skew (Bókony et al., 2019). The occurrence of sex reversal (temperature overriding an underlying GSD system) challenges this view. Sex reversal occurs in wild populations of two Australian reptiles: the central bearded dragon (Pogona vitticeps;Holleley et al., 2015), the eastern three-lined skink (Bassiana duperreyi; Holleley et al., 2016), and is suspected in six more species of lizard and turtle (Holleley et al., 2016;Wiggins et al., 2020). Sex reversal is also a frequent phenomenon in fish (Baroiller & D'Cotta, 2016).
Sex reversal provides a mechanism by which increasingly extreme conditions leads to the loss of the heterogametic chromosome (the Y or the W) and a transition from GSD to TSD (Holleley et al., 2015;Schwanz et al., 2020). Persistence of the population then requires one of two things; behavioural modifications that alter the nest environment, or evolution in the temperature threshold for sex reversal (Düsing, 1884;Edwards, 2000;Fisher, 1930). Though maternal behaviours can buffer against the effects of climate, this is not universally true and these behaviours are not always heritable (Doody et al., 2006;Ewert et al., 2005;Gutzke & Paukstis, 1983;Mitchell & Janzen, 2019;Refsnider et al., 2013;Refsnider & Janzen, 2016;Warner & Shine, 2008). Additionally, evolutionary responses take time and require existing heritable variation upon which to act (Schwanz et al., 2020). The time lag implicit in these evolutionary forces can leave thermosensitive species (including those with GSD and sex reversal) extremely vulnerable to local extinction under rapid or extreme climate change via demographic skews.
The central bearded dragon (P. vitticeps) displays a sex determination system characterized by a ZZ/ZW GSD system of female heterogamety (Ezaz et al., 2005) and sex reversal of the ZZ male genotype to phenotypically female under the influence of high temperature (Holleley et al., 2015;Quinn et al., 2007). Sex reversal occurs rarely at 31°C but rises in frequency with temperature until 36°C where almost 100% of ZZ genotype animals are feminized (Holleley et al., 2015). This Australian agamid is distributed across a range of climatic types spanning mesic to xeric climatic regions across central, south-eastern and northeastern Australia (Cogger, 2018;Rej & Joyner, 2018), providing a unique opportunity to study the dynamics of genetic and environmental influences on sex determination. Some populations of P. vitticeps may be in the early stages of transition from GSD to TSD via loss of the W chromosome (Holleley et al., 2015).
Here, we examine the distribution of sex reversal across the geographic range of P. vitticeps, using historical museum specimens and contemporary wild-caught individuals. Nest temperatures are predicted to vary across this range ( Figure 1). We predict that sex reversal frequency in the wild will be positively correlated with inferred incubation temperatures if there is a common and static threshold for temperature sex reversal. Conversely, if there is local adaptation in the threshold for sex reversal, we predict that the frequency of sex reversal will be independent of incubation temperature across the species range.

| Specimen collection and phenotypic sex identification
A total of 534 wild-caught dragons, collected from field trips (n = 337 dragons) and supplemented by vouchered museum specimens (n = 197) collected over a period of 38 years from 1980 to 2018, were K E Y W O R D S climate change, environmental sex determination, sex reversal, threshold evolution sampled in this study (File S1). Tissues from dragons collected during field trips were taken as either tail or toe clips, frozen blood or dried blood on DNA-preserving card storage systems (Whatman TM FTA TM Elute Cards, Macherey-Nagel Nucleocards, or PerkinElmer 226 Five Spot cards). Liver samples from museum specimens were sampled by museum staff at the time of collection, and either stored in ethanol at room temperature or frozen (−80°C). Specimens that were wildcaught and released in the field were phenotypically sexed by hemipenal eversion. In this approach, the hemipenes of adult males are reliably everted by running the thumb up the ventro-lateral surface of the base of the tail, while gently bending the tail to expose the cloacal aperture. Sex identification was confirmed secondarily by examining the animal for gravidity, observing nesting behaviour, and examining male secondary sexual characteristics, such as larger head size and dimorphism in femoral pores. The sex of museum specimen and road-killed animals was identified by dissection of the abdomen and inspection of the gonads. The few juvenile or hatchling animals captured or located during the study were only included if gonadal sex could be determined by dissection.

| Molecular sex identification
Genotypic sex (ZZ/ZW) was determined by extracting DNA from tissue and conducting a PCR sex test (Quinn et al., 2010, modified by Holleley et al., 2015. Three DNA extraction methods were performed as appropriate for each tissue type: (a) Qiagen's Gentra® Puregene® DNA purification kit for tail and toe clips, frozen liver, frozen blood, (b) Macherey-Nagel Nucleospin® Tissue kit for blood on frozen PerkinElmer 226 cards and Macherey-Nagel Nucleocards; (c) Whatman Elute quick extraction protocol for blood on Whatman TM FTA TM Elute Cards. Manufacturer's instructions were followed in all cases. Purified DNA was quantified using a NanoDrop 1,000 spectrophotometer (Thermo Scientific) and standardized to a concentration of 25 ng/µl prior to PCR. DNA extraction and PCR setup were automated using the epMotion 5,075 platform (Eppendorf, Germany), and included both positive and negative controls following Whiteley et al., (2017). PCR products were visualized on a 1.5% agarose gel. ZZ genotypic individuals displayed one band at 524 bp (both Z fragments) and ZW genotypic individuals displayed two bands, one at 524 bp (Z fragment) and the other at 357bp (W fragment). Internal ZZ and ZW control samples were run on every gel. Individuals with a ZZ genotype and a female phenotype were classified as sex-reversed. Chi-square tests were used to determine if there were sex biases in this dataset between males and females across the entire period of the study. Additionally, chi-square tests were conducted to determine if the proportion of concordant females and sex-reversed females varied across the years in which at least 15 phenotypic females were collected.

| Prediction generation: constant temperature equivalent
The rate of embryonic development of P. vitticeps is dependent on temperatures experienced by the egg within the nest, with embryos developing faster at hotter temperatures (Whiteley et al., 2017).
For embryos developing at higher temperatures, this means that We predict that sex reversal will occur most often in areas where incubation conditions are highest, and therefore expect to see a latitudinal gradient in the occurrence of sex reversal, with more cases of sex reversal occurring in the north of the species range [Correction added on 7 December 2020, after first online publication: the figure has been modified.] the period during which sex can be determined by temperature is shorter. Constant temperature equivalent (CTE) calculation allows the fluctuations in temperature experienced by a typical buried egg to be converted to a constant temperature which results in the same amount of embryonic development. First, the daily soil temperature was calculated over a single laying season (September to December in year y and January to February in year y + 1) using NicheMapR (M. Kearney, 2020; M. R. Kearney & Porter, 2017) and user-supplied weather input data (SILO climate gridded dataset). Secondly, the daily CTE was calculated using methods described in Georges (1989) and Georges et al. (1994). To generate range-wide long-term trends in CTE, the procedure was repeated over forty consecutive breeding seasons from 1975-2014, and for 1,000 locations randomly selected from within the range of P. vitticeps.
The microclimate model was driven by the Scientific Information for Land Owners (SILO) climate gridded dataset (Jeffrey et al., 2001).
For all the 1,000 randomly selected points, we extracted minimum and maximum temperature, relative humidity at time of minimum and maximum temperature, rainfall, vapour pressure, and solar radiation (used to estimate cloud cover, as outlined in NicheMapR documentation). In addition, we used a daily gridded mean wind speed dataset at both the minimum and maximum wind temperature (T. Viscarra Rossel et al., 2014aRossel et al., , 2014bRossel et al., , 2014cRossel et al., , 2014d. The proportion of clay, silt and sand at each location was used to determine the soil parameters based on table 9.1 from Campbell and Norman (1998), included in NicheMapR (M. Kearney, 2020). We used a custom R program (R Core Team, 2020) including the packages sp (Bivand et al., 2013;Pebesma & Bivand, 2005), Rfast (Papadakis et al., 2020), zoo (Zeileis & Grothendieck, 2005), timezone (Rundel, 2013) and raster (Hijmans, 2020) to collate climate variables into NicheMapR, run the microclimate model and extract the hourly soil temperature output for 10 soil depths between 0 and 200 cm. We assumed a nest depth of 15 cm (Pianka, 2005), and that nests were unshaded.
To calculate the daily CTE, we used the method described in Georges (1989) and Georges et al. (1994). The model assumes that daily soil temperature varies on a cycle about a stationary mean (equation 1), and as such provides a daily estimate of nest temperature corrected for fluctuations. Given the hourly soil temperatures for one day, T is the temperature at time t in degrees Celsius, M is the mean soil temperature, and R is the maximum deviation from M.
T 0 is the minimum temperature at which embryonic development is possible. In P. vitticeps, T 0 = 16.2°C (unpublished data). The CTE (T') occurs at time t'. Equation 2 is solved iteratively for t' until t'o ldt'n ew < 0.0001. If this did not occur within 10,000 iterations, we used the value of t' that occurred when the difference between iterations was smallest. The value for t' is then substituted into equation 1 and solved for T'.
If soil temperature dropped below T 0 , development ceased for part of that day. In that case, t 0 was calculated using equation 3 to estimate the time at which soil temperature drops to T 0 and development ceases. Then, t' was calculated iteratively using equation 4 until t'o ld -t'n ew < 0.0001. If this did not occur within 10,000 iterations, we used the value of t' that occurred when the difference between iterations was smallest. Finally, the CTE was calculated by substituting t' into equation 1. For further details on the derivation of these equations, please refer to Georges (1989).
For some days, there was little variation in soil temperature and M = T 0 , or M = T 0 and R = 0. On these days, the CTE could not be calculated using the method described in Georges (1989) andGeorges et al. (1994). When this occurred, we set CTE = T' = M. We calculated the daily CTE for each day in September to February. The mean CTE for a breeding season is the mean of all daily CTEs from September to February (n = 181 days) for that random location.
The high degree of correlation (Pearson's Correlation; R = 0.890, t = 61.739, df = 998, p < .001) between the average maximum temperature for each location from 1975 to 2014 and the corresponding CTE estimates (derived in part from maximum temperature) led us to use only the extracted SILO data for estimation of inferred incubation conditions for each collected specimen rather than the more derived CTE values.

| Geographic correlates of sex reversal
Multinomial logistic regression (R packages: nnet and tidy) was used to determine whether sex-reversed individuals displayed different geographical distributions when compared to the concordant sexes (where genetic sex matches phenotypic sex). A multinomial logistic regression was constructed using latitude, longitude and elevation as the predictor variables, with the sex of each specimen as the response variable, either concordant male (ZZm), concordant female (ZWf) or sex-reversed female (ZZf). The three categorical outcomes were then compared on a pairwise basis, to determine if there was an association between latitude, longitude, or elevation and sex.
Cluster analysis software SaTScan v9.6 (2018) was used to identify if sex reversal was clustered in the landscape. A Bernoulli probability model and moving window analysis was applied to identify significant clustering in the spatial distribution of sex-reversed compared to all concordant individuals. The maximum spatial cluster size (3) was set at 50% of the size of the population, a generally accepted standard (Kim & Jung, 2017), with high abundance clusters defined as containing at least ten instances of sex reversal.

| Environmental correlates of sex reversal: inferred incubation conditions
Inferred incubation conditions (Table 1) were determined for each individual using its estimated age, location of capture and parameters of the reproductive cycle of P. vitticeps. Temperature variables (maximum temperature and diel range) were selected for their direct relationship to sex reversal (Holleley et al., 2015) and rainfall for its effects on temperature. The deviation of these conditions during the inferred incubation period from the long-term average for each site was also included in analysis, to account for the possibility that it may not be absolute temperature which causes sex reversal, but rather deviation from long-term conditions, to which populations in different areas may have become adapted. Age was estimated qualitatively, and by using snout-vent length (SVL). Adult individuals were assumed to be between 3-5 years old, juveniles between 1 and 2 years old and hatchlings less than 6 months old. The breeding season (and therefore time during which incubation conditions could reasonably be experienced by each individual) was determined to be from September-February inclusive based on patterns of gravidity in dissected wild female specimens and patterns of egg laying in a breeding colony ( Figure S1). For adults, inferred incubation conditions during the third, fourth and fifth breeding season prior to collection were averaged. For juveniles, inferred incubation conditions during the first and second breeding season prior to collection were averaged. Hatchlings were excluded from analysis.
Environmental data during the inferred incubation period of each individual were extracted from the Scientific Information for Land Owners (SILO) database of Australian climate data, in gridded format (~25 km 2 ). We needed to estimate the home range size to ensure that we could reasonably assume that each specimen was incubated within the grid in which it was collected.

| Population structure
To determine whether population structure across the species geographic range may explain patterns in sex reversal, reduced representation sequencing was conducted on a sample of 218 individuals (13 ZZf, 73 ZWf and 132 ZZm; Figure S2). Each individual was genotyped by reduced representation Illumina short-read sequencing using commercial provider Diversity Arrays Technology (DArT Pty Ltd, Canberra, ACT, Australia). Briefly, SNP genotyping was performed using a combination of complexity reduction by using two restriction enzymes, implicit fragment size selection and next-generation sequencing (Kilian et al., 2012). A pair of restriction enzymes PstI (recognition sequence 5'-CTGCA|G-3') and SphI (5'-GCATG|C-3') was used for complexity reduction by double digestion. Sequences were processed using proprietary DArT analytical pipelines (Kilian et al., 2012) to yield SNP markers polymorphic within the set of samples. Calling quality was assured by high average read depth per locus (medium coverage, 10X). In addition, approximately one-third of the samples were processed twice from DNA to allelic calls as technical replicates. Scoring consistency (repeatability) was used as the main selection criteria for high-quality and low error-rate markers. Refer to Georges et al. (2018) for further detail. Additional filtering of loci and individuals and preliminary analysis was undertaken using the R package dartR (v.1.1.11; Georges et al., 2018). Loci with a call rate of less than 95%, repeatability of less than 99%, monomorphic loci and secondary SNPs were removed prior to analysis.
Individuals with a call rate of less than 75% (n = 2) were also removed from analysis. Structure across the landscape was visualized using principal coordinates analysis (PCoA) in the dartR package, and isolation by distance was assessed using Mantel's test in the R package adegenet (Jombart, 2008).

| Geographic distribution of sex reversal
Among the 534 wild-caught individuals sampled across the species range, 294 (55%) were ZZ males, and 240 (45%) were phenotypic females. Of those females, 28 (12%) possessed a ZZ genotype and were therefore sex-reversed. The bias towards males in our data was significant (χ 2 = 5.46, df = 1, p = .019) possibly due to capturebias and the higher visibility of males which display prominently.
Sex-reversed females were collected in 12 non-consecutive years of the 38 years sampled. The earliest case of sex reversal was in South Australia in 1983 and the most recent in New South Wales in 2018 ( Figure 2). There was no range-wide temporal trend in sex reversal. The proportion of females which were sex-reversed varied from 6% to 27% but without significant difference between years (χ 2 = 6.031, df = 5, p = .303) among the six years in which at least 15 phenotypic females were collected.
Sex reversal was detected across a substantial proportion of the range of P. vitticeps but was not randomly distributed ( Figure 2). Spatial cluster analysis using SaTScan identified a single large cluster in the south-eastern part of the sampling area where sex reversal occurred at a rate higher than would be expected if sex-reversed individuals were distributed randomly throughout the study range. The radius of the cluster was 410.83 km, containing 23 out of the 28 cases of sex reversal (Expected cases: 12.9, Log likelihood ratio: 8.182, p = .039). We found that the latitudinal, longitudinal and elevational distributions of concordant males (ZZm) and concordant females (ZWf) did not differ (Table 2).
Sex-reversed individuals (ZZf) were found not to differ in terms of latitudinal or elevational distribution compared to both ZZm and F I G U R E 2 (a) The distribution of sex reversal across the range of Pogona vitticeps. Cluster analysis identified a single cluster in which sex reversal occurred at a higher frequency than expected under the assumption of a uniform spatial distribution (Expected cases: 12.9, Log-Likelihood Ratio: 8.182, p = 0.039). (b) Temporal trends in sex reversal are difficult to determine owing to uneven sampling effort, but the earliest sex reversal event was detected in 1983 [Correction added on 7 December 2020, after first online publication: the figure has been modified.] ZWf, but were significantly different in their longitudinal distributions, being largely absent from lower longitudes (the western part of the species range; Table 2; Figure 2).
In a previous report detailing the extent of wild sex reversal (Holleley et al., 2015), the area examined spanned a minimum convex polygon (MCP) of ~ 260,000 km 2 (12.2% of the known species range), and sex-reversed females were found to occupy an MCP of ~ 26,000 km 2 (9.9% of the sampled range). Here we report a substantial increase in the sampling area, with an MCP of ~ 1.7 million km 2 (79.5% of the species range) sampled, of which sex-reversed females were found in a ~ 400,000 km 2 MCP representing approximately 24.2% of the sampled range.

| Environmental correlates of sex reversal
Overall, the binomial logistic regression comparing ZZ males and ZZ females using all six inferred incubation condition predictor variables had exceptionally poor predictive ability (R 2 = 0.032) ( Table 3). An unbiased model inference approach using MuMIn extracted nine "best" models to explain the binomial logistic regression. These models differed by less than two AICc values from the best model with the lowest AICc value and are therefore indistinguishable in terms of their explanatory power (Burnham & Anderson, 2002). Five of these were single-variable models, three contained two variables and one was the null model. That the null model was also included indicates that no combination of inferred incubation conditions was able to explain the binomial regression model of the frequency of sex reversal (ratio of ZZ females to ZZ males) more parsimoniously than the null model. This suggests that of the inferred incubation conditions included in the model, none are capable of adequately explaining the occurrence of sex reversal. We were unable to detect any broadscale geographic population structure. In a PCoA analysis, the first and second axes explain only 2.4% and 1.3% of variation respectively (Figure 3), indicating that despite the large sampling area, there is limited population differentiation. A relatedness matrix constructed using dartR identified a higher degree of relatedness in all specimens collected within a geographically isolated area south of the Murray River in South Australia. In the PCoA analysis, these individuals do cluster and appear to be genetically divergent from the rest of the sampled specimens, but low variation explained by the PCoA indicates that this divergence is minor. There was significant evidence of weak isolation by distance (Wright, 1943) across the entire range of the species (Mantel's test; r = .571, p = .001, Figure 3) suggesting a stepping-stone model of dispersal for this species (Kimura, 1953).

| D ISCUSS I ON
Here, we present the most detailed in situ account of naturally occurring sex reversal in a terrestrial vertebrate. Over a 38-year period encompassing a study area of 1.7 million km 2 (almost 80% of the species range), we show that sex reversal in P. vitticeps occurred across a quarter of the species range but was not directly explained by inferred incubation conditions. Sex-reversed ZZ females comprise 12% of all phenotypic females collected, demonstrating that the sexreversed phenotype is a substantial demographic in this species.
Counter to our predictions based on the established relationship between high temperature and sex reversal in a laboratory setting (Holleley et al., 2015;Quinn et al., 2007), we did not observe an association between environmental temperatures and sex in the wild. There is a noticeable lack of sex reversal in the northern and north-western part of the species range where nest CTE and ambient temperatures are the hottest (Figure 1). Instead, we observed a significant cluster of sex reversal in the south-eastern portion of the species range, where inferred incubation conditions are relatively cooler than in the north, and where long-term CTE's do not approach the 32°C sex reversal threshold. It must be noted that while we cannot conclude definitively that sex reversal is absent from the north of the species range, owing to lower sample size in this area, our clustering methods do account statistically for sampling heterogeneity and still identified a significant cluster. The occurrence of sex reversal in wild populations is clearly more complex than the simple relationship with temperature observed under controlled laboratory conditions, even after correcting for diel variation in temperature in natural nests (Georges, 1989;Georges et al., 1994), and could be explained by differences in the thermal threshold for sex reversal.
We argue that the geographically disjunct distribution of sex reversal is most likely explained by local genetic adaptation to environmental conditions and evolution in the thermal threshold for reversal. Such a mechanism is involved in transitions between GSD and TSD (Q uinn et al., 2011), and our population genetic data demonstrate that isolation by distance is the primary pattern of genetic differentiation in this species, providing a pattern of underlying genetic variation. In contrast to other reptiles with very large species distributions, P. vitticeps has no obvious barriers to dispersal and no significant population structure (Atkins et al., 2019;Melville et al., 2001Melville et al., , 2017Sovic et al., 2016). The Murray River in South Australia was the only weak and porous barrier to dispersal that we detected. We propose that the thermal threshold for sex reversal varies in a gradient F I G U R E 3 Genetic structure in the central bearded dragon (Pogona vitticeps). The structure displays a pattern of isolation by distance, but no strong genetic differentiation between populations throughout the species range. Slight genetic differentiation is apparent among specimens collected south of the Murray River, a geographical barrier to gene flow. (a) Principal coordinates analysis indicates that there is little genetic differentiation between individuals from across the species range. Axes 1 and 2 explain only 2.4% and 1.3% of variation, respectively, and this combined with the lack of clear delineation of groups of individuals, indicates that gene flow is occurring between groups of individuals across the species range, resulting in a geographical gradient of genetic differentiation. (b) Significant and relatively strong isolation by distance was detected (Mantel's test: r = 0.571, p = 0.001), indicating that across the species range, P. vitticeps genetic similarity decreases with increasing distance. PCoA and isolation by distance analysis were generated from a filtered SNP dataset of 11,128 binary SNPs and 216 individuals selected from across the species range [Correction added on 7 December 2020, after first online publication: the figure has been modified.] pattern across the landscape, matching the pattern of isolation by distance. We predict that threshold temperatures will be higher in the northern part of the species range, where incubation temperatures are the hottest.
For evolution of the thermal threshold to occur, heritable variation in the propensity to reverse must exist. Preliminary observations suggest that this is the case in P. vitticeps. Offspring of sex-reversed P. vitticeps have a lower threshold for sex reversal than the offspring of ZW females (Holleley et al., 2015). Additionally, individual females can produce offspring that are completely resistant to sex reversal at all temperatures (Holleley et al., 2015). There is also likely to be an epigenetic component to the sex reversal threshold. Epigenetic mechanisms already implicated in sex reversal include differential expression and temperature-specific splicing of chromatin modifying genes JARID2 and JMJD3/KDM6B (Deveson et al., 2017;Ge et al., 2018).
Compensatory maternal choice in nesting phenology, depth and location cannot fully explain the distribution of sex reversal in P. vitticeps. While we acknowledge that uncertainty in nest depth, location and specimen age may contribute to low power in this analysis, another study in a TSD reptile found significant environmental associations at this scale with a similar sampling and analytical approach (Ewert et al., 2005). Additionally, studies in   (Schwanz et al., 2020). Here, a poor relationship between temperature and the incidence of sex reversal suggests that local genetic adaptation to extreme conditions has already occurred as a sex ratio balancing mechanism. Populations that resist sex reversal then become a source of ZW immigrants, and even small rates of ZW immigration (1% p.a.) significantly decelerate the transition to TSD (Schwanz et al., 2020). Thus, local thermal adaptation and an influx of ZW immigration could allow P. vitticeps to remain stably as a mixed sex determination system in the long term.
Importantly, the existence of multiple evolutionary processes working in opposition means that the outcome of transitions between sex determining modes cannot be predicted based on climatic data alone (Schwanz et al., 2020).
The capacity of thermal thresholds to evolve is a critical component of the biological response to climate change, particularly in ectothermic animals. We demonstrate here that threshold evolution may have already occurred across the geographic range of a species with temperature sex reversal. Understanding this component of resilience to changing thermal regimes is essential to future conservation efforts.

ACK N OWLED G EM ENTS
For their assistance with museum collection access and tissue loans,

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The Scientific Information for Land Owners data used in this paper can be accessed at https://www.longp addock.qld.gov.au/silo/, and other sources of data were indicated. The associated collection and inferred incubation data for each specimen and CTE calculations for Stephen D. Sarre https://orcid.org/0000-0002-7158-2517 Clare E. Holleley https://orcid.org/0000-0002-5257-0019