Compensatory mortality explains rodent resilience to an invasive predator

Abstract

Invasive Burmese pythons (Python bivittatus) in the Everglades of Florida, United States, have drastically reduced populations of mammals, yet populations of some rodents appear unaffected by the invasion. To understand this pattern, we radio-tagged cotton rats (Sigmodon hispidus) in areas of high and low python occurrence densities (hereafter occurrence) and quantified the effects of python occurrence, seasonality, and sex on their survival and cause-specific mortality. Cotton rat survival was not influenced by difference in python occurrence (hazard ratio = 1.32, 95% CI = 0.77–2.26, P = 0.30). However, cotton rats were at greater risk from mortalities caused by mammals in areas of low python occurrence. In areas with elevated python occurrence, we attributed most cotton rat mortalities to birds of prey (48.6%) and reptiles (non-python = 24.3%, python = 16.2%). Where python occurrence was relatively low, we attributed cotton rat mortalities to native reptilian (28.6%), avian (35.7%), and mammalian predators (35.7%) with no python-related deaths. In total, pythons were responsible for 11.8% of all cotton rat mortalities. Finding no difference in the survival of cotton rats, despite differences in the causative agents of mortality, suggests that predation pressure from an invasive predator was compensatory for cotton rat population dynamics. This type of compensatory mortality is common for small mammals and helps explain why mammal communities in python-invaded portions of the Greater Everglades Ecosystem are increasingly dominated by cotton rats and other rodents.

Las pitones birmanas (Python bivittatus) invasoras han reducido drásticamente las poblaciones de mamíferos en los Everglades de Florida, United States, pero las poblaciones de algunos roedores parecen no verse afectadas por la invasión. Para comprender este patrón, marcamos con radio ratas algodoneras (Sigmodon hispidus) en áreas de alta y baja densidad de ocurrencia de pitones y cuantificamos los efectos de las densidades de ocurrencia de pitones, la estacionalidad y el sexo en su supervivencia y causa de mortalidad específica. La supervivencia de las ratas algodoneras no estuvo influenciada por la diferencia en las densidades de aparición de pitones (HR = 1.32, IC del 95% = 0.77–2.26, P = 0.30). Sin embargo, las ratas algodoneras corrían un riesgo mayor de mortalidad causada por mamíferos en áreas de baja densidad de ocurrencia de pitones. En áreas con elevadas densidades de ocurrencia de pitones, atribuimos la mayoría de los eventos de mortalidad de ratas algodoneras a aves rapaces (48.6%) y reptiles (no pitón = 24.3%, pitón = 16.2%). Donde las densidades de ocurrencia de pitones fueron relativamente bajas, atribuimos los eventos de mortalidad de ratas algodoneras a reptiles nativos (28.6%), aves (35.7%) y mamíferos depredadores (35.7%) sin muertes relacionadas con pitones. En total, las pitones fueron responsables del 11.8% de todos los eventos de mortalidad de ratas algodoneras. No encontrar diferencias en la supervivencia de las ratas algodoneras, a pesar de las diferencias en los agentes causantes de mortalidad, sugiere que la resiliencia de las ratas algodoneras a un depredador invasivo fue una función de la mortalidad compensatoria. Este tipo de mortalidad compensatoria es común para los pequeños mamíferos y ayuda a explicar por qué las comunidades de mamíferos en las porciones del ecosistema del Gran Everglades invadidas por pitones están cada vez más dominadas por ratas algodoneras y otros roedores.

Invasive species are one of the greatest threats to biodiversity across the globe (Vitousek 1990; Wilcove et al. 1998; Mooney and Cleland 2001; Levine and D’Antonio 2003; Salo et al. 2007). This is particularly true of invasive predators, which can attain elevated densities because of naive prey bases and release from their own native predators (Lockwood et al. 2013; Doherty et al. 2015). These predators can directly and indirectly alter populations and communities across trophic levels (Mooney and Cleland 2001; Arim et al. 2006; Sih et al. 2010; Gordon et al. 2015) with marked effects on biodiversity and ecosystem functioning (Simberloff et al. 2013; Doherty et al. 2016).

