Models of primary sex ratios at a major flatback turtle rookery show an anomalous masculinising trend

Abstract

Background

Quantifying primary sex ratios is essential for assessing how global warming will influence the population dynamics of species with temperature-dependent sex determination (TSD). Process-explicit (mechanistic) models can accurately estimate primary sex ratios but require the resolution of the key physiological parameters that influence sex determination and validation of the model by testing predictions against empirical data.

Results

To address these goals, we conducted incubation experiments on flatback sea turtle (Natator depressus) embryos from a large winter-nesting rookery at Cape Domett in the East Kimberley region of Western Australia. A TSD model fitted to laboratory and field nest data indicated that the pivotal temperature producing equal sex ratios was 29.4°C, with males produced below 27.7°C and females produced above 31.1°C. Back-switch experiments revealed that the thermosensitive period (TSP), when gonads differentiate into testes or ovaries, occurs between 43% and 66% of development to hatching. Integrating this new information with sand temperatures reconstructed from 23 years of historical climate data shows that male-biased sex ratios are likely if the TSP falls during the Austral winter. Annual variation in the simulated sand temperatures increased from 1990 to 2013, with cooler winters producing conditions that favoured male hatchlings for longer periods. The same model projected to 2030 and 2070 suggests that female-biased primary sex ratios will become more prevalent over time.

Conclusions

Our results show that accurate modelling of primary sex ratios depends on quantifying the thermal biology of embryos and on parameterising mechanistic models of sand temperatures with site-specific climate data.

1Background

In most animals, sex is determined at conception by inherited sex chromosomes [1] but in many reptiles, the sex of an embryo is determined by the temperatures experienced during incubation (temperature-dependent sex determination or TSD) [2]. Primary sex ratios are of particular interest in the context of climate change, as a warming climate will increase nest temperatures and could create or exacerbate existing sex ratio biases in hatchling cohorts. Persistent biases in the primary sex ratio can potentially lead to demographic collapse or localised extinction, as has recently been shown in an island population of tuatara [3].

Marine turtles are a major lineage of reptiles, and all extant species have TSD. In each species, females are produced at temperatures above a threshold and males below [4]. There is an increasing focus on modelling primary sex ratios in these taxa [5]. In part, this is due to the challenges of fieldwork at remote rookeries and also to difficulties in sexing hatchlings using non-invasive methods [6]. Moreover, primary sex ratios vary within and across populations [7], and measuring sex ratios at one site, or over a single nesting season, will not accurately reflect the primary sex ratio of a population over time. Methods for accurately predicting hatchling sex ratios in the complex thermal environment of natural nests [8] are needed to reduce the need for broad-scale collection of hatchings from natural nests and for predicting sex ratios under future climates.

Earlier attempts to predict hatchling sex ratios in reptiles with TSD were based on correlations between sex ratios and air or nest temperatures (e.g. [9]–[14]). Correlative models have a key constraint in that they should not be extrapolated to predict sex ratios when temperatures exceed the range used to generate the model [8], and this constraint can limit their utility for predicting the impacts of global warming. More recently, mechanistic (or ‘process-explicit’) models have been used to estimate nest temperatures for species with TSD [8],[11]. A further innovation has been to link predicted nest temperatures to a developmental model of an embryo, where key traits such as hatchling sex, development time and heat-induced mortality can be estimated after quantifying the physiological responses of an embryo to temperature [8]. An advantage of the mechanistic approach is that behavioural responses to counter the effects of global warming, or processes such as metabolic heating, can be simulated by adjusting model inputs or constructs.

Both correlative and mechanistic models require knowledge of the pivotal temperature (Tpiv—the temperature that produces a 1:1 sex ratio) and the transitional range of temperatures (TRT) where both sexes can be produced for the population under study. Mechanistic models further require estimation of the rate of embryonic development across the broad range of current and future environmental temperatures, as well as identification of the period of development where indifferent gonads are sensitive to temperature, known as the thermosensitive period (TSP). In the reptile species so far examined, the TSP broadly falls within the middle half to middle third of incubation, and in most species, it has not been resolved any further [2],[15]. Delineating the TSP more precisely allows the identification of the particular portion of the thermal profile of a nest that influences sex determination. In general, a subset of nest temperatures will be required for sex ratio estimation, and selecting the correct subset enhances the accuracy with which nest sex ratios can be predicted.

The flatback turtle (Natator depressus, Garman) is endemic to Australian waters and is listed as ‘data deficient’ by the International Union for the Conservation of Nature [16],[17]. While pivotal temperatures have been estimated for both Western and Eastern Australian populations of N. depressus[18]–[20], other TSD parameters are poorly resolved. Cape Domett in the East Kimberley region of Western Australia hosts one of Australia’s largest N. depressus rookeries where females have an extended nesting season from April to November that peaks in the Austral winter [21]. Most other populations of N. depressus show a nesting peak in summer; hence, a winter nesting population provides an important contrast when considering the relative vulnerability of flatback populations to climate change.

In this study, we designed laboratory experiments on N. depressus embryos to identify the timing and length of the TSP and to estimate the Tpiv and TRT for the Cape Domett population. With these key parameters resolved, we applied the modelling package NicheMapR [22] to reconstruct 23 years of sand temperatures at Cape Domett (Figure 1). We then applied a physiological model of embryonic development rate and sex based upon our new data (Figure 1) to estimate how primarily sex ratios would have varied within and between years. Finally, by forcing NicheMapR with higher air temperatures projected for 2030 and 2070, we expected to show that hatchling sex ratios at this major rookery will become increasingly feminised.

