Of the remaining 20% of the world's tropical forests still considered intact, many present cryptic and pervasive forms of degradation, particularly overhunting (Benítez-López et al. 2017). Even protected areas are threatened globally, with only 10% still free of intense human pressure (Jones et al. 2018). The rising number of people and their activities have imposed unprecedented threats on the vast majority of species. In particular, mammals have lost more than 50% of their geographic ranges, with most population extinctions concentrated in areas with high human densities, or with severe human impacts, such as intensive agriculture, grazing, and hunting (Ceballos and Ehrlich 2002). Large terrestrial carnivores have suffered significant population declines, and many are threatened with extinction. Their declines have been driven largely by the loss or degradation of their habitat and prey base, persecution by humans, and overexploitation (Ceballos and Ehrlich 2002).

Although biodiversity conservation generally benefits from protection of natural habitat across large contiguous areas (Primack 2008), many landscapes are already fragmented by increasing agriculture and/or livestock ranching (Reid et al. 2010). Small protected areas can remain valuable preserves of biodiversity in highly altered landscapes (Ricketts et al. 2005). They supplement larger tracts of natural habitat and may persist as high-quality remnants (Schwartz 1999). They usually cannot suffice alone to satisfy conservation goals, however, and are much more likely to suffer high levels of anthropogenic pressure from surrounding habitats (Hobbs 1993; Jones et al. 2018).

Large carnivores are particularly affected by all these issues. Each individual requires extensive suitable habitat within its home range (Newmark 1985), making their populations very prone to the effects of landscape fragmentation (Castilho et al. 2015), and consequently at risk of regional extinction due to a reduced potential for dispersal. Population persistence requires cooperative management of the public and private lands lying between protected-area boundaries that rarely encompass the biotic boundaries of large-carnivore populations (Newmark 1995).

Jaguars (Panthera onca) and pumas (Puma concolor) are the largest felids in the Neotropics. They coexist across almost the entire regional distribution of the jaguar (Sunquist and Sunquist 2014). Both felids remain susceptible to global extinction (IUCN Red List: jaguar, Near Threatened; puma, Least Concern, but populations decreasing; Nielsen et al. 2015; Quigley et al. 2017). Tropical forests in the Yucatán Peninsula contain one of the largest areas of jaguar habitat in Mexico. Suitable habitat, however, has suffered significant transformation and fragmentation by agricultural activities and by the expansion of the cattle frontier, especially during the 20th century (González-Iturbe et al. 2002; Conde et al. 2010). Two main initiatives have aimed to improve habitat connectivity for these large felids in this region. The Mesoamerican Biological Corridor, created in 1997, promotes transboundary landscape linkages between protected areas, to facilitate movement of larger mammals through all countries in Central America and five states of southern Mexico (Sarkar et al. 2009). Extending further south, Panthera's Jaguar Corridor Initiative is designed to connect and protect core jaguar populations in human landscapes, and to preserve their genetic integrity, in order to protect jaguars across their entire distributional range (from northern Mexico to Argentina; Panthera 2020). Despite these significant efforts, most of the habitat of the Selva Maya remains highly fragmented (Conde et al. 2010). Almost all of the ca. 50 protected areas found in this region (CONANP 2020) remain too small to safeguard the survival of large felids (Crooks 2002; Rodríguez-Soto et al. 2011). They have potential, nevertheless, to function as stepping stones for species requiring large home ranges (Boitani et al. 2007). Maintenance and restoration of connectivity among these protected areas is needed to offset the negative impacts of regional fragmentation and their small size (Crooks 2002).

Here, we describe a case study in the utility of small reserves designated for the conservation of large felids whose home ranges exceed reserve boundaries. We aimed to estimate the state of the regional jaguar population, the distribution of large felids, and their prey preferences, in and between two small private reserves of the Northern Yucatán Peninsula (NYP) of Mexico: El Edén Ecological Reserve (henceforth “EEER”) covering 25 km²; and El Zapotal Conservation Area (“EZCA”) covering 43 km². Our main goal was to evaluate the utility of small, protected areas in sustaining populations of top-predator felids in NYP. We had three objectives for addressing this goal: (1) to estimate jaguar population densities with the first application in NYP of a maximum-likelihood approach to spatially explicit and sex-dependent capture–recapture models applied to systematic camera trapping; (2) to test for spatial segregation among potentially competing jaguars and pumas; and (3) to identify prey within large felid scats, to test for overlap or segregation in the trophic niches of jaguars and pumas.

Materials and Methods

Study area.—Our study area (Fig. 1) comprised two small private reserves: (1) EEER (21°14′14″ to 21°09′48″N; 87°12′22″ to 87°09′08″W), containing medium-statured, semideciduous forest and secondary-growth semideciduous forest (Schultz 2003); and (2) EZCA (21°26′10″ to 21°19′41″N; 87°41′12″ to 87°32′51″W), containing semideciduous forest and secondary-growth semideciduous forest (Faller-Menéndez et al. 2005). Both reserves prohibit any habitation and all exploitation of natural resources. The study area also included an intermediate area (henceforth “IA”; 21°25′11″N to 87°28′12″W) located between EEER and EZCA, with medium-statured semideciduous forest and secondary-growth semideciduous forest. IA contained interspersed rural habitation, and “ejidos,” which are aggregations of land parcels entrusted by the Mexican government to rural communities for tenure as farmland. EEER is located close to the Southeastern corner of Yum Balam Flora and Fauna Protection Area (52,308 ha), and EZCA is close to the Southern limit of Ría Lagartos Biosphere Reserve (60,348 ha; CONANP 2020).

