Predicting the effects of climate change on incubation in reptiles: methodological advances and new directions

Predicting how biodiversity will respond to an increasingly warm and variable climate is one of the most important questions of contemporary ecology. Reptiles with temperature-dependent sex determination (TSD) are a model system for both empirical and analytical approaches to understanding impacts of climate change on vertebrate ectotherms. In this diverse group of oviparous species, which includes most turtles, some lizards, and all crocodilians and tuatara (Valenzuela and Lance, 2004), the thermal conditions that an embryo experiences direct the process of gonadal differentiation, i.e. whether it ultimately develops ovaries, testes, or – more rarely – an intermediate form. This discrete, measurable relationship between temperature and a physiological response that (1) has consequences not just for individual fitness (Warner and Shine, 2008) but also for a key population demographic parameter (Schwanz et al., 2010) and (2) can be exploited in a controlled laboratory setting (While et al., 2018) makes TSD a natural target for studying the potential ecological impacts of climate change and one that may have dire consequences for reptile populations (Janzen, 1994a), many of which are already subject to multiple anthropogenic stressors (Böhm et al., 2013; Stanford et al., 2020).

In this paper, we provide a brief overview of the TSD system and discuss some of the challenges in predicting the impacts of climate change on reptiles with TSD, with an explicit focus on how these issues are amplified as the biological (e.g. from embryo to clutch) and spatial (e.g. from nest to rookery) scales of our inquiry increase. Supported by a literature review, we also examine how research to date has modeled the relationship between temperature and gonadal differentiation and explore how recent advances in biophysical and microclimate modeling are placed to tackle these challenges. Finally, we discuss some steps that are necessary for research on TSD to maximize benefit from the latest modeling approaches, while calling for studies to remain grounded in observational data. Throughout this paper, we use the term ‘microclimate’ to refer both to (1) the microclimate layer of Earth's atmosphere, or the climate near the ground, within the context of modeling those conditions across a geographical extent and for a discrete period of time (Geiger et al., 2003) and to (2) the small-scale spaces that organisms inhabit and within which they can affect, and are affected by, local microclimate conditions. Within the context of reptile embryonic development, the relevant microclimate is that of a sub-surface nest.

In reptiles with TSD, gonadal differentiation occurs at the embryo scale but is modeled at the clutch or cohort scale as a curvilinear relationship between some measure of incubation temperature and hatchling sex ratio. The resultant reaction norm is then used to describe a population or, if sufficient population-specific data are unavailable, a species. While the prevailing terminology has evolved (Mrosovsky and Pieau, 1991), the definitions of the relevant parameters have changed little over the past three decades. A population's particular pattern of TSD is classified based on whether production of hatchlings with differentiated ovaries increases towards an asymptote at high temperatures (type Ia), low temperatures (type Ib), or both (type II) (Mrosovsky and Pieau, 1991). The type of TSD actually observed in a population is further described by two discrete parameters: the pivotal temperature (T_piv) and transitional range of temperatures (TRT). T_piv describes the inflection point(s) of the TSD reaction norm, where the ratio of female:male hatchlings would theoretically be 1:1 under constant incubation conditions (Mrosovsky and Pieau, 1991). The TRT is the range(s) of temperatures between the minimum and maximum female:male ratios (Mrosovsky and Pieau, 1991) (Fig. 1) and can be asymmetrical around T_piv. Because the rate at which the response asymptotes can vary on either side of T_piv, predicting hatchling sex ratios requires an accurate estimate of the TRT.

In addition to the parameters describing the shape of the TSD reaction norm, models benefit from knowledge of the thermosensitive period (TSP), the particular window(s) of development during which temperature affects hatchling sex ratio (note that the TSP is defined at the clutch scale or coarser) (Mrosovsky and Pieau, 1991). While the TSP itself is not included explicitly as a model parameter, only the thermal regime experienced during this period is typically used to define the ‘incubation temperature’ (Fuentes et al., 2017). The conditions embryos experience outside of the TSP affect the rate of development and, thus, the total amount of development completed at a given time point but do not have an impact on the outcome of gonadal differentiation. Probably because of the relative difficulty inherent in pinpointing the TSP via temperature-shift experiments (reviewed in Janzen and Paukstis, 1991), it has received little explicit attention in the TSD literature, relative to the TRT and T_piv. The TSP is often described as the ‘middle third of development’, an estimate based on results from shift experiments using two or more constant temperatures (but see Nelson et al., 2010). However, the range of developmental stages that correspond to the TSP can vary considerably within and among species, depending on the constant temperatures used in those experiments (Janzen and Paukstis, 1991; Mrosovsky and Pieau, 1991). In addition, ‘heatwave’ incubation experiments, in which embryos are exposed to short periods of relatively high or low temperatures, have shown that actual gonadal phenotypes can differ from those predicted when models assume equal effects of incubation temperature throughout the middle third of development (Woolgar et al., 2013).

