Feminization of hawksbill turtle hatchlings in the twenty-first century at an important regional nesting aggregation

Projected climate change is forecasted to have significant effects on biological systems worldwide. Marine turtles in particular may be vulnerable, as the sex of their offspring is determined by their incubating temperature, termed temperature-dependent sex determination. This study aimed to estimate historical, and forecast future, primary sex ratios of hawksbill turtle Eretmochelys imbricata hatchlings at an important nesting ground in northeastern Qatar. Incubation temperatures from the Arabian/Persian Gulf were measured over 2 nesting seasons. Climate data from same period were regressed with nest temperatures to estimate incubation temperatures and hatchling sex ratios for the site from 1993 to 2100. Future hatchling sex ratios were estimated for 2 climate forecasts, one mid-range (SSP245) and one extreme (SSP585). Historical climate data showed female-biased sex ratios of 73.2 ± 12.1% from 1993 to 2017. Female biases from 2018 to 2100 averaged 85.7% ± 6.7% under the mid-range scenario and 87.9% ± 5.4% under the high-range scenario. In addition, predicted female hatchling production was >90% from 2054 and 2052 for SSP245 and SSP585, respectively. These results show that hawksbill primary sex ratios in Qatar are at risk of significant feminization by the year 2100 and that hawksbill turtle incubation temperatures in an extreme, understudied environment are already comparable to those predicted in tropical rookeries during the latter half of the 21st century. These results can help conservationists predict primary sex ratios for hawksbill turtles in the region in the face of 21st-century climate change.