Fig. 1.

Camera traps deployed in Northern Yucatán Peninsula (NYP): El Edén Ecological Reserve (EEER), El Zapotal Conservation Area (EZCA), and intermediate area (IA), showing limits of EEER, EZCA, Yum Balam, Ría Lagartos, and study area (black rectangle in inset map). Black crosses represent large (>400 people) rural communities. The predominant habitat is forest, with increasing frequency of agricultural patches each <100 ha interspersed by small rural communities toward the southwest corner, to the south and west of EZCA.

Camera-trap sampling design.—We conducted camera-trap surveys during 2- to 5-month periods in 2015 and 2016. For EEER and EZCA, the sampling design followed the CENJAGUAR design (Chávez et al. 2011) for jaguars and prey > 6.8 kg. We deployed camera traps with a minimum separation of 1 km from any neighboring camera (see  Supplementary Data SD1 (gyac090_suppl_supplementary_data_sd1.docx) for a detailed description of CENJAGUAR). Camera traps were deployed in IA to bridge between EEER and EZCA, and to increase total sampling area of NYP (Fig. 1).

We used a combination of Cuddeback and Bushnell passive infrared digital cameras, activated by a heat-motion sensor. The different makes were interspersed across the survey area to avoid biases from minor differences in sensitivity (Quiroga et al. 2016). Each camera-trap station was set to function 24 h/day. Cameras were checked every 3–4 weeks. A total survey effort of 11,824 trap days was achieved in NYP (EEER, EZCA, and IA) from 2015 to 2016 ( Supplementary Data SD2, Table S1 (gyac090_suppl_supplementary_tables.docx)). Pumas could not be identified to individual due to the absence of coat spots or markings, and thus were excluded from analysis of density. Jaguar photographs that could not be reliably identified were also excluded from density analyses (27% of 359 jaguar images).

Models and methods used to estimate jaguar densities.—We employed two types of spatially explicit capture–recapture (SECR) models: a maximum-likelihood approach (R package secr; Efford et al. 2009); and a Bayesian approach of Markov and Monte Carlo simulations (R package SPACECAP; Gopalaswamy et al. 2012). We estimated jaguar density during 2015 and 2016 in EEER and EZCA, and across NYP (pooling EEER, EZCA, and IA). We restricted our camera-trap records to periods of ≤90 camera days, to minimize risk of violating the assumption of demographic closure, on which conventional SECR models depend for accurate density estimations (Otis et al. 1978; Royle et al. 2014). In order to obtain an insight into the closure state of our populations, we ran the population closure test developed by Otis et al. (1978;  Supplementary Data SD2, Table S2 (gyac090_suppl_supplementary_tables.docx)) for nonspatial capture–recapture models. To calibrate the above population estimates of jaguar densities in NYP against estimates that are more likely to violate the assumption of demographic closure, we conducted a second density estimation in NYP using all our camera-trap sampling records, taken from April to October (n = 179 days) in 2015 and from March to August (n = 144 days) in 2016. The magnitude of difference in density estimates from long- and short-period sessions should be informative about the influence of closure on our density estimations.

Model inputs and assumptions.—The maximum-likelihood and Bayesian models require the capture–recapture history of all identified individuals at each camera-trap station, in combination with the spatial distribution of the captures and recaptures (Noss et al. 2012). They assume: (1) closed model capture–recapture sampling (conventional SECR analysis); (2) independent activity centers for captured individuals; (3) fixed locations for activity centers during the sampling period; (4) a declining probability of detecting an individual at a camera station with increasing distance of the trap from the individual's activity center; and (5) independent capture events (Otis et al. 1978; Efford et al. 2009; Foster and Harmsen 2012; Gopalaswamy et al. 2012, 2014).

The maximum-likelihood approach provides a larger choice of possible models than the Bayesian approach. It allows population density estimations with three different models: (i) not accounting for sex differences in detection probability; (ii) incorporating sex-specific detection probabilities; and (iii) excluding all potential activity centers falling outside suitable habitat (Zimmermann and Foresti 2016). The Bayesian approach, on the other hand, deals well with issues presented by individual heterogeneity in capture probabilities. It also provides nonasymptotic inferences, which are more appropriate for the small samples of capture data that are typical of photo-capture studies (Royle et al. 2009; Gopalaswamy et al. 2012, 2014).