Nonetheless, the response of prey populations to invasive predators can be varied. Some populations have been rapidly extirpated (Savidge 1987; McCleery et al. 2015; McCreless et al. 2016) while others remain at reduced densities (Fritts and Rodda 1998; Blackburn et al. 2004; Campbell et al. 2012; Dorcas et al. 2012) or appear unaffected (Banks et al. 2018). The apparent resilience of unaffected populations may be a function of three different responses: (1) prey experience additive mortality and decreased survival from an invasive predator, but declines are offset by increased fecundity or compensatory reproduction (Abrams and Rowe 1996; Krebs 1999; Osmond et al. 2017); (2) prey survival increases because invasive predators release prey from their native predators and there are minimal direct impacts from invasive predators (Letnic et al. 2009; Gordon et al. 2015); and (3) for prey with generally high mortality rates, mortality from an invasive predator can be compensatory due to substitution of risk with density-dependent decreases in mortality from other predators or factors (Schmitz 2007; Creel 2011). Understanding how and why some species appear resilient is critical to predicting and mitigating the negative impacts of invasive predators as they present a mounting challenge to the conservation of native wildlife.

In the Greater Everglades Ecosystem of southern Florida (hereafter the Everglades), United States, the invasive Burmese Python (Python bivittatus) has greatly reduced the region’s native mammalian mesopredators, such as raccoons (Procyon lotor), bobcats (Lynx rufus), and opossums (Didelphis virginiana; Dorcas et al. 2012; Reichert et al. 2017; Taillie et al. 2021). However, small mammal populations, including cotton rats (Sigmodon hispidus) commonly found in python diets (Snow et al. 2007), appear unaltered, potentially even increasing where other mammals have declined (Dorcas et al. 2012; Hoyer et al. 2017; Soto-Shoender et al. 2020; Burkett-Cadena et al. 2021). While this pattern may be a function of the three mechanisms detailed above, it may also be a function of the unique seasonality of invasive predators and prey in the Everglades. Python feeding activity increases during the warm wet season in the Everglades (May–October) and decreases during cool dryer months (November–April; McCleery et al. 2015; Lodge 2016). Meanwhile, the abundance of cotton rats and other rodent species can double during the dry and cool season (Chapman 2019; Romañach et al. 2021). Accordingly, reductions in rodent survival during the warm wet season (Soto-Shoender et al. 2020) may be offset by compensatory reproduction during periods of reduced python activity.

Using the cotton rat as our model, we wanted to understand why some prey populations exposed to invasive predators appear to be resilient to population-level impacts. We investigated the effects of python activity and seasonality on cotton rat survival and causes of mortality (Table 1). Indicative of compensatory mortality, we expected that there would be no difference in survival between populations in areas of high and low densities of python occurrence (Dorcas et al. 2012; Reed et al. 2012; Taillie et al. 2021). In areas with high python occurrences, we expected that cotton rat mortalities would primarily be attributed to pythons, native reptiles, and birds of prey due to the loss of mammalian mesopredators from python predation (Dorcas et al. 2012; McCleery et al. 2015; Taillie et al. 2021). Lastly, we expected python predation of cotton rats to increase during the warm and wet season (McCleery et al. 2015).

Materials and methods

Study sites.

Our study sites were located within the Everglades, a 10,000 km² wetland complex located in southern Florida, United States. This large subtropical wetland (Richardson 2010) historically supported diverse plant and animal communities and includes 68 threatened or endangered species (Brown et al. 2006). The vegetation communities at our sites were predominantly grass prairies interspersed with tree islands and hardwood hammocks. The prairie communities were comprised of Sawgrass (Cladium jamaicense), Cattail (Typha spp.), and Coastal Plain Willow (Salix caroliniana). Additionally, there were several invasive plant species within our sites, including Wild Sugarcane (Saccharum spontaneum), Melaleuca (Melaleuca quinquenervia), and Burma reed (Neyraudia reynaudiana). The Everglades has a distinct, hot (30–35°C) and rainy season from mid-May to November (wet season about 114 cm), and a milder (12–25°C) dry season from December to mid-May with reduced rainfall (dry season about 38 cm; National Park Service 2021).