Figure 1

Framework for modelling hatching sex ratios, showing the relationships between the microclimate and physiological models. Bold text indicates parameters or functions developed in this study for N. depressus. AWAP refers to interpolated historical climate data sourced from the Australian Water Availability Project.

2Results

2.1 Incubation times and embryonic development rate

Of 296 eggs incubated in the laboratory under a range of thermal regimes (see ‘Methods’ section), 156 hatched and almost all eggs that did not hatch showed no signs of development (confirmed by dissection). All 156 hatchlings plus 138 late-stage embryos sampled from field nests were sexed using histological techniques, and the classification of sex was highly repeatable (binomial test; p < 0.0001, probability of the same classification = 0.97).

Incubation times varied between 43 and 60 days (Additional file 1: Figure S1). Eggs incubated predominantly at 32.5°C had a mean incubation duration of 45.4 ± 0.2 days and hatched earlier (unpaired t-test, p < 0.0001) than eggs incubated predominantly at 28°C (57.0 ± 0.4 days). A non-linear development rate function fitted to average incubation times had the parameters b₁ = 2.39, b₂ = 18.03, b₃ = 33.35, b₄ = 6 and b₅ = 0.4 (root mean sum of squares = 3.70; Figure 2).

Figure 2

A non-linear development rate function estimated for the Cape Domett population of N. depressus . The function is fitted to ten constant temperature incubation records at three similar temperatures (open circles). Slight differences in temperatures are due to variation between incubators and containers.

2.2 Thermosensitive period (TSP)

The reduced viability of eggs due to freight delays (see ‘Methods’ section) impacted our ability to interpret results of back-switch experiments designed to identify the thermosensitive period (Table 1). Consequently, data from the four containers in each treatment (two replicates from each of two incubators) were pooled to provide an adequate sample size for interpreting the effects of each back-switch treatment. Pooled data are described hereafter, but data from individual containers in the back-switch experiments are summarised in Additional file 2: Figure S1.

Table 1 Design of the back-switch experiments used to delineate the thermosensitive period

The controls at 28°C and 32.5°C (treatments 6 and 12) produced 92% male and 100% female hatchlings, respectively (Additional file 2: Figure S1). Female–male–female switches (32.5°C ➔ 28°C ➔ 32.5°C) produced 100% female hatchlings in most treatments (treatments 7, 8, 9 and 11), with the exception of treatment 10 where one of 12 hatchlings was male. Conversely, most male–female–male switches (28°C ➔ 32.5°C ➔ 28°C) produced mixed sexes. Treatments 1 and 3 were switched to female-producing temperatures at 43% of development, thereby inferring a lower boundary of the TSP (Figure 3). Treatment 5, where the switch commenced at 66% of development, produced all male hatchlings (Figure 3), indicating that this switch period fell after the TSP. Hence, the upper boundary of the TSP for a female trigger was 66% of development. As it was not possible to delineate the TSP for a male trigger, the feminising TSP of 43%–66% was hereafter used for estimating primary sex ratios.

Figure 3

Back-switch experiments that were useful in determining the timing and length of the thermosensitive period. Dashed black lines are the target switch temperatures and solid grey lines are the average actual hourly incubation temperature for the four containers in each treatment. The cumulative proportion of development completed is indicated by dotted lines (right vertical axis). %F indicates the proportion of female hatchlings in a sample of n eggs.

2.3 Pivotal temperature (Tpiv) and transitional range of temperatures (TRT)

Using all available sex-ratio data from the field and laboratory, the best fitting TSD function estimated Tpiv at 29.4°C (Figure 4A). When we instead fitted a TSD function to the laboratory data and subsets of field data (see ‘Methods’ section), it produced no differences in the mean residuals (one-way ANOVA; p = 0.997), indicating that fits were effectively the same (Table 2). However, a function fitted to laboratory data and data from a particular subset of field nests (1, 4, 6, 10 and 15) produced the lowest Akaike information criterion (AIC) of 25.9 and the lowest mean residual (Table 2). This function had the parameters P = 29.4, S = −0.02 and K = 0.1, which defines the TRT as between 27.7°C and 31.1°C (Figure 4B). This TSD function was applied for estimating sex ratios following the process outlined in Figure 1.

Figure 4

The relationship between temperature and the proportion of male hatchlings for the Cape Domett population. (A) TSD model fitted using all available data (Hill equation; red line) and the consensus model fitted using a subset of the field data (black line). (B) Shows the consensus model and only those data used for model fitting (laboratory data = open circles; field data = closed circles), with the pivotal temperature of 29.4°C indicated by the dotted line and the boundaries of the TRT by dashed lines.

Table 2 Variation in the Tpiv and TRT based on data from different subsets of field nests

2.4 Validation of the sand temperature model