For the maximum-likelihood approach, we started by generating a series of alternative areas of influence (masks) around the camera-trap arrays (Zimmermann and Foresti 2016), using increasing buffer widths of 1, 3, 5, 7, 9, 11, 13, 15, 17, 20, 25, and 30 km. We then ran a series of null nonspatial models (each one associated with each buffer), assuming constant values for the baseline encounter probability (g₀; i.e., encounter probability when the distance between the activity center of an individual and the camera trap is zero) and the spatial scale parameter (σ). Following Royle et al. (2014), we then selected the best model by choosing a buffer with a width of 2–3σ, which guarantees that individuals outside the area of influence have zero detection probability by the camera-trap array during the sampling period. Using the mask associated with the best model, we ran a mixture of alternative SECR equivalents to conventional capture–recapture models, in order to explore their effect on g₀, while keeping σ constant (Zimmermann and Foresti 2016). Selected models followed Zimmermann and Foresti (2016), and included a series of sources of variation (and various combinations of them): (1) M0—is the most basic capture–recapture model, where encounter probability is strictly constant for all individuals and occasions; (2) Mt, MT—all individuals have equal capture probabilities for a particular trapping occasion, but probabilities can differ among occasions (Mt), or there is a linear trend in baseline encounter probability over occasions (MT); (3) MB, Mb—individuals become either “trap happy” or “trap shy,” following either a learned response to the first capture (Mb) or a transient response to a preceding capture (MB); (4) Mk, MK—capture probabilities are specific to a particular site (local response: Mk), or to all sites (global response: MK); and (5) Mh2—a “finite mixture model,” in which individuals have heterogeneous capture probabilities across two undefined latent classes, and the likelihood uses a weighted sum over the classes. Our data set was too sparse to fit models Mt and Mbt, which were not included in any subsequent analysis. Models used either the Newton–Raphson algorithm (default) or the Nelder–Mead algorithm (less prone to settling in local minima; Otis et al. 1978; Royle et al. 2014; Zimmermann and Foresti 2016; Efford 2018a, 2018b). Model selection for estimating jaguar densities used Akaike's Information Criterion (AIC; Burnham and Anderson 2002; Royle et al. 2014). The best model had the lowest AIC value, and alternative models with ΔAIC < 2 also had strong support, where ΔAIC is the difference in AIC of the alternative from the best model (Zimmermann and Foresti 2016). Jaguar density was estimated from the best model, hereafter referred to as model “no_sex.”

We also accounted for possible sex-specific capture probabilities (hereafter referred to as model “sex”) with a hybrid-mixture model that accounted for missing information on sex for some individuals (Zimmermann and Foresti 2016; Efford 2018a), to minimize differences in captures depending on sex (males are captured substantially more on camera-trapping studies due to their preference for trails; Foster et al. 2010a; Foster and Harmsen 2012). We analyzed its effect on g₀ and σ by comparison to the null and the best spatial model found in the previous analysis with model “no_sex” (following Zimmermann and Foresti 2016). We compared the models with the AIC method described above, and averaged all best-fitting models with a ΔAIC value <2 to calculate jaguar densities (Burnham and Anderson 2002).

The above models had no restriction on the locations of jaguar activity centers. We compared these models to alternative sets that restricted potential jaguar activity centers to lie inside EEER and EZCA boundary limits (hereafter referred to as model “habitat”). The extreme assumption in this case was that only areas inside the reserve polygons were suitable jaguar habitat, and that everything outside was not. This model served to illustrate what density the jaguars would have if the habitat outside reserves was unsuitable for their survival. We created a habitat mask to eliminate all potential activity centers falling in nonjaguar habitat, using a mesh cell size of 1 km² ( Supplementary Data SD3, Fig. S1 (gyac090_suppl_supplementary_figures.docx)), and applied it to the models used for model “sex.” We estimated jaguar density by model selection using the AIC and ΔAIC method described above. All analyses used the R package “secr” (Efford 2018a), adapting R scripts in Zimmermann and Foresti (2016).

For the Bayesian approach, we followed recommendations by Noss et al. (2012) for parameter definition at estimating jaguar densities. We used the following parameter definitions: (1) “trap response absent,” which is the conventional, global behavioral response, with the same response to all traps everywhere after first capture; (2) “spatial capture–recapture” model; (3) “half-normal” detection function, which generally has a good fit to the data; and (4) “Bernoulli process” binary encounter model (the only encounter model available in SPACECAP), in which the probability of success is derived as the probability of a positive response under a Poisson encounter rate model (Gopalaswamy et al. 2014). For the same habitat grid that we generated with model “habitat,” we set Markov and Monte Carlo simulation parameters at: (1) “iterations” (total number of Markov chain Monte Carlo iterations) = 100,000; (2) “burn-in” (number of initial values to be discarded during the analysis) = 10,000; (3) “thinning” (proportion of iterations included in the output) = 1 (all iterations stored); (4) “data augmentation” (maximum number of uncaught animals in the whole state space) in EEER and EZCA = 100, and NYP = 200; and (5) “habitat mesh cell size” = 1 km².

Per capita space requirement by jaguars.—Average space requirement per jaguar was estimated as the inverse of jaguar density, in km²/individual. This is not the same as a home range, which may comprise both shared and exclusive space, and which would require radio-tracking data. It does, however, indicate the area that each individual requires for itself. For example, a density of 2 ind./100-km² inverts to 50 km²/individual, which could mean one individual making exclusive use of 50 km², or two individuals sharing 100 km², or three sharing 150 km², and so on. This is a useful estimate of jaguar space requirement insofar as it represents the pressure that each individual exerts on its prey resource.

For our density estimates, in EEER we used records from August to October (n = 71) in 2015 and from March to May (n = 65) in 2016. In EZCA we used records from June to September (n = 90) in 2015, and from April to July (n = 90) in 2016. In NYP, we used records from July to October (n = 75 days) in 2015, and March to June (n = 71) in 2016. In order to calibrate our closed-population estimates of jaguar densities in NYP against estimates that are more likely to violate the closure assumption, we conducted a second density estimation in NYP using all our camera-trap sampling records, taken from April to October in 2015 (n = 179 days) and from March to August (n = 144) in 2016. The magnitude of difference in density estimates from long- and short-period sessions indicated the influence of closure on density estimation.

