Introduction
Occupancy models have a long history of use in population status assessment and guiding management decisions for threatened species (MacKenzie et al., Reference MacKenzie, Nichols, Royle, Pollock, Bailey and Hines2017). Data are collected using either direct or indirect sampling of a species and used to account for imperfect detection and the proportion of sites where the species is actually present. Compared to standard or static occupancy models, dynamic occupancy models can be used to consider changes that may occur over time and space, which can be one of the main sources of error in monitoring studies (Yoccoz et al., Reference Yoccoz, Nichols and Boulinier2001; MacKenzie et al., Reference MacKenzie, Nichols, Lachman, Droege, Royle and Langtimm2002; Royle & Nichols, Reference Royle and Nichols2003). Failure to account for imperfect detection can result in sites appearing to be unoccupied that are actually occupied (i.e. false absences), which is common for cryptic species that may go undetected (Yoccoz et al., Reference Yoccoz, Nichols and Boulinier2001). Although occupancy alone is a useful state variable, extinction and colonization probabilities estimated over time can also be modelled in relation to site characteristics. These permit testing of hypotheses concerning the drivers of changes in occupancy and allow stronger inferences concerning the mechanisms underpinning occupancy dynamics.
The Bermuda skink Plestiodon (=Eumeces) longirostris is a cryptic species, categorized as Critically Endangered on the IUCN Red List (Cox & Wingate, Reference Cox and Wingate2021) and the only endemic extant terrestrial vertebrate on Bermuda (Davenport et al., Reference Davenport, Hills, Glasspool and Ward2001). Once considered abundant across Bermuda, the population continues to decline because of habitat loss, anthropogenic activities and invasive species. Declines were observed in the early 1900s when it was reported that skinks were rarely seen on the mainland but prevalent on offshore islands, particularly around coastal cliffs (Verrill, Reference Verrill1902). However, concern for the survival of the species only began nearly a century later, when a conservation campaign was undertaken through the Bermuda Zoological Society. Monitoring commenced shortly thereafter, with a focus on the sites with recent sightings, to elucidate population sizes and demographic characteristics (Raine, Reference Raine1998; Wingate, Reference Wingate1998; Davenport et al., Reference Davenport, Hills, Glasspool and Ward2001; Glasspool & Outerbridge, Reference Glasspool and Outerbridge2004). The island-wide skink population was estimated to be 2,300–3,500 individuals, with a suggestion that the true population size could be 5,000 or more given the likelihood of further sites being discovered (Edgar et al., Reference Edgar, Kitson, Glasspool and Sarkis2010). Determining the current distribution and population status is necessary to inform conservation management.
In this study we used 3 years of presence–absence surveys across Bermuda to model occupancy, colonization, local extinction and detectability. Inclusion of potential drivers of these parameters as covariates facilitated the identification of key areas for conservation of the species and appropriate management actions.
Study area
Forty locations were surveyed for the presence of skinks across Bermuda (Fig. 1) during 2015–2017. Surveys were undertaken during April–July when skinks are most active (Edgar et al., Reference Edgar, Kitson, Glasspool and Sarkis2010). Sites were selected based on historical records or were considered to potentially contain skinks based on existing knowledge of habitat preferences, with emphasis on sites that had not been surveyed for > 10 years (Glasspool & Outerbridge, Reference Glasspool and Outerbridge2004). Locations were isolated islets or islands, or situated at least 100 m apart on the mainland to ensure independence among sites.
Fig. 1 The 40 localities surveyed for the Bermuda skink Plestiodon longirostris during 2015–2017, indicating where it was and was not detected. There were three survey sites at Spittal Pond (Jefferies Hole, Checkboard and East End) and Inner Pear Rocks (Inner, Middle and Outer).
Methods
Surveys
The methodology followed that of Davenport et al. (Reference Davenport, Hills, Glasspool and Ward1997), with skinks captured using pitfall traps. As the Bermuda skink has a home range of c. 10 m2 (Davenport et al., Reference Davenport, Hills, Glasspool and Ward1997), traps were placed 5–20 m apart, with 10–72 traps at each site and a greater number of traps at larger sites to increase the chances of capturing skinks. Average trap density was 31 traps per site (0.008 traps/m2; Supplementary Table 1). The traps were opened during 11.00–16.00 and checked hourly to ensure any trapped skinks did not overheat. Skinks were measured, weighed and identified (e.g. Turner et al., Reference Turner, Griffiths, Outerbridge and Garcia2019) and then released at the point of capture. Each site was surveyed 2–15 times over the 3-year period (depending on weather and site access). To increase chances of detectability, surveys were not conducted during rain or high winds (> 40 km/h).
Statistical analyses
All statistical analyses were undertaken using R 3.3.2 (R Core Team, Reference R Core Team2016) with the function colext in the package unmarked (Fiske & Chandler, Reference Fiske and Chandler2011). Detection and non-detection data and a multi-season occupancy model were used to assess occupancy, colonization, extinction, detectability and probable distribution of P. longirostris across Bermuda. The parameters were ψ = probability of a site being occupied, p = probability of being detected given presence, γ = probability of colonization, and ɛ = probability of local extinction. Within each year, the four parameters were assumed to be constant, but changes in occupancy were modelled between years and could be a function of covariates (i.e. site-specific habitat features).
Model covariates
Occupancy probabilities may depend on covariates. Five site-specific covariates identified from previous studies (Davenport et al., Reference Davenport, Hills, Glasspool and Ward1997; Raine, Reference Raine1998; Wingate, Reference Wingate1998; Glasspool & Outerbridge, Reference Glasspool and Outerbridge2004) and discussions with experienced ecologists were included in the models (Table 1). These were constant for all visits to a site and thus permit testing hypotheses about drivers of occupancy, colonization and local extinction. Site-specific covariates included the site type (i.e. mainland or island), because anthropogenic disturbance on the mainland (coastal development, coastal and beach activities, and litter (i.e. discarded cans and bottles, which than be lethal) are known to threaten the skinks (Davenport et al., Reference Davenport, Hills, Glasspool and Ward1997, Reference Davenport, Hills, Glasspool and Ward2001; Raine, Reference Raine1998; Wingate, Reference Wingate1998). Habitat type was recorded as the most prevalent habitat type at each site (dense forest or coastal rock and scrub). Nesting seabird colonies and seasonal fruits provide the skinks with food (Davenport et al., Reference Davenport, Hills, Glasspool and Ward2001; Madeiros, Reference Madeiros2005), and therefore the presence of nesting seabirds (white-tailed tropicbird Phaethon lepturus catsbyii, Bermuda petrel Pterodroma cahow) and the prickly pear cactus Opuntia dillenii were included as site-specific covariates. As the number of traps varied between sites, this was included as a covariate.
Table 1 The five site-specific and five observation-level predictor covariates used in models of site use by the Bermuda skink Plestiodon longirostris. All covariates were dichotomous, except for number of traps, which was continuous.
1 White-tailed tropicbird Phaethon lepturus catsbyii, Bermuda petrel Pterodroma cahow.
Five introduced species (kiskadee flycatcher Pitangus sulphuratus, domestic cat Felis catus, yellow-crowned night heron Nyctanassa violacea, Jamaican anole Anolis grahami and black and brown rats Rattus spp.) were used as observation-level covariates as they are associated with known predatory threats. Although additional predators and competitors are present (Anolis leachi, Anolis sagrei, Anolis extremus, domestic chickens Gallus domesticus, and American crows Corvus brachyrhynchos; Wingate, Reference Wingate2011; Stroud et al., Reference Stroud, Giery and Outerbridge2017), these would have required more exhaustive surveys and therefore we did not collect data on their distribution sufficient for them to be included in the modelling. We visually confirmed the presence or absence of each of the five observation-level covariates during each trapping occasion.
Model selection and averaging
In our first analysis (the null model), we assumed that all four parameters (ψ, ɣ, ɛ, p) were constant across sites and surveys. This was denoted by (.) and no covariates were included. We used this model to provide a comparison with the unadjusted (i.e. naïve occupancy) proportion of sites where at least one skink was detected over the 3-year survey period.
In a second analysis, we included covariates (MacKenzie et al., Reference MacKenzie, Nichols, Lachman, Droege, Royle and Langtimm2002, Reference MacKenzie, Nichols, Royle, Pollock, Bailey and Hines2005) for each of the parameters hypothesized to affect occupancy, colonization, extinction or detection probabilities, using a maximum likelihood approach (MacKenzie et al., Reference MacKenzie, Nichols, Royle, Pollock, Bailey and Hines2005). All covariates were standardized (z-score) so that scores could be compared between different types of variable and that the logit scale coefficients 
$( \hat{\beta }) $ were not skewed by unevenly large ranges in the data (MacKenzie et al., Reference MacKenzie, Nichols, Royle, Pollock, Bailey and Hines2005). Where logit scale coefficient 
$( \hat{\beta }) $ estimates and associated standard errors (SE
$( \hat{\beta }) $) are reported, a positive number suggests a positive relationship between the covariate and the model parameter it is a predictor of, and vice versa for negative numbers.
Logit scale coefficient 
$( \hat{\beta }) $ estimates were back-transformed using the plogis function in R to give the model parameter estimates. Then, using the package AICcmodavg in R, we assessed the fit of our models using a goodness-of-fit test (MacKenzie–Bailey test) based on bootstrapping (10,000 iterations) and Pearson's χ 2. The level of significance was set at P = 0.05, with larger values indicating a poorer fit (MacKenzie & Bailey, Reference MacKenzie and Bailey2004; Wright et al., Reference Wright, Irvine and Rodhouse2016). To identify the most parsimonious and biologically plausible models for the observed data and assess which combination of covariates best explained the detection histories observed, we calculated and ranked Akaike's information criterion (AIC) values for model comparison (Akaike, Reference Akaike, Petrov and Csaki1973; Burnham & Anderson, Reference Burnham and Anderson2002). A set of candidate models was generated by selecting those that had a summed AIC weight of at least 0.95, indicating there was 95% confidence these models best explained the data. If there were multiple top ranked models, a weighted model averaging technique was applied (Burnham & Anderson, Reference Burnham and Anderson2002) to estimate the detection, occupancy, colonization and local extinction probabilities (with standard errors).
To optimize the survey design, the probability of detecting Bermuda skinks at least once if visiting a site K times was calculated from our best model to determine if the species was truly absent from a site, using the following expression (Pellet & Schmidt, Reference Pellet and Schmidt2005; Barata et al., Reference Barata, Griffiths and Ridout2017):
$$ K = \displaystyle{{{\rm log}\left( {1 -p^*} \right)} \over {{\rm log}\left( {1- p} \right)}} $$where p is the detection probability and p* is the desired probability of detecting the species at an occupied site on at least one of the K visits, set to 0.80, 0.90, 0.95 and 0.99.
Results
On average, we visited sites five times (range 2–15 visits) over the 3 years. Of the 40 sites surveyed, we confirmed the occurrence of the skink at 13 sites (Supplementary Table 1). The naïve occupancies (calculated as the proportion of sites where presence was confirmed) were 0.25 (8/32) in 2015, 0.48 (12/25) in 2016 and 0.11 (4/38) in 2017. The estimated probability of occupancy using the null model (AIC = 233.91) was 0.22 ± SE 0.08 across the 3 years.
The top three models were selected as they had a cumulative weight of 0.96, indicating these models best explained the data (Table 2). As the goodness-of-fit test had P < 0.05 we accept the hypothesis that the models adequately fit the data (model 1: χ 2 = 205.4, P = 0.02; model 2: χ 2 = 187.1, P = 0.03; model 3: χ 2 = 140.5, P = 0.01). The weighted model averaged estimates of occupancy (ψ), colonization (ɣ), local extinction (ɛ) and detection (p) probabilities were taken as the final estimates. The estimate of occupancy (0.25 ± 0.06) did not vary substantially from the naïve estimate because the detection probability was relatively high (0.45 ± SE 0.06), and therefore estimates are reasonably unbiased (MacKenzie et al., Reference MacKenzie, Nichols, Lachman, Droege, Royle and Langtimm2002; Table 3).
Table 2 Ranking of the best fitting models of occupancy (ψ), detection (p), colonization (γ) and extinction (ɛ) and their covariates, with AIC values, ΔAIC, AIC model weights (wi), cumulative AIC model weights (cw) and the number of parameters (K) in the top three ranked models of the occupancy of the Bermuda skink.
1 PR, presence of rats; PC, presence of cats; HT, habitat type; ST, site type; NT, no. of traps; PP, presence of prickly pear; SB, presence of nesting seabirds.
Table 3 Transformed parameter 
$( \hat{\beta }) $ estimates for the top three models listed in Table 2, with weighted model averages. All associated standard errors (SE) are included. The first, second and third set of coefficients are for model terms associated with the occupancy (ψ), colonization (γ) and local extinction (ɛ) parameters, respectively. The fourth set of coefficients explained heterogeneity in detection probabilities (p) associated with different surveys.
1 PR, presence of rats; PC, presence of cats; HT, habitat type; ST, site type; NT, no. of traps; PP, presence of prickly pear; SB, presence of nesting seabirds.
Parameter estimates using the top model (AIC value = 168.12) suggest there was a positive relationship between skink presence and coastal habitat type (
$\hat{\beta }$ 2.44 ± SE 0.61), the presence of prickly pear cacti (1.48 ± 0.55) and island sites (1.48 ± 0.55). The presence of cats negatively influenced colonization (−0.32 ± 0.65), whereas the presence of seabirds positively influenced colonization (3.10 ± 0.39). Local extinction was less likely in coastal habitat type (−3.10 ± 0.39) and where prickly pear cacti were present (−3.10 ± 0.39). Detection was positively influenced by the presence of seabirds (7.10 ± 0.21), the coastal habitat type (5.34 ± 0.22) and prickly pear cacti (3.54 ± 0.54), whereas skinks were unlikely to be detected where rats were present (−2.94 ± 0.45) and even less so with cats (−4.30 ± 0.25). As all three of the top models were influenced by the same covariates and in the same directions, this suggests there is considerable support for their inclusion. We found no support for any influence of the presence of kiskadees, yellow-crowned night herons or anoles, or the number of traps used, on ψ, ɣ, ɛ and p.
Minimum number of surveys
To improve the survey design for future monitoring, we calculated the number of surveys required to detect the Bermuda skink at a given site. Assuming skinks are imperfectly detected (detection probability < 1), detection probabilities of 0.80–0.99 were chosen to reflect a high chance of detecting skinks if present. We found that three (mean 2.74) surveys were needed for a 0.80 probability, four (3.92) for a 0.90 probability, five (5.10) for a 0.95 probability and eight (7.85) for a 0.99 probability that the Bermuda skink will be detected at a site.
Discussion
Occupancy probability
Occupancy was related to habitat type, site type and presence of prickly pear cacti. This provides evidence that rocky coastal habitats, particularly on offshore islands, are refuges that support and maintain skink populations (Glasspool & Outerbridge, Reference Glasspool and Outerbridge2004). Areas with greater levels of habitat degradation and loss (especially on the mainland) were less likely to be occupied than islands (skinks were detected on four mainland sites compared to nine island sites). This suggests that skinks require relatively undisturbed native habitats to thrive, and they could therefore act as important biological indicators of the condition of coastal habitats in Bermuda.
The presence of scrub vegetation and prickly pear cacti may provide a seasonal source of fruit, attract invertebrate prey species, provide shelter and refugia for skinks, and play a key role in erosion control (Le Houérou, Reference Le Houérou1996). Raine (Reference Raine1998) also noted an association between skinks and areas dominated by coastal vegetation such as sea ox-eye Borrichia arborescens and salt grass Spartina patens on Inner Pear, Charles Island and Spittal Pond. Therefore, clearing areas of invasive plants such as Brazil pepper Schinus terebinthifolia and asparagus fern Asparagus densiflorus in suitable coastal locations would be beneficial to skink survival. This would simplify migration between population fragments and increase the probability of population survival.
Although surveys conducted during 1998–2014 did not report occupancy or detectability of the Bermuda skink, our study indicates the species’ range is declining, as we found skinks at only 11 of the 26 sites occupied in the previous 20 years, and we detected a skink on more than one sampling occasion at only seven sites. As 67% of detections were in eastern Bermuda (i.e. within the Castle Harbour area), these subpopulations appear to comprise the majority of the population and would benefit from increased habitat management to control invasive vegetation and reduce the harmful effects of litter (i.e. discarded empty bottles and cans), which can be lethal to the skinks (Jones, Reference Jones2015).
Colonization and extinction probabilities
Seabirds were found to be important predictors of colonization. This confirms previous suggestions that skinks have a mutualistic relationship with Bermuda's nesting seabirds because skinks opportunistically forage in the seabirds' nests, feeding on failed eggs, carrion such as dead chicks and uneaten fish (Davenport et al., Reference Davenport, Hills, Glasspool and Ward1997). The installation of artificial nesting burrows has been an important component of the recovery of Bermuda's breeding white-tailed tropicbirds and Bermuda petrels (Madeiros, Reference Madeiros2008), so the continuation of this process may help to sustain the skink and encourage its colonization of suitable locations.
However, the probability of colonization by skinks in the presence of cats is low, as cats have been observed predating skinks on many occasions (Garber, Reference Garber1988). Domestic and feral cats are a global threat to many threatened species (Medina et al., Reference Medina, Bonnaud, Vidal, Tershey, Zavaleta and Donlan2011). Bermuda has a high number of domestic and feral cats relative to its size (McGrath, Reference McGrath2014), and a strategy for all aspects of cat management on Bermuda, including the creation of a legislative and regulatory framework, is needed.
Surprisingly, local extinction was not influenced by any predator covariates, but the absence of prickly pear cacti and rocky coastal habitat were key factors. The rocky remains of historical defence fortifications on some islands seem to provide good skink habitat. Changes in habitat are therefore most likely to be important indicators for predicting local extinction of the skink. Mean extinction probabilities tended to be higher than average colonization probabilities. The variability in average colonization and extinction probabilities suggests that the various subpopulations are going through temporal fluctuations in site occupancy.
Detection probability
We assessed the presence of multiple predator species alongside the detection of skinks at each site to determine which species pose threats to remnant skink populations. Although kiskadees, herons and anoles have previously been documented as predators (Davenport et al., Reference Davenport, Hills, Glasspool and Ward1997), these were not covariates retained in the top models (anoles: AIC > 234.93; herons: AIC > 234.04; kiskadees: AIC > 234.30). Birds are capable of accessing all sites and anole lizards are widespread across Bermuda (Macedonia et al., Reference Macedonia, Clark and Mcintosh2016). However, the recent arrival in 2011 and subsequent establishment of the Cuban brown anole Anolis sagrei has been identified as a potential threat to the Bermuda skink and warrants monitoring (Stroud et al., Reference Stroud, Giery and Outerbridge2017). The presence of rats and cats was found to negatively affect detection of the skinks. Continued monitoring of the prevalence of predators at each study site would determine the threat level they pose to the skink and whether the management of introduced predators or competitors at a site would facilitate recovery of the skink.
Improving parameter estimates
Precise estimates of occupancy require a large number of sites, but for threatened species there may be insufficient sites for rigorous replication. Because of occasional inclement weather and limited access, many of the sites in this study were visited only twice. Consequently, low precision, with high standard errors of estimates, is likely. To improve the precision of the estimated occupancy rate, increasing the number of surveys to at least five should result in a 0.95 detection probability, but there would be little gain in precision by undertaking more than five surveys.
Conclusion
Bermuda skinks exhibit high variation in occupancy and abundance between years (Turner et al., Reference Turner, Griffiths, Outerbridge and Garcia2019). The patchy distribution of this species and the high number of threats are constraints on dispersal and recruitment. Our findings demonstrate that the management of the remaining skink populations depends on (1) continuation of control of non-native predators such as cats and rats, and (2) restoring native coastal habitat on offshore islands. The apparently mutually beneficial relationship between skinks and seabirds and how this affects skink occupancy dynamics merits further research. Likewise, the impact on extinction–colonization dynamics of skink mortality in discarded bottles needs elucidation. New initiatives to explore opportunities for restoring populations on undisturbed offshore islands are needed to ensure the survival of the endemic Bermuda skink.
Acknowledgements
We thank the Government of Bermuda Department of Environment and Natural Resources for providing the necessary permissions and for their continued support; J. Maderios, A. Copeland, M. Meijas, N. Wellman, J. Labisko, S. Clayton-Green, D. Muldoon, N. Wright, J. Carney, D. D'Afflitto, M. Alonso, M. Shailer, P. Drew, R. Frith, K. Trott, R. Marirea, J.P. Rouja, P. Rouja and L. Thorne for assistance with fieldwork; R. McCrea and E. Matechou for assistance with statistical analyses; and the reviewers for their critiques. This research was funded by J. Summers Shaw, Chester Zoo, the State of Jersey, the British Herpetological Society and the Bermuda Zoological Society (Eric Clee Fund). This is contribution #295 from the Bermuda Biodiversity Project.
Author contributions
Analysis, writing: HT, RAG; experimental and statistical design: HT, RAG; field work: all authors; revision: RAG, GG, MEO.
Conflicts of interest
None.
Ethical standards
This research abided by the Oryx guidelines on ethical standards, was approved by Chester Zoo and the University of Kent Research and Ethics Committee and was conducted under permits issued by the Government of Bermuda's Department of Environment and Natural Resources. Capture and handling of skinks were undertaken in accordance with the conditions of the licence.
Open access
