The increasing frequency and severity of human-caused fires likely have deleterious effects on species distribution and persistence. In 2020, megafires in the Brazilian Pantanal burned 43% of the biome’s unburned area and resulted in mass mortality of wildlife. We investigated changes in habitat use or occupancy for an assemblage of eight mammal species in Serra do Amolar, Brazil, following the 2020 fires using a pre- and post-fire camera trap dataset. Additionally, we estimated density for two naturally marked species, jaguars Panthera onca and ocelots Leopardus pardalis. Of the eight species, six (ocelots, collared peccaries Dicotyles tajacu, giant armadillos Priodontes maximus, Azara’s agouti Dasyprocta azarae, red brocket deer Mazama americana, and tapirs Tapirus terrestris) had declining occupancy following fires, and one had stable habitat use (pumas Puma concolor). Giant armadillo experienced the most precipitous decline in occupancy from 0.431 ± 0.171 to 0.077 ± 0.044 after the fires. Jaguars were the only species with increasing habitat use, from 0.393 ± 0.127 to 0.753 ± 0.085. Jaguar density remained stable across years (2.8 ± 1.3, 3.7 ± 1.3, 2.6 ± 0.85 / 100km²), while ocelot density increased from 13.9 ± 3.2 to 16.1 ± 5.2 / 100km². However, the low number of both jaguars and ocelots recaptured after the fire period suggests that immigration may have sustained the population. Our results indicate that the megafires will have significant consequences for species occupancy and fitness in fire affected areas. The scale of megafires may inhibit successful recolonization, thus wider studies are needed to investigate population trends.

We used camera traps (Bushnell 119876, Panthera V4 and Cuddeback 1279) to survey the study area in December 2019 (session 1; year 1 – pre-fires) and December 2020 (session 2, year 2 – 2 months post-fires). Due to logistical constraints, we installed cameras in February 2022 (session 3, year 3 – 15 months post-fires) for an average duration of 53 trap nights (range 1-136, see SI_1 for complete details). All three surveys took place in the rainy season. Thirty-five stations were active in session 1, 43 stations in session 2, and 31 stations in session 3. Cameras were placed at a distance of 1.5 ± .5 km between stations and were located in different land covers (primary, secondary and gallery forest, savannah). Minimum convex polygons for each survey were 189.68 km2 in year 1, 272.26 km2 in year 2, and 245.95 km2 in year 3. We placed double stations to enable photographing both sides of each passing individual, thus enabling identification for naturally marked species like jaguars and ocelots. Each sampling station had 24-hour motion-triggered camera operation with a period of 30 seconds between photograph triggers. Geographic coordinates, camera serial number, date and time of camera installation, canopy cover, habitat, and whether the camera was on or off trail were recorded. 

Our survey design complied with methodological assumptions to estimate jaguar (Foster et al., 2020; Tobler et al., 2013) and ocelot densities (Boron et al., 2021; Satter, Augustine, Harmsen, Foster, Sanchez, et al., 2019; Wolff et al., 2019), and we kept a discrete distance between stations to obtain data for the wider mammal community (Boron et al., 2021; de Martins et al., 2006; Rovero et al., 2020; Rovero & Ahumada, 2017). Our survey was limited to less than 100 days per year, and fulfills overall capture-recapture model assumptions: a) the population needs to be considered closed and stable, and b) all individuals should have a chance of being captured (Otis et al., 1978; White, 1982).

Covariate selection and extraction

We selected a set of covariates to test our hypotheses related to pre- and post-fire habitat use/occupancy as well as density. Covariates were Normalized Difference Vegetation Index (NDVI), often used to assess habitat quality for mammals (Pettorelli et al., 2005; White et al., 2022), area burned derived from Normalized Burn Ratio (ΔNBR) which measures fire severity (Escuin et al., 2008), and distance from water (Boron et al., 2019). We additionally included effort as the total of trap nights per station; year, included to account for differences related to time variation as field staff and camera type on p (Gutiérrez-González et al., 2015; Kotze et al., 2012); and whether the camera station was on a trail or not for the probability of detection (p). Year or session was also used as a way to account for the heterogeneity of the detection probability, like seasonal activity of species and the possible loss of camera quality (Kotze et al., 2012; MacKenzie et al., 2003; Tobler et al., 2015).

NDVI was obtained for each study session from Copernicus-Sentinel-II sensors via Google Earth Engine (code here). NDVI calculates vegetation greenness on a normalized scale with denser vegetation approaching one and barren areas or water bodies closer to a value of zero (Pettorelli et al., 2005). Annual NDVI rasters were obtained on days with less than 10 percent cloud cover during the period of one month before camera installation with a grain size of 10 meters. We then extracted the mean NDVI value for a 500-meter buffer around each station.