Spatial distribution of jaguars and pumas in NYP.—We tested for evidence of spatial attraction or segregation between jaguars and pumas in the observed fractions of camera-trap stations that detected them both, or either one alone (plotted in  Supplementary Data SD3, Fig. S2 (gyac090_suppl_supplementary_figures.docx)). For each of 2015 and 2016, we generated 10,000 replicate random and independent redistributions of jaguar and puma detections among the actual set of stations. These replicates quantified the probabilities of noninteracting jaguars and pumas producing at least as many as the observed number of: (1) H_0A—stations with co-occurring jaguars and pumas; (2) H_0B—nearest-neighbor pairs detecting only jaguars at both stations; (3) H_0C—nearest-neighbor pairs detecting only pumas at both stations. Rejection of H_0A (P < 0.05) indicates evidence of nonrandom attraction between jaguars and pumas to the same stations. Rejection of H_0B (P < 0.05) indicates that jaguars create an enclave free of pumas. Rejection of H_0C (P < 0.05) indicates that pumas create an enclave free of jaguars. These tests assume that each species finds all stations equally attractive in the absence of the other species. This means the tests assume that none of the actually empty stations are inaccessible, and the only possible nonrandom influence on either species using a station is the presence there of the other species.

Prey consumption, prey exploitation, and niche overlap by large felids.—Dietary preferences by jaguars and pumas were analyzed from 147 scat samples. Felid scats were sought systematically, and collected opportunistically, along trails in EEER (n = 109 scats) and EZCA (n = 38 scats) from May to July in 2015, and from March to June in 2016 ( Supplementary Data SD3, Fig. S3) (gyac090_suppl_supplementary_figures.docx). Fresh or moist samples were immersed in a 95% ethanol solution for 24 h and subsequently dried in the sun prior to storage. All samples were stored with silica gel beads at room temperature (Wasser et al. 1997).

For the analysis of felid diets, scat samples were washed with a solution of water and a small amount of detergent, then sifted to isolate hair, bone, hooves, teeth, and undigested remnants of consumed prey. Bone, hooves, teeth, and other remnants were identified with a magnifying glass, a vertebrate bone collection, and reference guides by Schmid (1972) and Olsen (1973, 1982). Prey hair mounting was conducted using a modified protocol by Moore et al. (1974). Five to 10 hairs were fixed onto slides using the synthetic resin Entellan. A color light photomicroscope Nikon Eclipse E400 (20×) was used to examine medullar structure patterns. Species identifications were made by comparison to the reference hair collections of Harmsen (2006) and Valdes (2006) from mammals of Belize, and reference guides by Monroy-Vilchis and Rubio-Rodríguez (2003), Baca Ibarra and Sánchez-Cordero (2004), Debelica and Thies (2009), Palacio (2009), Pech-Canché et al. (2009), and Cornally and Lawton (2016).

Due to the very similar morphologies of jaguar and puma feces (Roques et al. 2011; Aranda Sánchez 2012), scat samples were identified to species of donor using molecular genetic techniques. DNA extraction followed the CTAB protocol of Doyle and Doyle (1990), designed to extract DNA from plant leaves and adapted for feces (Chapman M.A., unpublished). Primers were designed from published ND5 mitochondrial DNA sequences of jaguars and pumas (GenBank accessions  KF483864,  KM236783,  KP202264 [jaguar],  KX808222,  KP202261,  KX808229, and  KX808231 [puma]). Sequences were aligned in GeneDoc (Nicholas and Nicholas 1997). Three primer pairs were designed to amplify ca. 150- to 200-bp regions, from which was selected the one that provided the most reliable and clean PCR amplification product. PCR was carried out in 15 µl reaction volumes containing 1.5 µl 10× buffer (peqGold Taq buffer S), 0.2 mM dNTPs, 0.1 mM of each primer: ND5_1625F (TAACTATAAGCCAAAAATCCGCA), ND5_1771R (GCTAAGGGYTAAGGTGATTATGAA), and 0.5U Taq polymerase (peqGOLD Taq). Amplification success was determined on 1% agarose gels. PCR amplicons were prepared for sequencing by incubating with 4 units Exonuclease I and 0.8 units Shrimp Alkaline Phosphatase (Fisher Scientific, United Kingdom) at 37°C for 45 min, followed by enzyme denaturation at 80°C for 15 min. PCR sequencing was carried out using BigDye v3.1 (at 1/8th the recommended volume) and the BigDye protocol (Applied Biosystems, United Kingdom). Reactions were precipitated with ethanol and sent to the Department of Zoology at the University of Oxford for resolution on an ABI 3730xl (Applied Biosystems). Sequencing electropherograms were visualized using Chromas (Technelysium Pty Ltd 1998–2001) and diagnostic SNPs were used to identify whether the scat was from a puma or jaguar.