Temperature loggers (iButtons) distributed among 11 nests were retrieved from depths between 36 and 80 cm, and some were found clustered together and essentially provided replicate measurements (Additional file 1: Table S1). Thirty-two of 44 loggers deployed were found to contain a complete temperature record (other loggers failed or else recorded temperatures for only part of the incubation period). The diel temperature range decreased with depth, and this was reflected in both the empirical data and the modelled sand temperatures. On average, the difference in the mean measured nest temperatures relative to the mean of the modelled sand temperatures driven by weather data recorded at the rookery ranged from 0.35 to 1.30°C, with an average difference of 0.65°C ± 0.04°C. Sand temperatures modelled using data from the Australian Water Availability Project (AWAP)—a database of historical climate records for Australia [23]—differed from measured nest temperatures by 0.30°C–1.20°C, with an average difference of 0.66°C ± 0.05°C. In general, there was good agreement between actual and modelled sand temperatures (Figure 5D and E) but measured temperatures tended to be slightly lower than NicheMapR estimates near the start of incubation and higher than estimates towards the end of incubation (Figure 5A–C and Additional file 1: Table S1).

Figure 5

Examples of sand temperatures from September 18 to October 24, 2012 at Cape Domett. (A–C) NicheMapR predictions based on AWAP data (grey lines) or on weather station data (black lines) and the measured sand temperatures (red lines). Histograms of the r² values for the correlations between 32 iButton temperature traces and sand temperatures modelled using AWAP (D) and weather station; (E) data are also shown. Measured and modelled sand temperatures for an additional 29 nest/depth combinations are shown in Additional file 2: Figure S1.

2.5 Validation of the physiological model

Most (83%) of the 138 embryos sampled from field nests were female, but eight of the ten viable nests sampled contained embryos of both sexes. The sex ratios of embryos sampled from five field nests were predicted accurately when the physiological model (Figure 1) was based on measured nest temperatures, with the exception of nest 12 where mixed sexes were predicted but only female hatchlings were sampled from the nest (Table 3). Similarly, if sex ratios were predicted using sand temperatures generated using NicheMapR, they also generally agreed with the empirical data when NicheMapR was driven by weather station data (Table 3). However, the sex ratios predicted using sand temperatures generated when NicheMapR was driven by AWAP climate data agreed with the sampled sex ratios for only two of the five nests sampled (Table 3). In two of the three other cases, the physiological model estimated female hatchlings whereas hatchlings of both sexes were sampled from the nest.

Table 3 CTEs and sex ratios predicted for field nests based on measured and estimated sand temperatures

2.6 Historical sand temperatures and sex ratios

Driving NicheMapR with historical climate data (AWAP: 1990–2013) illustrates the annual variation in sand temperature at a typical nest depth of 50 cm (Figure 6). Annual maximum sand temperatures were stable over this period (linear regression analysis; p = 0.11 and r = 0.34), but minimum sand temperatures decreased (linear regression analysis; p = 0.017 and r = 0.49). Consequently, the annual range of sand temperatures increased between 1990 and 2013 (linear regression analysis; p = 0.004 and r = 0.57). In accordance with these trends, a greater portion of the year favoured the production of males in the most recent years (Figure 6).

Figure 6

Sand temperatures at 50 cm depth estimated for Cape Domett based on a 23-year hourly climate database. (A) 1990–1994, (B) 1995–1999, (C) 2000–2004 (D) 2005–2009 and (E) 2010–2013. Horizontal dashed lines show the bounds of the TRT and grey boxes indicate periods that would favour the production of males.

Contrasting a cooler year (2011) with a warmer year (2012) in detail (Figure 7) shows that in the former case, mixed sex ratios could be produced over a larger portion of the year. In the warmer year, sand temperatures mostly exceeded the upper boundary of the TRT (31.1°C); hence, females would be produced if the TSP fell within this period. The maximum sand temperature in each year was 34.0°C (December 2011 and November 2012; Figure 7). In both years, sand temperatures were below the lower boundary of the TRT (27.7°C) between May and September, indicating when male hatchlings would be produced.

Figure 7

Hourly sand temperatures at 50 cm depth predicted at Cape Domett for 2011 and 2012. Horizontal dashed lines show the bounds of the TRT for N. depressus, and grey or hatched boxes indicate periods that would favour the production of male or mixed-sex hatchlings, respectively. The upper limit of the shaded and hatched boxes indicates a potential upper limit of 35°C for embryonic viability. The red dashed bar shows the nesting season (the solid portion indicating peak nesting) [21], while the period when our fieldwork was conducted is shown by a horizontal black bar.

2.7 Sand temperatures and sex ratios under climate change

Adjusting NicheMapR inputs by increasing the 2007 AWAP air temperatures by the increments projected under 2030 and 2070 emissions scenarios shows that sand temperatures increase by a similar magnitude (Figure 8). Under a 2030 low emissions scenario, sand temperatures increase on average increase by 0.49°C relative to 2007 (Figure 8A). Sand temperatures predicted under the 2030 high emissions scenario and 2070 low emissions scenario are similar, with average increases of 1.22°C and 1.47°C relative to 2007. Under these latter scenarios, conditions that favour male hatchlings are limited to June–September. All other months favour female primary sex ratios (Figure 8B,C). Under a 2070 high emissions scenario, sand temperatures increased on average by 2.78°C, and the only months where temperatures fall below the TRT and thereby favour males are June–August (Figure 8D).

Table 4 Sand properties required by NicheMapR’s microclimate model

Figure 8

Sand temperature at 50 cm depth estimated for Cape Domett under four climate change scenarios: (A and B) 2030 low and high emissions, (C and D) 2070 low and high emissions (solid lines). In each case, sand temperatures in 2007 (dotted lines) are shown for comparison. Horizontal dashed lines show the bounds of the TRT, while grey or hatched boxes indicate the periods that would favour the production of male or mixed-sex hatchlings, respectively, under each scenario. The upper limit of the shaded and hatched boxes indicates a potential upper limit of 35°C for embryonic viability. The red dashed bar indicates the nesting season with the solid portion indicating peak nesting [21].