Burmese pythons are generalist apex predators capable of growing > 6 m in body length (Dorcas et al. 2017); they have been observed in the Everglades every year since 1995, but were likely established in the mid-1980s (Willson et al. 2011). The impacts from their introduction went undetected until the early 2000s when researchers noticed a decline in mesomammals that coincided in space and time with increasing Burmese python detections (Dorcas et al. 2012). Currently, in areas of the southern Everglades, once-common mesomammals such as marsh rabbits (Sylvilagus palustris), raccoons, opossums, river otters (Lontra canadensis), and white-tailed deer (Odocoileus virginianus) are now rare (Soto-Shoender et al. 2020; Taillie et al. 2021). These declines have been clearly linked to predation pressure from Burmese pythons (Dorcas et al. 2012; McCleery et al. 2015; Reichert et al. 2017). However, in contrast to patterns seen in mid-sized mammals, small mammals such as cotton rats still appear to be common (Dorcas et al. 2012; Hoyer et al. 2017; Gonzalez 2019; Taillie et al. 2021). At the time of this study, Burmese pythons were distributed across most of Miami-Dade, Monroe, Broward, and Collier counties (EDDMapS 2022). However, there appears to be a clear gradient of python detections decreasing from south to north across the region, with newer detections occurring in northern sites and consistent detections continuing in southern sites (Fig. 1; Hunter et al. 2019b; Burkett-Cadena et al. 2021; Taillie et al. 2021).

In order to understand how cotton rat exposure to different levels of python activity influenced their survival, we selected areas of relatively high and low levels of python occurrence densities (hereafter 'occurrence'). We categorized these areas using detections of python from the Early Detection and Distribution Mapping System (EDDMapS 2022; Fig. 1), selecting an area with a dense cluster of occurrences and an area with relatively few occurrences reported. The area of high python occurrence was a 4,100-ha strip of land (Frog Pond and Rocky Glades Public Small Game Hunting Areas; southern site on Fig. 1) adjacent to the eastern border of Everglades National Park. Pythons are regularly detected in this area (EDDMapS 2022), which had been used in previous research because of its relatively high python density (Reed et al. 2012). For an area with reduced python occurrence, we selected the Fran Reich Preserve (Fran Reich; northern site on Fig. 1), an easily accessible area with only two confirmed Burmese python sightings at the time of this study (EDDMapS 2022). Fran Reich is 647 ha and borders the Arthur R. Marshall Loxahatchee National Wildlife Refuge where mid-sized mammals were still present as of 2019 (Taillie et al. 2021). Additionally, given studies on mammal occurrence throughout the Everglades, we assumed that Fran Reich was representative of mammalian predator occurrence under low python conditions (Snow et al. 2007; Reichert et al. 2017; Taillie et al. 2021). Both sites were managed with prescribed fire and herbicide treatments to reduce invasive plant species and maintain grass-dominated vegetation structure.

Trapping.

