Catastrophic events, like hurricanes, bring lethal conditions that can have population-altering effects. The threatened Caribbean dry forest occurs in a region known for its high-intensity hurricane seasons and high species endemism, highlighting the need to better understand hurricane impacts in combination with habitat fragmentation and loss. However, such studies remain rare, and for reptiles are mostly restricted to Anolis. Here we used single-season occupancy modeling to infer the impact of the intense 2017 Atlantic hurricane season on the critically endangered Lesser Antillean Iguana, Iguana delicatissima. We surveyed 30 transects across eight habitats on St. Eustatius during 2017–2019, which resulted in 344 individual surveys and 98 iguana observations. Analyses of abundance and site occupancy indicated both measures for 2018 and 2019 were strongly reduced compared to the pre-hurricane 2017 state. Iguanas at higher elevations were affected more profoundly, likely due to higher wind speeds, tree damage and extensive defoliation. Overall, our results indicate a decrease in population estimates (23.3-26.5%) and abundance (22-23.8%) for 2018 and 2019, and a 75% reduction in opportunistic sightings of tagged iguanas between 2017–2018. As only small and isolated I. delicatissima populations remain, our study further demonstrates their vulnerability to stochastic events. Considering the frequency and intensity of hurricanes are projected to increase, our results stress the urgent need for population-increasing conservation actions in order to secure the long-term survival of I. delicatissima throughout its range. Given the projected increase and poleward shift of hurricanes, our study provides important insights from a non-model species.

To assess the impact of the 2017 hurricane season on the local iguana population, we compared data from 2018–2019 to the reference 2017 data collected before to the passing of both hurricanes. As our dataset contained a large number of non-detections and small counts of iguanas, we chose an N-mixture model approach over distance sampling (Royle, 2004). Changes in occupancy, defined as the probability of a site being occupied (MacKenzie et al., 2002), provide a robust proxy for population declines, particularly for sparsely populated and cryptic species (Beaudrot et al., 2016; Dénes et al., 2015). Models were based on four assumptions: (i) the iguana population in each site was closed to colonization and extinction during the sampling period; (ii) iguanas were either detected or not detected during surveys; (iii) iguanas were not falsely detected during surveys; and (iv) iguana detections at each site were independent (Murray and Sandercock, 2020). We followed the methods described by Madden et al. (2021). Briefly, we used likelihood-based single-season occupancy models (MacKenzie et al., 2017), using the R packages “wiqid” and “unmarked” (Fiske and Chandler, 2011; MacKenzie et al., 2002), to estimate site occupancy (ψ) in relation to habitat and elevation, while accounting for detectability (p). We modeled detection probability and occupancy probability to estimate iguana abundance (λ), and also tested the influence of independent variables on abundance estimates, including quadratic terms. Resulting models were ranked using Akaike’s Information Criterion, where the detection model with the lowest AICc is best-fitting (AICc; Burnham and Anderson, 2002), while correcting for small sample size. When there was not a clear single model (deltaAIC <2), we averaged across all models weighted by AICc to produce a model-averaged prediction (Burnham and Anderson, 2002). We checked the final model for multicollinearity using the package “VIF” (Lin, 2012), and inspected the residual plot. Finally, we tested for overdispersion by assessing the goodness-of-fit of the most parsimonious model in each year (MacKenzie and Bailey, 2004) with the “mb.gof.test” function in the package “AICcmodavg” using 1,000 simulations, which calculates a Pearson’s chi-square fit statistic from the observed and expected frequencies of detection histories for a given model. All analyses were performed in the R environment v3.5.1 (R Core Team, 2019), and our data are available at DataDryad (van den Burg et al., 2022).

For use with R packages wiqid and unmarked.