Relative occurrence (%) of each prey species in jaguar and puma diets was calculated as the number of prey items belonging to species x divided by the total number of prey items, and multiplying the result by 100. Number of prey items was defined as the number of scats containing each species (Foster et al. 2010a). Relative occurrence tends to overestimate smaller prey species when compared with larger ones; thus, relative biomass per species was also calculated, using mean body weights extracted from published data (Ceballos 2014 and Reid 2009 for mammals; del Hoyo et al. 2018 for birds). For prey species with mean body weights ≥2 kg, the correction factor of Ackerman et al. (1984) was applied to mean live weight to account for their incomplete presence in an individual scat. The mean live weight was calculated as a geometric mean following Jaksić and Braker (1983). Smaller species were assumed to have been eaten whole and therefore to have biomass given directly by their mean live weight (Foster et al. 2010a). The relative biomass consumed (RBC) of each species x was calculated as:

The observed food niche breadth index (B_obs) was calculated for jaguars and pumas in terms of their dietary diversity, following Levins (1968). In order to allow for comparisons of diet breadth between jaguars and pumas, the index was standardized for each species following Colwell and Futuyma (1971).

Dietary niche overlap between jaguars and pumas was calculated with Pianka's (1973) measure, using the R package “EcoSimR” (Gotelli et al. 2015). In order to determine whether the observed dietary niche overlap (O_jk) was higher or lower than expected from random, a null model was generated to simulate possible overlaps between jaguars and pumas, generating 1,000 Monte Carlo randomizations of proportions of different food items in the diets of jaguars and pumas. The algorithm retained the observed zero states (prey species not consumed in our data set were also not consumed in simulations) and randomized the dietary niche breadth (prey species with relative occurrence > 0 were replaced with random prey proportions in simulations). All prey species were assumed to have equal availability to both felid species (Gotelli et al. 2015).

Results

Per capita space requirements by jaguars.—We identified 21 individual jaguars ( Supplementary Data SD3, Fig. S4 (gyac090_suppl_supplementary_figures.docx)) in NYP (n = 361 captures), of which 10 appeared only in EEER (five males, two females, three undetermined), eight only in EZCA (four males, three females, one undetermined), and three only in IA (one male, one female, one undetermined). No individuals were detected in more than one of these three areas. Given the male bias in detection across all sites (10:6 male:female), we used the maximum-likelihood density estimates that accounted for sex differences in capture probabilities (model “sex” in Table 1).

Estimates of space requirements of jaguars for the NYP region were 44.64 km²/individual (from a density ± SE of 2.24 ± 1.32 ind./100-km²) for 2015 and 42.19 km²/individual (2.37 ± 0.79 ind./100-km²) for 2016. Little difference was found in estimates for NYP calculated over longer periods, at 36.23 km²/individual (2.76 ± 1.56 ind./100-km²; 179 days) for 2015, and 40.16 km²/individual (2.49 ± 0.77 ind./100-km²; 144 days) for 2016. Estimates of space requirement were lower for the reserves within NYP. Estimates for EEER were 32.57 km²/individual (3.07 ± 1.91 ind./100-km²) for 2015 and 28.25 km²/individual (3.54 ± 1.85 ind./100-km²) for 2016; estimates for EZCA were 32.47 km²/individual (3.08 ± 1.62 ind./100-km²) for 2016. We could not estimate the space requirements for EZCA in 2015, due to the small number of camera-trap records and recaptures for estimations of jaguar density. When we excluded potential jaguar activity centers outside reserve boundaries in EEER and EZCA, we obtained estimates of space requirement that were lower, and densities that were higher, by up to an order of magnitude (model “habitat” in  Supplementary Data SD2, Table S3 (gyac090_suppl_supplementary_tables.docx)).

Table 1.

Best estimates of population density and individual space requirement for jaguars in El Edén Ecological Reserve (EEER) and El Zapotal Conservation Area (EZCA), and Northern Yucatán Peninsula (NYP) (combining EEER and EZCA), using a maximum-likelihood spatially explicit capture–recapture model with differences depending on sex (model “sex”). The right-hand column gives the diffuse sampling area: the smallest area just encompassing all detectable individuals (detailed in Materials and Methods section). The full list of results can be found in  Supplementary Data SD2, Table S3 (gyac090_suppl_supplementary_tables.docx).

The Bayesian approach gave higher estimates than those from the maximum-likelihood approach that accounts for sex differences, but similar estimates to those that do not account for sex differences ( Supplementary Data SD2, Table S3 (gyac090_suppl_supplementary_tables.docx)). Estimates of space requirement of jaguars in NYP were 131.58 km²/individual (0.76 ± 0.32 ind./100-km²) for 2015 and 67.57 km²/individual (1.48 ± 0.41 ind./100-km²) for 2016. Again, they showed little difference in longer trapping periods, with estimates of 131.58 km²/individual (0.76 ± 0.26 ind./100-km²) in 2015, and 66.67 km²/individual (1.50 ± 0.41 ind./100-km²) in 2016. Estimates of space requirement of jaguars for EEER were 84.75 km²/individual (1.18 ± 0.53 ind./100-km²) for 2015 and 79.37 km²/individual (1.26 ± 0.52 ind./100-km²) for 2016, and for EZCA were 32.36 km²/individual (3.09 ± 1.16 ind./100-km²) for 2015 and 33.78 km²/individual (2.96 ± 1.21 ind./100-km²) for 2016.

