Forest cover mediates large and medium-sized mammal occurrence in a critical link of the Mesoamerican Biological Corridor

Connectivity of natural areas through biological corridors is essential for ecosystem resilience and biodiversity conservation. However, robust assessments of biodiversity in corridor areas are often hindered by logistical constraints and the statistical challenges of modeling data from multiple species. Herein, we used a hierarchical community occupancy model in a Bayesian framework to evaluate the status of medium and large-sized mammals in a critical link of the Mesoamerican Biological Corridor (MBC) in Costa Rica. We used camera traps deployed from 2013–2017 to detect 18 medium (1–15 kg) and 6 large (>15 kg) mammal species in a portion of two Jaguar Conservation Units (JCUs) and the Corridor linking them. Camera traps operated for 16,904 trap nights across 209 stations, covering an area of 880 km2. Forest cover was the most important driver of medium and large-sized mammal habitat use, with forest specialists such as jaguars (Panthera onca) and pumas (Puma concolor) strongly associated with high forest cover, while habitat generalists such as coyotes (Canis latrans) and raccoons (Procyon lotor) were associated with low forest cover. Medium and large-sized mammal species richness was lower in the Corridor area (x¯ = 9.78±1.84) than in the portions evaluated of the two JCUs (x¯ = 11.50±1.52). Puma and jaguar habitat use probabilities were strongly correlated with large prey species richness (jaguar, r = 0.59, p<0.001; puma, r = 0.72, p<0.001), and correlated to a lesser extent with medium prey species richness (jaguar, r = 0.36, p = 0.003; puma, r = 0.23, p = 0.064). Low estimated jaguar habitat use probability in one JCU (Central Volcanic Cordillera: x¯ = 0.15±0.11) suggests that this is not the jaguar stronghold previously assumed. In addition, the western half of the Corridor has low richness of large mammals, making it necessary to take urgent actions to secure habitat connectivity for mammal populations.