3Discussion

Cape Domett hosts the major winter nesting aggregation of N. depressus in Northern Australia [21]. Long-term data on nesting phenology does not exist for this rookery, but in the single year that this was assessed (April 2006–April 2007), peak nesting occurred in August and September [21]. Our sand temperature models show that from 1990 to 2013, the rookery experienced progressively cooler incubation temperatures in August and September (Figure 6), which are likely to have led to the production of predominantly male hatchings. This is contrary to a warming (and thereby feminising) trend that has been either documented for sea turtle beaches in recent years or else has been estimated from correlative models [9],[10],[24],[25]. Our anomalous result is explained by Cape Domett falling within a small region of north-western Australia where air temperatures have declined relative to historical benchmarks, due to large increases in rainfall and increase in atmospheric aerosols [26]. This trend is captured in the AWAP climate database that forced our model of historical sand temperatures (Figure 1). A change in peak nesting to spring months could compensate for the recent cooling of sand temperatures in winter. Once long-term data on nesting seasonality becomes available for this rookery, our mechanistic model could be used to assess whether phenological shifts would neutralise a male sex ratio bias (e.g. [8]).

Sand temperatures at 50 cm depth reconstructed over a 23-year period (1990–2013; Figure 6) showed similar seasonal patterns to the two recent years that we depict in greater detail (2011–2012; Figure 7). In most years, sand temperatures showed annual variation of approximately 10°C (Figure 6). As female N. depressus produce an average of 2.8 clutches per season and have a renesting interval of 12–16 days [16],[17], most females should produce offspring of both sexes in each year. Based on the new physiological data reported here, the TSP begins approximately 22 days after oviposition at the Tpiv of 29.4°C, suggesting that the peak nesting period identified in 2006–2007 [21] would have produced predominantly male or mixed sex hatchlings (Figure 6). As our field work occurred during September and October 2012, it was likely to have occurred towards the end of the peak nesting period. Sand temperatures are several degrees centigrade higher in these months than they are in August (Figure 7), and most embryos sampled from nests laid at this time were female (Table 3). Our sand temperatures and sex ratio predictions in 2012 showed good agreement with our empirical data (Table 3).

A trend evident in Figure 6 was an increase in the annual range of sand temperatures, driven primarily by the decreased winter temperatures. This pattern led to a greater proportion each nesting season where sand temperatures favour the production of male hatchlings. If these trends continue without a consequent shift in the timing of peak nesting, then primary sex ratios are likely to become more male-biased at this rookery, contrary to our original hypothesis. The increasing sand temperature range over the 23-year period also coincided with warmer summer temperatures in later years (e.g. 2005–2012; Figure 6D,E). These recent summer sand temperatures exceed the upper threshold of ~35°C for survival of marine turtle embryos [27], while temperatures predicted under future climates peak at 38.3°C in December (Figure 8D). Hence, year-round nesting would come under strong negative selection due to heat-induced mortality of embryos. Cooler winter and warmer summer temperatures could result in spring and autumn months becoming the most suitable for producing viable hatchlings of both sexes. Potentially, under future climates, this population could shift from predominantly winter nesting to bimodal nesting with peaks in both spring and autumn.

3.1 Validation of sand temperature predictions

The validity of the predicted sand temperatures is strengthened by comparison with actual measurements of nest temperatures at equivalent depths (Additional file 1: Figure S1 and Figure 5). NicheMapR’s microclimate model produced accurate estimates of nest temperatures (without any data-based calibration) across a wide range of depths, except towards the end of incubation when metabolic heating is likely to have increased nest temperatures above that of the surrounding sand. For example, in green (Chelonia mydas) and loggerhead turtle (Caretta caretta) nests, a temperature differential between 0.07°C and 2.86°C has been attributed to metabolic heating [28],[29]. A divergence in measured and predicted temperatures towards the end of incubation would not affect our sex ratio estimate, as metabolic heating should be most pronounced after the TSP. For example, in a direct comparison of sand and nest temperatures at identical depths in a study of N. depressus in the Pilbara, Box (2010) concluded that only two of seven nests showed any detectable sign of metabolic heating [20]. Nevertheless, it would be advisable to include a metabolic heating process (e.g. [30]) into a mechanistic model if traits such as hatchling mortality or incubation duration were of interest.

We did not examine the influence of fine-scale variation in nest orientation, shade or sand properties on nest temperatures in this study, as the Cape Domett rookery is fairly homogenous in these respects. The location of the nest within the beach profile would likely produce the most sand temperature variation due to variable sand moisture and potential shading from vegetation. These parameters can all be modified and mapped spatially using the NicheMapR system (e.g. [8]). Rapid changes in beach topography can make fine-scale modelling of nest shade and exposure challenging, but high-resolution, high frequency satellite and LIDAR imaging could ultimately allow the dynamic rendering of beach temperatures and sex ratios.

3.2 Delineation of the thermosensitive period (TSP)

