The study was conducted within a privately protected area in Las Piedras, Madre de Dios, Peru, centered around Hoja Nueva’s main research station (UTM 19L 442212E, 8667724N) ( Figure 1 ). The study area of 10 km² included both intact floodplain (seasonally flooded) and terra-firme (upland) forests, with a transitional zone in between. The land was primarily used for conservation and forest protection, with one incident of bushmeat hunting recorded during the study. Human activities were restricted to academic tourism, regular forest ranger patrols, and non-invasive wildlife research. The nearest mixed-use area with higher human impact, an agricultural association, was located more than 2 km south of the southernmost camera trap.

Camera traps were deployed at 62 sites over 80 nights: 33 sites within grid A from July to September 2022, and 29 sites within grid B from September to November 2022. Grid placement ensured representation of both floodplain and terra firme forests, accessibility within hiking distance from Hoja Nueva’s centers, and geographic proximity to represent the same forest type. Surveys were conducted sequentially to allow temporal comparison, all within the dry-wet transition season. Camera station locations were randomly generated for each grid, with an inter-site distance of approximately 450 meters to promote recaptures of both male and female margays. To minimize bias, all camera stations were strategically located away from human trails. Camera traps were deployed at semi-arboreal (2-5 meters) and terrestrial (40 cm) heights to investigate margay activity and behavior across the forest understory and forest floor, providing a comprehensive view of the species’ habitat use and movement patterns.

Upon reaching the predetermined coordinates, suitable margay trees were identified based on criteria such as an angle of 20 - 60°, minimum trunk diameter of 20 centimeters, and canopy connectivity ( Figure 2 ). Nearby locations for terrestrial cameras were then selected, ideally intersecting mammal trails where possible.

To assess occupancy and detection probabilities of margays across various habitat characteristics, data from all 62 sites were pooled. Arboreal and terrestrial camera traps were treated as part of the same physical site due to their proximity, often within a few meters of each other, focusing either on the ground level (terrestrial) or the lower canopy (semi-arboreal). This approach prevented artificially inflating the perceived survey area and avoided spatial autocorrelation errors in occupancy modeling. Differentiation between semi-arboreal and terrestrial cameras was achieved through a categorical observation-level covariate. This allowed robust comparisons within the same spatial site but across different environmental viewpoints.

Occupancy models, fitted using this configuration, recorded the presence or absence of margays for each camera setup during observation periods. The analysis considered not only environmental and habitat variables but also the specific impact of camera placement on detection probability. NDVI data from a May 2022 Sentinel-22 raster obtained from NASA’s (2024) Earth Observing System Data and Information System using the Harmonized Landsat Sentinel-2 (HLS) S30 version 2.0 was used to estimate vegetation density at each site. We also incorporated finer scale habitat characteristics including the average DBH within ten meters of each site, the number of trees larger than 15cm DBH at each site (TREES), and the SPI, which were all standardized and analyzed to determine their influence on margay occupancy. The SPI was calculated as the number of captures of prey under 5 kilos divided by the total number of trap nights per site. Other covariates, such as proximity to rivers and elevation, did not significantly influence margay occupancy in a preliminary analysis (unpublished data) and were therefore excluded from the models. None of the continuous covariates showed high correlation, allowing for independent inclusion in occupancy analyses.

The analysis was conducted using the ‘unmarked’ package (Fiske and Chandler, 2011) in R (R Core Team, 2022). An ‘unmarkedFrame’ object, incorporating standardized site-specific ecological covariates (NDVI, DBH, TREES, SPI) and an observation-specific covariate for detection (camera placement), served as the input for fitting the occupancy models. Occupancy and detection models were fitted using the ‘occu’ function, with sub-models for detection influenced by camera placement, and occupancy influenced by environmental covariates. Model selection, based on ecological relevance and statistical parsimony, was optimized using the ‘dredge’ function from the ‘MuMIn’ package (Barton, 2009), exploring all plausible covariate combinations and ranking them by Akaike Information Criterion (AIC-Akaike, 1973).