We calculated the area burned (AB) as the area within a 500-meter buffer of each camera station that presented moderate-low severity or higher according to ΔNBR. For AB, we included only ΔNBR values that represent moderate-low burn severity and higher (ΔNBR = 270+) (Keeley, 2009; Key & Benson, 2006) in order to differentiate from areas that may have had low affectation from the fires or that may have presented false-positive values where fires may have burned due to the lack of an NBR system for the region. We obtained surface water data from MapBiomas (https://brasil.mapbiomas.org) and calculated the Euclidean distance of each camera station from surface water. All geoprocessing was conducted in ArcMap Desktop 10.8 (ESRI Inc., 2020). We used Spearman’s correlation test to check for highly correlated covariates (>0.6) with the function ggcorr in the package “GGally v2.1.2” (Schloerke et al., 2022) in Program R v 4.2.2 (R Core Team, 2022). As covariates AB and ΔNBR were correlated (>|0.6|), we fit two global models , each including one of this covariates and selected the best model using Akaike Information Criteria corrected (AICc) for small sample sizes  (Burnham & Anderson, 2002). The model with AB covariate performed better than ΔNBR for most species and thus was used in the analysis. 

Dynamic occupancy

We determined the habitat use or initial occupancy probability of eight mammal species in the study area: jaguars (Panthera onca), ocelots (Leopardus pardalis), pumas (Puma concolor), giant armadillo (Priodontes maximus), lowland tapir (Tapirus terrestris), red brocket deer (Mazama americana), collared peccaries (Dicotyles tajacu), and Azara’s agouti (Dasyprocta azarae). We determined initial occupancy probability when we could assume closure (individual’s home range is less than the radius between camera trap stations) between camera trap locations (MacKenzie et al., 2002), for ocelots (Crawshaw & Quigley, 1989), red brocket deer (Varela et al., 2010), Azara’s agouti (Cid et al., 2013) and giant armadillo (Desbiez et al., 2020). And determined habitat use for jaguars (Kantek et al., 2021; Soisalo & Cavalcanti, 2006), pumas (Silveira, 2004), collared peccaries (Desbiez et al., 2009) and tapirs (Medici et al., 2022), whose home range surpassed the distance (1.5 km) between stations. Detection histories were created for each species, grouping camera data into a 7-21 days survey occasions based on the results of goodness of fit (GOF) tests (MacKenzie & Bailey, 2004) (SI_3,4). Dynamic occupancy models (DOM) estimate the probability of occupancy and detection and are particularly useful for monitoring changes in occupancy status over time (MacKenzie et al., 2018), allowing us to detect if certain variables were influencing the colonization (Ɣ) and extinction (Ɛ) trends. We scaled covariates before analysis for interpretability. We used “unmarked” package v 1.2.5 (Fiske & Chandler, 2011) in Program R v 4.2.2 (R Core Team, 2022) for all occupancy analysis.

The parameters used in DOM were Ψ = initial probability of a site being occupied; p = probability of a species being detected if it is present, Ɣ = probability of a new area to pass from unoccupied to occupied (or to unused to used) in the next year, Ɛ = probability that a species stops occupying an area, or to pass from used to unused.

We fit models for species individually, and selected models according to AICc (Burnham & Anderson, 2002). We used a stepwise method (Doherty et al., 2012) for model selection. We first fit models for detection (p) with all other parameters constant. We included survey effort, whether cameras were on a trail, and year as covariates for detection, and selected the best detection model based on AICc value. We proceeded with this best detection model and subsequently fit models using covariates describing occupancy, colonization, and, finally, extinction. We included NDVI and distance to water as covariates for occupancy, whereas area burned was applied to colonization and extinction.

We considered there to be satisfactory statistical evidence for an effect if the 95% confidence interval of logit scale coefficient estimates did not include zero (Muff et al., 2022). The β estimates were back-transformed to obtain model parameter estimates (MacKenzie & Bailey, 2004). We tested model fit by using a parametric bootstrap GOF test based on Pearson’s X2 where p>0.05 indicates adequate model fit (Fiske & Chandler, 2011) (SI_3). Finally, we derived annual probability for each year, and calculated standard errors for the derived values using a bootstrap method (Kéry & Chandler, 2012). To assess whether differences in occupancy were statistically significant, we used a nonparametric bootstrap with 1000 iterations and calculated the difference between 95% confidence intervals between years one and three. If the difference in confidence intervals did not contain zero, we considered it significant (Kéry & Royle, 2021).

Density

Jaguar and ocelot individuals were visually identified by their spots patterns and sex based on external genitalia (Boron et al., 2016). Individuals of unidentified sex were classified as unknown (SI_1). We generated three files to input in “secr” in R (M. G. Efford, 2020) for density estimation a) a capture history per species per session with trap number, individual ID, date and time occasion number of each record; b) a trap deployment file consisting of all traps locations and a binary record file of activity per occasion (defined here as a 24 hours period/day); c) a mask file to represent the state space habitat area, at a 500m spacing and within a buffer of 20 km using the suggest.buffer that corresponds to ~4 (Efford, 2018) that also matched the largest buffer suggested for jaguars previously (Greenspan et al., 2020; Noss et al., 2012). We used the half-normal (HN) detection function, which considers that the probability of capture (p) of any individual (i) decreases with the distance (d) from the activity center such that: pij = g0exp(−dij2/2σ2), g0 represents the probability of capture or detection when a trap (j) is located right at the center of the home range; and sigma (σ) is a proxy parameter related to home range size (Efford, 2004). Felid species have varying home ranges size according to sex (Massara et al., 2015; Satter, Augustine, Harmsen, Foster, & Kelly, 2019); however, as sex was not identified for all individuals, this was not included in the models. Felids are also known to use unpaved roads to varying degrees and jaguars may be more detectable on them (Sollmann et al., 2011). We therefore enabled detection (g0) to vary by whether the station was located on a trail, while sigma (σ) varied according to the different sessions to assess if the value is constant or varies across years. We additionally considered distance to water, NDVI and ΔNBR as covariates for density and selected models using AICc (Burnham & Anderson, 2002).