The precise delineation of the TSP is important for determining the sequence of incubation temperatures that contribute to sex determination [31]. The TSP is the least researched of all TSD parameters and most reptile studies delineate the TSP no more precisely than the middle half to the middle third of development [2],[32]. We showed that the TSP of N. depressus spans a period of about 23% of embryonic development (43%–66%), supporting our hypothesis that the TSP would be less than 33%–50% of the total development period. Hewavisenthi and Parmenter [19] provide the only other estimate of the TSP in N. depressus, concluding that it occurs between days 32 and 40 of incubation at 26°C. Based on incubation temperatures reported in the Hewavisenthi and Parmenter study, and the development rate model developed here (Figure 2), the proportion of development that would have occurred on days 32 and 40 would have been ~44% and ~54%, respectively. Both studies therefore infer that the TSP starts at about the same developmental stage. However, this study infers that the TSP is longer; encompassing approximately 23% rather than 10% of embryonic development shown in the Hewavisenthi and Parmenter study.

All hatchlings (except one) produced from female–male–female (FMF) treatments were female, indicating that either a masculinising switch does not fall within the TSP we identified, or else that our switch periods were too short or not sufficiently cool to trigger the development of males. The switch periods were designed to encompass the same proportions of development for switches in both directions. Instead, development of N. depressus embryos from Cape Domett occurred at slightly faster rates than for embryos collected from the Pilbara [20]. This caused a greater proportion of development to be completed by the start of the FMF switches than for the corresponding male–female–male (MFM) switches.

Switches to 32.5°C over as little as only 9% of development produced a feminising effect (Figure 3; treatment 5). Experiments on the pond turtle (Emys orbicularis) show that bipotential gonads develop into testes in the absence of a feminising temperature trigger during the TSP [33]. In contrast, when this trigger occurs, the enzyme aromatase acts on gonads to stimulate the production of oestrogen and development of ovaries [1],[33],[34]. In N. depressus, we hypothesise that the TSP is similar for triggering testicular or ovarian development, but that the proportion of development at feminising temperatures may be 9% (or less) to trigger ovarian development. Hence, if gonads are to develop into testes, the majority of the TSP needs to occur at a masculinising temperature. This is a possible reason for the lack of masculinisation seen in the FMF back-switch experiments (treatments 7–11), as experimental switch periods may have been too short for the majority of the TSP to occur at the masculinising temperature. Further experimentation is required to test this hypothesis, and its resolution will have important implications for estimating hatchling sex ratios in natural nests that fluctuate between male- and female-producing temperatures.

3.3 TSD parameters and embryonic development rate

Identifying the Tpiv, TRT and embryonic development rate were not major foci of this study, but the ability of our physiological model (Figure 1) to provide good estimates of hatchling sex ratios at Cape Domett (Table 3) suggests that these parameters and the boundaries of the TSP were reasonably accurate. The Tpiv estimated in this study of 29.4°C for the Cape Domett population is similar to the 29.3°C estimated by Limpus (1995, cited in [35]) and the 29.5°C estimated for Eastern Australian populations [18]. It is, however, lower than the 30.1°C estimated for a population of N. depressus from the Pilbara region of Western Australia [20]. Differences in the Tpiv estimates for the two Western Australian populations are concordant with the genetic differentiation reported between summer nesting populations from the Pilbara and winter nesting populations in Northern Australia [36].

We used sex ratio data from field nests to improve our estimates of the Tpiv and TRT, but field data is not an ideal substitute for carefully controlled laboratory experiments. In addition, we have not assessed sex ratios close to the estimated Tpiv (29.4°C). Consequently, constant temperature incubation experiments at temperatures at and near the Tpiv are needed to further refine the TSD parameters for the Cape Domett population. Incubation experiments at higher temperatures would also be highly relevant when considering impacts of global warming, as development rate increases with temperature only up to a certain point, after which it rapidly decreases due to enzyme deactivation [2]. High nest temperatures are a documented cause of mortality in marine turtle nests [12],[37] and well-resolved thermal response curves will be essential for predicting embryonic survival under climate change.

3.4 Implications for population viability and management

Marine turtles have adapted to past climate change and may also adapt to accelerated rates of climate change associated with anthropogenic global warming [38]. Adaptations to counter increasing temperatures include adjusting the timing of nesting, altering nest characters such as shade or depth and altering nesting latitude [13],[39]. As peak nesting at Cape Domett coincides with the coolest sand temperatures, it constrains opportunities for females to temporally adjust peak nesting times to counteract increases in incubation temperatures. Another behavioural characteristic—nest site choice with respect to shade cover—has shown some degree of plasticity in other reptiles, where females select sites of differing shade availability in different climates [40],[41]. The effectiveness of this response may be limited at Cape Domett, as shaded nest sites are scarce. Alternatively, females could dig deeper nests, where sand temperatures are usually cooler and more stable [14], but nest depth is limited primarily by limb length in turtles, meaning that nest depths are relatively fixed relative to those of other reptiles [40]. Similarly, cooler nest locations low on the beach risk flooding and suffocation of embryos.

The heritable variation in TSD traits presents an alternative avenue for adaptation to increasing temperatures. Microevolution of the Tpiv and/or TRT in response to climate change is possible [42],[43], and, unlike behavioural adaptations which alter temperatures experienced by embryos, changes in the Tpiv or TRT would alter the sex ratio produced at a given temperature. In general, data on the heritability of physiological traits in reptiles is scarce and as such provides an important avenue for future research [43]. Mechanistic models such as those developed here readily allow the effectiveness of microevolutionary change to be examined.