Even if the TSP could be adequately described as the middle third of development, we could not precisely identify the TSP using time (e.g. ‘days from oviposition’ or ‘duration’) as a proxy, unless incubation conditions were highly stable (Fuentes et al., 2017). In addition, climate change will probably shift the timing of the TSP via increasing incubation temperatures and, consequently, development rates or via a shift in nesting phenology, or both (Telemeco et al., 2013). Ideally, we would quantify the TSP as a range of percentage development, during which the outcome of gonadal differentiation remains indeterminate. Using the non-linear, direct relationship between temperature and embryonic development rate (Georges et al., 1994; Girondot and Kaska, 2014; Girondot et al., 2018) would reduce a substantial source of imprecision in predictions (Carter et al., 2019), potentially improving the accuracy of TSD models under fluctuating incubation conditions, generally, and as both the mean and variance of temperatures shift due to climate change. In practice, most estimates of the TSP remain too coarse to inform rigorous climate change-focused studies.

The various non-linear models used to predict development rates and hatchling sex ratios in reptiles with TSD have been reviewed most recently and comprehensively by Fuentes et al. (2017) and Abreu-Grobois et al., (2020) and implemented in the R package ‘embryogrowth’ (https://rdrr.io/cran/embryogrowth/). However, models of TSD describe processes that take place at spatial scales from local (populations) to fine (rookeries) to ultra-fine (nests) and at fine temporal scales, typically hourly or daily. To predict the impacts of a warming climate on embryonic development and hatchling sex ratios in reptiles with TSD, we must pay particular attention to how changes in scale affect the accuracy and precision of model predictions. The challenges of modeling TSD at the broad spatial and temporal scales typically sought by climate change ecology have rarely been explicitly addressed.

Two main problems of scale are inherent to research on TSD. The first, and the one that has received greater attention in the literature, is that we model an embryo's thermal developmental response at a coarser biological scale (i.e. clutch, cohort or population). Natural variation in abiotic conditions can lead to phenotypic outcomes that are difficult to predict, even if the number of possible outcomes is small. Conceptually, the relationship between an environmental variable and a phenotypic outcome at the organismal scale can take three forms: one-to-one, one-to-many, or many-to-one. We can safely assert that TSD in reptiles takes the latter form. The number of possible outcomes is very low (i.e. a single embryo may develop ovaries, testes or an intermediate form), but many patterns of natural incubation conditions appear to result in the same discrete outcome.

Both empirical and modeling approaches to TSD have recognized this many-to-one relationship implicitly by defining and estimating proxy variables (i.e. T_piv, TRT and TSP) that attempt, with varying degrees of success, to translate the phenotypic outcomes of natural incubation conditions to the laboratory environment, and vice versa. To predict the outcome of gonadal differentiation in reptiles with TSD, current models require three, interdependent pieces of information: (1) the period(s) of time – during incubation – during which gonadal differentiation occurs, (2) the incubation temperature during that period of time, and (3) the relationship between that incubation temperature and the sex ratio of the clutch or cohort. Each of these variables, derived from past constant-incubation experiments, introduces error into the model, which is then aggregated and potentially amplified in the resultant prediction. That modeled sex ratios are sometimes wildly inaccurate compared with those observed in natural nests is not particularly surprising.

Regardless of what proxies are used to describe the relationship between temperature and gonadal differentiation, we model TSD at the cohort scale, in terms of hatchling sex ratio: percentage female or percentage male, where hatchlings with intermediate gonads are typically excluded (Mrosovsky and Pieau, 1991). Sex ratios are often binned into ‘all female’, ‘all male’ or ‘mixed’ nests. However, coarsening the data resolution so that a nest containing 20% female hatchlings is binned into the ‘mixed’ category with a nest containing 80% female hatchlings would seem to serve little useful purpose, especially given that reptiles with TSD are a research focus within the context of climate change due specifically to the risk of demographic collapse and local extinction in populations with highly skewed sex ratios (Grayson et al., 2014; Janzen, 1994a; Mitchell et al., 2010; Schwanz and Janzen, 2008; Valenzuela et al., 2019). The potential consequences for a particular population in which 20% of hatchlings have ovaries versus 80% of hatchlings with ovaries would probably be very different (Grayson et al., 2014; Mitchell et al., 2010) and necessitate different types and degrees of management response.

Despite substantial research efforts since TSD was first described in a reptile by Charnier (1966), understanding how gonadal differentiation scales up from organisms to populations has remained challenging. Even though the number of possible outcomes of gonadal differentiation is few and finite, the complexity of thermal developmental plasticity does not translate exactly, or even well, between biological scales. This incongruity was perhaps inevitable. In the same paper that settled the definitions of the proxy variables now used widely in TSD research, Mrosovsky and Pieau (1991) explicitly cautioned against using clutch-scale estimates derived from constant-incubation laboratory experiments – often of only a few clutches – to represent wild populations. We should perhaps not be surprised that the predicted sex ratio differs depending on what proxy is used to estimate it (Fuentes et al., 2017). In addition, where data are available to validate predicted sex ratios, model performance varies considerably, depending on how the sex ratio itself is quantified in the validation dataset. For example, in addition to hiding important variation, binning the sex ratio of nests into discrete ‘all female’, ‘all male’ and ‘mixed’ categories (Telemeco et al., 2013) inflates model precision, compared with using stricter numerical criteria (Carter et al., 2019).

As climate change necessitates that research aims shift increasingly from acquiring an empirical understanding of TSD to predicting the impacts of increasing temperatures, greater caution is warranted when making decisions about how to model TSD, especially in populations for which few validation data are available. Even small differences between laboratory-derived proxies and the ‘real’ thermal responses of wild populations could lead to realized sex ratios that are vastly different from expected, based on modeling efforts. However, a perfect biophysical model of TSD, regardless of the biological scale it reflects, is useless without accurate and precise temperature data to drive it. Thanks to the ubiquitous iButton (Maxim Integrated) and other low-power, low-cost temperature data loggers, we can directly measure the abiotic environment inside nests as a proxy for the conditions experienced by individual embryos with relative ease. However, we are logistically constrained in the number of nests we can sample over both space and time. This limitation is a component of the second and, we argue, more critical issue that needs to be addressed in attempts to model the impacts of climate change on reptile sex ratios.

The second problem of scale, and the one that has received less attention in the TSD literature, is that understanding how climate change will affect gonadal differentiation requires not only scaling physiological responses up from organisms to populations but also scaling climatic processes down. Although incubation in natural nests is obviously not constant, our empirical understanding of TSD is based largely on experiments conducted under constant incubation conditions (Janzen and Paukstis, 1991; While et al., 2018), probably because the climate chamber technology needed to simulate more realistic conditions has not been accessible to most researchers. Thus, a substantial amount of effort has been dedicated to identifying or defining temperature variables that represent the stochastically fluctuating diel cycles of natural incubation and that reproduce the physiological effects of temperature fluctuations (Bowden and Paitz, 2018) with a single value (reviewed in Girondot et al., 2010).

Increasingly, both empirical and modeling approaches to TSD research are investigating methods that account for fluctuating diel temperature cycles, because summary variables (e.g. mean temperature) are poor correlates of hatchling sex ratio. In the predictive context, the most widely used temperature proxy is the constant-temperature equivalent (CTE), which was derived as a more accurate alternative to the mean (Georges, 1989) and represents incubation temperature as a single value based on development rate. The CTE is the temperature at the developmental median, i.e. the point at which 50% of development has occurred. However, the median can be referenced to the entire incubation period, to the middle third of incubation (time-based), or to the middle third of development (physiology-based) and, depending on the chosen definition, makes various assumptions about the effect of temperature on development rate and whether the growth rate itself affects gonadal differentiation. The predicted sex ratio can vary by more than 10%, depending simply on how the CTE is defined and weighted (Fuentes et al., 2017).

Another recently developed solution is the weighted cumulative temperature unit (CTU), which – similarly to a degree-time model – is the sum of hourly deviations from T_piv multiplied by the corresponding development rate and expressed as a daily value (Valenzuela et al., 2019). Both the CTE and CTU define development as a function of hourly or daily incubation temperature. Thus, both require temperature data on a fine temporal scale and, as predictions are being made for a single clutch, an ultra-fine spatial scale. This pre-requisite is where most analytical studies rely on thermal profiles recorded in actual nests or otherwise from within known nesting areas (e.g. at random points). Availability of these data has been identified as a key limitation on efforts to model the effects of fluctuating incubation conditions on embryonic development (Bowden and Paitz, 2018). However, we argue here that it is our reliance on these data as model inputs, rather than their availability per se, that limits our ability to predict hatchling sex ratios in the face of climate change.

Earlier research that sought to predict the effects of climate change on reptiles with TSD tended to use a proxy-based approach, where a single measure of local environment was associated with measured sex ratios. For example, Janzen (1994a) and Refsnider and Janzen (2016) used mean July air temperature from a proximate weather station as a stand-in for incubation conditions during the TSP to estimate hatchling sex ratios of painted turtles (Chrysemys picta), given an increase in that mean of 4°C. Schwanz and Janzen (2008) used cumulative heating degree-days of the prior winter to examine how climate-linked shifts in nesting phenology might affect incubation temperatures and, thus, sex ratios in painted turtles. The relationship between vegetation cover – as a proxy for incident solar radiation, which affects substrate temperature – and sex ratio has also been investigated (Doody et al., 2006; Janzen, 1994b; Schwanz and Janzen, 2008). These and other correlative approaches are limited in their ability to extrapolate to novel climate scenarios and, more importantly, do not capture the temporal heterogeneity that developing embryos actually experience – or in the case of whole nests, spatiotemporal heterogeneity – regardless of climate (Potter et al., 2013; Sears et al., 2011). However, studies that record more realistic nest temperatures (i.e. hourly substrate temperatures) do not necessarily relate them directly to embryonic development. Bevan et al. (2019) recorded hourly sand temperatures at 130–140 points across nine nesting beaches used by Kemp's ridley sea turtles (Lepidochelys kempii) annually for 4 years, a substantial and laudable effort. However, they used a much-coarsened comparison of contemporary mean temperatures across the entire nesting season to speculate about the impacts of climate change on incubation.

One strategy for using available nest thermal profiles to simulate climate change is to transform contemporary hourly data to reflect predicted increases in mean and variance (Stouffer and Wetherald, 2007). For example, Valenzuela et al. (2019) applied an affine transformation to hourly thermal profiles recorded in nests, increasing their variance by ±2, ±4 or ±6°C to simulate climate warming. This approach is a substantial improvement over extrapolating coarse measures of local environment to coarse developmental outcomes. However, it ignores the likelihood that, as the climate warms, maximum daytime temperatures during incubation should increase less than minimum night-time temperatures (Davy et al., 2017). That is, the pattern of thermal fluctuations at nest depth will probably look increasingly like a double sine pattern, with more extreme highs and dampened lows. Thus, symmetrical changes in variance are unlikely to accurately and precisely capture real nest thermal profiles in the future.

Regardless of the transformation applied, methods that rely on existing thermal profiles as model input are spatially constrained to predicting incubation temperatures for points at which temperature data have already been recorded. No method of direct sampling can capture future incubation conditions. We can see the implication of this limitation most readily in sea turtles, all seven species of which have TSD. Some nesting beaches are already warm enough to incubate embryos at temperatures approaching the high-temperature asymptote, producing 99% female hatchlings (Jensen et al., 2018). In addition, depending on the degree of climate warming actually experienced, beaches that are now identified as critical sea turtle nesting habitat may be inundated by sea-level rise over the coming decades (Fish et al., 2008; Fuentes et al., 2011; Mazaris et al., 2009). Similarly, painted turtles, common snapping turtles (Chelydra serpentina) and other aquatic species that nest terrestrially are experiencing loss of nesting grounds due to development (Lambert and Steen, 2019; Steen et al., 2012). In the case of both marine and freshwater turtles, we could spend substantial effort sampling rookeries that, in a few decades, will no longer exist. We need to be able to predict incubation conditions for locations that are not currently nesting areas. Thus, predictive models that do not rely on contemporary nest thermal profiles to predict sex ratios under warmer climates are key to identifying populations for which the risk of climate change to embryonic development is highest.

Thanks to continental and global scale research networks and remote sensing platforms, the scientific community has access to abiotic datasets at previously unattainable scales that, in some cases, have been extrapolated to relatively fine spatial and temporal resolutions (Abatzoglou, 2013; Maclean et al., 2017). Ideally, we could make realistic predictions about the impacts of climate change on embryonic development at these same scales and for any possible nesting location and climate scenario. However, remote sensing platforms that collect environmental data record primarily air and surface temperatures. In addition, microclimate conditions often vary substantially, over both space and time, from the data available in gridded macroclimate (Faye et al., 2016; Hannah et al., 2014) and microclimate (Kearney et al., 2014b) layers. Climate change is modeled anywhere on a global to regional spatial scale and on a temporal scale of seasons, years or decades. Hourly surface (Faye et al., 2016) and substrate (A.L.C., unpublished data) temperatures can vary on a spatial scale of a few centimeters. Although developing embryos are subject to the effects of their incubation environments instantaneously and at an ultra-fine spatial scale (i.e. the four-dimensional microclimate of an individual nest), the majority of research on TSD has not considered the effects of this spatial mismatch on predictions of embryonic development in wild populations.

Coupled with the high sensitivity of developmental outcomes, the implication of ultra-fine scale, natural variation in substrate temperatures is that the thermal differential among individual nests may be sufficient to shift hatchling sex ratios significantly from coarser (e.g. cohort scale) predictions. Thus, the spatial resolution of both the developmental and microclimate components of TSD models should, if at all possible, match that scale. However, adding in the logistical constraints on collecting real nest thermal profiles creates a trade-off between model resolution and model scale. The general solution for obtaining fluctuating incubation temperatures for unsampled locations or as-yet unrealized climates is to simulate them. Using a sinusoidal function to downscale daily air temperature minima and maxima to hourly values, a method borrowed from the field of building design (Chow and Levermore, 2007), then decreasing the mean and dampening the variance to account for the transfer of thermal energy through the substrate (van Wijk and de Vries, 1963) is one potential means of simulating different climate change scenarios. However, mathematical approaches to simulating hourly temperature fluctuations, even if they are simple to implement, would be difficult to validate in situ and would only be representative of a range of hypothetical situations.

A step farther than simulating substrate thermal profiles as a sinusoidal function is predicting microclimate conditions mechanistically (Kearney and Porter, 2009; Porter et al., 2010). Although more computationally intensive and, arguably, data-hungry than regression analyses or other methods of curve-fitting (Kearney and Porter, 2009), mechanistic models do not rely on existing microclimate data but on prevailing macroclimate conditions and local terrain to calculate microclimate conditions using the equations that describe atmospheric and soil physics (Barton et al., 2018; Kearney and Porter, 2017; Porter et al., 2010). Using mechanistic approaches, we can predict substrate temperatures at locations that have not previously been monitored for reptile nesting activity and for historical or hypothetical future climates (Fig. 2). In addition, models that do not rely on microclimate-scale inputs should not be constrained by scale. Thus, we can predict incubation temperatures at any spatial resolution.

So far, the development of mechanistic models for reptile-focused research has centered primarily on physiology, particularly thermal performance. However, a few studies have used mechanistic approaches to model substrate conditions, specifically within the context of TSD. Mitchell et al. (2008) used the microclimate module of Niche Mapper (McCullough and Porter, 1971; Porter, 2002), perhaps the most-recognized mechanistic modeling framework (authors’ personal observations), to simulate hourly incubation temperatures under climate change, then predicted hatchling sex ratios and hatching phenology using the Sharpe–DeMichele model (Schoolfield et al., 1981) for a population of tuatara (Sphenodon punctatus), one of the few species of reptiles with type Ib TSD (Mitchell et al., 2006). They ran multiple climate-change simulations for 52 actual nests but more than 27,000 grid cells across an entire island at a spatial resolution of 0.5 m, a scope of analysis that would be impossible with direct sampling.

The ability to simulate hourly incubation temperatures to within 2°C of ‘real’ values (Carter et al., 2015, 2018; Kearney et al., 2014b), with no constraints on spatial resolution, is invaluable for predicting hatchling sex ratios at the nest scale. However, although the subroutines were developed nearly 50 years ago, Niche Mapper has not seen wide adoption in TSD research, or in ecology more generally, perhaps due to the relative difficulty in implementing the original Fortran subroutines for researchers with little experience in computer science (authors’ personal observations). However, with the development of the ‘NicheMapR’ package for R (Kearney and Porter, 2017) and its integration with multiple physiological models (Kearney and Porter, 2020), we anticipate increased interest in both it and other mechanistic approaches. The main constraint on running mechanistic models at ultra-fine spatial resolutions is computational. To overcome that challenge, Carter et al. (2018) simulated substrate temperatures for three climate scenarios at five soil depths using the R-implementation of the Niche Mapper subroutines, at 500 random points (i.e. potential nests), a manageable sample size for a personal computer. They then used a geostatistical approach to interpolate CTEs, again estimated using the Sharpe–DeMichele development model, for tuatara across a 1 m grid containing more than 480,000 cells (Carter et al., 2018). Although the amount of error in interpolated CTEs was less than 2°C in most cases, that is high enough to drastically affect sex ratio predictions, especially for species with a relatively narrow TRT.

At this point, we are left with a mounting body of studies on thermal developmental plasticity in reptiles, much of it motivated by decades-old theory and experimental convention, now confronted (often awkwardly) by advanced technological and modeling developments. Our perspective is that future research efforts on this increasingly critical topic would be most productive if implemented with these advances at the forefront. We further this conversation by first providing an overview of research on thermal developmental plasticity, specifically TSD in reptiles, over the last 50 years, including a review of the literature, focusing on experimental approach.

We accessed the Reptile Development Database (RepDevo; https://github.com/RepDevo/ReptileDevelopmentDatabase?), a collection of research on thermally mediated developmental plasticity in oviparous reptiles (Noble et al., 2018), to identify relevant literature. Rather than relying on the database itself, which uses relatively strict inclusion/exclusion criteria, we filtered the associated ‘citations.csv’ file to identify studies from the past 50 years in which either (1) eggs were incubated in the laboratory through at least the middle-third of development, and multiple temperature treatments were reported, regardless of the trait(s) being measured, and/or (2) incubation temperatures were measured in nests in situ. We filtered the citations file, which includes papers published from 1969 to 2016, as well as a list of studies published from 2016 to 2018 that have been tagged for potential inclusion in a future release of RepDevo (Lisa Schwanz, personal communication). In addition, we conducted a search in Web of Science version 5.35 using the RepDevo search terms (Noble et al., 2018) for papers published from 2019 to mid-2020.

From those lists, we excluded papers that did not focus on an oviparous, non-avian reptile; were conference abstracts; were laboratory studies that did not report at least two temperature treatments; incubated cells or tissues, rather than whole eggs; or were unavailable or otherwise inaccessible or irrelevant (Noble et al., 2018). Although we excluded papers with unclear methods, we did not exclude studies that had used a flawed design or analysis (e.g. variable confounding or low sample sizes) or on the overall study aims or approach (e.g. simulations of climate change using temperature shift experiments or hormonal manipulations). These liberal inclusion criteria allowed us to create a timeline of research focus and intent that, even if some studies were not necessarily rigorously executed, provides a more complete picture of the field over time. Including papers in RepDevo, our final bibliography lists 504 papers published from 1974 to 2020, 247 of which measured or modeled hatchling sex ratios (see ‘Data availability’ section).

For each paper identified, we classified temperature treatments as ‘constant’, ‘fluctuating’, ‘irregular’ or ‘shift’ to differentiate among laboratory methods (Noble et al., 2018). Laboratory protocols including both constant and fluctuating treatments were labeled as ‘fluctuating’, as long as at least two fluctuating treatments were reported. We also tagged papers that modeled or directly measured hatchling sex, regardless of whether the study species is known to have TSD or whether sex was analysed explicitly as a dependent variable in the study. For example, a study on the effects of incubation temperature on size dimorphisms would still be required to record gonadal development, even though ‘sex’ would be an independent, rather than dependent, variable in the analysis. Finally, we searched all titles and abstracts for the keyword ‘climate change’ to identify studies that examined embryonic development and/or hatchling sex ratios specifically within the context of climate change.

Nearly five decades of research have sought to understand and predict the phenotypic outcomes of temperature-dependent embryonic development and sex determination in reptiles. However, the proportion of studies using non-constant temperature treatments, including those that have measured nest temperatures directly, has not increased appreciably (Fig. 3). Notably, 2019 was a bumper year for research on thermal developmental plasticity in reptiles, with around a third of studies that sampled hatchling sex using fluctuating temperature treatments. Overall, however, the field has not yet taken full advantage of the available incubation technologies and modeling approaches. Here, we offer some suggestions for how to push forward research on TSD and, more broadly, thermal developmental plasticity in reptiles, to tackle the challenges presented by a warmer and increasingly variable climate. Our observations and recommendations add to, and at times echo, those already made by Carter et al. (2019), Du et al. (2019), Refsnider et al. (2019) and Sears et al. (2019).

One of the main limitations on our understanding of TSD in reptiles is that key model parameters, such as T_piv, have been estimated using data from constant-incubation experiments that do not replicate the complex developmental outcomes of actual microclimate conditions. This constraint is at least partially technological in origin. The preceding 30 years have seen incubators become increasingly smaller, more accurate and precise, and less expensive. However, the large size, high power consumption and prohibitive cost of programmable environmental chambers that can reliably reproduce fluctuating incubation conditions has remained a logistical constraint for empirical work on TSD. Even when and where researchers have access to state-of-the-art climate chambers, their still-significant space requirements limit study design and opportunities for replication. As noted by While et al. (2018), empirical studies on thermal developmental plasticity in reptiles often spread a large number of eggs across very few incubation treatments (mean 4.3 per paper; While et al., 2018), suggesting that technological constraints are significantly impeding our ability to model thermal reaction norms, even at constant temperatures. As a potential alternative to unwieldly climate chambers, Greenspan et al. (2016) developed a method for converting a small, inexpensive constant-temperature incubator to a fluctuating incubator, capable of reproducing sinusoidal and stochastically fluctuating temperature profiles, using an open-source microcontroller. This or a similar approach to incubation in TSD research could greatly expand our ability to link nest thermal profiles and developmental outcomes in a statistically rigorous, reproducible way.

Another longstanding logistical challenge for research on TSD is the lack of a reliable, non-lethal and universal method for visually assessing gonadal development in very small hatchlings. Coelioscopy, a surgical technique in which a laparoscope is inserted through a small incision to inspect the gonads inside the body cavity, has been validated for marine turtles (Wibbels et al., 2000; Wyneken et al., 2007) and several species of juvenile tortoises that weigh as little as 10 g (Divers, 2015). Notably, this method has also been used successfully for visual inspection of gonads in juvenile red-eared sliders (Trachemysscripta), one of the most-studied freshwater turtles with TSD. However, gonadal histology in painted turtles and other freshwater species with hatchlings that are 5 g or smaller still requires euthanasia. This limitation presents an ethical and logistical challenge that limits sample sizes and slows research. In addition, we cannot validate TSD models for threatened species, nor should we destructively sample species that are cryptic or data deficient, as lethal sampling is antagonistic to conservation goals (Wyneken et al., 2007). Although minimally invasive immunoassays and blood hormone analyses have been developed for some species (Tezak et al., 2017, 2020; Wibbels et al., 2000), their application seems so far to be species-specific, and blood sampling is also size-limited (Owens, 1999; Owens et al., 1978) and potentially cost-prohibitive. Thus, the need to sex very small hatchlings and post-hatchlings remains a significant barrier to advancing research on TSD. Even so, just as this field has witnessed technological advances in nest temperature monitoring and incubator design, which have the potential to revolutionize and democratize research, substantial developments in sexing methods can be expected.

Although the wide distribution and relatively high abundance of some species with TSD have certainly been beneficial to study efforts, allowing researchers to amass invaluable, decades-spanning datasets on the same populations, sampling biases may have introduced unknown complications to the field. While et al. (2018) identified a substantial bias in the geographical foci of research on thermal developmental plasticity in reptiles, with the greatest number of studies – by far – published in the USA, followed by Australia. Not surprisingly, they also found a significant taxonomic bias for common North American species. Thus, much of our understanding of TSD in reptiles is based on studies with common aquatic turtles (painted turtles, red-eared sliders and common snapping turtles) and lizards (leopard geckos Eublepharis macularius and three-lined skinks Bassiana duperreyi) (While et al., 2018). We would also add marine turtles to the list of taxa that could be considered models for TSD research, because much of the most recent analytical work has used data from sea turtles, especially loggerheads (Caretta caretta) and green sea turtles (Chelonia mydas), to fit development models (Abreu-Grobois et al., 2020; Bentley et al., 2020; Girondot and Kaska, 2014; Girondot et al., 2018; Stubbs and Mitchell, 2018).

Regardless of geographical distribution, the model species bias in TSD research is probably due at least partly to political, economic and other logistical considerations. Research on TSD would benefit directly from supporting initiatives that seek to remove barriers to access for the global research and conservation-management communities. For example, paywalls, bundle subscriptions and other manifestations of the for-profit publication model substantially increase the difficulty of accessing the most recent papers for anyone who is not associated with a wealthy university or other institution. The English-language bias of science, generally, decreases accessibility and increases publication costs even further for researchers based outside of primarily English-speaking countries. Publishing in and accessing society journals, although they are less profit-motivated than large open-access journals, requires paying annual membership fees. Perhaps a no more relevant example of publication bias can be found than in the first known report of TSD, which was published by Madeleine Charnier – in French – in the local meeting records of a West African biological society (Charnier, 1966). An exciting research finding thus went largely unnoticed and was not even cited in the first study on TSD to receive widespread attention, published more than a decade later by researchers based at US institutions (Bull, 2004). Practices that support equitable access to data, empirical and analytical tools, funding and publication outlets are beneficial for global research on TSD.

Our call to use modeled microclimate temperatures to drive developmental models for reptiles with TSD does not imply that ‘real’ substrate temperatures have become less important to the field. In contrast, the ability to mechanistically model incubation conditions for any potential nesting location and any potential climate scenario frees us to use recorded thermal profiles as a means of model validation, rather than as model inputs. In addition to collecting thermal profiles from known nesting areas, we suggest identifying and monitoring ‘incubation conditions’ in potential nesting areas with a variety of habitat types, with the aim of determining whether populations of reptiles with TSD will have sufficient microrefugia from climate warming available. To the extent possible, we should leverage existing research networks, such as the Long-Term Ecological Research (LTER) network and National Ecological Observatory Network (NEON) to deploy thermal probe arrays at sites with existing weather and climate stations and flux towers. We should also continue to monitor known nesting sites, particularly those for which there are existing long-term datasets, with a focus on collecting robust data on fine-scale incubation temperatures and embryonic development both from nests and within broader nesting areas. Both are needed to identify whether warming incubation conditions are likely to increase populations’ risk of experiencing significant shifts in the timing, rate and/or developmental outcomes of embryonic development. Long-term data collection offers the best chance of responding proactively to climate change (Reinke et al., 2019).

Janzen and Paukstis (1991), in the most-cited paper on both thermal developmental plasticity in reptiles and on TSD specifically (While et al., 2018), called for the standardization of methods used for constant-incubation experiments and the associated statistical analyses for research on TSD within an eco-evolutionary framework. Nevertheless, issues with study design and data availability remain prevalent in the literature (Noble et al., 2018), a situation that is likely to be due in part to logistical constraints but nevertheless decreases opportunities for data sharing and other forms of collaboration that could move the field forward more rapidly. Thus, we find it necessary to renew the call for standardization of experimental approaches to understanding TSD, this time within an additional context of contemporary climate change. Technological advancements over the last three decades could eliminate our previous reliance on constant-incubation experiments, with the caveat that new technologies would have to be equitably distributed and adopted to have a global impact on TSD research. Given that, our aim should be to make constant-incubation experiments obsolete (Bowden et al., 2014), even as such studies might be more easily replicable and can enhance comparability of results among populations and species.

We should also move away from thinking about TSD in terms of discrete thermal parameters like T_piv and TRT. These parameters, which are analogous to the slope and intercept of the TSD reaction norm, are not feasibly identifiable parameters at the individual embryo scale but only at the clutch or coarser levels. Gonadal differentiation in reptiles with TSD is not differentiable from embryonic development. Biophysical and/or mechanistic developmental models are more likely to generate realistic predictions than correlative sex ratio models, particularly when paired with microclimate-scale thermal profiles (Mitchell et al., 2008). The decision to coarsen data on developmental outcomes, for example, from gonadal differentiation of individual hatchlings to clutch to nest cohort sex ratios, should be made cautiously and thoroughly justified. Similarly, decreasing the spatial and/or temporal resolution of nest thermal profiles allows us to compare a gross measure of local temperature across study sites. However, given the ultra-fine spatial scale and fine temporal scale over which thermal profiles vary, we should interrogate whether these thermal comparisons are indicative of ecologically relevant impacts on embryonic development or merely statistical differences. We should also question the value of using statistical tests alone to assess whether incubation temperatures or phenotypic outcomes are ‘different’ from each other. There is no a priori reason to assume that a statistically significant difference, or lack thereof, among populations or years is indicative of ecological importance or non-importance. Interpretations of statistical tests (i.e. their biological or ecological meaning) should be made explicit, never assumed based on statistical significance alone, and raw data should be deposited in an accessible format (Costello and Wieczorek, 2014).

Every new incubation technology or modeling approach has moved us farther down the path from purely associative studies to a fully mechanistic, integrated understanding of the relationship between temperature, embryonic development and developmental plasticity in reptiles with TSD. Research on TSD in reptiles is well on its way to achieving that understanding in the physiological and molecular dimensions, especially in model species (Czerwinski et al., 2016; Dorizzi et al., 1991; Weber et al., 2020; Yatsu et al., 2016). However, the rate at which the microclimatic dimension has been incorporated into those models has lagged behind. Although mechanistic models remain computationally intensive (Kearney et al., 2014a), high-performance computing – and the skills required to implement it – is becoming increasingly available (Carter et al., 2018). Ultimately, our ability to predict the impacts of climate change on embryonic development and, importantly, to scale those consequences up to populations and species will require integrating all of the diverse technological, empirical and theoretical expertise in the global herpetology community.

We thank the Journal of Experimental Biology for the invitation to contribute to this themed issue and Lisa Schwanz for valuable discussion about the RepDevo project. We are also grateful to the many students and colleagues who have helped shape and hone our perspectives on these increasingly pressing matters through the years. This is Kellogg Biological Station contribution number 2175.