Malayan kraits (Bungarus candidus) show affinity to anthropogenic structures in a human dominated landscape
Study siteThe study area covers the campus of Suranaree University of Technology (SUT) and its surrounding landscape in Muang, Nakhon Ratchasima, Thailand (14.879° N, 102.018° E; Fig. 1). The university campus covers about 11.2 km2, and comprises a matrix of human modified lands interspersed with mixed deciduous forest fragments (at the onset of this study we identified there were 37 mixed deciduous forest fragments on campus, mean = 7.36 ± 1.48 ha, range = 0.45–45.6 ha [note, “±” is used for standard error throughout the text]). More than 15,000 students are enrolled at SUT, and there are numerous multi-story classrooms, laboratory and workshop buildings, residential housing, parking areas, eating and sports facilities, an elementary school, and a large hospital on the university campus. During the first term of the 2019 school year, 7622 students, as well as numerous SUT staff, lived in on-campus residential areas. The landscape surrounding the university is primarily dominated by agriculture, though there are also patches of less-disturbed areas as well as several densely populated villages and suburban housing divisions among the monoculture plots of upland crops (e.g., cassava, maize, and eucalyptus).Figure 1Study site map illustrating the land-use types spanning the area where the Malayan kraits (Bungarus candidus) were tracked in Muang Nakhon Ratchasima, Nakhon Ratchasima province, Thailand. Map created using QGIS v.3.8.2 (https://qgis.org/) in combination with Inkscape v.1.1.0 (https://inkscape.org/).Full size imageThe study site is located within the Korat Plateau region with an altitude range of 205–285 m above sea level. Northeast Thailand has a tropical climate, and the average daily temperature from 1 January 2018 to 31 December 2020 in Muang Nakhon Ratchasima was 28.29 °C, with daily averages ranging from 19.3 to 34.1 °C38. The region receives an average annual rainfall ranging from 1270 to 2000 mm39. There are three distinct seasons in northeast Thailand: cold, wet, and hot, each are classified by annual changes in temperature and rainfall. Cold season is typically between mid-October and mid-February, hot season is generally from mid-February to May, while the highly unpredictable rainfall of the wet season is predominantly concentrated between the months May to October39,40.Due to the representation of agriculture, semi-urban, and suburban areas with patches of more natural areas all within a relatively small area, we determined the university campus provided an ideal setting to examine how land-use features and human activity influence the movements of B. candidus. Additionally, past studies have indicated northeast Thailand hosts the most bites by B. candidus in Thailand29,33, making sites like ours ideal.Study animalsWe opportunistically sampled Malayan kraits captured as a result of notifications from locals and ad-hoc encounters during transit due to low detectability in visual encounter surveys, in addition to those discovered through unstandardized visual encounter surveys. Upon capture, we collected morphometric data, including snout-vent length (SVL), tail length (TL), mass, and sex (Table 1, Supp. Table 1). We measured body lengths with a tape measure, measured body mass with a digital scale, and determined sex via cloacal probing, all while the snakes were anesthetized via inhaling vaporized isoflurane. We then housed individuals with an SVL > 645 mm and mass > 50 g in plastic boxes (with refugia and water) prior to surgical transmitter implantation by a veterinarian from the Nakhon Ratchasima Zoo. We attempted to minimize the time snakes were in captivity awaiting implantation; however, delays arose due to the veterinarian’s availability, the snake being mid-ecdysis, or the snake having a bolus that needed to pass through the digestive tract before implantation (n = 21 implantations, mean = 5.02 ± 0.61 days, range = 0.60–13.02 days). The Nakhon Ratchasima Zoo veterinarian implanted radio transmitters (1.8 g BD-2 or 3.6 g SB-2 Holohil Inc, Carp, Canada) into the coelomic cavity using procedures described by Reinert and Cundall41, while the snake was anesthetized. We assigned each individual an ID according to sex and individual detection number (e.g., M02 = a male was the second B. candidus individual documented during the study). We released the implanted individuals as close as possible to their capture locations (mean = 65.31 m ± 13.7 m, range = 0–226.42 m), though on six occasions we moved individuals ≥ 100 m because the individual came from either residential areas or a busy road (all but one were moved 800 mm; thus, nine of the males were adults and four were juveniles (though two of the males had an SVL > 720 mm, and therefore likely sub-adults). The single telemetered female was an adult.Individual tracking durations varied (mean = 106.46 ± 15.36 days, range = 28.5–222.77 days; Supp. Fig. 1), as many individuals were lost due to unexpected premature transmitter failures (n = 5) or unsuccessful recapture efforts due to individuals sheltering under large buildings as the transmitter reached the end of its battery life (n = 4). We only recorded one confirmed mortality in the study, M01, who was killed by a motorized vehicle when crossing a road (n = 1). Another three individuals were lost due to unknown reasons, which may have been due to premature transmitter failure, mortality, or the animal moving beyond radio signal despite extensive search efforts. Thus, we only successfully recaptured and re-implanted five individuals (M01 once, M02 twice, M07 once, M27 once, and M33 twice). Transmitter batteries generally lasted approximately 90–110 days, so we aimed to replace transmitters after ≥ 90 days of use. At the end of the study, only one individual was successfully recaptured to remove the transmitter.Data collectionWe used very high frequency radio-telemetry to locate each telemetered individual on average every 24.20 h (SE ± 0.41, 0.17–410.0 h; see Supp. Fig. 2 for distribution of tracking time lags). We aimed to locate each individual’s shelter locations once each day during the daylight (06:00–18:00 h); however, we were occasionally (n = 34 days) unable to locate a snake for several consecutive days when we were unable to obtain radio signal due to an individual having moved far away or deep underneath a large structure. There were also a few occasions where we were unable to track snakes due to prolonged and heavy rainfall (n = 4 days), as the moisture damages equipment, or other reasons (n = 4 days). We additionally located snakes nocturnally (18:00–06:00 h) ad hoc and in an attempt to observe nocturnal behaviors and movement pathways when animals were active. We defined fixes as any time a telemetered individual was located, and relocations (i.e., moves) as the occasions where we located an individual > 5 m from its previous known location.Each day we manually honed in on signal via a radio receiver to locate individuals (as described by Amelon et al.42, and recorded locations in Universal Transverse Mercator (UTM; 47 N World Geodetic System 84) coordinate reference system with a handheld global positioning system (GPS) unit (Garmin 64S GPS, Garmin International, Inc., Olathe, Kansas) directly above the sheltered snake. We generally approached within one meter of sheltering snakes during daylight to precisely record shelter locations and identify shelter type. Since we could not visually confirm snake locations, we methodically eliminated all possible locations where the snake could possibly be while at close range with the minimum possible gain on the radio receiver.Telemetered kraits tended to be inactive and sheltering underground during the daylight, thus we were confident that our diurnal location checks would not affect their movements. However, in some cases we resorted to determining an individual’s location via triangulation, where multiple lines cast from different vantage points towards the snake intersect on the snake’s location on the GPS, allowing us to determine the animal’s coordinate location from approximately 10–30 m away. This helped ensure that we recorded locations with greater accuracy when snakes sheltered underneath large buildings, as it allowed us to move away from large structures that hindered the GPS accuracy. This technique was also implemented during some nocturnal location checks when a snake was believed to be active among dense vegetation, in an attempt to prevent disturbance of the animals’ natural behavior. While we did hope to gain visual observations of active individuals during the night, we exercised more caution during nocturnal location checks, typically maintaining a minimum distance of approximately 5 m in attempt to lessen the chances of disturbing an active individual’s behavior. If the animal was active we recorded the animal’s observed behavioral state (i.e., moving, feeding, or foraging). When the radio signal was stable and the individual was not visible, we recorded the animal’s behavior as “sheltering”. We strived for an accuracy of 5 m difference), and land-use type (e.g., mixed deciduous forest, human-settlement, semi-natural area, agriculture, plantation; see Supp. Figs. 3 and 4 for photos of land-use types), behavior (e.g., sheltering, moving, foraging, or feeding), and shelter type (e.g., anthropogenic, burrow, or unknown, note we also recorded if we suspected the shelter to be part of a termite tunnel complex due to a close proximity to a visible termite mound; Supp. Fig. 5).During each location check we recorded the straight-line distance between the current and previous locations (distance moved/step length) with the GPS device. We then used step-lengths to summarize their movements by estimating the mean daily displacement (MDD; the total distance moved divided by the number of days the snake was located) and mean movement distance (MMD; the mean relocation distance, excludes distances ≤ 5). In order to limit biases due to some snakes being located multiple times within a given day/night, we limited our sample for estimating MMD and MDD to only include a single location per day. This was accomplished by manually removing “extra” nocturnal location checks that occurred within the same day, making sure to have all shelter relocations present within the dataset. When calculating MDD, we used the total number of daily location checks rather than the number of days between the individual’s tracking start and stop date since there were some days where individuals were not tracked. We also used the same one location check per day dataset to calculate movement/relocation probabilities and to examine each individual’s MMD, MDD, and relocation probability for the overall tracking duration as well as for each season.When feasible, we positioned a Bushnell (Bushnell Corporation, Overland Park, Kansas) time lapse field camera (Trophy Cam HD Essential E3, Model:119837) with infrared night capability on a tripod spaced 2–5 m from occupied shelter sites. We positioned the cameras so that we may gather photos of the focal snake as it exited the shelter site and/or behaviors exhibited near the shelter. We programmed the cameras using a combined setting, including field scan, which continuously captured one photo every minute, along with a motion sensor setting, which took photos upon movement trigger outside of the regular 1-min intervals.Space use and site fidelityAll analyses and most visualizations were done in R v.4.0.5 using RStudio v.1.4.1106 43,44. We attempted to estimate home ranges for the telemetered B. candidus individuals using autocorrelated kernel density estimates (AKDEs) using R package ctmm v.0.6.045,46 in order to better understand the spatial requirements of B. candidus. However, examination of the variograms revealed that the majority of the variograms had not fully stabilized (i.e., limited evidence of range stability in our sample), and many individuals had extremely low effective sample sizes (21.82 ± 9.75, range = 1.49–135.75; Supp. Table 4). Therefore, we do not report home ranges in this text, as the AKDE estimates would violate the assumption of range residency and either underestimate or misrepresent B. candidus spatial requirements. We also examined the speed estimates resulting from fitted movement models. Resulting variograms and tentative home range estimates are included in a supplementary file for viewing only (Supp. Fig. 6, Supp. Table 4). The original code is from Montaño et al.47.Since our data was not sufficient to estimate home range size for the telemetered B. candidus, we instead used Dynamic Brownian Bridge Movement Models (dBBMMs) with the R package move v.4.0.648 to estimate within study occurrence distributions. We caution readers that these are not home range estimates but instead modeling the potential movement pathways animals could have traversed49. Use of dBBMMs not only allows us to estimate occurrence distributions for each individual, thus helping us better understand the animal’s movement pathways and resource use, but it also allows us to examine movement patterns through dBBMM derived motion variance50,51. We selected a window size of 19 and margin size of 5, to catch short resting periods with the margin, while the window size of 19 is long enough to get a valid estimate of motion variance when the animals exhibit activity/movement. Contours however are somewhat arbitrary; therefore, we used three different contours levels (90%, 95%, 99%) to estimate dBBMM occurrence distributions (using R packages adehabitatHR v.0.4.19, and rgeos v.0.5.5), and show the sensitivity to contour choice52,53.All movement data, either including initial capture locations or beginning with the first location check ~ 24 h post release, was used for production of both the AKDEs and dBBMMs for each individual. We also estimated dBBMM occurrence distributions for each telemetered individual with the exception of M29, which only made three small moves within a burrow complex during the short time he was radio-tracked before transmitter failure.We compared space use estimates to two previously published B. candidus tracking datasets34,36, and one unpublished dataset shared on the Zenodo data repository54, all originating from the Sakaerat Biosphere Reserve (approximately 41 km to the south of our study site): two adult males from within the forested area of the reserve [one tracked every 27.8 ± 0.99 h over a period of 103 days, the other tracked every 38.63 ± 11.2 h over a period of 30.58 days]34,54, and a juvenile male from agriculture on a forest boundary [tracked every 50.19 ± h for 66.91 days]36. The previous studies on B. candidus only tracked the movements of a single individual each, had coarser tracking regimes, and used traditional—fundamentally flawed methods55,56—to estimate space use34,36. Therefore, we ran dBBMMs with these previous datasets using the same window (19) and margin size (5).To quantify site reuse and time spent at sites (residency time) we used recursive analysis with the R package recurse v.1.1.257. We defined each site as a circular area with a radius of 5 m around each unique location (matching the targeted GPS accuracy). Then we calculated each individual’s overall number of relocations, each individual’s total number of relocations to each site, and each individual’s site revisit frequency and residency time at each unique site. Then we plotted revisited locations on a land-use map with space use estimates (95% and 99% dBBMM) in an attempt to help identify and highlight activity centers for telemetered individuals (see Supp. Figs. 7–13). All maps were created using Quantum Geographic Information System (QGIS v.3.8.2).Habitat selectionWe used Integrated Step Selection Function models (ISSF) to examine the influence of land-use features on the movements of B. candidus at both the individual and population levels. We included movement data from all male individuals that used more than one habitat feature in our ISSF analysis. Therefore, we excluded F16 and M29 who both only used settlement habitat. Excluding M29 was justified by the individual having been tracked for the shortest duration (19 days) and had the fewest number of moves (n = 3), thus there were not enough relocations for ISSF models to work effectively. Using modified code from Smith et al.51 that used ISSF with Burmese python radio-telemetry data, we used the package amt v.0.1.458 to run ISSF for each individual, with Euclidean distance to particular land-use features (natural areas, agriculture, settlement, buildings, and roads) to determine association or avoidance of features. Cameron Hodges created all land-use shape files in QGIS by digitizing features from satellite imagery and verified all questionable satellite land-use types via on-ground investigation.The semi-natural areas, plantations, mixed deciduous forest and water bodies (such as irrigation canals and ponds which have densely vegetated edges) were all combined into a single layer of less-disturbed habitats which we refer to as “natural areas”. All feature raster layers were then converted into layers with a gradient of continuous values of Euclidean distances to the land-use features, and were inverted in order to avoid zero-inflation of distance to feature values and to make the resulting model directional effects easier to more intuitive. We were able to generate 200 random steps per each observed step (following Smith et al.51), due to the coarse temporal resolution of manually collected radio-telemetry data (i.e., we were not computational limited when deciding the number of random locations). Higher numbers of random steps are preferable as they can aid in detecting smaller effects and rarer landscape features59.To investigate individual selection, we created nine different models testing for association to habitat features, with one being a null model which solely incorporated step-length and turning angle to predict movement60, five examining land-use features individually (agriculture, buildings, settlement, natural areas, roads), and the other three being multi-factor models. Each model considers distance to a land-use variable, step-length, and turn-angle as an aspect of the model. After running each of the nine models for each individual, we then examined the AIC for each model, point estimates (with lower and upper confidence intervals), and p-values in order to identify the best models for each individual and determine the strongest relationships and trends among the samples. We considered models with ∆ AIC More