We captured cotton rats from March to August 2020. At each site we placed 220 traps (H.B. Sherman Traps, Inc., Tallahassee, Florida) in transects of 7-10 traps spaced about 12 m apart and baited with bird seed. We selected the number of traps per transect dependent on the amount of accessible cotton rat habitat. Upon capture, we recorded sex and weight, ear-tagged (Model 1005-1L1, National Band and Tag Company, Newport, Kentucky) each individual, and fitted adult cotton rats with Holohil (Carp, Ontario, Canada) RI-2DM (7.3 g) radio transmitters (collars) with mortality sensors (McCleery et al. 2022). These transmitters lasted about 24 weeks and had an effective range of about 120 m under field conditions. We tracked more cotton rats in areas with high python occurrence (low = 34, high = 81) to increase the probability of detecting any potential influence of pythons on cotton rat survival. This study was undertaken with the approval of the University of Florida Institutional Animal Care and Use Committee (IACUC #201910726) and followed protocols from Sikes et al. (2016).

Tracking and cause-specific mortality.

We tracked cotton rats every 48 h to check for a mortality signal and homed in on the rats at least once per week to ensure that they were alive. We obtained a visual observation if we detected a mortality signal, or if the rat did not move for consecutive tracking sessions. If we were unable to detect a signal from a cotton rat, we expanded our search to accessible areas within 5 km of the last location. We continued searching for lost collars for the remainder of the study. If the animal was dead, we evaluated the carcass and the surrounding area to determine the cause of death. We classified mortalities as Burmese python, reptilian (non-python), avian, mammalian, unknown/slipped collar based on visual, photographic, or genetic evidence (Henning et al. 2008; McCleery et al. 2015). We classified avian predations based on avian fecal sprays, feathers, and removal of fur and intestines near the carcass and/or beak marks on the collar (Conner et al. 2011). For mammalian predation events, we based our classification on scat and track evidence, skin pulled off the carcass, and bite marks on the rodent and/or collar (Henning et al. 2008; Conner et al. 2011). We categorized reptilian predation events based on regurgitation and scat evidence or if the radio signal emanated from the reptile itself (McCleery et al. 2015). In consultation with the Florida Fish and Wildlife Conservation Commission, we captured and removed any Burmese pythons encountered and identified all other reptiles to species. We recorded the weight, length, and sex of reptiles when possible.

When we found radio tags in snake scat and regurgitation (from unidentified species), we collected and froze samples until they could be analyzed for the presence of Burmese python DNA. We stored eight scat and regurgitation samples in freezers (−18°C) for 2–7 months before analysis. To test whether scat and regurgitation samples could be attributed to python, we evaluated samples for the presence of python DNA using species-specific polymerase chain reaction (PCR) primers (Hunter et al. 2015). In the U.S. Geological Survey (USGS) laboratory, DNA was extracted using ~100 mg of frozen sample placed in 50-ml tubes and completely covered with lysis buffer (1 M Tris-HCl, pH 8; 0.5 M EDTA, pH 8; 5 M NaCl; 1% SDS) and Proteinase K (1 mg/ml final concentration). Samples were vortexed and placed in a shaking water bath at 55°C for 24 h. DNA was isolated from the samples using a phenol–chloroform–isoamyl (25:24:1) protocol with a polyethylene glycol precipitation following Hunter et al. (2019a). All DNA extractions were assayed by amplifying and sequencing a 115-bp DNA fragment of the mitochondrial ND4 gene using PCR using primers described by Hunter et al. (2015). Each reaction was performed using the following 25 µl mixes: 0.2 µl JumpStart Taq DNA Polymerase without MgCl₂ (Sigma), 2.5 µl of 10× PCR buffer, 4.0 µl of 25 mM magnesium chloride solution, 4.0 µl of 1.25 mM dNTP solution, 2.5 µl of 4 mg/ml bovine serum albumin, 0.625 µl of each primer (10 µM stock concentration), 7.55 µl of water, and 3 µl of template DNA. The thermocycler (SimpliAmp; Life Technologies, Carlsbad, CA) cycling parameters were a single step at 94°C for 3 min, 39 cycles of 94°C for 1 min, 55°C for 1 min, 72°C for 1 min, and a final step at 72°C for 10 min. We used gel electrophoresis to visually asses PCR products for amplification of Burmese python DNA. The PCR products were purified using ExoSAP-IT for PCR clean-up (Affymetrix, Santa Clara, California) and sequenced on an ABI 3730XL using the Big Dye Terminator v.3.1 kit (Applied Biosystems, Foster City, California). To verify the presence of Burmese python DNA (GenBank accession number KF293729), the sequences were queried in NCBI BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) to assess the percentage of identity to the target organism (Burmese python).

Inundation and seasonal data.