Given the presence of multiple competitive models within two delta AIC units, we opted to employ model averaging to robustly estimate the effects of environmental covariates on margay occupancy. This approach allowed for the integration of uncertainty across a suite of plausible models, providing weighted average estimates of each parameter based on their relative support by the data. Specifically, we reported conditional averages, which incorporate only the models where the parameter appears, reflecting a more nuanced interpretation of influence when a parameter is not universally relevant across all considered models. This method ensures that our results maintain robustness against the risk of model selection bias, particularly important given the ecological complexity and potential collinearity among covariates in our study. The average probability of both detection and occupancy for margays were then recorded using the top model based on AIC. Lastly, as the home range diameter of margays exceeds the distance between our camera traps, and therefore violates the assumption of independence between sites, resulting parameters were interpreted as proportion of area used rather than proportion of area occupied (MacKenzie and Nichols, 2004).

Margay photo captures were analyzed to identify individual animals based on unique coat patterns ( Figure 3 ). Initial identification involved one reviewer examining images from each paired camera setup, focusing on distinctive markings on both sides of the animals. HotSpotter software (Crall et al., 2013) facilitated the selection of specific rectangular regions of interest (ROIs) on the images to match against a labeled database of all captured margays. Images lacking clear coat patterns were excluded to maintain identification accuracy. A second reviewer cross-verified each identification using the original images and HotSpotter results.

To estimate margay density and investigate potential sexual dimorphism in spatial ecology, SECR models were employed using the ‘secr’ package (Efford, 2023) in R. The analysis focused on a single grid (Grid A) with 33 sites, which gathered over 26 occasions spanning 78 trap nights. Data were imported and formatted for compatibility with the SECR modeling framework. Capture times were standardized to the local time zone (America/Lima), and sampling occasions were defined based on the study period from July to September 2022. Effective trap operation per occasion at each site was calculated and integrated into the capture histories. We assessed the appropriate buffer size and mask resolution for the SECR models by testing four different mask configurations with buffer sizes of 1000m and 2000m, and spacings of 200m, 250m, and 300m. The models were compared using AIC values, and the mask with a 2000m buffer and 200m spacing was selected due to its lowest AIC, indicating the best fit.

NDVI was incorporated into the selected mask by extracting NDVI values derived from a May 2022 Sentinel-22 raster layer. We used the extract function to obtain the mean NDVI values within a 150m buffer around each trap location, scaling these values to standardize them. These NDVI values were added as covariates to the mask, enabling us to test their effect on margay density. We fitted six different SECR models to test the effect of NDVI on density, as well as the effect of sex on lambda0 (encounter rate) and sigma (scale parameter), all using a hazard half-normal detection function with a Poisson observation model. Model comparisons were based on AICc values to select the best-fitting model based on small sample sizes.

To investigate margay activity patterns and compare them with ocelots and jaguarundis, camera trap data from 64 distinct margay captures were analyzed. Activity patterns were visualized using ‘activityDensity’ and ‘activityRadial’ functions from the ‘camtrapR’ package (Niedballa et al., 2016), centered around the diurnal midpoint (12 noon). The overlap package generated kernel density estimates and calculated coefficients of overlap (Dhat statistics) between species. Radial plots visually represented activity throughout the day and overlap analyses with bootstrapped confidence intervals assessed the significance of observed overlaps between margay, ocelot, and jaguarundi activity patterns.

An extensive dataset from eight grids encompassing 293 camera sites, collected over seven years (2015-2022), was then integrated to enhance robustness. This dataset included 147 margay captures, 1,249 ocelot captures, and 40 jaguarundi captures, enabling a more comprehensive activity analysis.

The conditional model averaging approach provided a nuanced understanding of the influences on margay occupancy across our study area. Table 1 presents the conditional averages for key parameters, reflecting the combined influence of the top models. From the suite of environmental covariates considered, DBH emerged as a significant predictor, with a negative effect on occupancy probabilities (Estimate = -1.289, p = 0.0489).

Contrarily, NDVI, a proxy for vegetation health and density, did not significantly influence margay occupancy (Estimate = -4.1517, p = 0.3243). Tree density, measured by the number of trees, showed a positive but non-significant association with margay occupancy (Estimate = 0.8594, p = 0.1483). SPI, while not included in the top model, appeared in several models close to the top model (Estimate = 0.4119, p = 0.3234).