Spatial distribution of jaguars and pumas in NYP.—Of the 65 camera-trap stations in 2015 across NYP, 23 (35%) detected both jaguars and pumas, 15 (23%) detected only jaguars, 11 (17%) detected only pumas, and 16 (25%) detected neither felid ( Supplementary Data SD3, Fig. S2 (gyac090_suppl_supplementary_figures.docx)). Random and independent redistributions of the 38 jaguar and 24 puma detections in 2015 had probability P = 0.09 of obtaining at least these many co-occurrences. Thus, the two felids tended toward mutual attraction in 2015. This suggestion of nonrandom attraction disappeared in 2016, however, among the 69 stations across NYP. Then, 24 stations (35%) detected both felids, 19 (28%) detected only jaguars, 14 (20%) detected only pumas, and 12 (17%) detected neither felid. Random and independent redistributions of the 43 jaguar and 38 puma detections in 2016 had P = 0.54 of obtaining at least these many co-occurrences. The difference in probability between years was due in part to the larger number of puma detections in 2016, and fewer empty stations.

Of the total of 46 pairs of nearest neighboring stations in 2015 across NYP, four pairs detected only jaguars, and two pairs detected only pumas. Random and independent redistributions of the 38 jaguar and 24 puma detections in 2015 had P = 0.47 of at least this many jaguar-only pairs, and P = 0.64 of at least this many puma-only pairs. There was thus no evidence of either species holding enclaves from the other during 2015. The evidence for enclaving increased in 2016, however, among the 48 pairs of nearest neighboring stations across NYP. Then, seven pairs detected only jaguars, and three pairs detected only pumas. Random and independent redistributions of the 43 jaguar and 38 puma detections in 2016 had P = 0.05 of at least this many jaguar-only pairs, and P = 0.32 of at least this many puma-only pairs. Thus, jaguars tended to have exclusive use of areas in and around EZCA in 2016, while pumas occupied more of EEER and its surroundings, although mostly sharing the space with jaguars ( Supplementary Data SD3, Fig. S2B (gyac090_suppl_supplementary_figures.docx)). Records of large felid occurrences per camera trap were higher inside than outside protected areas for EEER (mean inside = 4.4, outside = 3.4) and IA (inside = 3.8, outside = 1.2), but not EZCA (inside = 2.0, outside = 2.1), and were higher inside than out for the regional data of NYP (inside = 3.1, outside = 2.7).

Prey consumption, prey exploitation, and niche overlap by large felids.—A total of 102 scat samples (70% of the total analyzed) contained jaguar or puma DNA. A high proportion of the cervid remnants contained in jaguar and puma scats could not be distinguished between species (white-tailed deer Odocoileus virginianus, Yucatán brown brocket deer Mazama pandora, and Central American red brocket deer Mazama temama). We therefore grouped remnants belonging to any cervid species into a single “Cervids” category. We calculated an average of the median weights for these three species in the Yucatán Peninsula for analysis using body mass.

Across NYP, jaguar diet contained at least 15 species; 13 mammals, one bird, and one reptile (n = 53 scats; Table 2). We did not register consumption of domestic prey species by either felid, although four instances of canids consumed by pumas could have included domestic dog, if not coyote (Canis latrans). Puma diet contained at least 20 species (n = 49 scats), 17 mammals, two birds, and one fish. Jaguars most frequently ate two large-prey taxa (>10 kg body mass, with >5% relative occurrence and biomass); collared peccary (Dicotyles tajacu) and Cervids. These were followed in frequency by medium-sized prey (2–10 kg); white-nosed coati (Nasua narica) and nine-banded armadillo (Dasypus novemcinctus). Together, these large- and medium-sized prey constituted 56% of prey items consumed and 76% of total biomass eaten. Pumas also most frequently ate the same two large-prey taxa; collared peccary and Cervids, followed by medium-sized species, particularly unidentified canids. Together, they constituted 50% of prey items consumed and 70% of total biomass eaten. Niche breadths were very similar (B_sta = 0.277 and 0.280 for jaguars and pumas, respectively). The diet of jaguars was thus less species-rich, but similar in niche breadth, to that of pumas.

Table 2.

Body weights, relative occurrences, and consumed biomasses (absolute and relative) of prey species in jaguar and puma scats collected in Northern Yucatán Peninsula (NYP). Ordered by body weight class. *Grouped as “Cervids” in interpretations.

Jaguars and pumas from NYP showed a high dietary overlap (Pianka index = 0.86). This value was larger than the random expectation of 0.23 ± 0.02 (mean ± SE) obtained from 1,000 randomizations (P = 0.002). Jaguars and pumas both took collared peccary and Cervids more frequently than other species, with a suggestion of a higher preference for Cervids by pumas than jaguars ( Supplementary Data SD3, Fig. S5 (gyac090_suppl_supplementary_figures.docx)). Jaguars showed a strong preference for nine-banded armadillo and pumas for canids despite their low detection (although nine-banded armadillo was very likely underdetected by camera traps). Both felids avoided ocellated turkey (Meleagris ocellata), despite their great abundance in the area (this was the species with the largest number of camera-trap records).

Discussion

This study has aimed to compensate for information-poor observations by collating observational evidence from diverse sources on the population ecologies of jaguars and pumas in NYP. The Yucatán Peninsula is an exemplar of a mosaic habitat of small, protected forest patches surrounded by unprotected forest and an agricultural matrix. Historically, it has been subjected to intense human use, which has led to land-use changes and fragmentation of its landscape, especially during the last century (Vester et al. 2007). Our study reports the first usage to our knowledge of a maximum-likelihood SECR method for estimating densities, and space requirement, of jaguars in NYP. Best estimates suggested space requirement of jaguars of 28–45 km²/individual. We place most confidence in our space requirement estimations obtained from densities calculated with the maximum-likelihood model “sex” (including sex as a covariate; Sollmann et al. 2011).