Pythons are semiaquatic and it is possible that inundation of dry areas during the wet season allows pythons to more easily access prey congregated on slightly elevated patches of dryer terrain (McCleery et al. 2015; Smith et al. 2021). For our high python occurrence site (Frog Pond and Rocky Glades), we assessed groundwater inundation to delineate the effects of dry versus wet conditions on cotton rat survival. We compared daily groundwater level data from the Everglades Depth Estimation Network (EDEN) database (http://sofia.usgs.gov/eden, accessed March 2021; available only for our high python occurrence site) with our field-based assessments of inundation. We considered the site to be inundated when there was water on the surface for >4 days. Our sites were completely dry until mid-May when groundwater began to accumulate. Based on visual inspection of data (Supplementary Data SD1) and on ground observations of consistent inundation, we delineated the dry conditions as prior to 20 May 2020, and the wet conditions as after that date. Individuals whose survival spanned both seasons were censored at the end of the dry season, 20 May 2020, and reentered into the study for the wet season.

Analysis.

To visualize our data and estimate median and overall survival for all cotton rats, we first used a generalized nonparametric Kaplan–Meier model (Kaplan and Meier 1958; Murray 2006). We excluded mortalities in the first 5 days of the study to account for the stress and heightened risk of mortality from capture and radio-tagging (Morris et al. 2011; McCleery et al. 2015). We also excluded the first 5 days for all cotton rat survival times to eliminate upward bias of survival estimates for those remaining in the study (Murray 2006; Morris et al. 2011; McCleery et al. 2015) and adhered to the assumption that capture and tagging do not influence future survival of individuals (Pollock et al. 1989; Benson et al. 2018). Using the date of last known radiotelemetry signal when the animal was alive, we right-censored individuals whose fates we could not determine or who disappeared from the study areas (Israelsen et al. 2020; Baek et al. 2021). To determine if the proportion of right-censored individuals, we used a chi-squared test.

To determine whether variation in cotton rat survival could be explained by differences in our categorical variables (seasonal inundation and python occurrence), we used a Cox proportional hazards regression model (Therneau 2021). We also used this model to evaluate the potential for survival differences between sexes. The Cox proportional hazards model integrates staggered entry and right-censored data (Kaplan and Meier 1958; Fieberg and DelGiudice 2009; Goel et al. 2010) with the ability to evaluate multiple covariates (e.g., inundation and python occurrence). Using a semiparametric framework, the Cox proportional hazards model is a product of a nonparametric baseline hazard (i.e., instantaneous risk of death) for all individuals with effects of covariates modeled parametrically on the baseline hazard (Fieberg and DelGiudice 2009). To determine if the covariates helped explain variation in the survival of cotton rats, we assessed the hazard ratios (HRs). The HR was calculated using the hazard function of the two different groups (e.g., high and low python occurrence) and comparing them (Cox 1972; Singh and Mukhopadhyay 2011) to produce a ratio or effect size (Spruance et al. 2004). HRs can exceed 1 since they are the expected number of events per one unit of time. A ratio of 1 suggests the covariate has no explanatory power. HRs greater than 1 indicate an increase in risk from a predictor and those lower than 1 suggest a predictor variable reduces risk. First, we used our full set of data to determine if survival rates of cotton rats varied by sex. We pooled the data for further analysis if we found no difference. We conducted further analysis separately for males and females if survival varied between the sexes. Next, we developed an additive model with python occurrence (high and low) and inundation (wet and dry) to determine whether our data were suggestive of patterns of: (1) additive mortality leading to decreased survival; (2) increased survival from prey release (Letnic et al. 2009; Gordon et al. 2015); (3) comparable survival from compensatory mortality; or (4) seasonal shifts in survival that align with invasive predator activity patterns (Table 1).

Lastly, to determine whether any patterns seen in our data were a function of the increased sample size at the high python activity site, we ran a bootstrap resample function (DePatta Pillar 1998) to create 999 bootstrapped samples separately for high and low python activity data sets. We resampled cotton rat individuals to create resampled data sets in which the number of individuals from high and low python occurrence areas were equal (n_high = n_low = 34). We evaluated a single covariate Cox proportional hazards model using the data sets with n observations to determine if python occurrence explained variation in cotton rat survival when there were equal sample sizes.

We calculated HR using the cox.ph function in the survival package on the R platform (version 3.3.1; Therneau 2021) and considered HR of covariates with 95% intervals that did not overlap 1 to be significant. We also assessed the significance of P-values for each HR using the likelihood ratio test, which uses log-likelihood estimates to determine if the HR is significantly different from 1 (Kleinbaum and Klein 2012). For each model, we tested the proportional hazards assumption (HRs are constant) by plotting the Schoenfeld residuals and assessed model fit using the cox.zph function in the survival R package (Therneau 2021). Models do not violate the proportional hazards assumptions when Schoenfeld residual P-values were equal to or exceeded 0.05 (Hickey et al. 2019).

Next, to determine if patterns of survival were consistent with additive mortality, release from native predators, or compensatory mortality, we compared the causes of cotton rat mortality between areas with different levels of python occurrence. We did this using a nonparametric cumulative incidence function estimator as described by Heisey and Patterson (2006) and Fine and Gray (1999). We compared the differences in mortality risk based on python occurrence (high and low) using Gray’s test (Gray 1988; Zhang et al. 2008). Gray’s test determines if differences exist in cumulative incidence between groups for a particular event (e.g., death from pythons; Dignam and Kocherginsky 2008). Broadly, if areas of high python occurrence showed reduced survival and marked python predation, we considered this consistent with additive mortality. If areas of high python activity showed increased survival and measurable declines in predation from native mammals, we considered this to be consistent with prey release. If we found no difference in survival and a shift in causes of predation, we considered this to be consistent with compensatory mortality. Finally, if we found reduced survival after inundation, we considered this consistent with the prediction that environmental changes would make prey more accessible to pythons (Table 1). We conducted all survival and cause-specific mortality analyses in Program R (R Core Team 2021, version 4.0.3, https://www.R-project.org/) using ggplot2 (Wickham 2016), survminer (Kassambara et al. 2011), survival (Therneau 2021), and cmprsk (Gray 2020) packages.

Results

Survival.

We captured 126 total cotton rats between March and August 2020 and tracked them until October 2020. We analyzed the survival of 115 cotton rats (low python occurrence n = 34, high python occurrence n = 81) once we excluded 11 individuals who died or lost collars within 5 days of entering the study. Median survival time for all cotton rats in high and low python activity sites was 63 days, 95% CI 44–94 days (Supplementary Data SD2).

Using data from all 115 cotton rats, we found no indication that the sex of cotton rats explained variation in survival rates (HR = 0.87, 95% CI = 0.54–1.39, P = 0.60; Fig. 2A). Finding no difference, we pooled the sexes and found no influence of wet season inundation (HR = 1.07, 95% CI = 0.61–1.88, P = 0.80; Fig. 2B) which started on 20 May 2020 (Supplementary Data SD1). Additionally, we found no influence of python occurrence (HR = 1.32, 95% CI = 0.77–2.26, P = 0.30; Fig. 2C) on cotton rat survival. Using equal treatment sample sizes (n = 34) after running a bootstrap resample function on the data from the high and low python occurrence sites, we found no evidence of python occurrence influencing cotton rat survival (HR = 1.11, 95% CI = 0.59–2.10, P = 0.70; Fig. 2D). The proportional hazards assumption was met for all Cox proportional hazards models with P > 0.05.

Cause-specific mortality.

Of the 115 cotton rats in our study (low python occurrence n = 34, high occurrence n = 81), we classified the cause-specific mortality for 51 cotton rats—remaining cotton rats were either censored at the end of the study (n = 17), determined to have slipped their collar or cause of death could not be determined (n = 44), or their signal was lost (n = 3). Based on a chi-squared test we found no statistical difference between our sites in the proportion of cotton rats classified as cause of death undetermined or signal lost (low python = 11/34 [32%], high python = 36/81 [44%], X² = 0.23, P = 0.76). From the genetic identification analysis, two scat samples were positively confirmed as containing Burmese python DNA. The other samples likely originated from different species.

Of the 51 classified mortalities across all sites, most cotton rats were consumed by avian predators (45.1%; Table 2), followed by non-python reptiles (25.5% of known fates; Table 2). On the high python occurrence site, birds and non-python reptiles accounted for 48.6% and 24.3%, respectively, of classified mortalities for cotton rats, followed by Burmese pythons (16.2%) and mammalian predators (10.8%; Table 2). On the low python occurrence site, cotton rats were depredated by birds (35.7%), mammals (35.7%), and non-python reptiles (28.6%; Table 2). No cotton rats were consumed by pythons at the low python occurrence site.

Comparing cumulative incidence between high and low python occurrence sites for specific mortality risk (Fig. 3), we found no difference in the risk from avian predators (test statistic = 0.24, P = 0.62), non-python reptilian predators (test statistic = 0.16, P = 0.69), and pythons (test statistic = 1.52, P = 0.22). However, there was a greater risk of predation from mammals (test statistic = 4.16, P = 0.04) on our low python occurrence site compared to the high python occurrence site (Fig. 3).

Discussion

While most mammal populations in the Everglades are sensitive to invasive pythons (Dorcas et al. 2012; McCleery et al. 2015; Burkett-Cadena et al. 2021), our study helps explain the shift to communities dominated by cotton rats and other small mammals (Soto-Shoender et al. 2020; Taillie et al. 2021). Our finding of little difference in cotton rat rates of survival (Fig. 2), despite differences in the cause of mortality (Fig. 3), suggests that the apparent resilience of cotton rats to an invasive predator was a function of compensatory mortality (Table 1). Specifically, supporting our predictions, we found that birds and reptiles—including pythons—were compensating for the risks cotton rats faced from mammalian predators in areas with high python occurrence. This type of compensatory mortality has been shown for highly fecund prey species such as cotton rats (Creel 2011; Krebs 2013), where the substitution of predators can maintain overall density-dependent mortality risk (Creel 2011; Gordon et al. 2015). In fact, cotton rats have shown similar responses when their mesocarnivore predators have been removed or reduced in other systems (Wiegert 1972; Conner et al. 2011). This helps explain why, in contrast to other mammals, cotton rats and other rodents appear to be relatively resilient to invasive pythons (Dorcas et al. 2012; McCleery et al. 2015; Hoyer et al. 2019; Soto-Shoender et al. 2020).

Similar to previous research in Georgia, South Carolina, and Texas (Schnell 1968; Wiegert 1972; Kincaid and Cameron 1982; Conner et al. 2011), we found that the greatest mortality risk for cotton rats came from avian predators (Table 1; Supplementary Data SD3). This was particularly true on the site with more python activity where risk from mammalian carnivores was reduced (Fig. 3). Consistent with our prediction and a growing body of evidence from the Everglades, we found that cotton rats had little risk of predation from mammalian carnivores on our site with high python occurrence (Dorcas et al. 2012; McCleery et al. 2015; Soto-Shoender et al. 2020; Taillie et al. 2021; Fig. 3). In fact, on the site with reduced python occurrence density, the cumulative risk to cotton rats of predation from mammals was more than double that of risk of mammalian predation on the high python occurrence site (Fig. 3). Researchers have questioned whether the punctuated declines in mammalian predators in the Everglades was a function of direct mortality from pythons or the loss of their prey base (McCleery et al. 2015; Reichert et al. 2017). Our findings, along with the work of others (Dorcas et al. 2012; Hoyer et al. 2019; Soto-Shoender et al. 2020), suggest that for mammalian carnivores with diets dominated by rodents, including bobcats (Godbois et al. 2003), gray foxes (Urocyon cinereoargenteus; Deuel et al. 2017), and coyotes (Canis latrans; Grigione et al. 2011; Swingen et al. 2015) there is an available prey base. This in turn suggests that direct predation from pythons may be limiting mammalian predators (Dorcas et al. 2012; Reed et al. 2012).

While pythons are known to increase activity during the warm and wet season when the Everglades becomes saturated (McCleery et al. 2015; Mutascio et al. 2018), we found no evidence that seasonal inundation influenced cotton rat survival. Cotton rats usually experience higher mortality when food and cover resources diminish in the winter (Morris et al. 2011). In the Everglades, we would expect these resources to dimmish with seasonal inundation. Nonetheless, in this system cotton rats appear to be abundant year-round (Smith and Vrieze 1979; Soto-Shoender et al. 2020) and changes in inundation did not appear to have a marked influence on cotton rat survival.

One consequence of cotton rat resilience is that they may provide a consistent prey base for pythons (Soto-Shoender et al. 2020). In particular, gape-limited juvenile pythons require smaller prey like cotton rats (Shine 1991; Dorcas et al. 2017). Once they are adults, Burmese pythons can also consume larger prey including deer, coyotes, alligators (Alligator mississippiensis), wading birds, etc. (Snow et al. 2007; Orzechowski et al. 2019). If cotton rats provide a consistent food source, pythons may become an increased risk for sensitive or endangered species, which may be unable to withstand or recover from python-induced losses (Dorcas et al. 2017). Additionally, cotton rats are reservoirs for zoonotic viruses such as the Everglades Virus (EVEV) and hantaviruses (Coffey et al. 2004; Billings et al. 2010). With reductions in mammalian diversity, cotton rats and other rodents can amplify the prevalence of these diseases (Coffey et al. 2004; Hoyer et al. 2017; Burkett-Cadena et al. 2021). Alternatively, one benefit of a resilient cotton rat population is that they may provide a prey base for mesopredators if they eventually recover.

While our research suggests clear patterns of change in the cotton rat predator community but only negligible changes in survival with increased python occurrence, there are some limits to our inference. Specifically, with only a few accessible areas in the Everglades that have reduced python occurrence, our sample size in the low python occurrence site was smaller than that in the high python occurrence site. Additionally, methods for localized python population estimates have yet to be developed (Bonneau et al. 2016), so we used opportunistic presence data from EDDMapS to approximate differences in pythons across space. Over more than 15 years our low python occurrence site had markedly fewer reported detections of pythons compared with our high occurrence site (2 vs. >300 detections; EDDMapS, accessed January 2021).

Another constraint of our study was that it was less than a full year and thus did not include annual variation. Nonetheless, seasonality in the Everglades varies and our study captured a majority of the wet season of 2020 (May–October; Supplementary Data SD1) when pythons are most active (McCleery et al. 2015; Mutascio et al. 2018) and most likely to prey on cotton rats. Also, we only included adult cotton rats in our study due to radio-tag weight restrictions. Therefore, we do not have a background understanding of cotton rat population demographics for all age classes. However, our results are supported by recent research showing that cotton rat populations are abundant throughout the Everglades (Chapman 2019; Gonzalez 2019; Romañach et al. 2021). Finally, for a number of the cotton rats in the study we could not attribute a cause of death, or we lost signal. However, there was not a statistical difference between the portion of unknown and lost individuals on our sites and we attributed any minor differences in the proportion of unclassified animals to differences in accessibility between our sites. Accordingly, we have no reason to believe that our data were biased.

Invasive predators such as Burmese pythons are restructuring ecosystems around the world (Doherty et al. 2015, 2016). Their negative impacts are most often due to additive mortality of native prey, yet not all prey species experience the same consequences. Some prey experience release (Gordon et al. 2015) and others can persist relatively unperturbed due to density-dependent compensatory responses (Abrams and Rowe 1996; Rodriguez 2006; Abrams 2009; Osmond et al. 2017). As observed in this study, prey resilience can be a function of the compensatory responses that are most common in smaller mammals. In ecosystems devastated by invasive predators, this phenomenon can produce small mammal-dominated communities (Fritts and Rodda 1998; Hanna and Cardillo 2014). In the Everglades, cotton rats and other small mammals are clearly responding differently than large- and medium-sized mammals (Soto-Shoender et al. 2020). This pattern is likely to substantially impair restoration efforts in the Everglades as ecosystem processes (e.g., nutrient cycling, scavenging, seed dispersal) provided by once-common mammalian predators and herbivores are unlikely to be compensated for by rodents (McCleery et al. 2015). Our results also reflect a broader global trend of declines in larger species due to greater susceptibility to altered mortality risk from human-induced changes (Estes et al. 2011; Malhi et al. 2016; Loggins et al. 2019; Epperly et al. 2021). The resulting communities often have shorter food chains, coupled with altered disease dynamics and ecosystem functions (Mooney and Cleland 2001; Forys and Allen 2002; Hoyer et al. 2019).

Acknowledgments

We thank our collaborators and colleagues with the Florida Fish and Wildlife Conservation Commission, especially S. Funck, M. Spencer, J. Travers, and Officers J. Hart II and C. Ryan for their field and logistical support. We also thank our colleagues with the South Florida Water Management District, especially C. Stylianos, M. Porras, M. Kirkland, D. Cotter, J. Laing, and B. Garrett for permits, land access, and other logistical support. We are grateful to J. Hartfelder and L. Rittenburg for their dedicated assistance in the field, and M. Renshaw in the USGS laboratory. Funding for this project was provided by the Florida Fish and Wildlife Conservation Commission (grant #13416-A3044) as well as the United States Geological Survey (grant #G19AC00432). Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the United States Government.

Supplementary Data

Supplementary data are available at Journal of Mammalogy online.

Supplementary Data SD1.—Daily water level (by month) for high python activity site from mid-March through early October 2020. Water level data from 5 gages on or near the high python occurrence study site. Data accessed from DBHydro environmental database (https://www.sfwmd.gov/science-data/dbhydro, accessed April 2021). Vertical line denotes the 20 May date that partitions dates into wet and dry season.

Supplementary Data SD2.—150-day Kaplan–Meier survival curve with 95% confidence interval for all cotton rats captured in high and low python occurrence sites combined (n = 115) from March to July 2020. Intersecting dashed lines denote median survival time (63 days) for all individuals.

Supplementary Data SD3.—Cumulative cause-specific predation risk for all (n = 51) cotton rats with classified mortalities. Data collected from March to October 2020.

Literature Cited