Detection probability was significantly influenced by the intercept (Estimate = -2.7460, p < 2e-16), indicating a generally low detection probability across sites, which is typical for cryptic species like the margay. The effect of the detection covariate (detTer), representing terrestrial versus semi-arboreal camera trap placements, was positive but not significant (Estimate = 0.5036, p = 0.1549). 23 out of the total 47 detections were semi-arboreal, suggesting no substantial difference in detection probability between semi-arboreal and terrestrial camera setups ( Figure 4 ).

Predictions from the top-ranked model estimated the average occupancy probability across the study sites at approximately 53.82%, suggesting that nearly half of the sites likely support margay presence. However, the detection probability under this model remained low at about 6.57%, underscoring the inherent challenges in detecting these cryptic animals, even under presumably favorable conditions.

The SECR analysis provided comprehensive insights into the population dynamics and spatial behavior of margays within the study area. Over the 78-day study period, a total of 34 capture events were recorded, involving 11 unique individuals. The selected model, based on AICc, incorporated sex as a covariate for sigma but not for detection probability (lambda0). Model results comparing AICc values can be found in Table 2 .

Our results indicated that NDVI did not significantly affect margay density within the study area. The best-fitting model estimated a margay density of 71.46 individuals per 100 km², with a standard error (SE) of 26.85. Details for the best-fitting model can be found in Table 3 . The encounter rate (lambda0) was determined to be 0.025 (SE = 0.00748), indicating a low yet consistent detection probability across both sexes. Analysis showed that detection probability did not differ significantly between males and females, suggesting equal likelihood of capture when margays were present within the detection range of the traps.

Sigma, the parameter representing the spatial scale of animal movements, demonstrated notable sex-specific differences. For female margays, sigma was estimated at 448.78 meters (SE = 124.61). In contrast, male margays exhibited a significantly larger sigma of 846.72 meters (SE = 196.36), indicating a broader spatial range for males. To elucidate home range sizes, 95% home range areas were calculated based on sigma values. The estimated home range for females was 3.79 km², whereas for males, the home range was substantially larger at 13.50 km². These values highlight the larger movement patterns and potentially greater territorial needs of male margays compared to females. The resulting sex ratio was 0.45 males to 0.55 females, indicating a slight female bias in the sampled population. This sex ratio, combined with the observed differences in spatial behavior, provides important context for understanding margay population dynamics in the study area.

The activity analysis of margays based on 64 records indicated distinct nocturnal peaks at 12 am to 3 am, 5 am, and notably at 6 pm, as shown in the radial plot ( Figure 5 ). This pattern establishes the margay as primarily nocturnal, contrasting sharply with the diurnal activity pattern of jaguarundis, which was predominantly between 7 am and 1 pm. The overlap coefficient (Dhat1) between margays and jaguarundis was notably low at 0.21, suggesting minimal temporal overlap.

Conversely, the comparison between margays and ocelots revealed a higher degree of temporal overlap (Dhat4 = 0.79), with significant activity overlap during the early evening around 6 pm. Despite this, margays tended to avoid times when ocelots were most active, particularly during early morning (4 am-5 am) and late evening (11 pm-12 am) peaks. This pattern is visually represented in Figure 6 , which also shows a moderate midday activity overlap between margays and ocelots, further contributing to the higher overall overlap coefficient.

The bootstrap analysis supported these findings, providing a confidence interval for the overlap coefficient between margays and ocelots (0.706 to 0.909), indicating a statistically significant overlap despite their distinct peak activity times. In contrast, the bootstrap results for the overlap with jaguarundis were inconclusive due to the small sample size (n=3 for jaguarundis), resulting in a wide confidence interval (0.013 to 0.413) and highlighting the need for more data for robust comparative analysis.

The analysis of the expanded dataset provided new insights and likely more accurate overlap estimates based on the higher number of captures per species. Based on mean bootstrap estimates, the overlap between ocelots and margays increased from 0.79 to 0.88 (CI 0.87-0.96), and the overlap between jaguarundis and margays increased from 0.21 to 0.35 (0.20-0.39). This broader dataset not only strengthens our current findings but also provides a more comprehensive view of interspecies interactions and activity patterns within this vital ecological region.