Invasive snake causes massive reduction of all endemic herpetofauna on Gran Canaria

Invasive snakes represent a serious threat to island biodiversity, being responsible for far-reaching impacts that are noticeably understudied, particularly regarding native reptiles. We analysed the impact of the invasive California kingsnake, Lampropeltis californiae—recently introduced in the Canary Islands—on the abundance of all endemic herpetofauna of the island of Gran Canaria. We quantified the density in invaded and uninvaded sites for the Gran Canaria giant lizard, Gallotia stehlini, the Gran Canaria skink, Chalcides sexlineatus, and Boettger's wall gecko, Tarentola boettgeri. We used spatially explicit capture-recapture and distance-sampling methods for G. stehlini and active searches under rocks for the abundance of the other two reptiles. The abundance of all species was lower in invaded sites, with a reduction in the number of individuals greater than 90% for G. stehlini, greater than 80% for C. sexlineatus and greater than 50% for T. boettgeri in invaded sites. Our results illustrate the severe impact of L. californiae on the endemic herpetofauna of Gran Canaria and highlight the need for strengthened measures to manage this invasion. We also provide further evidence of the negative consequences of invasive snakes on island reptiles and emphasize the need for further research on this matter on islands worldwide.


Introduction
Invasive snakes have become a growing conservation concern as they are causing the decline or extinction of native and endemic species in various regions around the globe [1,2]. For example, the brown tree snake, Boiga irregularis, in Guam and the Burmese python, Python molurus bivitattus, in Florida, have severely reduced native bird and mammal communities [3][4][5][6]. Additionally, native reptiles became one of the most frequent prey of B. irregularis in Guam once native birds had been extirpated [5,7], which altered skink and gecko assemblages [8]. Similarly, the invasive horseshoe snake, Hemorrhois hippocrepis, affects native reptiles in the Balearic Islands [9,10], the common wolf snake, Lycodon capucinus, is presumed to be primarily responsible for the demise of endemic reptiles on Christmas Island [11], and the Indian wolf snake, L. aulicus, feeds on the endemic herpetofauna of La Réunion [12], though impacts per se have not yet been demonstrated, only presumed. These cases indicate that invasive snakes can potentially impact insular herpetofauna around the globe, even though such impacts remain unnoticed in most cases (see [1,2]).
Native and endemic reptiles play an essential role on islands, replacing the ecological functions normally carried out by other taxonomic groups in continental areas (e.g. [13,14]), and being central in island trophic webs (e.g. [15]). This ecological role is enhanced by the high densities that insular reptiles often attain [13]. Thus, the decline of insular reptile populations can be expected to generate far-reaching impacts on native communities and ecosystems (e.g. [15,16]). Furthermore, species abundance plays a major role in population dynamics as it governs reproduction and survival (e.g. [17,18]), determines extinction risk [18], shapes sexual and natural selection [19,20], and influences other ecological attributes (e.g. [21,22]). Thus, a notable density decrease can result in major disturbances to demographic dynamics of insular reptiles. Therefore, it is crucial to explore the impacts of invasive snakes on island reptile abundance as a means of addressing any negative effects on insular ecosystems.
Our study focuses on the invasive California kingsnake, Lampropeltis californiae, a colubrid snake native to western North America and Mexico [23], that was introduced via the pet trade to the island of Gran Canaria (Canary Islands, Spain) in 1998 [24]. Lampropeltis californiae is known to heavily predate the Gran Canaria giant lizard, Gallotia stehlini, the Gran Canaria skink, Chalcides sexlineatus, and Boettger's wall gecko, Tarentola boettgeri [25], which are the only three native (and endemic) reptiles on the island. However, nothing is known regarding the effects of this predation pressure on the abundances of the three species (but see [26][27][28]). In this context, we aim to explicitly evaluate snake impact upon the abundance of the three endemic reptiles of Gran Canaria based on the underlying hypothesis that invasive snakes will deplete their populations.

