Embryo deformities and nesting trends in Kemp’s ridley sea turtles Lepidochelys kempii before and after the Deepwater Horizon oil spill

Kemp’s ridley sea turtles Lepidochelys kempii were disproportionately affected by the Deepwater Horizon (DWH) oil spill, which began on 20 April 2010. Embryo deformities were documented in inviable L. kempii eggs before (2008−2010) and after (2011−2013) the DWH spill in 2 Texas (USA) nesting areas (Upper Texas Coast and Padre Island National Seashore). Additional nesting trends, including clutch size and hatching success, were also investigated. Total and late-stage embryo deformity prevalence were 1.5 times greater after 2010 than before, but low in all nesting seasons (mean ± SD: 0.7 ± 8.5% total; 0.6 ± 8.0% late-stage) and did not differ between locations. Craniofacial and carapace deformities were the most frequently observed deformity types. Documented nests in both areas declined in 2010 relative to previous years, ending an exponential increase observed beginning in 1995. Clutch size remained consistent before and after the spill. Hatching success averaged 87.0 ± 33.3% in all years, but no effects from DWH were determined. Collectively, these data represent useful benchmarks against which to judge impacts of future crude oil spills and other catastrophic events.


INTRODUCTION
On 20 April 2010, the Deepwater Horizon/Mississippi Canyon 252 (DWH) explosion and subsequent 87 d oil leak released 3.19 × 10 6 barrels (5.07 × 10 8 l) of oil into the northern Gulf of Mexico (GOM), resulting in the largest marine oil spill in history (DWH NRDA Trustees 2016). At its maximum extent, oil covered 112 115 km 2 of the GOM and ultimately contaminated the coasts of Louisiana, Mississippi, Alabama, and northern Florida (USA), including adjacent marine waters and wetlands (https:// response. restoration. noaa.gov/erma). Crude oil released into the marine environment has grave acute and chronic impacts on many taxa (Peterson et al. 2003), affecting growth rates, reproductive output, and overall fitness of individuals, as well as potentially impacting population dynamics. DWH oil directly and indirectly killed large numbers of organisms and contaminated important foraging and developmental habitats used by many species (Antonio et al. 2011, Henkel et al. 2012, White et al. 2012, Rozas et al. 2014, McDonald et al. 2017. Large amounts of fossil carbon from oil and gas released by the spill have been incorporated into the GOM ecosystem (Mitra et al. 2012, Cherrier et al. 2014, Quintana-Rizzo et al. 2015, Bonisoli-Alquati et al. 2016, Wilson et al. 2016, Shortly after DWH exposure, marine turtle mortality rates increased exponentially (Antonio et al. 2011) and marine turtle nest densities decreased significantly on beaches near the spill area (e.g. northern GOM, loggerhead turtles Caretta caretta; Lauritsen et al. 2017).
Turtles that survived the initial and acute impacts of the spill may have sustained permanent or temporary physiological injury, experienced difficulty finding food, ingested or inhaled oil, consumed oil-contaminated prey, or moved to other foraging areas (e.g. Milton et al. 2010, Moreno et al. 2013, Beyer et al. 2016). Furthermore, oil exposure in reproductive female turtles can lead to maternal transfer of the lipophilic components of oil and their associated metabolites to eggs, causing embryo deformities and death (Fritts & McGehee 1982, Van Meter et al. 2006, Shigenaka 2010, Field et al. 2017, Zychowski & Godard-Codding 2017. While L. kempii rarely experience complete developmental failure of clutches (Backof et al. 2019), oil and related components can have other deleterious effects on turtle nests. In C. caretta, crude oil experimentally added to sand above incubating eggs led to significant embryonic mortality and morphological changes (Fritts & McGehee 1982), while crude oil and oil-dispersant (Corexit) exposure led to clinicopathological abnormalities and failure to gain weight in hatchlings (Harms et al. 2019). Similarly, snapping turtle Chelydra serpentina eggs exhibited increased embryonic mortality and deformities following exposure to crude oil or oil components in vitro (Van Meter et al. 2006) and in situ, as did painted turtle eggs Chrysemys picta (Bell et al. 2006). This study investigated potential spill effects on embryo deformities, clutch size, and nest success in L. kempii based on statistical comparisons and quantification of data from Texas nesting seasons before (2008−2010) and after (2011−2013) the DWH spill. Research findings on nesting L. kempii in Texas are important to the entire nesting population be cause all individuals constitute 1 genetic stock and a single regional management unit (Wallace et al. 2010). This study was part of the US Department of the Interior Natural Resource Damage Assessment of DWH effects on marine turtles, and serves as a basis for future investigations to prepare for, and respond to, damages caused by future oil spills. Additional information on other aspects of the investigation are described elsewhere (Hooper & Schmitt 2015, Shaver et al. 2016a, Reich et al. 2017).

Study area
The study area consisted of Kemp's ridley Lepidochelys kempii nesting beaches along 2 portions of the Texas coast where long-term monitoring has been conducted: the 105 km length of Padre Island National Seashore (PAIS) and a 100 km section on the Upper Texas Coast (UTC) extending from the southern terminus of Surfside Beach at Brazos Inlet northeastward the length of Bolivar Peninsula to Rollover Pass (Fig. 1). These areas accounted for 56.5% of US L. kempii nests during the study years, 2008−2013 . Systematic patrols were conducted each year to identify, document, and protect L. kempii nests. Nests discovered at PAIS and on the UTC were excavated and packed into polystyrene boxes containing sand from the nest sites, then transported to the PAIS laboratory facility for protected incubation, following protection protocols for most nests found on the Texas coast (Shaver & Caillouet 2015). Incubated nests were categorized as those found at PAIS (PAIS-I) and those from the UTC for analyses. Eggs from some nests found at the southern end of PAIS were incubated in fenced beach corrals (PAIS-C). All hatchlings were released to the GOM at PAIS.

Data collection
Nesting females observed by biologists were examined for flipper tags and tagged with a passive integrated transponder (PIT) and 2 metal tags, if not already present. A biopsy sample was also collected from the trailing edge of one of the rear flippers, if possible, for genetic analyses. Each turtle was measured for maximum straight, maximum curved (cm; nuchal notch to posterior tip; SCL n−t , CCL n−t ), and/or minimum straight carapace length (cm; nuchal notch to posterior notch; SCL n−n ) at the first observation of the turtle in the study area each season. Conversion equations to determine SCL n−t from CCL n−t and SCL n−n for individuals missing SCL n−t were developed using linear regression: CCL n−t to SCL n−t : SCL n−t = 6.660 + 0.848 CCL n−t ; r 2 = 0.88, p < 0.01 (1) SCL n−n to SCL n−t : SCL n−t = 3.823 + 0.950 SCL n−n ; r 2 = 0.94, p < 0.01 (2) SCL n−t was recorded once per year for each female, because there is negligible within-season carapace growth, and measurements were applied to all clutches from each individual within a year. For clutches where a female was not observed nesting, genetic analyses were used to assign unknown clutches to known females (Frey et al. 2014), when possible, and SCL n−t measurements were assigned if they had been taken from the known individual that year. At termination of incubation, hatchlings were removed from the sand surface at the first signs of activity, and nest remains, including trapped live hatchlings, were exhumed. All hatchlings were counted, and live hatchlings were released at PAIS. To expedite release of large numbers of hatchlings, limit light exposure, and to follow pre-spill techniques for accurate comparison, live hatchlings were not examined for deformities. Unhatched eggs, including those that were pipped, were salvaged and refrigerated until examination, which occurred within 6 h. Embryos were examined for the presence of grossly visible abnormalities and assigned a developmental stage at mortality (Crastz 1982). Deformities were classified as one of the following: grossly visible craniofacial, carapace, plastron, limb, pigment, and other deformities (e.g. neck and tail, edema, multiple embryos). Tissue samples were collected from 5 embryos from nests where the female was unknown for genetic matching to identified females using the kinship analysis methods de scribed in Frey et al. (2014). The total number of embryos with each type of deformity, along with the total number of deformed embryos and the total number of late-stage (≥ stage 28; Crastz 1982) deformed embryos, were summed for each clutch for analyses and re ported as percentages.

Statistical analyses
We considered 2010 'pre-spill' in all comparisons because L. kempii nesting in Texas began on 24 April (4 d after the spill began), and it was assumed that most turtles departed foraging areas before oil spread to these regions (Shaver et al. 2016a, Reich et al. 2017. All analyses were conducted in R v.3.6.3 (R Core Team 2020). All comparative analyses were used to determine any differences pre-and post-DWH (2008−2010 vs. 2011−2013) and between incubation types/locations (PAIS-C, PAIS-I, UTC). Comparisons between specific years were not modeled because year and pre-/post-DWH variables were aliased. Preliminary data exploration was carried out following recommended protocols (Zuur et al. 2010).
A total of 34 irregular clutches (5.2%) ranging from 2 (2010) to 9 (2009) per year were excluded from analyses to eliminate confounding factors. Clutches were excluded if they had been disturbed by predators, inundated by seawater, incubated in situ (as opposed to corral or laboratory incubation), were surgically removed from a recently deceased female, or were produced in captivity by oxytocin administra-tion. For the 6 yr period, 650 clutches were included in analyses. Two additional clutches were removed from analysis of deformity rates because they had 100% developmental failure and thus no embryos were examined.
Deformity rates (total and late-stage) for 519 clutches, where the identity of the female was known through direct observation or genetic assignment, were evaluated using generalized linear mixed models (GLMMs). Deformity rates were modeled as the proportion of 'successes' (deformed embryos/ hatchlings) and the proportion of 'failures' (normal embryos/ hatchlings). Differences between specific deformity types were not analyzed, as some deformities occurred at such low frequencies that analyses could not be completed. Preliminary statistical testing indicated over-dispersion and zero-inflation of the data. Due to these factors, zero-inflated binomial GLMMs with a logit link function were used to analyze deformity rates, which incorporated a zero-model to account for excess zeroes and a binomial-model to determine the prevalence of the response. The proportions of deformities were modeled as functions of the fixed main covariates incubation type/location (categorical with 3 levels) and pre-/post-DWH (categorical with 2 levels) with nest day (continuous) as a nuisance variable. Maximum SCL n−t was not in cluded as a covariate as it was not correlated with any deformity types. Turtle identity was included as a random effect to account for non-independence of clutches from known females (i.e. multiple clutches as signed to a specific female within and between years).
Using linear regression, nesting data indicated the expected positive relationship between clutch size and SCL n−t (r 2 = 0.13, F 1, 464 = 70.65, p < 0.001) ( Fig. 2) (Witzell et al. 2005). A weaker negative relationship between clutch size and nest date (day of the nesting season, beginning with detection of the first nest in each study area each year) was also evident, with clutch size decreasing slightly as the nesting season progressed (r 2 = 0.02, F 1, 647 = 13.15, p < 0.001) (Fig. 2). There was no significant interaction between SCL n−t and nest date. Consequently, comparisons of clutch size between study locations and pre-/post-DWH required accounting for both nest date and SCL n−t . Clutch size was evaluated using a linear mixed model (LMM) and included data from 466 clutches from known females with associated SCL n−t measurements. Nest data without corresponding SCL n−t measurements (n = 184) were excluded from clutch size analyses, as SCL n−t was found to have a positive relationship with clutch size. Clutch size was modeled as a function of the fixed covariates pre-/post-DWH, incubation type/location, SCL n−t (continuous), and nest date (continuous) with turtle identity included as a random effect. Initial analyses indicated hatching success to be highly variable within and between years, time periods, and locations. Due to these factors, it was not possible to successfully model hatching success to fit the data to determine any effects post-DWH or between incubation types/locations. Therefore, hatching success is presented as observational data only.
LMMs were conducted using the 'lme4' package (Bates et al. 2015) and evaluated using standard methods (Zuur et al. 2010). GLMMs were conducted using the 'glmmTMB' package (Brooks et al. 2017) and were evaluated using the 'DHARMa' package (Hartig 2020). Models were fit to the data by testing the significance of each variable at each step and extracting non-significant terms (backwards stepwise) from the maximal models, using Akaike's information criterion for comparisons, until the most parsimonious (i.e. minimally adequate) model was obtained for each response variable (Crawley 2015). A significance level of α = 0.05 was used to judge the significance of all statistical tests. Interaction terms of fixed and nuisance covariates were included in maximal models (Crawley 2015), where appropriate (Zuur et al. 2010). Nakagawa's r 2 was used to determine the amount of variance explained by only the fixed (marginal r 2 ) and both the fixed and random (conditional r 2 ) effects of each final mixed model using the 'performance' package (Nakagawa & Schielzeth 2013, Nakagawa et al. 2017, Lüdecke et al. 2020. Only the main effect means (incubation type/location and pre-/post-DWH) are presented, with the focus on before/after comparisons. All data are reported as means ± SD. Data were analyzed after exclusion of irregular clutches. Females were directly examined for size and tags at 413 of the 650 nests included in the analyses (63.5%) during 2008−2013 (Table 1). UTC females were directly examined approximately half as often as PAIS females. An additional 108 clutches were assigned to known females using genetic analyses (Frey et al. 2014). Over the 6 yr study, using tag numbers and genetic analyses, 199 females were documented nesting only once during a nesting season, 122 were documented twice, and 26 were documented 3 times for a total of 521 clutches from 240 unique individuals (Table 1) assigned to a female. These unassigned nests were laid by unobserved nesting females who may have been first time nesters, turtles that nested previously in the study area but evaded detection, or turtles that have nested in the study area for the first time but nested elsewhere previously, thus identity could not be assigned. These clutches were not included in the mixed models. Unique females documented nesting in each year ranged from 47 (2010)  The total number of measured individuals is greater than 240, as known individuals returned to nest in more than one season, with new measurements obtained each year they were seen.

Embryo deformities
Two clutches from UTC, laid by the same female in 2010, had 100% developmental failure and were not included in deformity analyses. After the exclusion of these 2 clutches, 648 clutches were analyzed for general summary statistics and a total of 459 deformed dead embryos were observed in 241 clutches (37.2%) at an overall rate of 0.7 ± 8.5% clutch −1 ( Table 2). The majority of deformities were late-stage (developmental stage 28 or later; overall 0.6 ± 8.0% nest −1 ). Nine of the clutches with the greatest percentages of deformed embryos (> 5%, n = 11) were from 2011 or 2012, while 2 were from 2009. All clutches with elevated deformity rates, except one, were laid at PAIS.
From the 519 clutches laid by known females, excluding the 2 failed clutches, there were no significant interactions between variables, and incubation/ location type and nest date were also non-significant and were removed from the models to create the minimally adequate GLMMs which included pre-/ post DWH as the main effect, and turtle identity as a random effect. Final models indicated a 1.5 times higher rate of total number of deformities and of total late-stage deformities in L. kempii clutches post-DWH (total deformities: pre-DWH = 0.8 ± 0.1%, post-DWH = 1.15 ± 0.2%, p < 0.01; late-stage deformities: pre-DWH = 0.6 ± 0.1%, post-DWH = 0.9 ± 0.2%, p < 0.01) (Fig. 3, Table 3). Inclusion of the random effect increased the explanatory power of each GLMM (total deformities: marginal r 2 = 0.01, conditional r 2 = 0.14; late-stage deformities: marginal r 2 = 0.01, conditional r 2 = 0.18), but model explanation of the variability was low overall.

Clutch size and hatching success
The 650 clutches included in general summary statistical analyses contained 63 138 eggs during 2008− 2013, of which 55 140 hatched (87.3%) and from which 54 981 hatchlings (87.1%) emerged and were released (Table 4). Over the 6 yr study, 99.7 ± 5.4% of the turtles that hatched survived and were released. For the 466 clutches from known females with associated SCL n−t measurements, clutch size did not significantly differ pre-and post-DWH (p = 0.28) but did differ between locations, with PAIS-I and PAIS-C having lower average clutch sizes (PAIS-I: p < 0.01;    Table 3). Inclusion of the random effect increased the overall explanatory power of the LMM (marginal r 2 = 0.20, conditional r 2 = 0.41). Overall hatching success was 87.3 ± 33.3% for the 6 yr period, 87.8 ± 32.7% pre-DWH, and 86.9 ± 33.8% post-DWH. Hatching success was highly variable among nests from different locations within the study area, with higher rates of success for PAIS-I (88.4 ± 32.0%) and PAIS-C (91.1 ± 28.5%) than UTC (74.5 ± 43.6%).

DISCUSSION
International conservation efforts implemented since the mid-1970s resulted in an exponential increase of Lepidochelys kempii nests documented each year from 1999−2009 (Caillouet 2011, Gallaway et al. 2016a,b, Shaver et al. 2016b. However, nest numbers in Mexico and Texas, which are highly correlated and represent portions of the same nesting population ), decreased 37% in 2010 relative to 2009, as did the number of hatchlings released (Crowder & Heppell 2011), representing a major departure from trends in previous years. Nest numbers remained lower than predicted after the DWH spill , Gallaway et al. 2016b) and did not begin to increase again until 2016. Nesting has not resumed the exponential increase observed prior to DWH, causing further concern for the sta-   Table 3. Results of generalized and general linear mixed effects models (GLMM/LMM) of deformity rates and clutch size in Lepidochelys kempii clutches before (2008−2010) and after (2011−2013) the Deepwater Horizon (DWH) oil spill. Deformity rates were modeled using a zero-inflated binomial GLMM with logit link; clutch size was modeled using an LMM. PAIS-I: Padre Island National Seashore incubation facility; PAIS-C: Padre Island National Seashore corral; SCL n−t : straight carapace length, nuchal notch to posterior tip (cm) tus of this species (Gallaway et al. 2016b, Caillouet et al. 2018). Most L. kempii nesting occurs on the Gulf coast of Mexico in Tamaulipas and Veracruz and to a lesser extent in Texas . Nearshore waters of the northern Gulf, the same geographic area affected by the DWH oil spill, provide particularly important foraging and migratory habitat for L. kempii that migrate from nesting beaches in Texas and Mexico (Shaver & Rubio 2008, Shaver et al. 2013, Putman et al. 2015, Gredzens & Shaver 2020. Impacts from the DWH oil spill, if present, may have directly affected nesting females (oiling or disturbance) or indirectly affected them via ingestion of oil or related compounds from the water column or dietary sources (Zychowski & Godard-Codding 2017). Anthropogenic activities post-DWH, such as increased boat traffic during clean-up, may also have had detrimental effects on marine turtle populations (Lauritsen et al. 2017). It is unknown if oil contamination of the foraging grounds affected reproductive success and, ultimately, the observed slowed population recovery.
Deformity rates were higher in examined L. kempii embryos post-DWH. Statistical models explained a comparatively small proportion of the total variation in the data, thus indicating that factors not incorporated into the models are likely important. Possible confounding variables could be error introduced by embryo deformity examiners and/or a possible increased scrutiny of embryo examination after the oil spill due to awareness of the study. Future studies should work to mitigate these types of confounding factors. For this study, live hatchlings were not examined for deformities in order to protect them from exposure to white light and allow for quick release. Analysis of the data based only on females observed and examined while nesting or genetically matched to clutches indicates that differences between or among nests made by the same turtle account for substantially more of the total variation for some of the variables analyzed. Deformed embryo prevalence was low overall and highly variable (0.7 ± 8.5%), indicating probable low population-level impacts. In comparison, turtle embryo deformities from John Heinz National Wildlife Refuge, which is highly contaminated by multiple sources, ranged from 13−71% annually in 2 freshwater species, Chelydra serpentina and Chrysemys picta (Bell et al. 2006). Frequency of malformations for marine turtle eggs has also been reported in other species of marine turtles, for both in situ and relocated nests, ranging from 0.2−26.0%, with high variability between species and locations (Bárcenas-Ibarra et al. 2015). Embryo malformation rates for the most closely related species of marine turtle, the olive ridley Lepidochelys olivacea, are relatively high (1−2%) when compared to other species (Bárcenas-Ibarra et al. 2015).
Annual clutch size and frequency were monitored for changes after the spill, as fecundity is critical to , clutch size was not significantly different, but did differ between study locations after accounting for turtle size and nest date, with larger clutches found on the UTC. Differences in clutch size between locations are likely the result of more nesting at PAIS paired with increased patrol effort there late in the nesting season (August−September) . Analysis indicated that nesting day may affect clutch size, with smaller clutches later in the season (Fig. 2b). Other factors may also contribute, as the explanatory power of the model was low. In addition, there may be other unexplained factors as to why UTC clutch sizes are larger, such as mean individual age or health or the closer proximity of the UTC to important foraging areas (Shaver et al. 2013(Shaver et al. , 2016a. Over 36 yr (1978−2014), annual mean clutch size did not vary greatly for L. kempii nesting in Texas . Considered together, the number of nests detected per known female was the same in 2011−2013 (1.5 ± 0.6) and 2008−2010 (1.5 ± 0.6). However, identified turtles undoubtedly made some of the nests laid by unidentified females, indicating that the number of nests per female is lower than the estimated total at 1.8 nests yr −1 (Frey et al. 2014).
Hatching success for L. kempii nests in the study areas was monitored for changes after DWH. However, no significant differences in hatching success between 2008−2010 and 2011−2013 were found due to high variability in the data. Hatch success is very complex and is impacted by many factors (e.g. nest temperature, moisture, sand quality, genetics, health  . For these reasons, this parameter may be of limited benefit in future studies examining shortterm impacts of oil spills, unless the impacts are profound. Long-term impacts of persistent environmental exposure (Peterson et al. 2003), ontogenetic ef fects (Bell et al. 2006), or latency factors may delay detectability for many years. The timing of the spill likely mitigated impacts to L. kempii during 2010, as females nesting during that season migrated away from the oiled habitat to nesting grounds prior to the spill (Shaver et al. 2016a). Data from satellite-tracked adult nesting females and stable isotope samples indicated that 51.5% of the L. kempii in the northern GOM may have been exposed to DWH oil in 2010(Reich et al. 2017, consistent with other reports of marine turtle interactions with oil in the northern GOM (McDonald et al. 2017). These assumptions are further supported by the fact that chemical analyses of eggs and embryos salvaged from 2010−2012 nests contained no detectable, or extremely low, concentrations of oil metabolites (Hooper & Schmitt 2015). These observations do not preclude the possible involvement of the spill in the lower than expected nest numbers in Texas and Mexico during 2011−2015 (Crowder & Heppell 2011) (i.e. in addition to direct mortalities of adult turtles, oilaffected turtles that survived the spill may not have migrated and/or nested).
Observations of nesting turtles and nest success represent only the beginning and the end of the marine turtle reproductive cycle. Tracking data from 2010−2013 indicate that many nesting female L. kempii traveled in the oiled area after nesting in Texas (Shaver et al. 2016a, Gredzens & Shaver 2020, and immature turtles exposed to oil from the DWH spill and the progeny of exposed adults will not mature for many years. Impacts observed in this study remain low but could be drastically different if a similar environmental disaster occurred closer to the main nesting areas (Lauritsen et al. 2017), where significant numbers of oil and gas platforms overlap with critical habitat ). There is significant potential for future catastrophic oil releases within critical foraging habitat and L. kempii nesting beaches . Continued evaluation of data from nesting beaches (i.e. nest numbers, nest success, remigration intervals, growth rates, and re -cruitment) is necessary to provide evidence of whether oil contamination at the foraging grounds affects reproductive success and, ultimately, population recovery. Data and methods reported here represent useful benchmarks against which to judge impacts to marine turtles by future oil spills and can serve as a basis for future investigations to prepare for, and respond to, damages caused.