Invasive species are a major threat to global biodiversity (Simberloff et al. 2013). Understanding the ecology and life history of an invasive species is integral to advancing management strategies. Because spatiotemporal patterns can directly influence population and community dynamics (Legendre and Fortin 1989, Collinge 2001), understanding the spatial ecology of an invasive species is a critical precursor to understanding opportunities for control and management (With 2002, Meier et al. 2014).
An invasive population of Burmese pythons has been established in southern Florida for several decades (Willson et al. 2011). They have caused major declines in meso-mammal populations in Everglades National Park (ENP; Dorcas et al. 2012, McCleery et al. 2015) and consequently have become a major focus of management and control efforts. However, their cryptic habits lead to extremely low detection probabilities (Nafus et al. 2020) and this may bias perceptions of habitat use toward locations where pythons are most likely to be detected by humans.
The field of spatial ecology explicitly acknowledges the crucial role played by spatiotemporal patterns in ecological processes and how they influence population and community dynamics (Legendre and Fortin 1989, Collinge 2001). Understanding the ecology and life history of Burmese pythons is integral to advancing management strategies for this species within its invaded range (With 2002, Meier et al. 2014). The first research study using radiotelemetry for Burmese pythons within the Florida Everglades (ENP) commenced in 2005, and tracking efforts continued through 2012 (Hart et al. 2015). This study revealed the aspects of python spatial ecology, including their ability to navigate (Pittman et al. 2014), home range sizes and movement patterns (Hart et al. 2015), and details of habitat use (Hart et al. 2015, Walters et al. 2016). Further field evaluations of whether spatial ecology may inform the development of control tools followed (i.e., the Judas technique; Smith et al. 2015, 2016), while other efforts relied on theoretical modeling to characterize movement ecology (Bonneau et al. 2016, Mutascio et al. 2017a) or on presence-only data (Mutascio et al. 2017b), rather than telemetry. Most of these studies were focused within ENP. However, the Greater Everglades ecosystem comprises multiple habitat types extending well beyond ENP, inclusive of area outside the boundaries of ENP, and it was not understood whether findings from ENP were applicable elsewhere in the area occupied by pythons.
Our objective was to describe broad patterns of the spatial ecology and habitat use of invasive Burmese pythons in Collier County, southwestern Florida. To do this, we analyzed a VHF radio-tracking data set of 25 pythons to quantify their home range size, movement rates, and second- and third-order habitat selection. We compared our findings to those from in ENP and discuss the implications of our findings for areas of peninsular Florida outside of the Greater Everglades Ecosystem.
Materials and Methods Study siteThe study area in eastern Collier County, Florida, USA, is within the Big Cypress Basin Watershed, a western component of the Greater Everglades Ecosystem (Duever et al. 2005). The study was conducted on public conservation lands in Rookery Bay National Estuarine Research Reserve (RBNERR, 445.2 km2), Collier Seminole State Park (CSSP, 26.0 km2), Picayune Strand State Forest (PSSF, 295.4 km2), and on adjacent private lands (Fig. 1).
Enlarge this image.Fig. 1. Study area for radio tracking 45 Burmese pythons. The study area is northwest of Everglades National Park on Florida’s southwest coast. It encompasses the city of Naples, as well as Rookery Bay National Estuarine Research Reserve, Collier Seminole State Park, and Picayune Strand State Forest.
The Big Cypress bioregion is similar to the freshwater Everglades system in habitat diversity, although natural communities are primarily forested and form a mosaic of habitat types, as opposed to vast expanses of herbaceous communities (Duever et al. 2005). Rookery Bay National Estuarine Research Reserve includes nearly 100,000 acres of coastal habitat in southwestern Florida, roughly half of which is subtidal. Although most of the area is aquatic or marine, it includes a variety of ecosystems from mangrove forests and freshwater wetlands to upland pine flatwoods and rare xeric oak habitats in relict dune ridges (Barry et al. 2013). The developed portion of the study area includes both urban land cover and agricultural land cover. Historically, rain-driven water in this region primarily flowed south to the Gulf of Mexico (Klein et al. 1970). The hydrologic cycle can be described in simplest terms as bimodal, with a summer wet season and a winter dry season, and the associated rise and fall of a shallow, expansive, and ephemeral sheet-flow water regime. This surficial aquifer is above ground during the wet season and below ground during the dry season.
Python tracking and capturesWe encountered Burmese pythons during both active searching and conducting radiotelemetric activities. A two- or three-person team searched along canal banks, agricultural levees, or roadsides for Burmese pythons basking or hidden in the vegetation. Searches were primarily conducted from the back of a truck but were conducted on foot in places inaccessible to vehicles. Locating radiotelemetry subjects during the breeding season (December–March; Smith et al. 2016) often led to additional adult pythons of both sexes in healthy body condition, some of which were selected for use in the study.
Snakes selected for the radiotelemetric study were captured by hand, measured and weighed, and surgically implanted with a Holohil model AI-2 radiotransmitter (25 g; 158–170 MHz; Holohil Systems, Ontario, Canada) using standard surgical procedures (Reinert and Cundall 1982). We recorded mass, lengths (total length and snout–vent length), ventral scale count, and pattern photographs while under anesthesia. After implantation, we allowed pythons at least 24 h of recovery time before they were released at the original capture location.
We radio-tracked pythons using a RA-23K VHF two-element antenna (Telonics, Mesa, Arizona, USA), Yagi three-element antenna (Titley Scientific, Columbia, Missouri, USA), and truck-mounted RA-5A omni-directional antenna (Telonics, Mesa, Arizona, USA) attached to a R-1000 telemetry receiver (Communication Specialists, Orange, California, USA). We recorded visual observations or obtained the approximate location of a snake within five meters whenever possible. In general, we located all pythons on foot approximately once a week during the breeding season (December–March), unless the individual was inaccessible (e.g., in dense mangrove swamps or aquatic habitats with water depths greater than 1 m), and approximately twice per month for the remainder of the year (April–November). When a study animal was located, we recorded its location (UTM coordinates, NAD83, obtained with handheld GPS unit), habitat type, and the time of observation. Unless a physical measurement was required, we attempted to minimize disturbance to the animal. To monitor long-range movements, we used a Cessna Skyhawk 172 with wing-mounted RA-2AHS antennae (Telonics, Mesa, Arizona, USA) to execute aerial surveys twice per month.
Data analysisTo define the specific extent of the study area for spatial analyses, we fit a 100% minimum convex polygon (MCP) around all locations of all pythons. We then added a 1-km buffer around the entire polygon to capture areas near the extreme locations. Finally, we removed the sliver of the resulting polygon that extended out into the Gulf of Mexico by clipping it with the digital elevation model (see Second-order habitat selection). We radio-tracked pythons until mortality occurred or the study ended. For this spatial analysis, we cut off data collection in April 2018, and we retained only individuals with at least 15 locations over at least 250 d. After filtering individuals, we summarized each individual’s movement history. We defined December–March as breeding and April–November as nonbreeding seasons. All geographic data processing and analyses were done in either QGIS version 3.0.0 (
An animal’s home range is often defined as the area an animal uses during the course of its normal daily activities, such as foraging, mating, and raising young (Burt 1943), and a wide range of methods exist for estimating this space from tracking data (Kie et al. 2010). For comparison to previous python tracking studies (Hart et al. 2015), we fit 95% minimum convex polygons (MCPs), 95% kernel density estimates (KDEs) or home ranges, and 50% KDEs or core use areas to each python’s entire data set. The MCP is the simplest method and has been in use the longest (Blair 1940). A 95% MCP removes the 5% of locations farthest from the centroid of all locations, and then fits a convex polygon by connecting the outermost locations. We fit 95% MCPs for each python using the function mcp() in the package adehabitatHR (Calenge 2006). The utilization distribution is the probability distribution that describes an animal’s position in two-dimensional space. It can be estimated using a nonparametric kernel method on a discrete raster grid, which results in the kernel utilization distribution (KUD). This KUD can be used to estimate the home range, resulting in the kernel density estimate (KDE) of the home range (Worton 1989). A 95% KDE is the smallest area that contains raster pixels summing to 0.95, and a 50% KDE is the smallest area that contains raster pixels summing to 0.5. We interpreted the 95% KDE as representing the home range, while we interpreted the 50% KDE as representing the core area.
We fit the estimated KUD in R using the function kernelUD(), estimated the size of the 95% and 50% KDEs using the function kernel.area(), and created the contours of the 95% and 50% KDEs for mapping using the function getverticeshr(), all from the package adehabitatHR (Calenge 2006). We used least-squares cross-validation (LSCV) to select the smoothing parameter (h), as recommended by Worton (1995). However, this LSCV parameter could not be minimized for many cases (which is a common problem), and so we took the average value of h (255) from those cases that did converge, and set the smoothing parameter for all individuals to this value.
VHF telemetry locations are usually irregularly sampled, as opposed to GPS data, which are often collected on a set schedule (Cagnacci et al. 2010, Smith et al. 2018). Our VHF technology resulted in temporally irregular sampling, which presents an issue for step-based analyses where the step (the link between consecutive locations) is treated as the sampling unit (Fortin et al. 2005). Given that step length distributions are typically skewed, with many small steps and only occasionally much longer steps, summary statistics (i.e., mean or median step length) do not accurately represent the true distribution of steps over a given time period. Despite these challenges, we were still interested in reporting mean daily movements, as have been reported in previous studies of python movement (Hart et al. 2015). Following the methodology of Hart et al. (2015), we connected successive locations, calculated the straight-line distance between them, and divided by the total number of days between locations. To better represent the distribution of steps, we report the mean and each quartile (minimum, 25th percentile, median, 75th percentile, maximum) for each animal. We also report the mean, minimum, and maximum number of days between locations for each animal to better represent the sampling regime.
Habitat selectionResource selection is key to understanding an animal’s spatial ecology and occurs when an animal uses a resource in greater proportion to that resource’s availability (Johnson 1980, Manly et al. 2002). It is therefore critical to define what is meant by availability, and Johnson (1980) recognized that selection occurs at multiple scales. First-order selection is placement of the species’ geographic range on Earth. Second-order selection is the placement of an individual’s home range within the species’ range. Third-order selection is the use of particular sites (e.g., feeding or bedding) within that home range. Fourth-order selection is the use of particular items (e.g., a particular food resource) at those particular sites. Understanding first-order habitat selection by an invasive species in its new range is often difficult, as their ranges often have yet to reach equilibrium (Phillips et al. 2008), and applying insights from that species’ native range might better approximate the realized niche than the fundamental niche (Rodda et al. 2011).
To understand how pythons select habitats in southwestern Florida, we focused here on the second and third scales of habitat selection, which are comparable to findings from the Everglades. For second-order habitat selection, we combined all pythons and all seasons to estimate home range placement. For third-order habitat selection, we separated the seasons into breeding and nonbreeding periods. During the breeding season, males and females are often found together, so we lumped together both sexes for that analysis. During the nonbreeding season, males and nongravid females should be selecting habitat in which they can find prey, whereas gravid females should be selecting nesting and brooding habitat. We therefore grouped males and nongravid females together for an analysis of feeding pythons, and we separately analyzed gravid female nest placement. We considered that nest placement could occur anywhere in the study area, and thus treated that subset as a second-order analysis. We used resource selection functions (RSFs; Boyce and McDonald 1999, Manly et al. 2002) to estimate habitat selection at the two scales, which are defined to be any function, which is proportional to the probability of use (Manly et al. 2002).
Second-order habitat selectionWe used a use–availability design to estimate habitat selection, where the observed telemetry points were considered used, and we randomly sampled from the entire study area to create points considered available. This random sampling was done by randomly placing 10× as many random points within the study area as there were observed points, so that the null probability of selection was ~9% (i.e., 10/11 points are background and 1/11 are observed). We randomly placed these locations in R using the function spsample() from the package sp (Pebesma and Bivand 2005, Bivand et al. 2013). Then, we extracted the underlying spatial covariates for both the observed and random points using the function extract from the package raster (Hijmans 2017). Using the extracted spatial data, we fit a set of RSFs by using a generalized linear mixed model (GLMM) with a binomial distribution (i.e., mixed logistic regression) to estimate the relative probability of selection. We included individual python identity as a random intercept to account for unequal sampling and individual variation (Hebblewhite and Merrill 2008). We included habitat variables in the model as fixed effects, and we represented habitat using three variables: elevation, land cover, and distance to urban areas. Elevation came from the digital elevation models (DEMs) created by the U.S. Geological Survey (USGS) as part of the 3D Elevation Program (3DEP;
Previous work has shown that Burmese pythons may avoid urbanization (Reichert et al. 2017), so we included distance to urban areas as a habitat covariate. We calculated the distance to urban areas in QGIS by using the GDAL tool Proximity (raster distance) to calculate the distance from each cell of our raster (again, matching the elevation raster) to the nearest urban feature in the CLC map. Elevation and distance to urban were continuous numeric variables, so we included quadratic transformations of each as predictors in our model to allow for a possible parabolic relationship (i.e., a peak or valley in some intermediate values). We also evaluated collinearity between these continuous predictors prior to model fitting to ensure this was not an issue (Pearson’s r = 0.002). Land cover was a categorical variable, and so we treated it as a factor in R. To build our candidate model set, we began with a full model and then used backward selection to remove insignificant variables and grow the model set. The full model estimated 10 fixed effects: an intercept, five parameters for the land cover classes (with the sixth being captured by the intercept), two parameters for elevation (the slopes of the linear and quadratic terms), and two parameters for distance to urban (the slopes of the linear and quadratic terms). If the quadratic term was found to be a significant predictor in our model, we always forced the inclusion of the linear term for that variable. We selected the model that best fit our data using AICc (AIC corrected for small sample sizes). AICc tables were constructed using the function model.sel() from the package MuMIn (Barton 2016). We used the resulting top model to draw inferences about how pythons select habitat. We also evaluated goodness of fit for that model by using the pseudo-R2 method of Nakagawa and Schielzeth (2013), and we reported both the marginal R2 (R2GLMM(m); fixed-effects only) and conditional R2 (R2GLMM(c); fixed and random effects). From the top model, we created model prediction figures to illustrate how each variable was affecting the relative probability of selection. We used the predict() function in R to generate model predictions, and used the function bootMer() from the package lme4 to estimate the 95% confidence intervals around each prediction via bootstrapping. Lastly, we plotted the predicted relative selection on a map of the study area to visualize where hot spots of high selection would occur.
Third-order habitat selectionWe again used a use–availability design to estimate third-order habitat selection, but this time we separated each animal’s trajectory by season, and we defined availability as those habitats found within that animal’s home range (i.e., their 95% KDE) for that season. As above, we placed 10× as many random background locations as we had observed locations in that season. We included all animals of both sexes in the breeding season data set, and only males and nongravid females in the nonbreeding season data set. We extracted elevation, land cover, and distance to urban data from the underlying rasters using the same procedure as for second-order habitat selection. We built our candidate model set using backward selection from a full model, as with the second-order analysis, but this time each term of our full model contained an interaction with the categorical season variable (two levels, breeding or nonbreeding), thus doubling our number of parameters from 10 for the second-order full model to 20 for the third-order full model (i.e., all the parameters from the second-order model were fit for each of the two seasons). We again selected a top model using AICc, drew inference from that model, and we calculated the pseudo-R2 statistics to measure goodness of fit. We created prediction figures with 95% confidence intervals using the same method as above, and then, we mapped the predicted selection across the entire study area.
Nest site selectionWe excluded any gravid females from the nonbreeding season analysis of third-order habitat selection. Instead, we determined nest sites of all radio-tracked gravid females and combined this with information on old python nests (as evidenced by eggshells) and recent clutches that were exposed by heavy equipment operations. To consider that females could place their nest anywhere in the study area (not necessarily within their home range) and include non-tracked individuals, we treated this as a second-order habitat selection analysis (i.e., the whole study area was sampled to estimate availability). We followed procedures very similar to the second-order habitat analysis described above, including generating 10× as many background locations as we had observed locations. However, for this analysis, we dropped the random effect for the individual level (because we had only a single observation for almost all animals), thus making our model an ordinary GLM as opposed to a GLMM. We again selected a top model using AICc, drew inference from that model, and we calculated the pseudo-R2 statistics to measure goodness of fit. For this GLM, we calculated pseudo-R2 using the method of Nagelkerke (1991). We produced prediction figures as for the previous analyses; however, estimates of standard error for a GLM are computed by the predict() function in R, and therefore, bootstrapping was not required to generate 95% confidence intervals.
ResultsBetween January 2013 and June 2018, we tracked 45 adult Burmese pythons (30 male and 15 female). After filtering out pythons with insufficient data points, we were left with 25 individuals for our analysis (Table 1). Of these 25, nine were female (36%) and 16 were male (64%). The average total length (TL) of females was 360 cm (SD = 68), and the average TL of males was 321 cm (SD = 35). The average mass of females was 31.6 kg (SD = 13.3), and the average mass of males was 18.8 kg (SD = 7.2). Females were tracked for an average of 848.8 d (SD = 493.3), during which time they were located an average of 51.6 times (SD = 22.1). Males were tracked for an average of 794.1 d (SD = 370.5), during which they were located an average of 62.6 times (SD = 22.1).
Table 1 Summary of tracked Burmese pythons included in the analysis.
TL is total length.
Home rangesThe areas of home range metrics were variable but within the range of previously reported ranges (Table 2) for other pythons tracked in south Florida. The 95% MCPs ranged in size from 0.9 km2 to 18.3 km2, and average size of the 95% MCP for all pythons was 5.3 km2 (SD = 4.4). Average size of the 95% MCP for male pythons was 6.75 km2 (SD = 4.83), whereas average size of the 95% MCP for female pythons was smaller (2.82 km2, SD = 1.72). The 95% KDEs ranged in size from 3.5 km2 to 14.8 km2 (Fig. 2), and average size of the 95% KDE for all pythons was 7.5 km2 (SD = 2.9). Similarly, average size of the 95% KDE for male pythons was 8.49 km2 (SD = 2.96) and average size of the 95% KDE for female pythons was smaller (5.62 km2, SD = 1.86).
Table 2 Home range sizes for all tracked pythons.
| Python ID | 95% MCP (km2) | 95% KDE (km2) | 50% KDE (km2) |
| 17-AGA | 4.0 | 6.9 | 1.6 |
| 14-ALE | 1.9 | 5.5 | 1.0 |
| 15-BRA | 2.6 | 6.2 | 1.4 |
| 01-ELV | 15.8 | 14.8 | 2.7 |
| 35-END | 3.0 | 5.8 | 1.6 |
| 26-FIO | 2.7 | 4.7 | 0.7 |
| 36- FRE | 10.8 | 10.1 | 2.4 |
| 04-GEN | 4.5 | 6.6 | 1.2 |
| 22-GRE | 6.3 | 8.1 | 1.1 |
| 24-JAE | 7.5 | 10.3 | 2.8 |
| 21-JOH | 18.3 | 13.7 | 2.3 |
| 13-KIR | 7.4 | 7.1 | 1.6 |
| 19-LUT | 8.6 | 11.0 | 2.1 |
| 27-MAL | 7.8 | 9.4 | 2.5 |
| 06-MAR | 0.9 | 4.0 | 0.8 |
| 03-NOO | 1.7 | 4.9 | 0.9 |
| 20-ORR | 1.9 | 4.9 | 0.9 |
| 05-PAU | 1.1 | 3.5 | 0.8 |
| 28-SHR | 4.1 | 7.3 | 1.4 |
| 10-STE | 5.7 | 8.3 | 1.3 |
| 08-SUG | 4.4 | 7.1 | 0.8 |
| 02-SWE | 4.9 | 8.5 | 2.3 |
| 30-TRE | 1.6 | 4.2 | 0.9 |
| 09-VAL | 1.4 | 4.3 | 1.2 |
| 29-DOS | 4.6 | 9.2 | 1.6 |
MCP is minimum convex polygon; KDE is kernel density estimate.
Enlarge this image.Fig. 2. Kernel density estimate (KDE) home ranges. Maps show the locations of 95% KDE (A) and 50% KDE (B) home ranges in the study area for female (left) and male (right) Burmese pythons.
Pythons moved on average 0.07 km/d (Table 3), with a mean range from 0.02 to 0.14 km/d. The maximum movement observed corresponded to a daily movement rate of 1.94 km/d. The mean maximum daily movement rate for all pythons was 0.52 km/d, showing that while pythons have low average movement rates (10 s of m/d), they are capable of much longer steps (0.5–1.0+ km/d) between consecutive locations. Pythons often exhibited no movement over long periods of time, and only one male python (#36) exhibited a nonzero minimum movement rate (0.01 km/d, i.e., he was never found in exactly the same place twice consecutively).
Table 3 Movement rate summary statistics along with tracking frequency summary statistics.
| Python ID | Movement rate (km/d) | Sampling regime (Δ days) | |||||||
| Mean | Min | 25th % | Median | 75th % | Max | Min | Mean | Max | |
| 0.05 | 0.00 | 0.01 | 0.03 | 0.08 | 0.22 | 1 | 16.9 | 261 | |
| 14-ALE | 0.08 | 0.00 | 0.01 | 0.06 | 0.10 | 0.31 | 2 | 18.1 | 262 |
| 15-BRA | 0.05 | 0.00 | 0.01 | 0.04 | 0.07 | 0.20 | 6 | 24.3 | 258 |
| 01-ELV | 0.06 | 0.00 | 0.01 | 0.03 | 0.08 | 0.45 | 2 | 18.0 | 263 |
| 35-END | 0.08 | 0.00 | 0.01 | 0.05 | 0.10 | 0.62 | 1 | 12.9 | 265 |
| 26-FIO | 0.04 | 0.00 | 0.01 | 0.02 | 0.07 | 0.22 | 1 | 14.4 | 250 |
| 36- FRE | 0.14 | 0.01 | 0.02 | 0.07 | 0.16 | 1.02 | 1 | 10.9 | 52 |
| 04-GEN | 0.05 | 0.00 | 0.01 | 0.02 | 0.08 | 0.21 | 2 | 16.1 | 252 |
| 22-GRE | 0.07 | 0.00 | 0.01 | 0.04 | 0.10 | 0.84 | 1 | 14.4 | 257 |
| 24-JAE | 0.06 | 0.00 | 0.02 | 0.04 | 0.06 | 0.32 | 1 | 21.4 | 259 |
| 21-JOH | 0.09 | 0.00 | 0.02 | 0.04 | 0.11 | 0.83 | 1 | 12.9 | 274 |
| 13-KIR | 0.05 | 0.00 | 0.02 | 0.03 | 0.06 | 0.30 | 1 | 16.3 | 259 |
| 19-LUT | 0.11 | 0.00 | 0.01 | 0.04 | 0.12 | 0.99 | 1 | 18.3 | 279 |
| 27-MAL | 0.11 | 0.00 | 0.02 | 0.06 | 0.10 | 1.94 | 1 | 14.6 | 292 |
| 06-MAR | 0.11 | 0.00 | 0.02 | 0.08 | 0.13 | 0.38 | 1 | 20.6 | 251 |
| 03-NOO | 0.02 | 0.00 | 0.00 | 0.00 | 0.02 | 0.18 | 4 | 16.9 | 257 |
| 20-ORR | 0.05 | 0.00 | 0.01 | 0.03 | 0.06 | 0.52 | 1 | 11.6 | 41 |
| 05-PAU | 0.06 | 0.00 | 0.00 | 0.02 | 0.08 | 0.51 | 2 | 9.7 | 30 |
| 28-SHR | 0.10 | 0.00 | 0.03 | 0.06 | 0.13 | 0.72 | 1 | 16.4 | 257 |
| 10-STE | 0.05 | 0.00 | 0.01 | 0.03 | 0.07 | 0.21 | 1 | 14.5 | 54 |
| 08-SUG | 0.06 | 0.00 | 0.01 | 0.03 | 0.07 | 0.29 | 1 | 12.2 | 39 |
| 02-SWE | 0.04 | 0.00 | 0.01 | 0.02 | 0.04 | 0.51 | 1 | 22.5 | 264 |
| 30-TRE | 0.03 | 0.00 | 0.00 | 0.01 | 0.03 | 0.25 | 1 | 13.0 | 266 |
| 09-VAL | 0.05 | 0.00 | 0.00 | 0.03 | 0.06 | 0.25 | 2 | 15.7 | 252 |
| 29-DOS | 0.08 | 0.00 | 0.01 | 0.04 | 0.11 | 0.60 | 1 | 13.1 | 253 |
We had 1466 observed python locations from 25 individuals available for the second-order habitat selection analysis. We generated 14,660 background points, resulting in 16,126 points for the GLMM. The full model containing the CLC categories and quadratic terms for both elevation and distance to urban (Table 4) was the top model, receiving 100% of the model weight. The next closest model was 27.8 ΔAICc from the top model, indicating very strong support for the top model. The top model had R2GLMM(m) = 0.843 and R2GLMM(c) = 0.843, indicating a high proportion of variance explained, and also demonstrating that the random effect (individual ID) had no discernable impact on the intercept. The top model indicated that pythons selected positively for agriculture, freshwater wetland, saline wetland, and upland land cover classes, while avoiding open water and urban land cover classes (Fig. 3A). Relative selection was highest for pythons at an elevation of 0.58 m (Fig. 3B). Relative selection was highest for pythons at a distance to urban of 515 m (Fig. 3C). The map of relative selection shows hotspots of python relative habitat selection that occur primarily between the urban development and the wetter mangrove fringe (Fig. 4).
Table 4 Model selection table for second-order habitat selection analysis.
| Model | No. of parameters | AICc | ΔAICc | ω |
| CLC + Elev2 + Urb2 | 11 | 8254.90 | 0.00 | 1.00 |
| CLC + Elev + Urb2 | 10 | 8282.74 | 27.83 | 0.00 |
| CLC + Elev2 + Urb | 10 | 8437.83 | 182.92 | 0.00 |
| Elev2 + Urb2 | 6 | 8883.10 | 628.20 | 0.00 |
| Elev + Urb | 4 | 9492.17 | 1237.27 | 0.00 |
| Null | 2 | 9829.12 | 1574.22 | 0.00 |
All models were mixed models and thus contained an additional parameter for the variance of the random effect. Possible fixed effects considered in the model set were the cooperative land cover categories (CLC), elevation (Elev), and distance to urban (Urb). A superscript 2 indicates that a quadratic term was included for that variable. All models with a quadratic term also included the lower-order linear term. The null model contained only an intercept and the random effect variance. AIC is Akaike information criterion.
Enlarge this image.Fig. 3. Second-order habitat selection model predictions. The top model from our model set contained the cooperative land cover categories; (A) categories and a quadratic term for elevation (B); and distance to urban (C). Pythons were found to select for agriculture, freshwater wetland, saline wetland, and upland habitats, while avoiding open water and urban habitats. They selected maximally for habitats near 0.5 m in elevation and habitats around 600 m from the closest urban area. Vertical bars or dashed lines show the 95% confidence interval for the prediction, and the solid red line represents the null probability of selection (1/11 or ~9%).
Enlarge this image.Fig. 4. Map of second-order habitat selection throughout the study area. Hotspots of python relative habitat selection occur primarily between the urban development and the wetter mangrove fringe. Null selection was 0.09, so locations with estimated values lower than this (darker colors) are being avoided, while values higher than this (brighter colors) are being selected for by pythons.
We had 1,307 observed python locations from 24 individuals available for the third-order habitat selection analysis. Of these observed locations, 621 were during the breeding season, and 686 were during the nonbreeding season. We generated 13,070 background points and were left with a data set of 14,377 points for the GLMM. The top model in this model set was the full model, with an interaction between season and CLC, an interaction between season and a quadratic term for elevation, and an interaction between season and a quadratic term for distance to urban (Table 5). The top model was strongly supported, with 98.2% of the model weight, and with the second-best model 8.8 ΔAICc below the top model. The top model had R2GLMM(m) = 0.161 and R2GLMM(c) = 0.173, a low-to-moderate proportion of variance explained. The random effect (individual ID) had a larger impact on the intercept than in the second-order analysis, reflecting a small amount of individual variation in habitat selection.
Table 5 Model selection table for third-order habitat selection analysis.
| Model | No. of parameters | AICc | ΔAICc | ω |
| Seas × (CLC + Elev2 + Urb2) | 21 | 8445.32 | 0.00 | 0.98 |
| Seas × (CLC + Urb2) + Elev | 18 | 8454.11 | 8.79 | 0.01 |
| Seas × (CLC + Urb2) + Elev2 | 19 | 8455.48 | 10.16 | 0.01 |
| Seas × (Elev2 + Urb2) + CLC | 16 | 8463.63 | 18.31 | 0.00 |
| Seas × (CLC + Elev2) + Urb2 | 19 | 8470.80 | 25.48 | 0.00 |
| Seas × (CLC + Elev2) + Urb | 18 | 8475.06 | 29.74 | 0.00 |
| Seas × (Urb2) + CLC + Elev2 | 14 | 8475.64 | 30.32 | 0.00 |
| Seas × (CLC) + Urb2 + Elev2 | 17 | 8479.59 | 34.27 | 0.00 |
| Seas × (Elev2) + CLC + Urb2 | 14 | 8491.18 | 45.86 | 0.00 |
| CLC + Elev + Urb2 | 10 | 8504.61 | 59.29 | 0.00 |
| CLC + Elev2 + Urb2 | 11 | 8506.05 | 60.73 | 0.00 |
| Seas + CLC + Elev2 + Urb2 | 12 | 8508.00 | 62.68 | 0.00 |
| CLC + Elev + Urb | 9 | 8508.95 | 63.63 | 0.00 |
| CLC + Elev2 + Urb | 10 | 8510.19 | 64.86 | 0.00 |
| Elev2 + Urb2 | 6 | 8640.54 | 195.22 | 0.00 |
| Elev + Urb | 4 | 8664.81 | 219.49 | 0.00 |
| Null | 2 | 8763.51 | 318.19 | 0.00 |
All models were mixed models and thus contained an additional parameter for the variance of the random effect. Possible fixed effects considered in the model set were season (Seas), the cooperative land cover categories (CLC), elevation (Elev), and distance to urban (Urb). A superscript 2 indicates that a quadratic term was included for that variable. Many of the models propose an interaction between season and one or more other covariates—those interactions are shown in parentheses. All models with a quadratic term also included the lower-order linear term. The null model contained only an intercept and the random effect variance. AIC is Akaike information criterion.
The top model indicated that during the breeding season, pythons selected agriculture and upland land cover types, while strongly avoiding urban areas (Fig. 5A). It indicated that during the nonbreeding season, pythons selected freshwater wetland land cover, while avoiding open water and urban (Fig. 5A). We found strong selection of higher elevations, with this selection being more pronounced during the breeding season (Fig. 5B). We found a weak relationship between distance to urban and relative selection, with the predominant pattern being an avoidance of distances far from urban during the breeding season (Fig. 5C). The map of relative selection shows that selection on the landscape is primarily driven by elevation, with the caveat that heavily developed areas are avoided (Fig. 6).
Enlarge this image.Fig. 5. Third-order habitat selection model predictions. The top model from our model set contained interactions between season and (A) land cover categories, (B) a quadratic term for elevation, and (C) a quadratic term for distance to urban. Vertical bars or dashed lines show the 95% confidence interval for the prediction, and the solid red line represents the null probability of selection (1/11 or ~9%).
Enlarge this image.Fig. 6. Map of third-order habitat selection throughout the study area. Panel A shows breeding season relative selection, and Panel B shows nonbreeding season relative selection. Selection on the landscape is primarily driven by elevation, with the caveat that developed areas are avoided. Null selection was 0.09, so locations with estimated values lower than this (darker colors) are being avoided, while values higher than this (brighter colors) are being selected for by pythons.
We had 16 observed python nest locations available for the nesting habitat selection analysis. We generated 160 background points and were left with a data set of 176 points for the GLMM. The top model in this model set was the full model, containing the CLC categories and quadratic terms for both elevation and distance to urban (Table 6). The second-best model was just 1.2 ΔAICc from the top model and differed only in that it had only a linear term for distance to urban (it was only missing the quadratic term). Together, these top two models received 80.9% of the model weight. The next closest model was 2.01 ΔAICc from the top model, indicating weak support for separation of these top two models. That third model was missing the CLC categories, so the support for keeping land cover in the model is weak. The top model had R2 = 0.59, indicating a moderate proportion of variance explained. The top model indicated no significant selection or avoidance of any CLC class (Fig. 7A), although this could be due to low sample size, so mean patterns should be considered. The elevation variable showed the strongest pattern, and relative selection was highest for pythons at an elevation of 1.7 m (Fig. 7B). The relationship with distance to urban showed a large amount of uncertainty (as illustrated by the large confidence intervals), but was highest for pythons at a distance to urban of 394 m (Fig. 7C). The map of relative selection shows several hotspots (Fig. 8). These hotspots tend to be concentrated on the fringe of the urban development on the border of agricultural fields, or else in the sandy uplands.
Table 6 Model selection table for nest site habitat selection analysis.
| Model | No. of parameters | AICc | ΔAICc | ω |
| CLC + Elev2 + Urb2 | 10 | 73.34 | 0.00 | 0.52 |
| CLC + Elev2 + Urb | 9 | 74.51 | 1.17 | 0.29 |
| Elev2 + Urb2 | 5 | 75.35 | 2.01 | 0.19 |
| CLC + Elev2 | 8 | 88.54 | 15.21 | 0.00 |
| Elev2 + Urb | 4 | 90.34 | 17.00 | 0.00 |
| CLC + Elev + Urb2 | 9 | 91.51 | 18.17 | 0.00 |
| CLC + Urb2 | 8 | 93.77 | 20.43 | 0.00 |
| Null | 1 | 109.25 | 35.92 | 0.00 |
Possible fixed effects considered in the model set were the cooperative land cover categories (CLC), elevation (Elev), and distance to urban (Urb). A superscript 2 indicates that a quadratic term was included for that variable. All models with a quadratic term also included the lower order linear term. The null model contained only an intercept. AIC is Akaike information criterion.
Enlarge this image.Fig. 7. Nesting habitat selection model predictions. The top model from our model set contained the cooperative land cover categories; (A) categories and a quadratic term for elevation (B); and distance to urban (C). Pythons were not found to select for a particular land cover type. They selected maximally for habitats near 1.5 m in elevation and habitats around 375 m from the closest urban area. Vertical bars or dashed lines show the 95% confidence interval for the prediction, and the solid red line represents the null probability of selection (1/11 or ~9%).
Enlarge this image.Fig. 8. Map of nesting habitat selection throughout the study area. Hotspots of python relative habitat selection occur primarily between the urban development and the high and dry levees along the agricultural areas. Null selection was 0.09, so locations with estimated values lower than this (darker colors) are being avoided, while values higher than this (brighter colors) are being selected for as nest sites by pythons.
This long-term radiotelemetric study provides summaries of home range, movements, and habitat preferences for adult Burmese pythons inhabiting a coastal zone in southwestern Florida. Spatial patterns of habitat use for 25 pythons revealed selection for a mosaic of habitat types. This information can be incorporated into designing control strategies throughout the python’s invaded range.
Habitat selectionPrevious telemetry studies of Burmese python habitat selection at finer scales (second- and third-order habitat selection) were based on data collected in ENP (Hart et al. 2015, Smith et al. 2016, Walters et al. 2016). These studies found, for example, that pythons select for higher elevations in the marsh, such as tree islands. Through our analyses of this long-term data set, we were able to explain a large proportion of the variation in second-order habitat selection (R2GLMM(m) = 0.843). However, we found that pythons tend to place their home ranges at intermediate elevations, where they have access to both water and upland features. We found that they generally avoid urban areas, consistent with previous findings that the impacts of pythons on mammal communities are reduced near urban development (Reichert et al. 2017). However, pythons also selected for habitat that was somewhat near urban areas, with selection peaking at a distance of just 515 m from urban development. Specific land management practices in the urban environment such as horticultural operations that produce and store large volumes of plant waste material are attractive to small- and medium-sized mammals, potentially elevating prey levels for pythons along this gradient. Hobby farms on urban fringes often raise chickens and goats, which may provide prey cues for pythons. Similarly, housing developments bordering natural areas often include artificial lakes that attract aquatic birds, and in which we occasionally located telemetered pythons. Pythons therefore have multiple reasons for remaining in proximity to the urban interface.
The map of relative selection shows several hotspots of python activity, which are concentrated between urban development and low-lying mangrove fringe, especially where upland features are near water features (Fig. 4). The land cover class most strongly selected for was upland. Indeed, the closely related Indian python (Python molurus) often uses dry, sandy uplands in its native range (Schleip and O’Shea 2010), and even Python bivittatus needs higher ground for nesting and reproduction. Hunter et al. (2018) showed evidence that the python population in Florida may represent a hybrid of these two closely related species. In any case, the animals in this study exhibited strong second-order selection for upland habitats, which indicates that these are an important feature of python home ranges.
We were not able to capture the variability in third-order habitat selection and for second-order habitat selection (R2GLMM(m) = 0.161), likely because pythons select different dimensions of habitat such as prey availability or game trails, for which we lack data. We did find strong support for an effect of season on all of our covariates at this scale, especially with elevation. Whereas in the second-order analysis we found selection peaked at an intermediate value (~0.58 m), within their home ranges, pythons are selecting for the highest elevations available. This effect is even more pronounced in the breeding season, a period during which female pythons select dry ground for oviposition and nest attendance.
Our analysis captured an intermediate amount of variance in nesting habitat selection (R2 = 0.59); however, note that when comparing these values the methods for estimating pseudo-R2 for GLMMs (Nakagawa and Schielzeth 2013) and GLMs (Nagelkerke 1991) differ slightly. We again found that elevation was a strong predictor, and similarly to the second-order habitat selection analysis, the strongest selection occurred at an intermediate value (1.7 m). This analysis had the smallest sample size and so has larger uncertainties. Nonetheless, we were able to identify potential hotspots of nesting habitat in the study area (Fig. 8).
Home range sizeHart et al. (2015) reported 95% MCPs from ENP ranging from 1.7 to 87.4 km2, with a mean size of 22.5 km2. The average 95% MCP in the Everglades was thus 4.2 times bigger than those we report here. Hart et al. (2015) only calculated KDE home ranges for the five pythons for which they could minimize the smoothing parameter with least-squares cross-validation (LSCV). From those five pythons, they reported a mean 95% KDE of 7.3 km2 and a mean 50% KDE of 1.5 km2. These sizes were very similar to the home range sizes we reported here—the 95% KDEs in ENP were just 0.2 km2 smaller than ours (7.3 km2 vs. 7.5 km2) and the 50% KDEs in ENP were the same size as ours (i.e., 1.5 km2).
The Greater Everglades is an extensive interconnected wilderness environment, and this discrepancy in MCP estimates of home range size might indicate that higher habitat quality for invasive pythons in southwestern Florida—with a greater availability of habitat types and higher prey abundances—could provide these animals all they need in a more restricted area compared with ENP. However, the concordance in KDE estimates, especially for the core home ranges, might suggest that the core ranges are actually similar in size, with ENP pythons simply making longer occasional forays but otherwise occupying a similarly sized core area. Alternatively, the five pythons in Hart et al. (2015) for which the LSCV algorithm converged may not be representative of the entire data set, so interpretation of these patterns in KDEs should proceed with caution. Disentangling potential drivers of home range size is difficult in general and typically requires long tracking durations across a range of potential spatial and temporal drivers of space use (Börger et al. 2006). Additional tracking studies on Burmese pythons in portions of their invaded range that include a mix of habitat types away from the influence of the coast and the urban interface could provide information useful for further comparisons. The current study could also be extended or replicated to identify the relationship between home range size and prey availability or the effects of removal programs on home range size on Burmese pythons Florida in southwestern Florida over time. Further analyses utilizing dynamic Brownian bridge movement models (dBBMMs) may be useful for comparing spatial use patterns to Burmese pythons within the native range (Silva et al. 2020, Smith et al. 2021).
Movement ratesHart et al. (2015) reported mean daily movement rates of 0.04–0.18 km/d, very close to our reported range of 0.02–0.14 km/d. Hart et al. (2015) also reported that their maximum movement rates ranged from 0.11 to 1.40 km/d, similar to our reported range of 0.18–1.94 km/d. Thus, pythons in ENP and southwestern Florida exhibit similar movement rates, despite occupying different habitat types. The urban environment includes a network of canals and drainage ditches that provide movement opportunities for pythons. The maximum movement rate reported in southwestern Florida was from a male python (#27) during the breeding season and likely occurred while swimming along a canal at the urban interface.
Management implicationsDeveloping an effective control strategy for the nonnative Burmese python is challenging due to their extremely low detection probabilities (Nafus et al. 2020). This radiotelemetry study provides new summaries of python home range, movements, and habitat preferences for adult pythons inhabiting a habitat mosaic in southwestern Florida. Patterns observed here can be useful for designing search protocols and studies both within and outside of the immediate study area. For example, searches could be designed in winter months across areas the size of home ranges and within habitat types selected by pythons in this study. There are large-scale habitat differences between the expansive, herbaceous wetland system of southeastern Florida in ENP and the mosaic forested and coastal mangrove marsh systems in the southwestern portion of the State. Tracking data in this study indicated locations of increased python activity during the breeding season on or immediately adjacent to elevated habitat types or features. Within coastal conservation lands of southwestern Florida, this elevation corresponds with scrubby oak or xeric pine habitat features that are often concentrated and interspersed with long stretches of suboptimal breeding habitat such as mangrove forests and brackish marshes. Within the urban environment, increased python activity can be found along the elevated levee features adjacent to canals and associated with large-scale agricultural operations. These findings suggest that searching for adult pythons during the breeding season (November–April) in elevated habitat with appropriate cover that may include subterranean refugia (Metzger 2013, Bartoszek et al. 2018) may increase detection probability and removal rates for reproductive female pythons. Locating these potential features in the landscape could be facilitated using detailed Lidar (light detection and ranging) elevation maps.
Refinement of a management tool targeting adult, reproductive, female pythons can have strategic application for localized control efforts across the invaded range. The Judas technique (Taylor and Katahira 1988, Carrick et al. 1990, McIlroy and Gifford 1997) is based on the idea that a radio-tagged animal subject will seek out and associate with conspecifics during the course of its routine social and reproductive behaviors, thus allowing humans to detect and remove those associated individuals. Smith et al. (2015, 2016) reported on the potential for this control tool to yield breeding females, which for demographic reasons are a primary target for control. We located an additional 108 adult pythons (male, n = 53; female, n = 55) in direct proximity to our radio-tagged pythons. Average size for additional females was 3.22 m SVL and 30.4 kg and for additional males was 2.47 m SVL and 12.0 kg. We removed >2300 kg of reproductive python biomass during the study. Continued investment in the use of these scout snakes may help to identify additional areas where reproductive activity occurs in our study site, and could be applicable to other areas invaded by invasive boas and pythons with similar predictable annual associations among individuals.
This investigation is the first to quantify the scope and scale of invasive python movements and habitat use in southwestern Florida. Continued tracking of pythons across south Florida during the breeding season may contribute to additional removals of reproductively active snakes. Future tracking efforts focused on collecting data to estimate demographic parameters of the python population could be optimized by using this movement and habitat selection information.
AcknowledgmentsDr. Paul Andreadis was instrumental in study design and inception of this project. We thank Ian Easterling, Anthony Flanagan, Cailin Prokop-Ervin, Jaimie Kittle, Monica Lasky, Michelle Bassis, Katherine King, and Alex Furst for invaluable assistance in the field. We thank staff at Rookery Bay National Estuarine Research Reserve, Collier Seminole State Park and the Florida Fish and Wildlife Conservation Commission for logistical support. Research was permitted under animal care protocol USGS-FORT ACUC-2013-1. Radiotelemetry activities were permitted under Florida Fish and Wildlife Conservation Commission permits EXOT-15-29 EXOT-15-29a EXOT-15-35 and EXOT-17-06. Funding and project support was provided in part by the USGS Greater Everglades Priority Ecosystem Science Program, USGS Invasive Species Program, the Naples Zoo Conservation Fund, and Private Philanthropy through the Conservancy of Southwest Florida. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. The authors have no conflict of interest.
© 2021. This work is published under http://creativecommons.org/licenses/by/3.0/ (the “License”). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.