Male-biased capture probabilities have been found in a number of studies of jaguar density (e.g., Silver et al. 2004; Harmsen 2006; Foster 2008; Tobler and Powell 2013). Extreme heterogeneity in capture probability of a subset of the population will lead to lower reliability in abundance estimates (Harmsen et al. 2010). The inclusion of sex as a covariate has been shown to improve estimates (Sollmann et al. 2011). Estimates using alternative models that did not account for sex differences in capture probabilities (i.e., Bayesian approach and maximum-likelihood model “no_sex”) generally provided lower values, consistent with underestimating an undetected female component of density. Where the best-fitting model “no_sex” included heterogeneous capture, however, its estimates tended to be more similar to those of the model “sex.” The capture heterogeneity was likely responding to the sex difference that is built into the model “sex.”

A choice of camera-trap locations based on optimization of capture probabilities could result in biased sampling if the camera locations are only optimal for a subset of the sampled population (Foster and Harmsen 2012; Wearn and Glover-Kapfer 2017). Extrapolations of density estimates across areas encompassing lower-quality habitat should be done with caution (Foster and Harmsen 2012), and correction for this heterogeneity can ultimately influence density estimates in SECR models (Sollmann et al. 2011). Our study sites of EEER and EZCA, as private reserves, are likely to maintain good-quality habitat inside their limits; however, surrounding areas are owned by communal “ejidal” lands, which are likely to have lower-quality agricultural land. At EEER and EZCA, our camera-trap stations were placed both inside and outside reserve limits (which was a consequence of the area requirements needed for the CENJAGUAR design). Likewise, some stations in IA were located in communal “ejidal” lands. Thus, our sampling design allowed us to cover a wide range of habitat types, and ultimately to minimize this source of bias.

Due to the difficulty of identifying individual pumas from camera-trap photographs, it was not possible to estimate puma densities in our study. Harmsen et al. (2010) found that pumas in Belize tended to use forest trails more than jaguars, with jaguars more likely to use the forest matrix. If this holds true also for our study sites, then pumas will have had higher capture probabilities than jaguars in our sampling sites, due to our deployment of camera-trap stations mainly along trails. In this event, we have reasonable grounds for concluding that pumas may be present at somewhat lower densities than jaguars in EZCA, and at commensurate densities in EEER.

Jaguars and pumas in NYP showed some indication of spatial aggregation, although jaguars appeared to dominate some unshared areas. The two felids also showed a high dietary overlap with each other, prioritizing the same prey species, mainly collared peccary and deer, which were abundant and widespread in NYP.

The camera-trap stations between EEER and EZCA, and in IA, that we used to estimate jaguar densities in NYP (the pooled data), were spaced with larger distances than those recommended for studies of large felids (Wegge et al. 2004; Tobler and Powell 2013). Cameras spaced too widely may fail to detect individuals if they occupy home ranges that fall between trap locations, breaking the assumption of nonzero probabilities of capture for each individual, necessary for conventional capture–recapture models (Foster and Harmsen 2012). Nevertheless, predictions are possible outside the range of the data by making inferences from the sample to individuals that live in these holes, based on the explicit declaration that SECR models apply to any area within the state space, even to unsampled areas (Royle et al. 2014). We therefore consider our density estimates for EEER and EZCA to be more accurate than those for NYP. We interpret the NYP estimates with more caution, given the sparse distribution of camera traps between EEER and EZCA (including IA), which will likely result in underestimations of density, and overestimations of space requirements. Our NYP estimates nevertheless provide a useful insight into jaguar density in this area, in terms of contiguous occupancy with little mixing of individuals between EEER and EZCA.

Although our study lasted only 2 years, our estimations of jaguar population densities in NYP varied little, either between years or spatially between EEER and EZCA. By contrast, findings in a previous study (Ávila-Nájera et al. 2015) in EEER, using a Bayesian approach to estimate jaguar densities (not accounting for sex differences in detection), showed considerable variations in a 4-year study (0.7 ± 0.05 to 3.65 ± 1.39 ind./100-km²). Faller et al. (2007) also found highly variable estimates in an area within the NYP (including EZCA), using a nonspatial approach in a 3-year study (1.82 ± 0.17 to 6.18 ± 0.33 ind./100-km²). The temporal consistency in our study could reflect consistent environmental conditions throughout the 2 years, with no natural catastrophes impacting on abundance (Michalski and Peres 2007). In contrast, during the study of Faller et al. (2007), two hurricanes (Emily and Wilma, both category 5) and forest fires (during 2006) severely affected the region. During the study of Ávila-Nájera et al. (2015), a forest fire (May 2011) affected EEER.

Some ejidos surrounding private reserves of EEER and EZCA have joined the Mexican program of Payment for Ecosystem Services (SEMARNAT 2013) to promote the conservation of their natural resources. The jaguar space requirement estimated for EEER and EZCA that assumed suitable habitat only inside reserve limits (model “habitat”) showed evidence of the important hold that these ejidos have on the survival of jaguar populations in NYP. If they are lost, jaguar populations would not persist within the remaining protected fragments of EEER and EZCA. Although the movements of individual jaguars can overlap considerably within and between sexes (Harmsen et al. 2009), the space requirement of 2–7 km² per individual that would result from constraining the current population to only EEER and EZCA (model “habitat” in  Supplementary Data SD2, Table S3 (gyac090_suppl_supplementary_tables.docx)) would be an order of magnitude smaller than those found in any free-ranging populations. This unrealistically high packing exposes the vital contribution of the surrounding communal forested areas to the maintenance of jaguar populations in NYP.