The overwhelming consensus that anthropogenic climate change will increase global temperatures means that female biases in primary sex ratios could become more prevalent in marine turtles [44],[45]. Consequently, population viability may decline and in extreme cases, localised extinctions could occur [10]. Projected global warming in Queensland suggests that within 50 years, ratios of one male to four females are possible for many marine turtle rookeries in the Great Barrier Reef [38], and such biases may not be sustainable when combined with other anthropogenic stressors such as habitat modification and loss of rookeries [39],[42]. However, polyandrous mating systems may allow marine turtles to persist under female-biased operational sex ratios (e.g. [46]). Understanding responses to operational sex ratio biases and accurately predicting primary sex ratios will both be important when assessing population viability under a warming climate.

4Conclusions

Our reconstruction of cooling in sand temperatures at Cape Domett in recent decades (Figure 6) contrasts with our projections of warmer sand temperatures under future climates (Figure 8). While our projections of future sand temperatures could be further refined by applying downscaled global climate models that anticipate changes in rainfall and cloud cover for this region [26], it is clear that warmer sand temperatures would lead to more female hatchlings being produced at Cape Domett relative to the primary sex ratios produced historically. Furthermore, the periods of year when sand temperatures remain below an upper limit for embryonic viability are reduced appreciably. If the trend of an increasing range in annual sand temperature continues (as seen in the 23-year reconstruction in Figure 6), then the winter sand temperatures that currently coincide with peak nesting could decrease further, potentially favouring male hatchlings. The fact that Cape Domett falls within a small region of Australia where annual minimum air temperatures have decreased [26] means that monitoring and conservation of this rookery should be a key priority, as rookeries that produce predominantly males can balance a feminising trend elsewhere in the nesting distribution.

5Methods

5.1 Ethics Statement

All procedures were reviewed and approved by the University of Western Australia animal ethics committee (UWA animal ethics permit RA/3/100/1145) and the Department of Parks and Wildlife (scientific licence SF008844).

5.2 Study site and egg collection

The N. depressus rookery at Cape Domett occurs on a gently sloping, 1.9 km-long, north-west facing beach in the East Kimberley region of Western Australia, Australia (14.798° S, 128.415° E). Freshly oviposited eggs were collected for transport to a laboratory in Perth, Western Australia on September 20 2012. Approximately, whole clutches (52 eggs) were collected from each of six nesting females (total 296 eggs) and stored in damp vermiculite within a portable refrigerator (model: Engel MT45FP) set to 8°C to suspend embryonic development thereby facilitating egg viability during transport (Harry & Limpus 1989). Eggs were transported by helicopter to Kununurra where they awaited further transport via a commercial airline to Perth. Inadvertently, eggs were held over in Kununurra due to a lack of cargo space, and consequently eggs reached the University of Western Australia about 100 h after oviposition. This time period fell outside the 72-h window of viability defined by Harry and Limpus [47].

In situ incubation temperatures were recorded in 11 nests on the nights of September 16 and 17 in 2012. Four Thermochron® iButtons (Maxim Integrated Products; DS1921H; accuracy ±1°C, precision 0.125°C) were placed at intervals into each nest chamber during oviposition, with the goal being to record temperatures at different depths within the nest chamber. The nests were situated at various locations within the ‘pink sand’ sections of the beach [21] and were marked with a GPS device (Garmin eTrex Vista HCx) and photographed for ease of relocation. A weather station (WeatherHawk, Signature Series 232) was erected midway along the beach on a small dune above the high tide mark. Weather data (air temperature, humidity, wind speed, solar radiation) were recorded at 30-min intervals.

The rookery was revisited 40 days later (25–26 October 2012) and all nests were relocated and excavated. Between 12 and 16 eggs that were closest to each of the four iButtons were removed from each nest (total of 138 from all nests). The depths of iButtons below the sand surface were measured to the nearest 10 mm before removal. All remaining eggs were reburied. The embryos inside the sampled eggs were euthanised by chilling for approximately 24 h. Chilling was achieved on the beach using ice, and eggs were relocated to a portable refrigerator set to 2–4°C within 2–3 h of collection. Chilling was maintained during transport by boat and road back to Kununurra. There, dead embryos were removed from their shells and fixed in 10% buffered neutral formalin for road transport to Perth. Ten of the eleven nests (nest IDs 1, 2, 3, 4, 5, 6, 10, 12, 13 and 15) contained viable embryos; eggs from the remaining nest showed no signs of development.

5.3 Incubation experiments

The six live clutches returned to Perth were either used in ‘back-switch’ experiments designed to delineate the TSP, or else were incubated at a constant temperature close to the estimated pivotal temperature. Ten eggs (2 each from 5 clutches) that were not used in back-switch experiments were incubated at a constant temperature of 30.1°C (Model i180), which was the estimated pivotal temperature for a population of N. depressus from the Pilbara region of Western Australia [20].

For the back-switch experiments, one egg from each clutch was buried in damp sand within one of 48 partially sealed containers. Containers were placed in one of four Steridium incubators (models: i180 and i500) set to either a male-producing temperature (28°C; 24 containers) or a female-producing temperature (32.5°C; 24 containers). During development, each container of 5–6 eggs was switched to the opposite temperature treatment for either 10% or 15% of embryonic development (96–240 h, depending on the time and temperature), after which time they were switched back to their original temperature (Table 1). Each back-switch treatment was replicated within a pair of incubators (one pair of 180 L incubators and one pair of 500 L incubators), and iButtons were placed next to an egg when each container underwent a switch to record the temperatures experienced by the embryos.

