Effects of natural nest temperatures on sex reversal and sex ratios in an Australian alpine skink

Altered climate regimes have the capacity to affect the physiology, development, ecology and behaviour of organisms dramatically, with consequential changes in individual fitness and so the ability of populations to persist under climatic change. More directly, extreme temperatures can directly skew the population sex ratio in some species, with substantial demographic consequences that influence the rate of population decline and recovery rates. In contrast, this is particularly true for species whose sex is determined entirely by temperature (TSD). The recent discovery of sex reversal in species with genotypic sex determination (GSD) due to extreme environmental temperatures in the wild broadens the range of species vulnerable to changing environmental temperatures through an influence on primary sex ratio. Here we document the levels of sex reversal in nests of the Australian alpine three-lined skink (Bassiana duperreyi), a species with sex chromosomes and sex reversal at temperatures below 20 °C and variation in rates of sex reversal with elevation. The frequency of sex reversal in nests of B. duperreyi ranged from 28.6% at the highest, coolest locations to zero at the lowest, warmest locations. Sex reversal in this alpine skink makes it a sensitive indicator of climate change, both in terms of changes in average temperatures and in terms of climatic variability.

Climatic data. Thermal data-loggers (iButton® model DS1921G, Maxim Integrated, San Jose, CA, USA; accuracy ± 1 °C from − 30 °C to + 70 °C; dia 17 mm, height 6.2 mm, mass 3.2 g) were placed to measure the soil temperature at the surface, 10 cm and 20 cm depth at each field location. For each location, daily projections of climatic data were obtained from a publicly available spatial data repository (SILO, Scientific Information for Landowners 26 constructed from records provided by the Australian Bureau of Meteorology, Canberra. The data were for the period 1889 to present at a grid resolution of 0.05° latitude by 0.05° longitude (approximately 5 km × 5 km) 26 .
Timing of breeding season and fieldwork. Female B. duperreyi lay clutches of eggs in communal nests (Fig. 1c). in late November to early January each year depending upon seasonal conditions. Nests are typically found under rocks in exposed areas subject to high solar radiation [27][28][29] . Therefore, we started our fieldwork in the first week of November to locate nests as soon as possible after laying (typically 3-7 days incubation before discovery) in both open and forested areas in both seasons. We searched thoroughly for nests from November to mid-January in each year.
Nest characteristics and nest temperature. Each nest was marked with a plastic flag and had the GPS locations recorded. The flag was temporarily removed, and we recorded minimum nest depth (to top of the www.nature.com/scientificreports/ shallowest egg) and maximum nest depth (to the bottom of the deepest egg). If nest depth was less than 10 cm, one thermal data-logger (iButton® model DS1921G, Maxim Integrated, San Jose, CA, USA; accuracy ± 1 °C from − 30 °C to + 70 °C; diameter 17 mm, height 6.2 mm, mass 3.2 g) was placed in the core of the nest; if nest depth was greater than 10 cm, two iButtons were placed immediately above and below the egg mass of each nest. The iButtons were factory calibrated and were set to record the temperature at hourly intervals throughout the incubation period. We monitored nests weekly (9 weeks) for their condition, except at Dartmouth (7 weeks Egg collection, incubation and sample collection. After 9-10 weeks of development in the field (approximately 90% of the incubation period), we removed 415 eggs from 42 randomly selected nests and transferred them to plastic boxes in which they were buried in moist vermiculite (4 parts water to 5 parts vermiculate by weight). Each egg was separated by plastic partitions. The eggs were transported in a portable incubator set to 23 °C with high but unmeasured humidity. All experimental protocols were conducted with permission and in accordance with the procedures of Animal Ethics Committees at the University of Canberra and the CSIRO. Eggs were weighed with an electronic balance (± 0.01 g), and egg lengths and widths were measured using digital vernier callipers (± 0.01 mm). Eggs were incubated 23 °C (± 0.5 °C), which typically produces a balanced sex ratio 15 ; this was a precaution only, because the eggs were harvested after the temperature is likely to exert an influence on offspring sex 15 . The boxes were gently rotated inside the incubator every couple of days, and hatchlings (hatching success was 95.2%) were removed as soon as they emerged from the egg. Sex was identified by manually everting the hemipenes of males 15,31 . 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 90% ethanol. This the study was conducted and is reported in accordance with ARRIVE guidelines (https:// arriv eguid elines. org).
Molecular detection of sex reversal. DNA was extracted from tail tips using a Gentra Puregene commercial kit (Qiagen Science, Maryland, U.S.A.) following manufacturer protocols; DNA was extracted from blood samples following manufacturer 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, N.S.W., Australia). The genotypic sex was identified using a PCR test based on seven Y-specific markers 32 . Briefly, in applying the test we used 1 × MyTaq™ HS Red mix (Bioline U.S.A. 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. This was followed by 30 cycles of 95 °C denaturing (20 s), 65 °C annealing and 72 °C extension (10 min). PCR products were visualised 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 recognised as XY individuals, whereas as the samples for which a band was not amplified in all seven markers were recognised as XX individuals. The seven markers always concurred in their identification of genotypic sex, as did they in the original study published by Dissanayake et al. 32 . False negatives arising from recombination events are thus highly unlikely as they would present as some but not all markers detecting the presence of a Y chromosome. No XY females were observed, another indication that recombination and/or mutation involving these loci is negligible and has not affected the accuracy of genotypic sex assignment. Phenotypic male lizards showing genotype-phenotype discordance were classified as sex-reversed 21,32 . All molecular sex tests were conducted blind to the phenotypic sex of the individuals.

Analysis.
Results are presented as means ± standard errors unless otherwise indicated. Correlation and regression analysis were used to describe relationships between variables, and the Student's t-test used to compare mean weekly air and nest temperatures. To analyse the relationship between natural thermal regimes and sex ratio of B. duperreyi, we corrected nest temperatures to constant temperature equivalents (CTE) 33,34 using the reaction norm for development rate against temperature estimated using the developmental model developed by Dallwitz and Higgins 35 ( Fig. 2; Figure S1A). Sexual fate is sensitive to temperature in a wide range of reptiles during the middle third of incubation 33,34 . This has not been established for Bassiana duperreyi, but we nevertheless restricted our modelling to the middle third of development. We estimated the middle third of development by summing development as a function of temperature in small but finite increments, back from the point of hatching. The likelihood of reversal was based on two criteria: raw temperature and temperature corrected for developmental rate (CTE). Data analyses were performed in R (R Core Team, 2017) and GraphPad Prism version 9 for Windows (GraphPad Software, La Jolla California USA).
Simple linear regression was used to determine the relationship between the frequency of sex-reversed hatchlings in each population with respect to elevation and to characterise the trends in mean nests temperature, yearly mean maximum (Tmax), and yearly mean minimum temperature (Tmin) against elevation in each field location. 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 vapor pressure (hPa)] including elevation recorded at each field site. We used a t-test to compare the sex reversal frequency of adult males (Dissanayake et al. 21 ) with the sex reversal frequency our nests.

Results
Climate data. Thermal profiles showed considerable diel variation at each field location ( Circus: F 1,517 = 539.6, P < 0.0001, R 2 = 0.51; Cooma: F 1,517 = 537.9, P < 0.0001, R 2 = 0.50; Dartmouth: F 1,517 = 627.5, P < 0.0001, R 2 = 0.54). Air temperatures were consistently cooler at higher elevational locations than the lower elevations, which inversely correlated with elevation when B. duperreyi eggs were incubation (F 1,22 = 4.795, P < 0.039, R 2 = 0.17). The weekly mean Tmax and Tmin showed that temperatures fluctuated substantially during the 9-week eggs incubation period (Fig. 4). The highest mean rainfall events (17.0 ± 23.87 mm) were recorded at the Piccadilly Circus in the first week of egg incubation period in the first season. The highest mean rainfall events (5 ± 11.85 mm) were recorded at Mt Ginini during the third week of B. duperreyi egg incubation. The highest rainfalls were recorded during the fourth week of the egg incubation period at Cooma (5 ± 6.3 mm) and Dartmouth ( (Fig. 4). Nest depth was inversely related to elevation (F 1,86 = 39.80, P < 0.0001, R 2 = 0.32; Figure S4). . The nests showed significant differences in mean nest temperature maxima in all locations except Dartmouth, but mean nest temperature minima shows a significant difference in all locations (Table 1).
Mean weekly nest temperature was correlated with mean weekly air Tmax (R 2 = 0.42-0.74, P < 0.05) and Tmin (R 2 = 0.70-0.93, P < 0001) at all field locations. The significant warming trend was recorded during the incubation weeks in first season at the Piccadilly Circus (F 1,300 = 31.74, P < 0.001, R 2 = 0.09), but a significant cooling trend was recorded in the second season (F 1,124 = 52.63, P < 0.001, R 2 = 0.29). In the second season a significant cooling trend was also recorded at Mt Ginini (F 1,79 = 132.3, P < 0.001, R 2 = 0.62). Whereas at Cooma and Dartmouth showed no significant trend as the season progressed.
Mean nest temperatures during the incubation period were inversely correlated with elevation (F 1,2 = 41.71, P < 0.05, R 2 = 0.95) (Fig. 6a). The highest and lowest mean daily CTE were recorded at Dartmouth (30.35 ± 0.12 °C) and Mt Ginini (26.2 ± 1.98 °C), respectively. The mean daily CTE was significantly inversely correlated with elevation (F 1,74 = 11.39, P < 0.001, R 2 = 0.17). When the CTE dropped below the 20 °C (the threshold for reversal, Shine et al. 15 ) during the thermosensitive period, for even for a short time during the incubation period, sex reversal was observed (Fig. 7 and Figure S1 B-F).  www.nature.com/scientificreports/ Nest depth was not significantly associated with mean nest temperature (F 1,84 = 0.081, P = 0.77) at any location. Examination of mean nest temperatures suggested a slight warming trend at Piccadilly Circus during last 22 years, but this was not statistically significant (F 1,8 = 3.49, P = 0.09). Nests that produced sex-reversed hatchlings had mean daily nest temperatures cooler than nests that did not produce sex reversed hatchlings. However, mean weekly temperature regimes for these two categories of nests were not significantly different at any field location (Table S1).
Genotypic sex identification and the frequency of sex reversal. During the study period, we collected a total of 415 eggs from the natural nests. The hatching success rate was 95.2%. Therefore, a total of 395 hatchlings were able to be phenotypically identified. Of them, 262 (66.3%) were phenotypic males and 133 (33.6%) were phenotypic females. The sex ratio (  21 (chi-squared test with Yates correction, χ2 = 4.05, df = 1, P < 0.05). A total of 59 (7.03%) phenotypically male hatchlings were sex-reversed. The highest frequency of sex reversal in hatchlings was recorded at the highest elevation site (28.6%, Mt Ginini, 1,640 m a.s.l.), and zero sex reversal was observed at the lowest elevation (Dartmouth, 380 m a.s.l) (Fig. 1a). This observation concorded with the previous study has been contacted for adult individuals (see Dissanayake et al. 21 ). The frequency of sex reversal was positively correlated with elevation (F 1,3 = 41.71, P < 0.05; R 2 = 0.95) (Fig. 6b) and each population showed a negative correlation with mean nest temperature (Pearson's correlation coefficient r = − 0.95, P < 0.05); Tmax (R2 = 0.99, P < 0.05) was significantly negatively correlated with sex reversal frequency (Table S2). When we compared the current study with Dissanayake et al. 21 for adult sex reversal frequency, the frequency of sex reversal in the nest was higher than the sex-reversed adult in the same field locations, but not significantly (P = 0.36). However, both hatchlings and adults (Dissanayake et al. 21 ) rates of sex reversal frequency positively correlate with their respective elevation (F 1,2 = 15.17, P < 0.005; R 2 = 0.77) ( Figure S5) (see also Dissanayake et al. 21 ).

Discussion
Over the last two decades, many field studies and laboratory experiments have been conducted to understand Bassiana duperreyi nesting ecology and phenotypic plasticity 29,30,[37][38][39][40] . However, the consequences of cold nest temperatures on sex reversal in B. duperreyi have been little studied in the wild (but see Holleley et al. 18 ). We show that sex reversal in natural nests of this species is common at elevations above 900 m a.s.l. The frequency of sex reversal in nests ranges from 28.6% at the highest, coolest location (Mt Ginini, 1,640 m a.s.l., mean nest www.nature.com/scientificreports/ temperature 19.7 ± 0.62 °C) to zero percent at the lowest, warmest mean nest temperature recorded location (Dartmouth, 380 m a.s.l., 25.6 ± 1.25 °C). Thus, both hatchlings (this study) and adults as observed by Dissanayake et al. 21 exhibit a comparable relationship between the frequency of sex reversal and elevational variation in ambient temperatures. These substantial levels of sex reversal in nests and the adult population, derived from the rates of sex reversal we have observed in the nest, indicate that the sex-reversed phenotype is an important component of the demography of this species. How sex reversal will come to influence demographic processes that govern the persistence of local populations in the Australian high country requires additional information on the fertility of sex-reversed individuals. The presence of sex-reversed adults 21 at a similar rate that observed in the nest (this study) establishes the viability of sex-reversed individuals. But are they fertile? If the B. duperreyi sex reversal yields infertile male individuals, then males of the species are subject to latent mortality-a significant component of the population could comprise viable but infertile males. As these are generated by reversal of XX individuals otherwise destined to be females, the effective population size, which depends on reproductive female number, is drawn down numerically. It is also drawn down because many of the remaining XX females will potentially mate with the infertile males to no effect, and this proportion of females will not contribute to effective population size. Together, this could lead to a rapid population decline in effective population size and local extinction should there be a sequence of episodic cold seasons or cold spells at a time when the embryos in the nests are thermosensitive. In contrast, if the sex-reversed males are both viable and fertile, then sex reversal under cooler temperatures will lead to an overproduction of males and, in the context of other evolutionary and phenotypic responses, Fisher's frequency-dependent selection 41,42 may be invoked in support of traits that bring the population sex ratio back to equilibrium. This will manifest as a reduction in the sensitivity of sex to temperature or reduction of the frequency of the Y chromosome in the population or both and, potentially, its loss 7,21,43,44 . Suppose for example, a future change in climate or annual weather is sudden, such that the species is subject to substantial sex reversal but unable to evolutionarily adjust the male-biased sex ratio. In that case, effective population size will decline, again putting local populations at greater risk of extinction. Clearly, under either scenario, sex reversal has profound implications for the demography of B. duperreyi at the high elevation distributional limits of its range whether or not sex-reversed males are fertile.
The above impacts of sex reversal can be ameliorated if the frequency of sex reversal can be constrained by phenotypic responses. In particular, nest site choice and nesting phenology are potentially both important in ameliorating the frequency of sex reversal. We recorded that B. duperreyi selected open grassland and 98% of rock substrate at all field locations for laying their eggs, similar to the findings of previous studies 29, 30 . We show that the consequence of egg laying in exposed sites is a high diel variance in nest temperature. B. duperreyi typically select nesting sites under rocks compared to other alpine skinks that live at the same field location that typically nest in logs and rotting organic material 29 . Furthermore, B. duperreyi appear to select nest sites on both thermal averages and high thermal variability both of which have measurable effects on hatchling phenotype 40 . Adjusting for diel fluctuations, we found that when the CTE dropped below the 20 °C (the threshold for reversal, Shine et al. 15 ) during the thermosensitive period, for even for a short time during the incubation period, sex reversal occurred (Fig. 7). Therefore despite the observed temperature fluctuations in nests, extreme cold temperature events could lead to 100% sex reversal of the XX genotype in the alpine populations at the highest elevations.
When considering the alpine populations, B. duperreyi embryo development occurred under high thermal fluctuations and some nests reached beyond the embryo physiological limits for brief periods. The high level of diel fluctuation of nest temperature leads some nests reaching both upper and lower limits in development temperatures. At Mt Ginini (i.e., highest elevational location for this study), four nests show that their temperature exceeded both the maximum and minimum temperature for embryo survival during the incubation period (see Table 1). However, we observed that these nests survived and successfully completed their natural embryonic development to hatchling. Therefore, a short period of exposure to extreme temperature events did not affect egg survivorship but may well have influenced the frequency of sex reversal.
The embryonic survivorship under the current increasing global temperature trend is a significant challenge for oviparous species globally. The current trend in global warming inescapably will increase the nest temperatures of many oviparous reptiles 30,46,47 . Therefore, natural selection should favour potential nesting sites, nest phenology and other maternal behaviors that enhance embryonic survivorship. Egg laying of B. duperreyi shifted temporally, with the nesting season starting from mid-November to early December, confirming a similar observation by Telemeco et al. 30 . However, in our second season (2018/19), females initiated egg laying in early January, as they did in 1968 27 . B. duperreyi appear to shift their egg laying from year to year, presumably in response to variation in natural conditions conducive to successful incubation and emergence (see also Pengilley 27 ). Incubation conditions, and so timing or egg laying, also have strong effects on offspring phenotypes 15,48 and sex determination, including sex-reversal. In addition to altering the timing of nesting, B. duperreyi has progressively been digging deeper nests over the last two decades (Fig. 5 and see Telemeco et al. 30 ), likely to reduce exposure to extreme daily temperatures and have a complex impact on average nest temperatures depending upon soil composition and structure.
We have shown a significant relationship between the frequency of sex reversal in hatchlings and nest temperatures during their natural incubation period. Our work establishes and quantifies current elevational trends in the frequency of sex reversal in response to cold temperatures both in the nests and in the adult population 21 . This forms a baseline for examining the effects of climate change, and in particular climate warming. On the basis of the trends in primary (egg) and operational (adult) sex ratio with elevation, presumably driven by associated variation in thermal conditions, we expect to see populations without sex reversal, that is, governed by GSD, to increase in elevation as the climate warms. At the highest altitudes, if climate change is accompanied by increased variability, this will potentially disrupt the demographic processes on a local scale, leading to local population crashes and local extinction events, driven by fluctuating frequency of sex reversal. Temperature induced sex www.nature.com/scientificreports/ reversal in B. duperreyi and its demographic consequences likely make this species a very sensitive indicator of climate change. Monitoring changes in the frequency of sex reversal through time along our elevational gradient and at the species highest distributional limits is recommended as one option for assessing the impacts of climate change on the biota of the Australian high country.

Data availability
All data and materials are presented in the main paper and supplementary information. The Scientific Information for Land-Owners' data used in this paper can be accessed at https:// www. longp addock. qld. gov. au/ silo/. The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.