Jaguar home ranges in the Yucatán Peninsula are typically an order of magnitude larger than the areas of our reserves, and show low levels of overlap between their territories (Chávez 2010; Cruz et al. 2021). With reserve areas of, respectively, 25 km² and 43 km², each reserve would barely suffice to support a single individual at the space requirement of 28–45 km²/individual estimated from our study. These small reserves can, however, provide refuge for large felids, given the effective deterrence of human hunting at least inside reserve limits. The observed preference of large felids for areas inside reserve limits suggests that the reserves serve this refuge function to some extent.

Our analysis of prey remains in jaguar and puma scats revealed broad diets for both species and substantial overlap between them, consistent with considerable opportunism in their carnivorous diet (Foster et al. 2010c). Ungulates were the favored prey of jaguars and pumas, whereas collared peccary had slightly higher biomass in jaguar diet, and deer in puma diet. With high dietary overlap and similar niche breadths, the principal difference between the two felid diets is in the longer list of rarely consumed species for the puma. Wider evidence for omni-carnivory in these felids is found in the notable geographic variations in diet according to prey availability. Previous studies have shown different degrees of dietary breadth and overlap, with jaguars often selecting for larger species, particularly peccaries (Polisar et al. 2003; Scognamillo et al. 2003; Weckel et al. 2006; Foster et al. 2010a). If we consider jaguars and pumas to be opportunistic carnivores, then their limiting resource is less likely to be any specific prey than the availability of enough habitat to support prey of any kind. These felids are efficient users of limiting resources in this respect, tolerant of disturbance to the extent that little suitable habitat is left unused by their populations and by the populations of their prey (Doncaster et al. 1996). Although the sample size of jaguar and puma scats available to estimate jaguar and puma diets in NYP was smaller than the minimum recommended to make accurate estimates of diet, it sufficed to detect the most common prey species (≥5% occurrence; Foster et al. 2010b).

Stable coexistence of carnivore species requires some level of segregation along trophic, spatiotemporal, or behavioral niche dimensions (Schoener 1974; Karanth and Sunquist 2000). Jaguars and pumas are thought to coexist by differential use of their habitats (Scognamillo et al. 2003; Harmsen et al. 2009; Sollmann et al. 2012). The felids have few options for spatial segregation in NYP. The highly disturbed landscape of flat karst limestone with practically nonexistent surface water (Bauer-Gottwein et al. 2011) presents a largely homogeneous secondary tropical dry forest containing all their prey, which are also sought by hunters, and inhospitable land converted to human uses of agriculture or settlement (Primack et al. 1998; González-Iturbe et al. 2002). In our study, jaguars and pumas showed possible evidence of mutual attraction (2015), which in the light of their high dietary overlap is consistent with a forcing together of the two species. Jaguars moreover made exclusive use of some habitat in EZCA (2016), which is consistent with an outcome of competitive dominance in otherwise shared areas. Coexistence of jaguars and pumas in NYP could be sustained by a constriction of niche separation between the felids, enforced by the combination of extreme daytime temperatures (which in the area typically exceed 30°C and frequently exceed 40°C) and lack of shelter, limiting their ability to thermoregulate, and the high levels of human disturbance outside the small reserves (Astete et al. 2017).

In conclusion, our results are consistent with spatial, but not dietary, partitioning in habitat use playing an important role for coexistence of felids in NYP, in combination with competitive dominance. Here, we identified a possible component of spatial segregation that may fluctuate between years. Our findings, nevertheless, indicate the low levels of segregation that can only be sustained by dominance–subordinate relationships in a disturbed environment (Doncaster 2009). We recommend that the spatial niche dimension be further assessed in future studies using GPS or satellite collars on individuals, to determine its role in facilitating coexistence of these two species in NYP. Sex-dependent temporal differences in habitat use have also been identified as a mechanism that can drive coexistence between these felids (Azevedo et al. 2018). We recommend that future studies assess intraspecific variation in activity patterns of jaguars and pumas in NYP.

The small reserves of EEER and EZCA cannot provide an entire home range for even a single jaguar or puma; however, they act as stepping stones for these large felids in NYP, connecting unprotected ejidal lands with large, protected areas of Yum Balam and Ría Lagartos reserves (Fig. 1). Overall, these four protected areas represent the northernmost boundary of continuous habitat for large felids in the Yucatán Peninsula, and play a crucial role in facilitating movement of jaguars across highly modified landscapes (Luja et al. 2017).

NYP still maintains vast areas of jaguar habitat with healthy conservation status (Sanderson et al. 2002), and it constitutes one of the areas in Mexico where jaguar conservation has good chances of long term of success (Carrillo et al. 2007). A particular concern, however, is the effect on these felids of subsistence hunting by people for game species that are the main prey in the diets of large felids. Long-term game hunting may force a higher foraging effort on large felids, which will eventually lead to an increased likelihood of encounters with humans and their livestock, resulting in lower felid carrying capacities and threats to their long-term persistence (Novack et al. 2005).