Detection has been a long-standing challenge to monitoring populations of cryptic herpetofauna, which often have detection probabilities that are closer to zero than to one. Burmese Pythons (Python bivittatus [=Python molurus bivittatus]), a recent invader in the Greater Everglades Ecosystem of Florida, are cryptic snakes that have long periods of inactivity. In addition, management actions such as removal of every python encountered create challenges for estimating population size and quantifying effects of management using traditional statistical approaches. We used Bayesian analysis of data collected from 59 visual surveys (144 person-surveys) covering a total distance of 485.6 km (1,185.1 person-km) and radiotelemetry to estimate detection probability for Burmese Pythons. These estimates can improve interpretation of encounter and removal data. We found that detection probability ranged from 0.0001 to 0.0146 depending on whether effort units accounted for total human effort across multiple surveyors and statistical method used. On the basis of our surveys, detection probabilities for Burmese Pythons are therefore likely <0.05, but factors such as the number of searchers or time of day may improve detection probability. Traditional capture–recapture or visual surveys are, however, unlikely to yield accurate information on Burmese Python population size or trends across time without cost-prohibitive effort. Consequently, novel methods development to monitor or measure Burmese Python populations, including techniques better equipped to handle very low detection, are critically needed for informative and reliable inferences about population size or the management effects of python removal.
Population management of reptiles can be complex, in part because traditional field methods and data analysis techniques are poorly suited to handle species with low detectability (Steen, 2010; Durso et al., 2011). Within reptiles, cryptic snakes can be particularly challenging because of coloration and behavior, long periods of inactivity, and their use of fossorial and arboreal habitats (Durso et al., 2011). Despite this challenge, an important component of population estimation and species management includes calculating species-specific detection probabilities. Detection and capture probability provide quantitative information used to estimate the population size, number of survey visits required for a species to be detected at a given level of certainty, and the effort required to detect changes in population size (Kéry, 2002; Sewell et al., 2012). Methods for calculating detection probabilities that accommodate very few detections are a vital research need for reptiles, especially those of management concern.
Burmese Pythons (Python bivittatus [=Python molurus bivittatus]) are native to Southeast Asia and were once considered a subspecies of P. molurus, but were elevated to full species status because of distinct morphological characters, co-occurrence in parts of their respective ranges, and an apparent absence of natural interbreeding (Jacobs et al., 2009; Schleip and O'Shea, 2010). They are one of the largest snake species in the world, with adults exceeding 5 m in length. Despite their large sizes, they remain difficult to detect (Reed and Rodda, 2011). Invasive populations of Burmese Pythons have been established in southern Florida for over 2 decades (Meshaka et al., 2000; Willson et al., 2011a). Their presence has been linked to a decline of native mammals in Florida's Everglades ecosystem (Dorcas et al., 2012; McCleery et al., 2015), with potentially cascading ecological implications (Hoyer et al., 2017; Reichert et al., 2017; Willson, 2017). In the Everglades, they are known to consume a diverse range of vertebrate prey species, most of which are native and a subset of which are also classified as endangered under the U.S. Endangered Species Act (Snow et al., 2007; Dove et al., 2011; McCleery et al., 2015). Estimates of site occupancy, rate of spread, and population sizes are thus priority topics for the management of Burmese Pythons in Florida and the native species on which they feed.
Because of their cryptic coloration and behavior and the habitat Burmese Pythons tend to occupy in Florida (Hart et al., 2015; Walters et al., 2016), surveys are physically difficult and costly to implement. Efforts to quantify presence, population size, and effects of management actions have thus been hampered by an inability to effectively detect individuals (Falk et al., 2016). Alternative detection methods such as environmental deoxyribonucleic acid (eDNA) and “Judas snakes”—radio-tagged individuals that lead searchers to untagged animals—have been explored as methods to improve detection and removal rates (Piaggio et al., 2014; Smith et al., 2016). Although novel field methods for detection are invaluable, integrating their information into traditional statistical approaches to infer effort required for detection, occupancy, and population size estimation requires continued development. Better estimates of detection probabilities for Burmese Pythons are critically needed to monitor rate of spread, site occupancy, population size, and to subsequently determine whether removal efforts are contributing to population declines or reduced site occupancy.
This study was designed to estimate detection probabilities for Burmese Pythons from a combination of transect-based visual surveys and radiotelemetry, because traditional mark–release–recapture methods do not align with resource managers' goals to remove every python encountered. We used radiotelemetry on a subset of individuals to estimate detection probabilities for Burmese Pythons when they were known to be available for detection but did not attempt to estimate capture probability because of zero recaptures. We used detection probabilities to estimate survey effort (distance and number of site visits) required to detect at least one individual or declare their absence with 95% confidence (McArdle, 1990). We also evaluated environmental and survey conditions that correlated with detection.
Materials and Methods
We conducted our research at the southernmost tip of Florida, in Southern Glades Wildlife and Environmental Area (WEA). The Southern Glades WEA is in Miami-Dade County, owned by the South Florida Water Management District, and managed by the Florida Fish and Wildlife Conservation Commission to help protect wildlife habitat as part of the Greater Everglades Ecosystem restoration goals. The habitat consists of sawgrass marsh, marl prairie, and tree islands intermixed with artificial canals. We surveyed a single canal (C110) such that a survey consisted of walking south on the western side of the canal (4.1 km), crossing the canal (0.03 km), and walking north on the opposite side of the canal (4.1 km) and ending at the starting point on the opposite side of the canal (8.23 km total transect line; transect start and end located at approximately Universal Transverse Mercator (datum WGS84: 17R 551001, 2806523). We a priori selected the survey site because the canal was a known python hot spot on the basis of previous captures mapped on the Early Detection Distribution Mapping System (EDDMapS; https://www.eddmaps.org), had good visibility, was easily accessible, and thus optimized detection conditions.
We completed 59 surveys from 17 December 2012 to 10 April 2013 from 0822 h to 1452 h. We elected to complete surveys during diurnal hours to take advantage of increased basking behavior during cold months (Dorcas et al., 2011) and to better synchronize with hours when telemetry could occur. Surveys were typically completed by a two-person team, although as many as five individuals participated on some days. We captured three male pythons during 17 December 2012–19 January 2013 and implanted transmitters in them. Those individuals were removed for periods of 5–13 d to complete surgery. During their absence, up to five visual surveys occurred. We removed every female and all male pythons found after January 2013. Additionally, during this study participants in the 2013 Python Challenge® removed four pythons from C110 (Mazzotti et al., 2016), which may have reduced detection probability estimates. At the onset of each of our surveys we recorded time, temperature, humidity, and wind speed. We also recorded weather in the categorical format of rain, overcast, partly cloudy, or clear. We downloaded nocturnal minimum temperatures from the nearest National Oceanic and Atmospheric Administration weather station (Homestead GEN Aviation Airport, Florida US USC00084095) for the dates on which surveys occurred.
We transported three male pythons to the University of Florida, Fort Lauderdale Research and Education Center, Davie, Florida, and intracoelomically implanted two radio transmitters in each of them (frequency 167.000–167.999), approximately one-third of the body length anterior to the cloaca. The second transmitter served as a backup to minimize loss of snakes (Walters et al., 2016). We used two transmitter sizes: Holohil models SI-2 (11 g) and AI-2 (25 g) (Holohil Systems, Ltd., Ontario, Canada). We held snakes up to 11 d before transmitter implantation and released them at the point of capture within 48 h of surgery. After their release, we tracked pythons using handheld radio receivers on most days that a visual survey was conducted, resulting in a tracking interval of 2–8 d. Tracking occurred from 16 January to 11 April 2013. When we located an individual, we recorded location (±3 m), whether the python was visible, their behavior when visible, and a description of the habitat. At the completion of this study all pythons were removed and humanely euthanized using a captive bolt gun through their cranium, to adhere to Institutional Animal Care and Use Committee and permit requirements.
Detections were too rare to quantify the effect of observer on detection. Additionally, because our goal was to maximize total detection, we did not hold observer number constant by survey day; instead we increased the number of observers when additional labor was available. We were thus uncertain of the degree to which individual observers or having multiple observers conducting surveys on the same transect on a given survey day affected detection. Therefore, we computed survey effort in two formats: 1) cumulative kilometers surveyed across all survey days (km) and 2) cumulative kilometers surveyed across all survey days multiplied by number of observers (person-km). In other words, effort for a two-person survey team for a single day was calculated as either 8.23 km or 16.46 person-km. The latter calculation increased the per-unit effort and thus generated a more conservative estimate of python detectability. We estimated detection probabilities from two python data sets: total python encounters and those exclusive to pythons after radio transmitters were implanted and pythons were known to be in the survey area on the day of a survey. We used telemetry points to confirm that telemetered pythons were at C110 during visual surveys and thus available to detect. For calculations of survey effort specific to radiotelemetered animals, we only included surveys during which telemetered snakes were known to be within 10 m of C110 because they could be radio tracked. For telemetry-based detection calculations, we omitted detections before transmitter implantation because outside of the detection event presence on C110 was unknown, yielding zero detections for telemetered pythons. Because surgical implantation can have behavioral effects on snakes (Weatherhead and Blouin-Demers, 2004) and telemetered pythons may have had abnormal detection probabilities, using total captures despite their unknown availability during nondetections was desirable.
Because management agencies do not allow traditional mark–release–recapture methods for pythons and because of extremely low detection probabilities for cryptic snakes (Kéry, 2002; Christy et al., 2010; Durso et al., 2011), mark–recapture methods for estimating cryptic reptile population parameters may poorly fit the data and generate biased estimates (MacKenzie et al., 2002; Steen, 2010; Dorcas and Willson, 2013; Rodda et al., 2015). Pythons were also removed temporarily—for transmitter implantation—and permanently upon capture during these surveys and by Python Challenge participants, causing sampling to occur with and without replacement. We therefore used two Bayesian estimation approaches more traditionally used during risk analysis to estimate the probability of occurrence and associated 95% binomial confidence intervals (CIs) when zero to few events occur. Bayesian approaches allowed for removal and nonremoval sampling, as well as large numbers of surveys. The methods were thus suited to handle intensive but uneven sampling effort and extremely low detections (Tuyl et al., 2008). Bayesian probability of occurrence was interpreted as an estimate of Burmese Python detection probability (p) that accounted for effort, but not time to detection.
The first approach we applied was Jeffreys prior, which is a compromise between Laplace binomial estimation and the more traditional frequentist approach of maximum likelihood estimation (Tuyl et al., 2008). The Jeffreys interval does not bias the beta coefficient (p) to be centered close to 0.5. Jeffreys prior has default hyperparameters that specify an equal-tailed beta distribution of α = 0.5, β = 0.5, which we a posteriori modified to a posterior beta distribution (α = 1.5, β = 4). The modified distribution reflected the assumption that p ≥ 0 and unlikely to be p > 0.50 (Fig. 1). We made this assumption because detection cannot be negative and we detected pythons (p > 0), but studies consistently generate low detection coefficients for cryptic snakes (Kéry, 2002; Christy et al., 2010; Durso et al., 2011). Secondarily, we applied the method of variance estimates recovery (MOVER) developed from the Newcombe method using Jeffreys interval in the ratesci package (Laud, 2018). The MOVER method provided flexibility to account for temporary and permanent python removals during the survey period. For consistency, we modified the distribution to the a posteriori assumption of a reasonable detection distribution for Burmese Pythons defined above. All calculations were completed in program R v3.4.1 (R Core Team, 2013). We report p and 95% CIs for eight models on the basis of alternate effort calculation techniques for 1) only telemetered animals postimplantation and 2) all survey and detection events.
From the p estimates, we calculated the number of site visits required to generate confidence at a specified level for presence and absence. To calculate the number of site visits that would be required to achieve up to 100% confidence that a python would be detected when present (P), we used the formula P = 1 − [1 − p]K (where K = number of site visits) described by MacKenzie and Royle (2005) for each of the four p estimates generated for the Jeffreys method. To calculate the number of site visits that would be required to interpret lack of detection as true absence with 95% confidence for each of our p estimates, we used the formula described by Kéry (2002): N = log(0.05)/log(1 − p).
To evaluate if search conditions could improve the probability of detecting a python, we ran a binomial linear regression with a dependent variable of python detected (1) or not detected (0) during a given day. We included survey start time, temperature, humidity, and wind speed at the start of the survey as well as minimum temperature the previous night, observer number, and elapsed time (total time to complete an 8.23-km transect by all observers minus any time spent capturing a python) as continuous predictors with linear and quadratic effects on p. We also included season (breeding [December–March] = 1; nonbreeding = 0) and cloud cover (present = 1; absent = 0) as binomial predictors. Because of small sample sizes and comparatively large numbers of potential predictor variables collected, we used an information-theoretic approach to evaluate model support (Akaike's information criterion adjusted for small sample size [AICC]). We used a forward and backward stepwise approach to identify the top five candidate models using the step function and AICcmodavg package in program R (R Core Team, 2013). We selected the top model using the model with the lowest ΔAICC to report potential relationships between predictors and python detection from the five best-fit models (Burnham and Anderson, 2002). Only one of the top five models had a ΔAICC that was greater than 2, which is the minimum value recommended to differentiate between candidate models (Burnham and Anderson, 2002). Therefore, the frequency with which parameters appeared in the top five models is reported as a measure of effect size. We also report P values for the top model to provide interpretative information, although Burnham and Anderson (2002) recommend caution in interpreting P values for information-theoretic modeling approaches.
The completed 59 survey days (144 person-surveys) cumulatively yielded 485.6 km surveyed or 1,185.1 person-km after accounting for survey team size. Within those 59 survey days, 38 (90 person-surveys) were conducted on days in which one or more pythons were available for detection on the basis of the locations of telemetered snakes (Table 1). Six detections of six unique pythons occurred before any transmitter implantations (0.10 pythons/survey day or 0.005 pythons/person-km surveyed; Table 2). After implanting transmitters into the three males there were no detections of telemetered pythons during visual surveys, yielding zero detection events for telemetry-only p estimations (Table 1). Of 62 total tracking events across all three pythons, they were fully exposed only 14% ± 7% (range: 4–27%) of tracking events on average.
Detection probabilities for Burmese Pythons were low. Depending on statistical method and effort unit, p ranged from 0.0001–0.0146 (Table 3). The most liberal calculation for python detection, which ignored unequal availability across time and the effect of multiple observers, generated an upper 95% confidence limit of p = 0.03 (Table 3). Incorporating all python detections improved detection probabilities, as compared with only using telemetered animals (Table 3). Jeffreys methods, which did not account for animal-specific effort, generated more optimistic detection probabilities, such that incorporating variation in effort by animal using the MOVER method decreased detection probabilities by 8- to 10-fold (Table 3). The most optimistic detection probabilities suggested that at least 35–149 km of canal line need to be searched (Fig. 2), ideally broken across repeated site visits (Fig. 3; at least five visits to complete 35 km given our study design) to ensure detection where pythons were present.
Model selection yielded a top model that included observer number and a quadratic term for start time as predictor variables (Table 4). Observer number had a positive correlation with detection, such that for each extra observer up to five, probability of detection increased onefold (z = 2.5, P = 0.01, 95% CI = 0.3–2.3). Observer number was in all five of the top candidate models. Start time as a quadratic term was in two of the five top models and increased for later start times (z = 1.8, P = 0.06, 95% CI = −0.3 to 14.7) until roughly noon and then declined again (z = −2.0, P = 0.04, 95% CI = −17.9 to −0.9). Peak probabilities for detecting a python occurred when surveys were initiated from ∼1030 h to 1315 h (Fig. 4). Elapsed time was selected in two of the five top models and breeding season was selected in one (Table 4).
Even the least conservative models suggested very low probabilities that Burmese Pythons would be detected using diurnal visual surveys with one to five observers. The least conservative scenario generated similar estimates to capture probabilities for captive Burmese Pythons in a large mesocosm (Dorcas and Willson, 2013), whereas the more conservative scenarios suggested detection probabilities close to zero. Monitoring population trends for cryptic reptiles has been suggested to require cost-prohibitive effort because of their low detectability (Ward et al., 2017). Additionally, in an occupancy framework—which we did not model because of data collection methods—species with low detection probabilities require numerous site visits to confirm absence with reliability (Durso et al., 2011; Sewell et al., 2012; Steen et al., 2012). The number of visits required to declare absence is inversely proportional to detection probability, such that species with a detection probability < 0.05 may require >60 visits to confirm absence on the basis of occupancy models (Durso et al., 2011). In contrast to our detection estimates, field application of eDNA to Burmese Pythons has yielded a detection probability of 0.91 on the basis of occupancy estimation (Hunter et al., 2015), which suggests that management concerns such as range expansion may be more accurately monitored through eDNA sampling than via visual surveys. Although survey methods across the two studies are not directly comparable, average per survey detection across the six pythons detected was 0.03 ± 0.02 in this study, which was similar to estimates of 0.01 ± 0.02 python detected per day in the field when using Judas snake methods (Smith et al., 2016). Thus, although alternative detection techniques may improve occupancy estimation, Burmese Pythons are still likely among those species for which monitoring population trends requires extensive survey efforts and even still such methods would likely yield large confidence intervals. A study on box turtles, with a similarly low detection probability (0.03) as estimated here, yielded a population abundance estimate with a 95% CI that ranged from 28 to 1,360 individuals (Refsnider et al., 2011), a range that represents vastly different management goals and challenges. As a result, monitoring population trends or estimating total population size remains an elusive goal until better detection and population estimation methods are developed for cryptic species.
The least conservative calculations of p assumed that the number of observers did not improve detection probability. The regression analysis, however, indicated that observer number had a positive effect on detection probability. Detection probabilities estimated using person-km as a measure of search effort may thus have encompassed important variation because of observer effects, observer number, or total effort (Christy et al., 2010; Ward et al., 2017). Differences in human searcher ability is a known source of detection error, which can be mitigated for but not necessarily eliminated by searcher experience (Lardner et al., 2016). As such, the models that used person-km as the effort unit may have generated more realistic detection probability estimates for Burmese Pythons. Alternatively, survey effort was often bolstered through the assistance of volunteer observers. The perceived positive effect of observer number may, therefore, reflect that more people opted to participate when conditions were perceived as good rather than that greater numbers of people improved the probability of detection. Thus, the effect of observer number may instead reflect variance in environmental conditions that increased the odds that Burmese Pythons were detectable.
With only six total detections we anticipated that the model evaluating survey conditions that correlated with detection would have limited power to detect significant relationships. Therefore, the lack of an effect of any environmental variable that was measured in this study was probably not a reliable indication of the effects of environmental conditions on python detection. Factors such as minimum and maximum temperatures have been correlated with python movement (Hart et al., 2015), and thus environmental conditions almost certainly can affect detection. Given low detection probabilities, telemetry studies may be a more cost-effective and powerful method to understand environmental conditions that support python activity and detection. Our models did, however, support that survey conditions affected detection, with start time identified as a predictor of detection. Start time is likely correlated with air temperature effects on python behavior or diurnal activity patterns. Regardless, methods used during surveys and environmental conditions during the survey are likely to alter the probability of detecting Burmese Pythons.
During this study the field methods implemented could not fully control for known sources of bias. For instance, most of the pythons detected were male, as were all the telemetered animals. Male Burmese Pythons may have different movement patterns than females (Hart et al., 2015; Smith et al., 2016). Movement patterns in turn affect detection (Keiter et al., 2017; Siers et al., 2018). Consequently, male snakes can have higher detection coefficients than females (Christy et al., 2010). Additionally, the bulk of these diurnal surveys were conducted during the breeding season. Many snakes have seasonal or temporal changes in behavior that can affect detection and parameter estimation (Brown et al., 2002; Willson et al., 2011b). Habitat also has known effects on detection (Melbourne, 1999; Kéry, 2002). Because these surveys focused on a single transect line to minimize habitat heterogeneity, the detection probabilities reflect those expected along canals. In stark contrast to our highly accessible transect along a canal, most of the Everglades is remote, difficult to navigate, and seasonally inundated, which are effects that would be anticipated to reduce detection probabilities. Detection inferences for Burmese Pythons throughout the Everglades thus cannot be made, although overall detection probabilities are unlikely to be greater elsewhere.
A cold snap in 2010 resulted in an unknown proportion of Burmese Pythons dying in the Everglades (Mazzotti et al., 2016), such that Burmese pythons may have been less abundant during our study than otherwise (Falk et al., 2016). Population density should also affect detection such that denser populations should have greater detection probabilities (Keiter et al., 2017). Variance in python density across their invasive range will thus contribute to variance in detection probabilities in specific locations, particularly when considering the invasion front where populations exist at low densities. Although we did not know the true density of pythons at the survey site, >100 pythons have been detected or removed from the vicinity of the C110 canal over the last decade on the basis of EDDMapS points. Continued removal of numerous pythons and low detection probability argue against a small local population.
Results from this study suggest that Burmese Pythons in Florida have extremely low detection probabilities, much like other cryptic snakes (Kéry, 2002; Durso et al., 2011). Although sampling using attractive lures, such as artificial refuges, can provide detection enhancement tools for cryptic herpetofauna (Engelstoft and Ovaska, 2000; Michael et al., 2018), they have not been documented to function for large constrictors. Other methods such as line distance or N-mixture models are suggested to be ineffective for snakes (Rodda and Campbell, 2002; Ward et al., 2017). Novel tools, such as removal models, for invasive species population monitoring are viable alternatives to mark–recapture approaches, but tend to require relatively high removal rates to generate accurate population estimates (Davis et al., 2016). Consequently, novel method development of population estimation techniques that are better equipped to handle detection probabilities that are closer to zero than to one, low rates of removal, or improve detection for cryptic reptiles such as Burmese Pythons are critical needs for making informative and reliable inferences about population size or the effect of removal efforts on Burmese Pythons and other cryptic invasive species.
The benefits of obtaining reliable population estimates from ongoing control practices are the ability to monitor spatial and temporal changes in populations, evaluate the effectiveness of management strategies, and determine the cost/benefits of different management actions.
We thank L. Barraco, M. Collier, S. Connor, K. Creely, M. Curtis, J. Duquesnel, J. Edwards, S. Featherstone, C. Gillette, R. Hill, E. Metzger, M. Rochford, B. Russ, C. Speroterra, M. Squires, A. Staniewicz, J. Vinci, J. Wasilewski, N. Wasilewski, S. Williams, and others who contributed their time to complete python surveys. We also thank L. Bonewell for project management. This work was done in accordance with specifications listed under Florida Fish and Wildlife Conservation Commission (FWC) Permit # EXOT-12-107 and EXOT-13-05 and Animal Resource Center protocol # 012-08FTL. Work was supported by the U.S. Geological Survey (USGS) Greater Everglades Priority Ecosystem Science Program, South Florida Water Management District, USGS Invasive Species Program, British Broadcasting Corporation, and the University of Florida. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.