We used a development rate function developed for N. depressus by Box (2010) to calculate the day that embryos should have completed 35% of development to hatching when incubated at either 28°C or 32.5°C. Switches then began on this day (or when embryos were further developed—refer Table 1) and concluded when a further 10% or 15% of development to hatching should have been completed at the switch temperature (28°C or 32.5°C). Procedural controls were also included for the earliest switches, in which case embryos incubated at 28°C were switched to a different incubator at 28°C for the appropriate period, and then back again to the original incubator. A procedural control was also applied to eggs incubated at 32.5°C.

Eggs in both the back-switch and Tpiv experiments were weighed at the start of incubation and thereafter every seven days to assess embryonic viability. Unviable eggs (those that lost weight and became discoloured) were removed from incubation boxes. On each occasion, containers were repositioned randomly in the incubators to minimise any impact of subtle temperature gradients within incubators. Containers were checked daily after 35 days of incubation (one week before the earliest expected hatch date) and incubation time was recorded as the time from the start of incubation to when pipping commenced. The four days between oviposition and the start of incubation were not included in the incubation time, as eggs were cooled during this period and embryonic development was suspended [47]. Once pipping commenced, hatchlings were weighed and euthanised by an intracoelomic injection of 0.4 mL sodium pentabarbitone at a concentration of 160 mg/kg. Hatchlings were labelled and stored in 10% buffered neutral formalin for later dissection.

5.4 Identification of hatchling sex

The gonads of marine turtle hatchlings are small (<500 μm) and attached to the kidney [48]. Consequently, the left kidney of each preserved individual was removed. Kidneys were then prepared as paraffin-embedded sections and stained with haematoxylin and eosin for light microscopy [31]. The sex of each specimen was identified based on the differentiation of gonadal medulla and cortex, where males have seminiferous tubules in the medulla and a regressed cortex, and females have a disorganised medulla and a thick, well-developed cortex [49],[50]. Each specimen was classified as male, female or unknown on three separate occasions, without reference to previous assessments, and a repeatability analysis was conducted to determine the reliability of the classification. Any specimens where a gonad was not visible or was unable to be classified were resectioned and reexamined until each individual could be classified as male or female.

5.5 Calculation of a non-linear function for embryonic development rate