Introduction
Biodiversity is essential to maintain the resilience of ecosystems and the stability of their functions [1,2]. It is also critical for supporting a range of ecosystem services, reducing the risk of spread of infectious diseases, and maintaining productivity of several agricultural systems (e.g. facilitating pollination) [1,3,4]. While forests represent one of the most biodiverse ecosystems [5], global forest cover is being lost at a rate of 0.6% per year, largely due to conversion to agroindustrial land uses [6]. The increasing isolation of intact forest areas may not be enough to guarantee the long-term survival of species that require large spatial extents to sustain viable populations [7,8]. Thus, connectivity between separated forested areas through dispersal or biological corridors in human-dominated landscapes is crucial for the conservation of wideranging species [9,10].
With only 0.50% of the world's land area, the Mesoamerica region is recognized as a global hotspot, holding~7% of the world's biodiversity [11]. Because of its relatively small size and its geographic position between North and South America, this isthmus has functioned for millennia as a natural bridge for wildlife, becoming arguably the most critical region for habitat connectivity in the Americas. Consequently, in 1997 the governments in the region created the Mesoamerican Biological Corridor (MBC), an initiative to preserve biodiversity and connect protected and other natural areas from southern Mexico to Panama [12]. Nevertheless, forest is being lost in several parts of the corridor, with approximately 271,600 ha lost across the region in the last ten years [13], and some researchers have already highlighted possible areas where connectivity may be broken or close to broken [9,[14][15][16][17].
Medium and large-sized mammals play an important role in ecosystem dynamics, frequently representing the bulk of wildlife biomass in a given area [18], serving as predators exerting topdown control over other vertebrates [19], or performing ecological roles as prey species [20][21][22], seed consumers [23-25], or seed dispersers [24,26,27]. Medium and large-sized mammals are also among the most hunted animals by humans, and may comprise the main source of protein for some communities [18]. Importantly, large mammals incur substantial energetic costs and require large areas to maintain viable populations, placing many of them at higher risk of extinction and greater vulnerability to habitat alterations through human-related pressures [28][29][30].
The idea of large-sized mammals as "sentinel" species, decreasing in abundance or disappearing from areas where perturbations occur, has been a major driver of research into mammal ecology and distribution [29,[31][32][33]. Much work has focused on single species under the "umbrella" concept, meaning that the conservation of one wide-ranging species will conserve other sympatric species [30,[34][35][36][37]. However, some investigators suggest that the "umbrella" concept has not been tested appropriately [38][39][40], and single-species work is frequently confronted with small sample sizes that limit the analytical toolkit and scope of inference that can be applied to a specific area [41][42][43].
To overcome the problems of low sample sizes and/or detection rates, researchers often turn to occupancy models. Modern occupancy models account for imperfect detection by incorporating detection probability [43,44], producing less biased estimates compared to naive occupancy estimation [45]. The extension of occupancy models to a community or multi-species framework [46][47][48] allows the estimation of species and community-level occupancy by drawing species-level estimates from community-level hyperparameters [45,[47][48][49][50], which can improve estimates for rare or elusive species. For these reasons, a multiple species approach to assess biodiversity, evaluate the effects of impacts or management actions, and evaluate connectivity has proven to be a useful alternative to the aforementioned methods [16,[46][47][48]50].
In this paper, we use a multi-species community occupancy approach [46,50] to establish baseline information for medium and large-sized mammals, including two large predators, in a critical wildlife corridor and the adjacent protected areas in Costa Rica. This corridor was primarily created to secure jaguar (Panthera onca; IUCN: near threatened; [51]) population connectivity. Thus, this area is considered a crucial corridor for the MBC and the Jaguar Corridor Initiative (JCI), the largest-scale carnivore conservation effort to date [9,15,52]. The JCI aims to preserve jaguar populations and range-wide habitat connectivity from Mexico to northern Argentina by identifying and securing dispersal corridors between core populations, also known as Jaguar Conservation Units or JCUs [9]. Jaguars, alongside pumas (Puma concolor; IUCN: least concern; [53]), are the largest predators in the region and their presence elsewhere has been associated with prey biomass and availability [22,[54][55][56]. However, some studies have also found that prey richness can be related to large carnivore richness [57], occupancy [15], or dietary niche breadth [58].
Our objectives were to: (1) determine the environmental and human-related factors driving the occurrence of 24 medium (N = 18) and large (N = 6) native mammal species, plus the nonnative domestic pig (an important prey item for jaguars and pumas in the study area), and (2) evaluate differences between the corridor and adjoining JCUs in terms of (a) species richness of medium and large-sized native mammal species and (b) habitat use probability for jaguars and pumas.
We expected that medium and large mammal species richness would be lower in the corridor than in the JCUs, given that the latter have higher forest cover and are classified as protected areas or indigenous territories. We also hypothesized that habitat use probability of jaguars and pumas would be related to large prey species richness given carnivores' high energetic requirements. We discuss the implications of our findings and the practicality of our model for this critical link of the MBC and the JCI, and for initiatives directed at preserving and monitoring biodiversity.

Study area
The study area comprises the Barbilla-Destierro Biological Corridor (hereafter "the Corridor") and a portion of two adjacent JCUs: the Central Volcanic Cordillera (CVC) JCU and the Talamanca-Cordillera Central (TC) JCU (Fig 1). The JCUs are expert-defined areas that are believed to have resident jaguar populations, an adequate prey base, and high habitat quality for this species [8,59].
The CVC JCU comprises a continuous block of protected areas in central Costa Rica (IUCN categories II and VI;~1,153 km 2 ) to the west of the Corridor, while the second JCU (TC JCU) is shared between Costa Rica and Panama and connects to the south-eastern side of the Corridor. The TC JCU encompasses a continuous group of protected areas (IUCN categories II and VI) and indigenous territories (~7,002 km 2 ) on the Costa Rican side. Primary and secondary forests represent~75% of the JCU area, with the remaining area largely dominated by pasture [60].
The Corridor between those two JCUs is approximately 362 km 2 , of which 58% is protected (IUCN category VI; Forest Reserves and a Protected Zone) or indigenous territory with a certain level of protection. No other suitable connections for jaguars have been identified between the CVC and TC JCUs, and more broadly between Nicaragua and Panama [15]. Almost all the area within the corridor is privately owned, comprising small farms. About 64% of the Corridor is covered by primary and secondary forest (2015 Aster Image; AIST https://gbank.gsj.jp/ madas/), with the rest of the area dominated by pasturelands for livestock (20%), and agriculture (14%). There are two two-lane paved roads within the Corridor that connect the towns of Turrialba and Siquirres.

Study design and data collection
We conducted camera trap surveys to assess the presence of medium (1-15 kg) and large (>15 kg) mammal species over the study area. The study area was divided in four different blocks (CVC JCU, TC JCU, and Corridor Blocks 1 and 2). We created a grid system of 63 16 km 2 cells over the entire Corridor and portions of the CVC and TC JCUs (Fig 2). Cell size represented the approximate home range size of jaguars in Central America, one of the target species in our analyses, and presumably the species with the largest home range size [61][62][63].
The 16 km 2 grid cells were subdivided into four sub-cells of 4 km 2 each (Fig 3). We sampled two sub-cells per grid cell with two stations (one camera trap per station) in each sub-cell for a period of approximately three months. Selection of these sub-cells was random, but had to be adjusted when there was no forest, permits were not given, or access was difficult (e.g. presence of very steep slopes). Blocks of the study area were sampled at different times due to limitations in the number of cameras available and logistical considerations, as well as to allow for a higher . No major land use or management changes occurred during these periods, though we account for potential baseline differences in mammal habitat use across blocks via random intercepts as described later.
We placed motion sensitive camera traps (Panthera1 V3, V4, V5 and V6) in forested areas, strapped to trees at approximately 0.4-0.5 m parallel to the ground. Cameras were set to function continuously and to take three shots in every event during the day and one shot during the night. One camera of each 4 km 2 surveyed sub-cell was placed off a trail and the other one was placed on a human-made trail (when available), in an attempt to detect species that may avoid or use human trails (Fig 3).

PLOS ONE
Camera location (±30 m) was recorded using a GPS device (Garmin1). Cameras were checked approximately every six weeks to download pictures and perform camera maintenance. Data were processed on Panthera IDS (Integrated Data Systems; version 1.13.786), where the species, date, time, and number of individuals in each photograph were recorded. For the occupancy-based analyses, we grouped data from all stations in each grid cell to make detection histories per week ("1" = detected, "0" = not detected, "NA" = inactive camera) for each of the 24 medium (N = 18 species) and large (N = 6) native mammal species, plus domestic pig (Total N = 25) ( Table 1). We included domestic pig as they were frequently recorded and were consumed by jaguars and pumas in the study area, especially in indigenous areas where they roam freely in the forest (R. Salom-Pérez, personal observation).
Site covariates were selected a priori and hypothesized to be the main drivers of occupancy of medium and large sized-mammals in the study area [9,15,16,59,64,65] (S1 Table in S1 File). Site covariates expected to have a positive relationship with occupancy of medium and large mammals were enhanced vegetation index (EVI) [66], forest cover [65], distance to a primary road [65], distance from any human settlement, and distance from major human settlements (see detail in S1 File) [65]. Site covariates expected to have a negative relationship with occupancy of medium and large mammals were elevation [67], terrain ruggedness [68], distance to strictly protected areas (IUCN Ia & II categories) [65], distance to JCUs (Panthera unpub. data) and human presence (calculated as the number of human detections per 1,000 trap nights per stations; S1 File). All covariates were normalized prior to inclusion in the models such that the magnitude of regression coefficients could be compared within and among models [69]. To calculate effort, a covariate on detection, we calculated and standardized the sum of all trap nights on each occasion for every grid cell.
Permits for data collection were granted by the Costa Rican National System of Conservation Areas (SINAC; Spanish).

Multi-level community occupancy model
Given that our survey occurred over a long period of time (~3.5 years), the closure assumption (i.e., there are no changes in the occupancy status of grid cells during the survey period [44]) is likely violated. Thus, occupancy (C) is interpreted as probability of habitat use [70], and we assume that any changes in the occupancy status of grid cells over our survey period were random and that there were no major changes in the area throughout the study period.

PLOS ONE
We used a multi-level (or hierarchical) community occupancy model in a Bayesian framework to calculate jaguar and puma habitat use probability and species richness [46- 48,50]. Unlike traditional occupancy models, the Bayesian framework can account for unobserved heterogeneity across species, space or time through random effects [15]. This type of heterogeneity could be expected in the current investigation given the relatively large spatial extent and the fact that the blocks were surveyed at different time periods.
To determine the covariates included in the global community occupancy model, we first ran single-species, single-season occupancy models for each medium and large-sized mammal species using R package RPresence [71]. We included effort as a covariate on detection in all models. On the occupancy side, we used all possible combinations of our 11 site covariates (additive only, no interactions), taking care to exclude variables from the same model if correlated at |r| > 0. 60 ([72]; S2 Table in S1 File). In total, there were 79 models for each species (S2 File). Following Fieberg et al. [73], we then summed the log-likelihoods of each model across our 25 species and calculated Akaike Information Criterion with correction for small samples (AIC c ), with n = number of species [71]. This approach assumed independence between species observations [44].
Our final hierarchical occupancy model included forest cover, human presence, EVI, distance to strictly-protected area, and terrain ruggedness on habitat use, and effort on detection (S1 and S2 Tables in S1 File). Habitat use probability for each species at each site was estimated Table 1. Medium and large-sized mammal species included in the global analysis (N = 25). Those included in the habitat use analysis for jaguars and pumas are indicated by a check mark. Barbilla-Destierro Biological Corridor and portions of the Central Volcanic Cordillera (CVC) and Talamanca-Cordillera Central (TC) Jaguar Conservation Units (JCUs), surveyed with camera traps from 2013-2017.

PLOS ONE
as: where ξ il * Normal(μ ξ ,τ ξ ) is the random intercept for each species i at each block l (CVC JCU, Corridor Blocks 1 and 2, and TC JCU), μ ξ is the community-level or hyperparameter mean for the intercept on habitat use, and τ ξ is its precision. The random intercept was used to account for potential spatial and temporal heterogeneity in habitat use due to differences by species and survey block [47,48,50]. In addition, α i are estimated beta coefficients on habitat use for species i, where α i * Normal(μα,τα); μα is the community-level or hyperparameter mean for each beta coefficient, and τα is its precision. D j are the standardized values for each covariate at grid cell j. Detection probability for each species i at each site j and in each week k, was estimated as where v il * Normal(μ v ,τ v ) is the random intercept for each species i at each block l, μ v is the community-level mean for the intercept on detection, and τ v is its precision. In addition, β i is the estimated beta coefficient for effort for species i, where β i * Normal(μ β ,τ β ), μ β is the community-level mean for effort, τ β is its precision, and effort(j,k) are the standardized values for effort at grid cell j on occasion k. We defined true occurrence z( i,j ) as a binary variable in which z( i,j ) = 1 if species i occurred in grid cell j and = 0 otherwise. We modeled occurrence from a Bernoulli random variable, where z i,j~B ern(C i,j ), where C i,j is the probability that species i occurs at grid cell j. Importantly, species richness per grid cell was calculated as a derived parameter from the summation of z i,j values.
To account for imperfect detection, we modeled observed data y(i,j,k) as Bern(p i,j,k � z i,j ), where p i,j,k is the detection probability of species i in grid cell j in the survey occasion (week) k. The model accounted for the effect of species abundance on detection probabilities via species correlation parameter rho (ρ), a correlation between habitat use and detection probability [46].
We fit the Bayesian models in R 3.5.1 (R Core Team 2018) using package jagsUI [74], specifying three MCMC chains of 30,000 iterations, a burn-in of 5,000, and a thinning rate of three.
We estimated jaguar and puma habitat use probability as a function of species richness from the community model. We calculated species richness (again, via summation of z i,j values) for (1) all medium (N = 18) and large-sized (N = 6) native mammal species, including jaguar and puma and excluding the domestic pig; (2) jaguar and puma large prey species, including domestic pig (N = 4); (3) jaguar medium prey species (N = 10); and (4) puma medium prey species (N = 11 species). In order to select the medium and large prey species for jaguars and pumas, we conducted a literature review of publications on jaguar and puma diet and predation reports from Mexico to Panama (S4 Table in S3 File).
Lastly, we estimated the correlation between jaguar and puma habitat use probabilities and (1) medium prey richness and (2) large prey species richness, both of which were derived parameters from the prey community occupancy model as stated above.

Results
Camera traps operated for 16,904 total trap nights across 209 stations. Fifty-five out of 63 total cells were surveyed, covering 87.30% of the study area (880 km 2 ). We registered 2,946 independent records ("independent" defined as records separated by 24 hours or occurring at different camera sites) of medium (N = 18) and large-sized (N = 7) mammal species, including domestic pig (S5 Table in S4 File).
All variables informing habitat use had a 95% CI overlapping zero, suggesting imprecision in parameter estimation and likely high variability of covariate effects among species (Table 2 and Fig 4 and S1-S4 Figs in S5 File) [46]. Percent forest had the greatest influence on species richness, with 89.26% of its posterior distribution above 0 and an effect size more than three times that of any other covariate (Table 2). At the species level, this relationship was clear (i.e., 95% CI did not overlap 0) and positive for collared peccary, jaguar, domestic pig, paca, puma, agouti and ocelot, while it was negative for coyote, nine-banded armadillo and raccoon (Fig 4). There was no clear effect of forest for the other species.
While the two JCUs performed similarly on most other metrics, they greatly differed with respect to jaguar habitat use, with jaguars having very low probability of habitat use in the CVC JCU (� x = 0.15 ± SD 0.11) in comparison to TC JCU (� x = 0.58 ± SD 0.16) (Fig 6A).

Discussion
To our knowledge this is the most intensive camera trap study on a continuous area in Costa Rica [75], and the first to take place in this critical link between two JCUs. We found that forest cover was the main driver of medium and large-sized mammal habitat use at the community level, with evidence of widely differing relationships at the level of individual species. The hierarchical community approach allowed for the estimation of species and community-level effects, and for the incorporation of data from rare species. The model also accounted for heterogeneity in the sampling process through the incorporation of random effects.
While the medium and large-sized mammal species richness in the Corridor is slightly lower (� x = 9.78±1.84 spp.) in comparison to the JCUs (� x = 11.50±1.52), the number of records of large-sized species, is considerably lower in the former (Corridor: 38 independent records in 8,952 trap nights vs JCUs: 200 independent records in 7,952 trap nights). Additional research on the CVC JCU (where jaguar habitat use was low) is necessary to evaluate the viability of this critical link of the MBC and the JCI for jaguars, as jaguar records were entirely on the east side of the study area (Corridor Block 2 and TC JCU). This baseline information will be of paramount value to measure the outcome of conservation initiatives in the years to Forest cover was the most important covariate in our model related to medium and largesized mammal habitat use, having an effect that was at least three times higher than any other covariate. This covariate seemed to be especially important for certain species, including the collared peccary, jaguar, domestic pig, puma, agouti and ocelot. These species, with exception of the domestic pig, are known to depend on, or at least be associated with, vegetation cover [51,53,[76][77][78]. Almost all domestic pigs we detected were in the eastern side of the corridor, specifically in or near indigenous territories with high forest cover. These animals belong to the indigenous people and roam freely in the forest. On the other hand, coyote, nine-banded armadillo and raccoon seem to avoid areas with high forest cover in the study area. This was not surprising, as these are adaptable species that can be found in open and/or disturbed areas PLOS ONE [76,[79][80][81]. Thus, our model results were in line with disturbance-related species ecology and can help inform management decisions.
Contrary to our predictions, human presence and EVI were positively and negatively associated with overall species habitat use, respectively, though both effects were weak. While human presence is often used as a proxy for disturbance and has found to be negatively associated with habitat use at the community level [48], in our case the presence of humans may have a different interpretation. For example, people in our study area (e.g. tourists, hunters, indigenous people) may be actively looking for these animals and selectively walking on the same trails or in the same areas where they occur, leading to a positive association with mammal habitat use. This potentially important result requires additional study.
As we anticipated, medium and large-sized native mammal species richness was lower in the Corridor area than in the two JCUs. Nonetheless, these differences are subtle, indicating that there are some areas of the Corridor that are still in good condition, especially in the eastern portion close to the TC JCU. The lowest values for richness for these mammals, especially for large prey species, as well as for jaguar and puma habitat use probabilities, were in the western half of the corridor.
Low jaguar habitat use probability in the western JCU (CVC) was unexpected, given its apparent suitability for this species based on forest cover and presence of prey [52,65]. However, jaguar presence is likely not supported by other characteristics of the area, as most of the surveyed area in this JCU is not strictly protected (IUCN category VI), the terrain is very rugged and some areas have high elevation (2,000-3,300 m.a.s.l.) [51,82]. The highest habitat use probability for jaguars was located inside or near a strictly protected area, Barbilla National Park (IUCN category II) in the eastern JCU (TC). In contrast, pumas have high habitat use probability in both JCUs. Pumas are known to occur more frequently in higher elevations than jaguars, and could be benefitting from the apparent absence of a direct competitor in the CVC JCU [53,82,83]. A recent investigation further west into the CVC JCU found no sign of jaguar, adding extra support to the hypothesis that this area is not the jaguar stronghold previously assumed (Velado et al. unpub. data). It remains to be established whether low jaguar presence

PLOS ONE
in this JCU is explained by the site conditions mentioned above or if it is the result of more historical pressures (e.g. hunting, isolation).
Low richness of medium and large-sized prey within the Corridor was largely found between the two main roads. Although the proximity to road covariate did not make our global habitat use model, this and other types of infrastructure are known to have a negative effect on the presence of certain mammal species [32,84,85]. We detected the presence of large prey species, jaguars and pumas near these roads, but they were absent from the area between them. This area has less forest cover and higher settlement density, factors highly related to the presence of roads. Additionally, they are the Corridor grid cells located at greatest distance from the strictly protected areas and JCUs. We recommend further investigation into the potential barrier effect of these roads (and not just evaluate road proximity or density), such as the use of telemetry and non-invasive genetic tools to investigate movement and gene flow [14,84,86].
We found that puma and jaguar habitat use probabilities were strongly correlated with large prey species richness and, to a lesser degree, with medium prey species richness. Our results are consistent with other research on trophic interactions [57] and energetic constraints [87] showing a stronger link between large-bodied mammal predators and large mammal prey species in comparison to smaller prey species. The relationship between jaguar and puma occupancy with prey richness should be explored further in future studies, controlling for potentially confounding factors that can be correlated with prey richness, such as forest cover and biomass.

Conclusions
Our results highlight the importance of generating on-the-ground information on the status of multiple species within population source sites and corridor areas, as well as using a hierarchical modeling framework for robust parameter estimation at the community and individual species levels to inform management decisions. The ability to account for heterogeneity in the sampling process (e.g. data taken over several sites and different years) makes this model versatile and easily adapted to different species and study systems.
Urgent actions are needed to secure connectivity of mammal populations in our study area and within the greater Mesoamerican Biological Corridor, and should focus on (1) increasing forest cover in the western half of the Corridor, notably between the two main roads, and (2) increasing habitat quality and conditions for prey species, with a particular emphasis on large species.