INTRODUCTION
Climate change has been shown to have a significant effect on and pose a substantial risk to ecosystems, taxa and ecological processes around the world (Walther et al. 2002, Hays et al. 2003, Fuentes et al. 2010, Muhling et al. 2011, Riegl & Purkis 2012, Alongi 2015, van Hooidonk et al. 2015 Marine turtles are at particular risk from climate change since hatchling sex is determined by their incubation temperature (temperature-dependent sex determination, TSD; Yntema & Mrosovsky 1980, Mrosovsky et al. 1992, Ackerman 1997, Godfrey et al. 1999, Godley et al. 2002, LeBlanc et al. 2012b). The pivotal temperature (PT), approximately 29°C in most marine turtle species, is the incubating temperature at which a clutch will produce a 1:1 male: female sex ratio (Gross et al. 1995, Godley et al. 2002. However, recent research has shown that the PT may vary considerably between species, populations and even between clutches (Howard et al. 2015, Bentley et al. 2020. PTs between different flatback turtle Natator depressus populations in western Australia varied up to 1.5°C and were as high as 31. 1°C (Bentley et al. 2020). While sea turtles show resilience to climate variability and heightened temperatures, the magnitude of forecasted climate change in the 21 st century is likely to induce significant in creases in incubation temperatures, causing primary sex ratios to skew (Booth & Freeman 2006, Tezak et al. 2018. While previous research has predicted hatchling sex ratios using climate data (Hays et al. 1999, 2003, Hawkes et al. 2007, Fuentes et al. 2010, Askari Hesni et al. 2016, Laloë et al. 2016, Monsinjon et al. 2019b, Patrício et al. 2019), relatively few studies exist from the extremities of global marine turtle distribution. The Arabian/Persian Gulf (hereafter referred to as 'the Gulf') is an extreme environment in which hawksbill turtles Eretmochelys imbricata nest (Pilcher et al. 2014a, Pazira et al. 2016, Chatting et al. 2018, Rees et al. 2019. The sea is a shallow basin with mean and maximum depths of ~35 and 160 m, respectively, and regularly experiences summer sea surface temperatures (SSTs) in excess of 36°C (Price et al. 1993, Sheppard et al. 2010. The region also experiences extreme seasonal heating and cooling cycles (Price et al. 1993), and many species persist on the threshold of their physiological tolerances (Hofmann & Somero 1995, Pörtner & Knust 2007, Riegl & Purkis 2012. Hawksbills in the Gulf display reduced fecundity compared to tropical populations; for example, it has been reported that they lay fewer eggs per clutch and may nest for a shorter period each year (Chatting et al. 2018). Despite many studies reporting on the nesting ecology of hawksbills in the region (Pilcher 1999, Pilcher et al. 2014b, Pazira et al. 2016, Chatting et al. 2018, Rees et al. 2019, incubation temperatures have received little attention, and it is still not known how incubation temperatures in this extreme environment compare with other regions. In addition, the Gulf is expected to experience significant climate warming in the 21 st century (Pal & Eltahir 2016). As such, extreme temperatures may already be skewing primary sex ratios, and climate change is likely to pose a serious threat to marine turtles in the region throughout the 21 st century.
The current study represents the first published account of hawksbill incubation temperatures in the Gulf. This study also aimed to predict incubation temperatures from air temperatures and use this relationship to statistically model primary sex ratios from 1993 to 2100. The results reported here can add to global understanding of hawksbills and can inform conservation and management practices to help mitigate against the potential effects of climate change.

Study site
Nesting site incubation temperatures were measured along the northeastern coast of Qatar in the 2016 and 2017 nesting seasons. This area of Qatar has a well-documented history of hawksbill turtle nesting and represents a significant nesting aggregation for the species in the Gulf (Pilcher 1999, Tayab & Quiton 2003, Pilcher et al. 2014a, Chatting et al. 2018. As this work was part of a wider conservation program, 22 of the total number of clutches in which incubation temperatures were recorded (n = 26, see Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ n044 p149 _ supp. pdf) were relocated to protect against tidal flooding or predation. Clutches were laid anywhere from 26.1275°N, 51.3181°E to 26.0047°N, 51.4005°E along the northeastern coast of Qatar (Fig. 1a) and were relocated to a protected hatchery site on Fuwairit Beach (26. 0311°N,51.3752°E,Fig. 1b).
Incubation temperatures were measured during the 2016 (n = 2) and 2017 (n = 24) nesting seasons using TidbiT v2 loggers with a resolution to 0.02°C and an accuracy of ± 0.21°C. Relocation efforts were only made if clutches were thought to have been deposited <12 h earlier to minimize potential damage to developing embryos (Ahles & Milton 2016). When relocating clutches, eggs were placed in relocation boxes to minimize disturbance. Data loggers were deployed into the middle of the clutch immediately post relocation. The middle of the clutch was determined by counting the eggs upon excavation, replacing half in the nest, deploying the temperature logger and then replacing the second half of the eggs. If clutches were left in situ, eggs were excavated, counted and then placed back into the nest.
Loggers were deployed during this process after half of the eggs had been returned to the nest. All contact with clutches was minimized and monitored by the Qatar Ministry of Municipalities and Environment to ensure ethical treatment of clutches. Prior to engaging in this work, all participants were trained in the handling and relocating of marine turtle eggs. Mean daily incubation temperatures for each of the 26 nests in which loggers were deployed were calculated. Nests where loggers were deployed spanned all 4 months of the annual nesting period in Qatar (April to July) (Chatting et al. 2018) and the loggers recorded temperatures every 4 h.

Nest temperature reconstruction
Recorded incubation temperatures were regressed with climate data to reconstruct an hourly time series of past and future nest temperatures based on a previously published method (Girondot & Kaska 2015, Monsinjon et al. 2017, 2019b, Laloë et al. 2020. Generalised linear models (GLMs) with a Gaussian distribution were used to predict daily mean nest temperature and daily nest temperature amplitude sep a rately. Increasing lag between successive nest temperature measurements and predictors (from 1 to 15 d) were contrasted to determine the best model for daily mean nest temperature (Girondot & Kaska 2015, Monsinjon et al. 2017). Predictors used for daily mean nest temperature were daily mean air temperature at 2 m height (TAS) and site-adjacent daily SST as reported in previous studies (Girondot & Kaska 2015, Monsinjon et al. 2017. The best predictors of TAS, SST and daily maximum (T max ) and minimum (T min ) temperatures were used to model daily nest temperature amplitudes. Multicollinearity was addressed by removing explanatory variables with a variance inflation factor > 3.3 (Kock & Lynn 2012). Models that best fit the data were selected by the lowest Akaike's information criterion (AIC) value (Zuur et al. 2009). Daily nest temperature maximum and minimum (daily mean ± amplitude/2) as well as time of day they were recorded (minimum: 02:00 h and maximum: 14:00 h) were then used to model a sinusoidal function of hourly time-series nest temperatures using the 'HelpersMG' package in R (Girondot 2017). Daily mean TAS and SST directly seaward of the nesting beach have previously been used to predict daily nest temperatures (Monsinjon et al. 2017). TAS, T min and T max datasets were obtained from the National Center for Environmental Prediction (NCEP, www. cpc. ncep. noaa. gov) (Kanamitsu et al. 2002), and SST was ob tained from the National Oceanic and Atmospheric Administration National Climatic Data Centre (NOAA NCDC) (https:// psl. noaa.gov) (Reynolds et al. 2007). Climate datasets (daily TAS, T min , T max and SST) were extracted from NCEP and NOAA data for the relevant location in Qatar from 1993 to 2017 to reconstruct past incubation temperatures. The entire NCEP dataset ranges from 1979 to 2017; however, missing data in the years 1979−1992 prevented extraction from these years. Precipitation, which has previously been reported as an influential factor in nest temperatures (Wyneken & Lolavar 2015), was excluded, as Qatar is a desert nation that experiences a general lack of rainfall (Sheppard et al. 2010). A mechanistic microclimate model, which uses weather data, local topography and soil thermal properties, has previously been shown not to have greater predictive power than correlative models using TAS and SST to infer incubation temperatures and so was not used in the current study (Laloë et al. 2020). To predict future incubation temperatures from 2018 to 2100, forecasted climate data in Fuwairit Beach, northeastern Qatar, were extracted from the Coupled Model Intercomparison Project phase 6 (CMIP 6) climate forecasts. Shared socio-economic pathways (SSPs) 245 and 585 (https:// esgf-node. llnl. gov/ projects/ cmip6/) were selected. SSP245 and SSP 585 represent mid-and high-level climate scenarios, respectively. Climate datasets (daily TAS, T min , T max and SST) were extracted from an ensemble of 6 climate projections ( Table S2 in the Supplement) for both scenarios. The climate ensemble was bias-corrected for each location using historic observations and hindcast climate data from the same model ensemble using the equation (Luo et al. 2018): where Cor T min (d) represents corrected T min for the d th day. Hist T min and Obs T min refer to hindcast historical and observed T min data, and μ is the mean.
Weighting coefficients for each bias-corrected ensemble member were then calculated following Muhling et al. (2011) (Table S2). From weighting coefficients, a mean weighted climate forecast dataset was calculated for daily TAS, T min , T max and SST. Incubation temperatures and 95% CIs were predicted on an hourly scale up to the year 2100.

Modelling embryonic development of hatchlings
Sex of marine turtle hatchlings is determined in the middle third of development, which is termed the thermosensitive period (TSP) (Girondot 1999, Godley et al. 2002, Wyneken et al. 2007, LeBlanc et al. 2012b, Girondot & Kaska 2014, Girondot et al. 2018. The timing of the TSP is not constant and changes depending on incubation temperature. It can be estimated from the thermal reaction norm model (Girondot & Kaska 2014, Fuentes et al. 2017, Girondot et al. 2018) using the 'embryogrowth' (Girondot 2019) package in R. The thermal reaction norm model was fit using maximum likelihood estimation from incubation data collected by loggers deployed in the field. Confidence intervals of the thermal reaction norm were generated by using Monte Carlo simulations over 10 000 iterations. From the thermal reaction norm, hatchling growth throughout the incubation process can then be modelled using the equation (Fuentes et al. 2017): where X(t) represents the size of a developing hatchling at time t, and X (0) where T i is the temperature at any time i during incubation, W i is the growth from time i to i + 1, and N represents the length of incubation in the same time unit as i. A metabolic heating coefficient was added to CTE w values to account for the increase in incubation temperatures during the incubation process with respect to the surrounding sand (Monsinjon et al. 2017). The metabolic heating value was calculated by assuming a final difference of nest and surrounding sand temperature of + 2.05°C at the time of hatching (Monsinjon et al. 2017) and that metabolic heating only affected sex ratios up to the end of the TSP. Hence metabolic heating follows the developmental curve modelled by Eq.

Annual sex ratio estimation
It was assumed that the timing and distribution of the annual nesting period remained the same each year from 1993 to 2100 and lasted from the start of April to the middle of June and peaked in the first 2 wk in May (Chatting et al. 2018). CTE w and 95% CI were calculated from reconstructed hourly incubation temperatures for each day of the nesting season starting on 1 April until 15 June for each year. This approach captured 76 unique incubation periods annually and their corresponding CTE w . A logistic regression model using a quasibinomial distribution predicted sex ratio from each unique CTE w and its associated CIs to determine the proportion of each clutch that was female from 1993 to 2100. As no previous efforts have investigated PTs in hawksbills in the Indian Ocean, the logistic regression model was developed by combining observations from Godfrey et al. (1999) and Mrosovsky et al. (1992), the only 2 previous published attempts to develop such a model for hawksbill turtles (Fig. S1 in the Supplement). The PT (mean ± SE) for this model was 29.8 ± 0.1°C while the transitional range of temperature was from 28.6 to 30.8°C. Previously, incubation temperatures > 35°C have been assumed to be the lethal limit for marine turtle clutches (Valverde et al. 2010, Howard et al. 2014). This assumption was not applied in this study, as there was no decrease in hatching success in 7 out of the 26 clutches monitored (M. Chatting unpubl. data) where incubation temperatures regularly exceeded 35°C. To calculate total annual sex ratios, a normally distributed frequency of nesting over each nesting season was assumed from 1 April to 15 June. The results of the logistic regression were then weighted depending on when in the season they occurred, and the sum of all of these weighted sex ratio estimates was the total sex ratio for the year (Fig. S2 in the Supplement). This was repeated for upper and lower CI bounds to propagate uncertainty levels. All statistical and numerical analyses were performed in R version 4.0.2.

RESULTS
Mean ± SE measured daily clutch temperature from April to June was 31.1 ± 0.1°C and ranged from a minimum to maximum mean daily temperature of 25.9 and 34.9°C, respectively. Mean monthly incubation temperatures showed an increase throughout the season from April (26.9 ± 0.2°C) to May (29.7 ± 0.1°C), June (32.1 ± 0.05°C) and July (33.8 ± 0.1°C). The thermal reaction norm exhibited by hawksbills in Qatar showed an increase in embryonic growth rate up to 35.1°C, after which the rate decreased (Fig. 2). The logistic model used to predict proportion female from TSP incubation temperature showed a PT of 29.6°C and a transitional temperature zone (the range of temperatures in which a mixed sex ratio occurs) from 28.6 to 30.8°C (Fig. S1).
Annual mean TAS and SST showed an increasing trend from 1993 to 2100 (Fig. 3). TAS and SST averaged 26.8 and 26.4°C, respectively, in 1993. By 2100, they increased to 29.7°C (TAS) and 28.5°C (SST) under the mid-range scenario (SSP245, Fig. 3a) and to 32.4°C (TAS) and 31.6°C (SST) under the highrange scenario (SSP585, Fig. 3b). Historical (1993− 2017) reconstructed hourly incubation temperatures showed slight increases in maximum temperatures from highs of 36.9°C in 1993 to highs of 38.3°C in 2017. Post-2017 increases were more significant; highs > 39°C were detected every year from 2062 onwards under SSP245 and even earlier (2039) under SSP585. The maximum incubation temperature predicted under SSP245 was 40.5°C in 2095, while under SSP585, hourly incubation temperatures reached a maximum of 43.7°C in 2098. Mean annual CTE w prior to 2018 was 30.9°C (Fig. 3), while forecasted CTE w values increased by 1.8°C (32.7°C) and 2.7°C (33.7°C) under mid-(SSP245) and high-range (SSP585) scenarios from 2018 to 2100. A significant female bias was detected throughout the entire study period (1993−2100) (Fig. 4). Two years out of the 107 yr study period were male biased (< 50% female hatchling production). Data also showed a steady increase in female hatchling production from 1993 to 2100 (Fig. 4). Historical data revealed that annual female hatchling proportions ranged from a minimum of 49.0 ± 22.6% to a maxi- Fig. 2. Thermal reaction norm used for modelling hawksbill turtle embryo growth throughout the incubation process. The solid line is the mean growth rate and dashed lines are 95% confidence intervals mum of 89.5 ± 8.1% prior to 2017 (Fig. 4). This resulted in a female bias of historical primary sex ratios of 73.2 ± 12.1% for the first 25 yr of the study period (1993). The following 25 yr of the study period (2018−2043) remained comparable under both scenarios (SSP245: 73.5 ± 9.9%, Fig. 4a; SSP585: 75.2 ± 10.3%, Fig. 4b), while the final 25 yr of the study period (2075−2100), resulted in the greatest increase in female hatchling production (SSP245: 95.1 ± 2.9%; SSP585: 98.6 ± 0.3%). Over 90% annual female hatchling production occurred regularly by 2054 and was predicted to occur in 41 out of 107 years (38.3%) under the mid-range scenario (SSP245), while it occurred earlier under the high-range scenario (2052) and was predicted in 47 years (43.9%) between 2018 and 2100.

DISCUSSION
Results from this study predict a significant increase in the proportion of female hatchlings pro-duced in Qatar in the 21 st century, an important nesting ground in the Gulf (Tayab & Quiton 2003, Pilcher et al. 2014a, Chatting et al. 2018). These results support other predictive studies forecasting significant feminization of primary sex ratios in marine turtle species around the world (Fuentes et al. 2009, Laloë et al. 2016, Reneker & Kamel 2016, Monsinjon et al. 2019b, Patrício et al. 2019). An almost 'complete feminization' of green turtles Chelonia mydas in the northern Great Barrier Reef has been predicted (Jensen et al. 2018, p. 154). In the northeast Caribbean, 0.4% male green turtle hatchling production has been predicted by 2090 (Laloë et al. 2016), and 76− 93% female hatchling production is predicted in West Africa (Patricio et al. 2019). Our results showed current incubation temperatures at the end of the nesting season in the Gulf (July mean: 33.8 ± 0.1°C) to be comparable to the 33.0°C forecasted in the Caribbean by the year 2060 (Laloë et al. 2016). In addition, maximum incubation temperatures to which clutches were exposed in the present study (36.8°C) were higher than incubation temperature predictions for the Caribbean by the year 2090 (34.2°C) (Laloë et Fig. 3. Annual constant temperature equivalents weighted by hawksbill turtle embryo growth (CTE w ) and 95% CIs (thick black lines with dashed lines), annual air temperature (red lines) and annual mean sea surface temperature (SST, blue lines) directly seaward of the nesting beach for (a) the mid-range scenario SSP245 and (b) the high-range scenario SSP585 Fig. 4. Forecasted annual mean and 95% CIs for the proportion of female hawksbill turtle hatchlings from 1993 to 2100 for (a) the mid-range scenario SSP245 and (b) the highrange scenario SSP585 al. 2016). However, unlike in the Gulf, some recent work has found male-biased or balanced sex ratios from nesting sites at higher latitudes and in tropical areas where thick vegetation shaded nests, providing cooler incubating temperatures (Steckenreuter et al. 2010, LeBlanc et al. 2012a, Kamel 2013, Patrício et al. 2017, Jensen et al. 2018. Recent research of green turtles in the Great Barrier Reef traced adults back to the beach on which they incubated and found that adult males originated from sites with cooler temperature profiles (Jensen et al. 2018). These regional male-producing nesting grounds may become increasingly important for sea turtle resilience to climate change.
Annual female hatchling production in the first half of the 21 st century was comparable to historical estimates, whilst the latter half of the 21 st century showed a significant increase in female sex ratios. Prior to 2040, annual female hatchling production was ~75% over both scenarios, which is comparable to hawksbill juvenile ratios found in the Dominican Republic (2.71:1, F:M) (Leon & Diez 1999). In contrast, adult sex ratios in Qatar have been reported to be 1:4 (F:M) (Pilcher et al. 2015). It is still unclear what the ideal operational sex ratio is for marine turtle species. Recent work has shown that males breed more frequently than females; subsequently, fewer males are needed to reproduce in order to sustain the population (Hays et al. 2010(Hays et al. , 2014. It has also been suggested that small increases in female hatchling sex ratios would in crease the number of nesting females and subsequently could boost population numbers (Hays et al. 2017, Patrício et al. 2019. The majority of years in the second half of the 21 st century produced a > 90% female ratio. Since hawksbill nesting in Qatar has historically been from spring into summer, this period has provided a wide enough range in temperatures suitable for producing males. However, such high female proportions by the year 2100 could result in insufficient numbers of males to sustain the population. Moreover, the high incubating temperatures we forecast would likely impede clutch development and significantly reduce hatching success (Fisher et al. 2014, Santidrián Tomillo et al. 2015, Rafferty et al. 2017. Recent tracking work has also shown that hawksbills in the Gulf migrate relatively short distances compared to other marine turtle species (Pilcher et al. 2014a, Rees et al. 2019). Significant feminization of hatchlings in Qatar, coupled with limited migration distances of adults, would mean nesting grounds within the region would need to produce enough males to sustain the entire Gulf hawksbill population.
In estimating historical and future primary sex ratios, numerous assumptions had to be made from previous work. It was assumed that the PT is within the range of those reported by Mrosovsky et al. (1992) and Godfrey et al. (1999). However, other marine turtle species have shown PTs to vary within species, and population-specific PTs are likely to exist (Wibbels 2003, Howard et al. 2015, Bentley et al. 2020. Local adaptations of the TSD reaction norm may be present in hawksbills of the Gulf, which is possible, as 7 of the 26 nests where temperatures were recorded experienced incubation temperatures > 35°C, the previously considered upper limit for thermal stress in incubating marine turtle eggs, with no detriment to hatching success (Valverde et al. 2010, Howard et al. 2014). Likewise, under higher temperatures, the TSP may be shorter than the those estimated in the current study, which would affect sex ratio estimates (Fuentes et al. 2017). Future research should investigate the PT of hawksbills in the region to reduce the uncertainty on estimates provided in this study.
We also assumed that the annual timing and distribution of nesting in Qatar would remain constant throughout the 21 st century, but this may not be the case, as marine turtles may change their nesting phenology in response to environmental variations and climate change (Mazaris et al. 2013, Patel et al. 2016, Monsinjon et al. 2019a). However, as the region experiences a huge change in climate throughout the year, many biological processes are driven by environmental change (Price et al. 1993). As other tropical populations may rely on spatial heterogeneity of nesting grounds to produce male and female hatchlings, temporal plasticity in nesting periodicity may be a more important driver for sex ratios in the Gulf. However, long-term annual nest timing datasets are required to confirm this. In addition, the underlying assumptions of the approach used assume that the relationship between climate data (SST and TAS) and incubation temperature will remain the same throughout the 21 st century. Climate change is forecasted to occur at an unprecedented rate in the region (Pal & Eltahir 2016), and this relationship may not remain constant. Mechanistic models that include microclimate and topography data could offer better predictive power, but this is not always the case (Laloë et al. 2020). The approach used in the current study has previously been used to forecast marine turtle incubation temperatures (Patrício et al. 2019) and can still inform conservationists today of the potential impacts of climate change. In addition to these assumptions regarding modelling the in -creasing feminization of primary sex ratios, sea level rise may also be a significant threat to marine turtle rookeries in the region. Clutches in Qatar are already being relocated to avoid inundation as part of a wider conservation program.
Despite these assumptions, these results can be important for conservationists throughout the region and globally. Our study provides information on hawksbill turtles in an extreme environment that already experiences incubation temperatures comparable to forecasted temperatures elsewhere in the world. This study can help conservationists in the region mitigate against potential impacts of climate change by monitoring incubation temperatures and primary sex ratios going forward into the 21 st century, as well as any possible impacts from climate change, for example, reductions in hatching success or temporal shifts in nesting phenology. Our results also highlight important directions of research in the region. For example, determining a specific PT for hawksbills within the Gulf would help reduce uncertainty in sex ratio estimates provided in this study. Recent work has shown that hatchlings can be reliably sexed with just a small blood sample (Tezak et al. 2020). This new approach may allow conservationists to track future primary sex ratios and may help detect extreme sex ratio biases, to further inform conservation actions. Investigations into the thermal properties on more nesting grounds across the region would help identify potential key male-producing sites that would need to be protected to ensure future population viability. Continuation of long-term reproductive monitoring projects can also help identify any timing shifts in annual nesting behaviour.