Development rate, expressed as a function of temperature, was calculated from data on incubation duration using the program DEVARA [http://delta-intkey.com/devar/] [31],[51]. Hourly temperature records (from iButtons) and the average incubation time (in days) of eggs in each container provided the inputs for DEVARA. Data from eggs incubated only at constant temperatures (Tpiv experiment and two back-switch controls) were used, and temperatures were sub-sampled to six values per day to reduce the amount of information read by the program. DEVARA fits a non-linear model expressing development rate (r _a ) as a percentage per day, as a function of temperature:

$$r_a=b_{1}10^{‐v^2\left(1‐b_5+b_5{v}^2\right)}$$

where

$$u=\left(\left(T‐b_3\right)÷\left(b_3‐b_2\right)\right)‐c_1$$

$$v=\left(u+e^{b_{4}u}\right)/c_2$$

$$c_1=1÷\left(1+0.28b_4+0.72ln\left(1+b_4\right)\right)$$

$$c_2=1+\left(b_4÷\left(1+1.5b_4+0.39{b_4}^2\right)\right)$$

The parameters fitted by DEVARA describe the maximum development rate (b₁), its corresponding temperature (b₃) and the temperature at which development rate approaches zero (b₂). Parameters b₄ and b₅ control the asymmetry and steepness of the curve and were fixed at 6 and 0.4, respectively, as is recommended when development rates at extreme incubation temperatures are unknown [51].

5.6 Delineation of the thermosensitive period (TSP)

Using the parameters estimated for the non-linear development rate function defined above, hourly temperature records from the back-switch experiments were converted into developmental increments. These increments were integrated to determine the cumulative proportion of development that was completed on each day. The length of each switch window was thereafter described as a developmental proportion (e.g. 12% of development). If any back-switch regime produced a sex ratio different to that expected from the dominant incubation temperature, then the TSP was assumed to fall within the portion of development where the temperature was switched.

5.7 Determination of the pivotal temperature (Tpiv) and transitional range of temperatures (TRT)

The relationship between incubation temperature and sex ratio was calculated using TSD software version 4.0.3 [http://max2.ese.u-psud.fr/epc/conservation/TSD/index.html] developed by [52]. This program is widely used to compare the fit of up to five threshold models using maximum likelihood (Richards/a-logistic, Weibull, Weibull*, Hill and Hill*) where * indicates the a-logistic or ‘asymmetrical’ version of the model. Data from constant temperature incubation were intended to be used for this analysis, but the program requires that three or more temperatures produce mixed sexes in order for the fit of different models to be compared [52]. Of the five constant incubation temperature regimes applied in this study (Tpiv experiment, male back-switch controls and female back-switch controls), only one produced mixed sexes. Hence, sex ratio and temperature data from field nests were added to the laboratory data to produce a more rigorous estimate of a sex ratio function.

As the TSD software is designed for constant temperature data, temperature records from field nests were converted to constant temperature equivalents (CTEs; [53]). Effectively, the CTE is the temperature above which half of development occurs (i.e. a developmental median) and thereby the relationship between development rate and temperature must first be established. Furthermore, when incubation temperatures are variable, the CTE during the TSP (and not the entire incubation period) is most relevant for fitting the TSD function [31]. Hence, we used our development rate function (Figure 2) to convert the nest temperatures that fell within the TSP into hourly developmental increments. These increments were ranked and integrated to determine the temperature above which half of development occurs (the CTE).

The CTE parameter and sex ratios measured from all ten viable field nests were used in combination with laboratory data, and the Hill equation produced the lowest AIC value. The Hill equation has the form:

$$\mathrm{sr}\left(t\right)=\frac{1}{1+e^{\frac{1}{s}\left(ln\left(P+K\right)‐ln\left(t+K\right)\right)}}$$

where sr(t) is the sex ratio at a given temperature t, P is the pivotal temperature, S describes the steepness of the transition from male to female producing temperatures and K is a parameter that describes the asymmetrical shape of the function. The fitted values of P, S and K and the corresponding TRT were noted and used as reference points for a further ten models, which were fitted from random subsets of five field nests (Table 2). By using only five of the ten field data points available, we were able to meet our goal of retaining an independent data set to test our mechanistic framework (Figure 1). For each of the ten subsets, there was more than one equation that produced an equally good fit (∆AIC <2); however, the Hill equation consistently produced the lowest AIC value across all subsets. Hence, we used the Hill equation to make a standardised comparison of the TSD equation generated by each of the ten datasets, and used the mean residuals to select a consensus model that could be used to estimate the Tpiv and TRT of the Cape Domett population.

5.8 Sand temperature reconstruction and projection

NicheMapR is an R version of the mechanistic modelling program NicheMapper^TM; [54] and was used to predict hourly sand temperatures by simultaneously solving heat and mass balance equations based on physical properties of beach sand and on regional climate data (e.g. [11],[22]). The parameters of the microclimate model were not ‘tuned’ to the observed sand temperatures. Rather, the sand properties for the Cape Domett rookery were estimated from relevant literature (e.g. [55],[56]) or were measured empirically (Table 5). Percent sand moisture was measured by determining the wet and dry weights of sand samples collected from a range of depths from each field nest during excavation in October. Solar reflectance of sand sampled from two nests in the wavelength range 300–2,100 nm was measured using two spectrometers (Ocean Optics USB2000 for the UV-visible range and NIRQuest for the NIR range) and two light sources (Ocean Optics PX-2 pulsed xenon light for the UV-visible range and HL-2000 tungsten halogen light for the visible-NIR range), all connected with a quadrifurcated fibre optic. The probe on the end of the fibre optic was held within an ocean optics RPH-1 probe holder at a constant angle (45°) and distance from the surface, and each measurement was expressed relative to a Spectralon 99% white reflectance standard (Labsphere, Inc., North Sutton, NH, USA), and weighted by solar irradiance. There was no significant difference in the solar reflectance of the two samples (unpaired t-test; p = 0.959).

Daily maximum and minimum temperatures, relative humidity and wind speeds were collated from data generated by the weather station we installed at the rookery between September 18 and October 24, 2012. Rainfall data during this period were obtained from the nearest Bureau of Meteorology weather station at Wyndham (http://www.bom.gov.au). These weather values were used as inputs into NicheMapR, and sand temperatures were predicted for user-specified depths (range 36–80 cm) using the parameters in Table 2 and assuming 0% shade. Predicted temperatures at a particular depth were then compared to actual nest temperatures measured at the same depth, using the r² statistic.

Historical (1990–2013) climate data for Cape Domett (daily maximum and minimum temperatures, relative humidity, rainfall and solar radiation) were obtained from the Australian Water Availability Project (AWAP). Gridded long-term average wind speed data were obtained from Australian National University Climate software package (ANUCLIM; [57]). Both the AWAP and ANUCLIM databases are derived from interpolated data records from weather stations across Australia. These historical climate data were interpolated at a point approximately 25 km south of the Cape Domett rookery (15.003° S, 128.383° E).

The AWAP climate data, combined with the shade and sand parameters defined earlier, were used as inputs into the microclimate model within NicheMapR to estimate sand temperatures at 50 cm depth for the 23 years from 1990 to 2013. To investigate the influence of climate change on sand temperatures, the 2007 air temperatures from the AWAP database were adjusted in accordance with the Commonwealth Scientific and Industrial Research Organisation’s (CSIRO) projections of future air temperatures for Australia [58]. Under a low emissions scenario, AWAP air temperatures were increased by 0.6°C or 1.8°C for the years 2030 and 2070, respectively. Under a high emissions scenario, air temperatures were increased by 1.5°C or 3.4°C for the same years.

5.9 Model validation

The DEVARA function was used to convert the hourly sand temperatures predicted by NicheMapR (or the actual temperatures measured in nests) into developmental increments for the five nests not used to fit the TSD function. Development was assumed to have started on the date of oviposition; hence, the first developmental increment we calculated was for the hour immediately following oviposition. All the hourly development increments calculated for a particular temperature record were integrated to determine the proportion of development completed on each day and thereby identify the dates that formed the boundaries of the TSP. We calculated a CTE for all temperatures that fell within the TSP (as described earlier) and CTEs were inputted into the TSD function (Figure 4B) to estimate the sex ratio of the embryos sampled at the various depths in each nest. Sex ratios were classified as female (<5% males), mixed (5–95% males) or male (>95% males). Predicted sex ratios were then compared to the actual sex ratios of the sampled embryos, which allowed us to assess how effectively our physiological model estimated sex ratios from hourly temperature records (Figure 1).