Methods (a) Study zones
We measured the impact of L. californiae on the abundance of the endemic herpetofauna of Gran Canaria-one of the biggest islands of the Canary Islands-in all three zones invaded by 2017: North, East and South zones (figure 1). The East zone invasion started in 1998, although establishment was confirmed only in 2007, whereas the North and South zones were confirmed as incipient populations in 2010 and 2015, respectively [28]. All zones are covered by native scrublands and formerly cultivated areas [29] and presented comparable climatic conditions [30].
To define the invaded area in each zone, we first analysed all records of snake captures, observations or snake presence signals-skin shedding or excrements-from 2009 to 2017 (compiled and provided by GESPLAN S.A.) on QGIS Essen v. 2.14.2 (QGIS Development Team, https://www.qgis.org/en/site/). We calculated the minimum convex polygon containing the highest number of records no more than 300 m apart, including a 200 m circular buffer around each record to account for the average species home range [31,32]. Second, we determined for the biggest invaded zone (the East Zone) the number of sampling sites that could be undertaken according to the available resources. We then calculated the number of sampling sites we should adopt for the North and South zones to ensure equal sampling effort, following the same proportion of sites per km 2 as in the East zone. We located sampling sites in accessible places, less than 200 m from the closest snake record, and roughly 100 m or more from the closest sampling site. Third, we selected for each zone the same number of uninvaded sites as invaded ones, so that they presented similar elevational range, vegetation and habitat features, and placed them at least 450 m apart from any snake record from 2015 onwards.

(b) Endemic herpetofaunal abundance
We quantified G. stehlini density (sensu [33]) using two complementary sampling methods-spatially explicit capturerecapture (hereafter SECR [34]) and distance sampling [35]and accounted for C. sexlineatus and T. boettgeri abundance through active searches under rocks [36]. We sampled on warm and sunny days between 09.00 and 19.00 to match G. stehlini activity above ground, and from 08.00 to 20.00 to find C. sexlineatus and T. boettgeri under rocks. Finally, we identified the abiotic variables influencing our field data (for each of the three species separately) and used them as independent variables in each of our subsequent models (see electronic supplementary material, S1 for the analyses to correct for the influence of abiotic variables on our sampling results).
In order to apply SECR for the calculation of G. stehlini density, we set 8, 4 and 2 sampling sites in the East, North and South zones, respectively, half in invaded and half in uninvaded areas. Invaded and uninvaded capture sites were separated from each other by greater than 3 km in all zones. At each sampling site, we performed two 1 h trapping surveys per month between May and September 2018, except for one invaded and one uninvaded site in the East zone that we first visited in June. For later analysis, we grouped visits into five sampling periods (one per month). Within months, we interposed an average of 4.03 ± 1.47 days between trapping surveys to minimize the likelihood of learned responses to trapping. In each trapping survey, we used 15 pitfall traps (30 × 40 × 50 cm, L × W × D) baited with tomatoes and sardines, unevenly placed along with surface features (e.g. rocks, rock walls) a minimum of ca. 4 m apart, and minimizing sun exposure to avoid lizards' overheating. We registered trap coordinates on the Real-Time map from UTM Geo Map app and recorded the amount of time each trap was active (trap sampling time, min). We set this time to 0 when a trap was inactive at the end of the trapping survey. We marked with a non-toxic paint the abdomen of all captured G. stehlini to quickly identify recaptures between two consecutive surveys and photographed the tops of their heads to individually photoidentify recaptures throughout the whole period (see electronic supplementary material, S2). We used SECR-derived data to build a capture-history dataset on the secr package [36] with separate sessions for each capture site and sampling period combination (each session included two occasions); we set the trap argument to 'multi' and included trap sampling time as a measure of sampling effort in each sampling site. We used the resulting dataset to estimate G. stehlini density following Borchers & Efford [38] (see electronic supplementary material, S3 for further details).
To quantify G. stehlini density through distance-sampling, we selected 120, 60 and 24 transects of distance sampling (approx. 60 m each) in the East, North and South zones, respectively, equally divided between invaded and uninvaded sites. The lowest distance between invaded and uninvaded sites ranged  royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20211939 from 1.9 km in the South zone to 3.0 km in the North zone. We repeated each transect a mean of 4.15 ± 0.61 times from May to September 2018. While walking slowly and constantly (approximately equal to 1 km h −1 ) along with each transect, a single observer (J.C.P.) counted all G. stehlini observed within a 40 m band and estimated their distance from the observer in 1 m classes. We registered sampling duration time and recorded lengths and paths of each transect using IGN Mapas de España v. 1.0 (http://www.ign.es/web/dir-aplicaciones-moviles), to later calculate the surveyed surface for each transect. We built separate datasets for each month (n = 5) to ensure model convergence of subsequent analyses using the unmarked package [39]. To increase the accuracy of density estimations, we truncated all datasets at 8 m distance (i.e. we only considered lizard sightings within 8 m of the transect line for the abundance estimation), following the 5% cut-off criteria [40,41]. We used the resulting data to derive G. stehlini densities following Kéry & Royle [42] (see electronic supplementary material, S3 for further details).
To account for the cryptic habits during daylight of C. sexlineatus and T. boettgeri [42,43], we quantified species abundances through active searches under rocks [42] in a total of 18, 36 and 8 sampling sites in the North, East and South zones, respectively, equally distributed among invaded and uninvaded sites and separated by at least 1 km in all zones. From March to September 2019, a team of 2-5 people placed in parallel nonoverlapping transects sampled each site every ca. five weeks for a time summing to 40 min among all the observers. We counted the number of individuals found of C. sexlineatus and T. boettgeri and wrote down the total number of rocks lifted in each session as a proxy for refuge availability. We ran a linear regression for each species with the number of individuals per site as a dependent variable (transformed using Yeo-Johnson power transformation due to the existence of 0 values [44]) and previously detected abiotic factors as independent variables. We extracted the residuals from each linear regression as a proxy for each species' abundance in all following analyses.
(c) Impact of Lampropeltis californiae on endemic reptile abundances We analysed the impact of L. californiae upon the three reptile species by comparing the density or abundance scores obtained between invaded and uninvaded sites. We ran generalized linear mixed models on glmmTMB [45] and set the estimated density or abundance of each species per method as the dependent variable, invaded versus uninvaded sites as a fixed factor, and zone, sampling period and sampling site as random factors.
Owing to spatial and temporal heterogeneity in the density or abundance of the three species, we fitted the dispersion parameters as a function of the zone and sampling period in all models. Data of SECR and distance sampling followed a zeroinflated Gaussian distribution, whereas the remaining models assumed a Gaussian distribution. Since zero-inflation GLMs are two-step analyses that first model the probability of encountering zeros as a binomial process (zero-inflation model) and subsequently model the data with reduced zero-inflation (conditional model) [46], we referred our results to both parts of the model. We modelled the occurrence of zeros in SECR and data of distance sampling as a function of snake presence. We checked model residuals using the DHARMa package to ensure we met model assumptions [47]. We quantified model main effects (i.e. snake impact) for both the zero-inflation and the conditional parts of the models by calculating type-II Wald chi-square tests. Finally, we used QGIS Essen v. 2.14.2 to plot the density or abundance of endemic herpetofauna against snake records in each study zone. To accentuate the distribution and occurrence of snake records, we divided each study zone with a 250 × 250 m grid and calculated the number of snake records per cell.
We performed all statistical analyses in R v. 4.0.2 [48] (see electronic supplementary material, S4 for our R code file). All results are expressed as mean ± s.d. unless stated otherwise.

Results
We devoted 206.7 h of trapping to SECR (104.8 h in invaded and 101.8 h in uninvaded sites), capturing 1 individual of G. stehlini in invaded sites and 269 in uninvaded sites (227 distinct individuals); that is, 99.6% of captures occurred in uninvaded sites. The zero-inflation model for SECR showed that the probability of capturing individuals was significantly higher in uninvaded sites (Wald chi-square: The zero-inflation part of the model showed that the probability of G. stehlini density being greater than 0 lizards ha −1 was significantly higher for uninvaded than invaded sites (Wald χ 2 : x 2 1 = 101.19, p < 0.001). Gallotia stehlini densities using this sampling method were 8.80 ± 37.29 individuals ha −1 on invaded sites versus 112.03 ± 144.28 individuals ha −1 on uninvaded sites (figure 2). The conditional part of the model showed invaded and uninvaded sites presented similar densities after removing zero-inflation (Wald χ 2 : x 2 1 = 3.65, p < 0.056). We did not detect any lizard in 86.3% of invaded sites, whereas that only occurred in 13.7% of uninvaded sites; however, most of the invaded sites where G. stehlini was still present were located on the outer edge of the snake's invaded range (electronic supplementary material, S6 and figure S6.1). When considering each zone separately, density was 93.8% lower on average in invaded sites versus uninvaded sites. This difference was particularly notable in the North zone, with a density of 0 in invaded sites but 78.97 ± 106.54 individuals ha −1 in uninvaded sites.
Active searches accumulated a total sampling effort of 206.7 h for C. sexlineatus

Discussion
Our study indicates that L. californiae is highly likely to be responsible for a severe abundance reduction of the three endemic species of Gran Canaria's herpetofauna, with a more pronounced decrease in those zones invaded earlier.
This impact showed the lowest magnitude for T. boettgeri and the highest for G. stehlini, which is driven to extinction in most areas. We detected a reduction in the number of observed individuals of 99.6% and 93.5% of G. stehlini (SECR and distance sampling, respectively), 82.8% for C. sexlineatus and 52.1% for T. boettgeri in invaded sites compared to uninvaded sites.
Our results are consistent with the fact that G. stehlini and C. sexlineatus are the most consumed [25] and most preferred prey [49] of L. californiae on the island. The observed differences between the impact that the new predator is causing upon the three endemic species could also be related to potential differences in prey behaviour (e.g. activity patterns, refuge selection) and adaptation to the invasive predator. However, considering that no snakes have ever naturally occurred on Gran Canaria [50,51], all endemic reptiles in this island are likely lacking antipredator responses against L. californiae, as no co-evolutionary history exists between these taxa [52][53][54]. Alternatively, their vulnerability could vary according to each species's ecological habits and their interaction with the new predator habits. Lampropeltis californiae is a generalist and opportunistic predator that in its native range chiefly consumes lizard species that overlap with its diurnal activity cycles and microhabitat selection [55]. In this sense, the higher impact on G. stehlini could result from this species being a diurnal and surface-dwelling lizard that exhibits similar ecological habits as the native prey of L. californiae [55]. The slightly royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20211939 lower impact we detected in the case of distance sampling could potentially be explained by methodological differences between this method and SECR [56]. Alternatively, it could also relate to the spatial distribution of each type of sampling site, as distance-sampling allowed detecting lizard populations at the edge of the invaded range of L. californiae (see electronic supplementary material, S6), whereas SECR samples were located in the core of the invaded area, where lizards were mostly absent. In the case of C. sexlineatus, although the impact is also noticeable, it is lower with respect to G. stehlini, which could be potentially linked to the fossorial habits and the small size of this species [42] that may lie in crevices that are unreachable to snakes. Finally, since T. boettgeri is also a small-sized and nocturnal gecko that perches upside-down under rocks to take refuge during daylight [57], the lowest impact detected for this species abundance could be related to its different refuge selection and nocturnal activity to that of L. californiae [58]. Determining population density is always a challenging endeavour [59,60] that often involves a certain amount of limitations and uncertainty [56,60,61]. In our estimates, our main uncertainty derives from the SECR method, as the low within-session recaptures produced estimates with broad dispersion measures-in some cases exceeding the threshold proposed by Efford [34] (see electronic supplementary material, S5). Nonetheless, the absence of G. stehlini in most invaded sites and the consistency of our results irrespective of the sampling method constitute solid evidence of drastic lizard depletion. For C. sexlineatus and T. boettgeri, we did not estimate densities, but abundance indices are recognized as a cost-effective and reliable tool, so long as the magnitude of spatial and temporal variation in species detectability is lower than the specific trend that researchers aim to detect [60]. Thus, considering the magnitude of the abundance reduction for both species, our results are also a conclusive indication of the negative impacts of L. californiae predation pressure on them.
The existence of a high-magnitude impact in all invaded areas, regardless of the invasion time, matches expectations from predator-prey theory. Accordingly, a severe depletion of native species is expected from the early stages of a toppredator invasion [62]. The greater reduction of most endemic prey in those areas first invaded could relate both to long-term exposure to predation pressure and the endemic prey's inability to cope with the new predator. Also, the magnitude of the depletion could indirectly impact on the species fitness through Allee effects [63,64], which may further complicate species recovery even after snake removal [65]. Additionally, the decline of reptile populations in Gran Canaria could potentially lead to complex demographic processes, such as sourcesink dynamics between the invaded range and adjacent populations or even the degradation of subpopulation networks, as detected for other invasive predators [66,67]. As we currently lack empirical information about reptile population connectivity in Gran Canaria, further research is needed to test these hypotheses. On the other hand, the endemic herpetofauna plays a keystone ecological role as plant pollinators [13,68], seed dispersers [69,70] and consumers of invertebrates [42,71] in native ecosystems on the Canary Islands. Therefore, L. californiae is likely to trigger cascading ecological events, such as the disruption of pollination and seed-dispersal processes (e.g. [72][73][74]) or the alteration of food webs and the increase of invertebrate abundance (e.g. [15]).
From a global perspective, our study is one of the few to provide evidence of the capacity of invasive snakes to decimate endemic fauna on islands and drive native communities to collapse. Considering that invasive snakes have been successfully introduced to numerous archipelagos around the globe [1,75], there is an urgent need to evaluate additional impacts of these predators. Invasive snakes are often only narrowly distributed on islands or mainland locations [1,75], which might decrease our perception of their importance for the conservation of global biodiversity compared to more ubiquitous invaders [76]. On the other hand, the use of functional responses to predict the impact of invasive predators [77][78][79] might mislead our predictions regarding invasive snakes due to their low energetic requirements [80]. Nonetheless, despite being less widely distributed than other invasive predators, previous examples have shown that invasive snakes have the potential to cause irreversible biodiversity loss in the ecosystems they invade (e.g. [3,6,73]). Moreover, functional responses are highly context-dependent [81], and their use to predict the risk associated with snake invasion on islands should be taken with care. Regardless of the reasons behind the current lack of information on most invasive snakes, our study clearly underlines that research on their impacts is urgently needed to understand the threat they pose and to prevent further biodiversity loss on islands worldwide.

Conclusion and management implications
Our study provides compelling evidence that L. californiae is causing a massive reduction of the endemic herpetofauna on Gran Canaria, which may trigger cascading events that could increase the species's impact and eventually cause major ecological disturbances, as observed for other high-profile invasive snakes [72]. Sharing our results with the regional and insular administrations has increased political engagement, leading to a rise in funding for snake management and reinforced in-site control actions. The present study has led to the implementation of the Strategic Control Plan to Fight against the California kingsnake 2019-2022, which anticipates an increase in the number of workers involved in the control of the snake and the implementation of several actions to improve biosecurity between islands in the near future. This demonstrates the importance of raising awareness among local and regional stakeholders for any control action to be successful [82,83]. We emphasize the urgent need for efficient management measures to prevent L. californiae spreading further on Gran Canaria and the rest of the archipelago and to minimize its impacts. Considering experience from other invasive snakes [84,85], the use of acetaminophen-treated aerial baits could be promising in Gran Canaria [86]. Predation pressure on endemic reptiles could be further reduced by increasing snake capturability through the use of synthetic odour-based lures matching snake prey preference [49] or sexual pheromones [84]. Complementarily, endemic reptiles could be trained to avoid snake predation [87], following previous experiences with other species (including closely related ones [88]), which would enable species translocation to reinforce local populations and prevent their total extirpation. From a broader perspective, we underline the need to fill the current information gap regarding invasive snakes and their negative ecological consequences, so as to raise awareness and promote the implementation of effective management measures royalsocietypublishing.org/journal/rspb Proc. R. Soc. B 288: 20211939 to counter the ecological impact they cause on island biodiversity worldwide.
Data accessibility. Raw data on endemic herpetofauna abundance are available from the Figshare repository: https://doi.org/10.6084/ m9.figshare.14443727. R code is available in the electronic supplementary material, S4. The data are provided in the electronic supplementary material [89].