Introduction
Compared to other vertebrate taxa, studies of snake population trends within the Anthropocene have been severely lacking (Mullin and Seigel Reference Mullin and Seigel2011), although several regional and landscape level studies have identified severe declines (Winne et al. Reference Winne, Willson, Todd, Andrews and Gibbons2007; Reading et al. Reference Reading, Luiselli, Akani, Bonnet, Amori, Ballouard, Filippi, Naulleau, Pearson and Rugiero2010). These declines are particularly troubling given the lack of sufficient data to quantify anthropogenic effects on most snake populations. This data paucity is attributed to snakes’ elusive behaviour, sporadic activity and the lack of suitable methods for reliable detection (Durso and Seigel Reference Durso and Seigel2015). Additionally, snake detectability is considerably lower than other co-occurring vertebrates (Durso et al. Reference Durso, Willson and Winne2011). In temperate regions, cover boards are used to detect snakes that must actively thermoregulate (Sewell et al. Reference Sewell, Guillera-Arroita, Griffiths and Beebee2012; Ward et al. Reference Ward, Griffiths, Wilkinson and Cornish2017). In the tropics however, high stable temperatures and structurally complex habitats confound this method. Surveys of tropical snakes predominately rely on active searches via encounter surveys defined by time/distance spent searching, or transect sampling of pre-defined sites (Rödel and Ernst Reference Rödel and Ernst2004; McDiarmid et al. Reference McDiarmid, Foster, Guyer, Chernoff and Gibbons2012). Maximising survey effort, regardless of strategy, is the baseline for increasing detections of species with low detection probability (MacKenzie and Royle Reference MacKenzie and Royle2005). However, care should be taken when increasing specific aspects of survey effort. Previous research identified a decrease in Boiga irregularis detections per unit distance during surveys with quicker walking speeds, yet increases in the total number of detections per unit time compared to slower surveying speeds (Lardner et al. Reference Lardner, Yackel Adams, Savidge and Reed2019). Maximising snake detections by increasing survey effort will increase the number of snake detections, thus improving statistical power of subsequent population models. However, abiotic and biotic covariates associated with snake activity patterns/microhabitat use may affect snake detection probability, and thus the reliability of such models.
Activity and detectability of tropical snake species are often associated with multiple climatic, temporal and biotic drivers (Durso and Seigel Reference Durso and Seigel2015). These poorly understood activity associations include short-term climatic fluctuations such as rainfall, temperature and humidity (Daltry et al. Reference Daltry, Ross, Thorpe and Wüster1998; Oliveira and Martins Reference Oliveira and Martins2001; Brown and Shine Reference Brown and Shine2002). However, longer-term fluctuations, including seasonality, increased yearly temperatures and cyclical weather events such as El-Niño, are also associated with snake detectability (Brown and Shine Reference Brown and Shine2002; Zamora-Camacho et al. Reference Zamora-Camacho, Moreno-Rueda and Pleguezuelos2010; Chaves et al. Reference Chaves, Chuang, Sasa and Gutierrez2015). El-Niño results in both increases and decreases in rainfall, mast fruiting/flowering events and in many cases increased small mammal density (Curran et al. Reference Curran, Trigg, McDonald, Astiani, Hardiono, Siregar, Caniago and Kasischke2004; Nakagawa et al. Reference Nakagawa, Miguchi, Sato, Shoko and Nakashizuka2007; Wells et al. Reference Wells, Kalko, Lakim and Pfeiffer2007). Shifts in snake detectability and activity during these El-Niño periods are poorly understood due to a lack of baseline detection data. Therefore, determining snake species’ detection probability and their associations with short-term (climatological, temporal and survey effort) and long-term (El-Niño) factors will be useful for informing future studies which seek to elucidate reliable population estimates, particularly where mark-recapture is impossible.
Here, we provide the first estimates of detection probability and detectability associations for Bornean riparian snakes. Specifically, we addressed the following questions: (1) How are climatological, temporal and survey effort factors associated with the detectability of river-associated snakes and snakes in general? (2) How do detectability associations vary across these snake taxa? and (3) How does the detectability of these snakes, and Bornean snakes in general, shift as a result of El-Niño?
Methods
Study site
We conducted this study within the Deramakot (5°14–28’N, 117°19–36’E) and Tangkulap (5°17–30’N, 117°11–21’E) forest reserves in central Sabah, Malaysian Borneo (Supplementary material, Figure 1). The reserves encompass an 82,600 ha area of hilly dipterocarp forest at varying stages of regeneration following conventional and reduced impact logging. Within the reserves, we collected snake observation data from 39 stream transects (100 m in length), using Standardized Visual Transect Survey (SVTS) methods, and eight river sites (900–3100 m in length) using Visual Encounter Survey (VES) methods. The SVTS sites were surveyed at 3.3 m per minute, whereas VES sampling involved considerably greater survey speeds (20–40 m per minute). Surveys were conducted between March and July in 2017, 2018 and 2019 for stream transects and between March 2019 and July 2019 for VES sites. In both survey types, two surveyors walked the length of the transect/VES site at night (between 18:30 and 23:30 hrs), recording all visually detected snakes. All detected snakes were identified to species level following Stuebing et al. (Reference Stuebing, Inger and Lardner2014).
Coefficient estimate sizes of five covariates on the detection probability of the five snake species, including the “pooled species” species category: humidity, survey length, month, time since sunset (TSS) and year (divided by sampling year 2017–2019). Thinnest line indicates upper confidence intervals (2.75–97.5%), and thicker line indicates lower confidence intervals (25–75%). Significant effect sizes (those not overlapping zero) marked in red, non-significant effects sizes marked in black and covariates not included within the species optimum models marked in grey.
Detection covariate collection
Climatological and temporal covariates assumed to affect detectability in Bornean snake species and that have direct associations with snake detectability and activity in other regions (Daltry et al. Reference Daltry, Ross, Thorpe and Wüster1998; Brown and Shine Reference Brown and Shine2002; Durso et al. Reference Durso, Willson and Winne2011; Sewell et al. Reference Sewell, Guillera-Arroita, Griffiths and Beebee2012) were collected for each survey. These included temperature, humidity, maximum daily rainfall, total rainfall over 60 days, time since sunset (TSS) (0–4 hours), month (represented as ordinals), sampling year, survey distance (100–3100 m) and survey speed (3.3–40 m per minute). Temperature and humidity daily averages and maximum daily rainfall data were collected from a weather station located within 0.4–13.5 km of all survey sites. Although this station may not represent the subtle variations of temperature and humidity in each stream during surveys, we previously identified significant correlation between temperature/humidity values recorded at a subset of stream sites and this weather station. As such, the weather station appears to provide a suitable proxy for local climatological patterns. Rainfall data were also summed to represent the total rainfall over a 60-day period prior to each survey. Stream width, which is associated with the occurrence of multiple Bornean snake species (Stuebing et al. Reference Stuebing, Inger and Lardner2014; Asad et al. Reference Asad, Wilting, Siku and Rödel2020b), was recorded at each site. Stream width was recorded every 10 m of the stream transects and the mean taken. Stream width of VES sites was recorded roughly every 150 m along the site, with the mean taken.
Statistical analysis
We analysed detection data for the most frequently detected snake species (>10 detections) and a combination of all snakes detected within our study, hereafter referred to as “pooled snake.” To determine species/pooled snake detection associations, data were incorporated into single-species, single-season occupancy models using the package unmarked ver 0.12-3 (Fiske and Chandler Reference Fiske and Chandler2011). This method allows the modelling of a species’ detection (p) and occupancy (psi) probability as a function of covariates, without requiring individual recognition (MacKenzie et al. Reference MacKenzie, Nichols, Lachman, Droege, Royle and Langtimm2002). Prior to modelling, all covariates were scaled to have a mean of zero and a variance of one. We tested for collinearity between covariates using Spearman’s rank correlation in the package Hmsic ver 4.2-0, with substantially correlated covariates (coefficient >0.7) excluded from subsequent analysis (Dormann et al. Reference Dormann, Elith, Bacher, Buchmann, Carl, Carre, Marquez, Gruber, Lafourcade, Leitao, Munkemuller, McClean, Osborne, Reineking, Schroder, Skidmore, Zurell and Lautenbach2013). All covariate combinations including a null model (no covariate effects on detection/occupancy) were included in the model selection process. The covariate combinations resulting in the optimum model for each species/pooled data were determined via Akaike’s information criterion (AIC) model selection in unmarked (Burnham and Anderson Reference Burnham and Anderson2004). We conducted all analysis in R 3.6.1 (R Core Team 2019).
Results
During this study, we conducted 274 stream transect surveys (73, 117 and 84 surveys in 2017, 2018 and 2019, respectively) and 21 VES (2019 only), totalling 89.4 km. We detected 23 snake species, with only five species detected on ten or more occasions: Psammodynastes pictus, Opisthotropis typica, Boiga dendrophilia, Psammodynastes pulverulentus and Xenochrophis trianguligerus (Supplementary material, Table 1). Of these species, null detection probability (assuming no effects of covariates on detection) ranged from 0.05 ± 0.01 to 0.14 ± 0.05 (Table 1). The only significantly correlated covariates were survey distance and survey speed. As survey distance is a measure of survey effort and speed is a methodological covariate, we chose to include survey distance rather than speed in the subsequent model selection process. Optimum detection models for each species yielded a range of detection covariates with humidity the most represented (three species) followed by survey distance, month and TSS (two species each), with temperature, year and 60-day rainfall represented in one model each (Table 1, Figure 1). Humidity significantly improved detection probability for two species (P. pulverulentus and X. trianguligerus), whilst survey distance exhibited a significant positive association with P. pictus detectability. Later survey months exhibited a significant positive and negative association with the detectability of P. pulverulentus and O. typica, respectively. Opisthotropis typica detectability was significantly positively associated with 60-day rainfall. Stream width proved to be a poor predictor of occupancy, represented in only three optimum models and exhibiting large, unreliable confidence intervals. This is likely due to the low number of occupied sites by the five species and the lack of sufficient explanatory variables to determine patterns in site occupancy.
Results of the optimum single-species occupancy models for five snake species and the “Pooled snake” category (all 23 species) generated from the combined stream transect/VES dataset. The coefficient estimates ± the standard errors and significance (*) of covariates within the optimum models for each species are included. Covariates not included within the species/pooled snake optimum models are marked with a dash. Null detection and occupancy probabilities (models with no covariates included) for each species/pooled snake data are also provided.
Pooled snake data yielded a 0.29 ± 0.02 probability of detecting at least one of the 23 species within the null model (Table 1). Humidity, survey distance and year were associated with pooled snake detections, with survey distance showing a significant positive association. The 2017 sampling year was significantly negatively associated with pooled snake and O. typica detectability, with both significantly positively associated with the 2019 sampling period (Figure 1). Additionally, the total number of snakes encountered within the stream transects (surveyed in each of the three years), increased in the final year (2017 = 2.05, 2018 = 2.47 and 2019 = 4.53 snakes per km) (Supplementary material, Table 1). This increase coincides with an El-Niño event, which occurred throughout Malaysia between February and May 2019 (Payus et al. Reference Payus, Huey, Adnan, Rimba, Mohan, Chapagain, Roder, Gasparatos and Fukushi2020).
Discussion
Detection probability for species in this study was low, highly variable and exhibited species-specific associations with climatic, temporal and survey effort covariates.
Whilst increasing surveyed area has resulted in increased snake detections elsewhere (Lind et al. Reference Lind, Welsh and Tallmon2005; Ward et al. Reference Ward, Griffiths, Wilkinson and Cornish2017), survey distance was only significantly positively associated with the detectability of Psammodynastes pictus and the pooled snake records in our study. This lack of a uniform positive association across taxa between effort and detection probability may simply be due to a lack of statistical power, as we detected the majority of species (except P. pictus) on less than 16 occasions. Additionally, the low predictive power of survey distance at our site may be due to the dynamics of temperate compared to tropical forest ecosystems. In tropical forests, high, stable thermal conditions make the use of cover boards and checking optimal thermal regulation sites obsolete (Shine and Madsen Reference Shine and Madsen1996), most importantly the greater habitat heterogeneity and complexity increases the difficulties in visually detecting snakes (Lardner et al. Reference Lardner, Yackel Adams, Savidge and Reed2019). Tropical ecosystems also include large bodied predators (primates and pigs), which may exploit forest floor refugia and potential basking sites for reptiles (Jolley et al. Reference Jolley, Ditchkoff, Sparklin, Hanson, Mitchell and Grand2010; Headland and Greene Reference Headland and Greene2011). Regardless, increasing survey distance did increase the total number of snakes detected. Detection probability may therefore be higher during shorter surveys under optimal detection conditions rather than during longer surveys under sub-optimal conditions. As such, determining the associations between a species’ detection probability and climatological/temporal factors should be included in future studies.
We detected no uniform associations between climatological/temporal covariates and detectability across species. Only two species (Psammodynastes pulverulentus and Xenochrophis trianguligerus) exhibited similar, significant positive associations with humidity. Whilst P. pulverulentus and X. trianguligerus were detected resting nocturnally on riparian vegetation, both are active diurnal hunters (Stuebing et al. Reference Stuebing, Inger and Lardner2014; Asad et al. Reference Asad, Wilting, Siku and Rödel2020b; Charlton Reference Charlton2020). Higher humidity has been linked to increased snake activity in West Malaysia (Calloselasma rhodostoma (Daltry et al. Reference Daltry, Ross, Thorpe and Wüster1998)), Brazil (Bothrops atrox (Oliveira and Martins Reference Oliveira and Martins2001)) and the USA (Pantherophis alleghaniensis (Carlson et al. Reference Carlson, Williams and Langshaw2014)). This humidity-associated activity could be to reduce water loss during movements (Daltry et al. Reference Daltry, Ross, Thorpe and Wüster1998), influence ecdysis (Carlson et al. Reference Carlson, Williams and Langshaw2014) or to coincide with increased prey activity (Brown and Shine Reference Brown and Shine2002).
Later months of our survey seasons (ranging from March to July) were linked to a significant increase and decrease in detection probability for Psammodynastes pulverulentus and Opisthotropis typica, respectively. This could indicate some form of seasonality in the activity patterns of these species as these later months coincide with Borneo’s subtle dry season (Kitayama et al. Reference Kitayama, Ushio and Aiba2021). A study of Australian water-associated snakes identified significant variation in activity between months, with these shifts purportedly linked to either increased prey activity or unmeasured proximate factors (Brown and Shine Reference Brown and Shine2002). As a semi-aquatic, nocturnal snake specialised in small anuran and fish predation, O. typica ecology differs markedly from that of P. pulverulentus (a diurnal small lizard predator) (Stuebing et al. Reference Stuebing, Inger and Lardner2014; Charlton Reference Charlton2020). Additionally, O. typica was the only species with a significant (positive) association with long-term (60 day) rainfall patterns. Higher rainfall would also affect aquatic habitat dynamics (increased stream size, higher siltation and increased stream speed) and associated prey (changes in density and capturability). These may increase the available foraging habitat and prey availability for O. typica.
While year was only significantly associated with the detection of O. typica, it exhibited a significant association with the pooled snake records. Combined with the increase in snake detections per km in stream transects during 2019, this suggests a possible association between snake detectability and the El-Niño event of the same year. Increased snake-bite prevalence during El-Niño events has been recorded in Costa Rica (Chaves et al. Reference Chaves, Chuang, Sasa and Gutierrez2015) and Ecuador (Gonzalez-Andrade and Chippaux Reference Gonzalez-Andrade and Chippaux2010), with the authors suggesting increased snake activity as the likely cause. Although these increases were speculated to be associated with increased rainfall, the cause of possible El-Niño associations with snake detectability at our study site is highly debatable. In Sabah, the El-Niño event involved a significant drought period from February to May (Payus et al. Reference Payus, Huey, Adnan, Rimba, Mohan, Chapagain, Roder, Gasparatos and Fukushi2020), yet the 60-day rainfall covariate yielded poor support in pooled and species-specific models. However, the El-Niño event also triggered major shifts in aquatic habitats within our site. During the dry period, many water bodies either changed in composition (lower water volume, slower speeds and heavy vegetation build up) or dried up completely, which also resulted in shifts in amphibian occurrence. Loss and shift in water and potential prey availability led to increased activity in a terrestrial snake species (Elaphe obsoleta) in the USA (Sperry and Weatherhead Reference Sperry and Weatherhead2008). Additionally, the number of species associations with humidity suggests that tropical snakes may prefer more humid habitats (such as riparian forest) during climatic extremes. In our site, this may explain the increased detection of species not normally associated with aquatic habitats (pooled snake data), and the lack of support for year in the detectability of our riparian-associated species which may prefer larger more drought-resistant aquatic habitats.
Caveats
The associations outlined in this study provide insights into the detectability and activity of these species and the detection of snakes at our study site. Unfortunately, the analysis of data collected over three successive years did not account for shifts in site-level occupancy due to reduced impact logging (RIL) activities within eight stream transects during the second year. Stream and habitat measurements obtained from a long running amphibian-monitoring project identified associations between recent RIL activity and aboveground carbon density, stream siltation and stream width variance (Asad et al. Reference Asad, Abrams, Guharajan, Lagan, Kissing, Sikui, Wilting and Rödel2021). These covariates, whilst important for localized site level occurrence of anurans (Asad et al. Reference Asad, Abrams, Guharajan, Sikui, Wilting and Rödel2020a; Asad et al. Reference Asad, Abrams, Guharajan, Lagan, Kissing, Sikui, Wilting and Rödel2021), likely have less effect on the localized occurrence of snakes due to this taxa’s greater resistance to desiccation, generally larger home ranges and independence from aquatic habitats for egg laying (Gibbons et al. Reference Gibbons, Scott, Ryan, Buhlmann, Tuberville, Metts, Greene, Mills, Leiden and Poppy2000). However, another confounding factor due to the analysis of data collected over three years is unmeasured yearly population fluctuations. Snake population variance between years resulting from stochastic processes varies in magnitude between species and regions (Weatherhead et al. Reference Weatherhead, Blouin-Demers and Prior2002; Winne et al. Reference Winne, Willson, Todd, Andrews and Gibbons2007; Reading et al. Reference Reading, Luiselli, Akani, Bonnet, Amori, Ballouard, Filippi, Naulleau, Pearson and Rugiero2010). Although snake populations are susceptible to yearly climatic fluctuations, these usually have a delayed effect, resulting in population increases/decreases in subsequent years (Moreno-Rueda and Pleguezuelos Reference Moreno-Rueda and Pleguezuelos2007; Reading et al. Reference Reading, Luiselli, Akani, Bonnet, Amori, Ballouard, Filippi, Naulleau, Pearson and Rugiero2010). As the El-Niño event of this study occurred on the final year of sampling, we believe its influence on snake population dynamics will be negligible until the following year. As such, due to the limited impacts of RIL on stream habitats and the relatively stable climatic conditions over the course of the study (with the exception of the final year), we believe that population fluctuations had only a minor impact on our occupancy model results. A further caveat is the addition of eight sites with different survey methods in the final year, which may complicate our conclusions regarding an increase in snake detectability during El-Niño. Although this is certainly a confounding issue, we also recorded an increase in snakes detected per km within the stream transects (surveyed in each year) during the El-Niño period. As such, whilst we would caution against considering our models an indication of the strength of detectability increases during El-Niño, we still believe these results suggest a general increase in snake detections during this period.
Conclusions
This study contributes to our understanding of the detectability associations of these five poorly known species and identifies increased snake detection (of pooled species records) during an El-Niño year. We tentatively suggest that this increase in snake detectability along stream habitats could be linked to El-Niño induced water/prey availability and site-level humidity. Additionally, considering the mixed support for survey distance, climatological and temporal covariates within our models, we suggest that shorter surveys under optimal detection conditions in some cases may increase detectability compared to longer surveys under sub-optimal conditions. As such, studies of tropical snake species should include considerations for short-term (climatological and temporal) and long-term (El-Niño) effects on snake detection.
Acknowledgements
The authors would like to thank Vitor Vitalis, who provided invaluable assistance as the research assistant, and Roshan Guharajan and Andreas Wilting, whose assistance with data collection and the analysis made this project possible. We thank the Sabah Forestry Department, particularly Johnny Kissing and Peter Lagan. We also thank the Sabah Biodiversity Council (SaBC) for granting us the access licence to allow research at the site; research permit number: JKM/MBS.1000–2/2 JLD.7 (63). This manuscript was greatly improved by the comments of two anonymous reviewers.
Financial support
This research was funding by Columbus Zoo, Auckland Zoo, the German Society of Herpetology and Herpetoculture (DGHT), the Museum für Naturkunde Berlin and the Leibniz Institute for Zoo and Wildlife Research.
Supplementary material
For supplementary material accompanying this paper visit https://doi.org/10.1017/S0266467421000468
