Calibrating occupancy to density estimations to assess abundance and vulnerability of a threatened primate in Tanzania

The current decline of mammals worldwide makes quantitative population assessments crucial, especially for range-restricted and threatened species. However, robust abundance estimations are challenging for elusive or otherwise difficult to detect species. Alternative metrics requiring only presence/absence data, that is, occupancy, are possible but calibration with independent density estimates should be foreseen, although rarely performed. Here, we calibrated density estimates from acoustic surveys to occupancy estimates from camera-trapping detections to derive the abundance of the endangered Sanje mangabey ( Cercocebus sanjei ) across its entire range in the Udzungwa Mountains of Tanzania. We found marked occupancy – density relationships for the two forest blocks where this primate occurs and used them to derive spatially explicit density estimates. Occupancy increased in montane forest zones at mid-elevation but decreased slightly


INTRODUCTION
In the Anthropocene, wildlife species are experiencing an unprecedented decline, with 26% of mammals currently threatened with extinction (International Union for Conservation of Nature [IUCN], 2021).While the highest concentration of threatened and range-restricted species are found in tropical forests, these represent one of the most affected biomes, with habitat loss and harvesting that disproportionally impact mammals (Roberts et al., 2021;Schipper et al., 2008).Nonhuman primates are among the most threatened order of mammals, with ~65% of the 504 existing species being threatened with extinction and ~75% of the populations facing decline (Estrada et al., 2017;Fern andez et al., 2022).By functioning as seed dispersers, these species serve a critical ecological role in forest regeneration, affecting plant species diversity and demography (e.g., Andresen et al., 2018;Holbrook & Loiselle, 2009;Lambert, 2010).For these reasons, nonhuman primates are considered excellent ecological indicators and are highly sensitive to anthropogenic disturbance (Marsh, 2003;Rodrìguez-Luna et al., 2013).Therefore, research and conservation efforts focused on primates may have the added benefits of preserving important ecological functions as well as other species (Gippoliti & Sousa, 2004;Lambert, 2010;Martins & Valladares-Padua, 2005).However, the quantitative knowledge on abundance, conservation status, and vulnerability of threatened primates needed to promote sound conservation remains limited (Estrada et al., 2017), and rarely does this cover the entire range of species (e.g., Davenport et al., 2022).
A key limiting issue of such assessments is the pervasive difficulty of estimating population abundance and its spatial variation, especially for rare and elusive free-ranging animals.Hence, accounting for imperfect detection is important, especially for species that are difficult to detect, and methods that include it are highly recommended (e.g., Cavada et al., 2019;Spehar et al., 2015).In place of direct counts of individuals or groups, a practical approach is the use of detection/nondetection data (Joseph et al., 2006;Kühl et al., 2008).The latter are suitable for estimating occupancy, which is broadly considered a surrogate for abundance, and can be modeled with spatial covariates to assess species response to drivers of change (Kéry & Royle, 2016;Linden et al., 2017;MacKenzie et al., 2002MacKenzie et al., , 2006)).However, the assumption of a positive relationship between the two metrics is often untested and ideally requires the calibration of occupancy with independently derived density estimates (Linden et al., 2017).
Here, we present the results of an approach that calibrates occupancy from camera-trapping detections to density from acoustic surveys to estimate population abundance and spatial variation of an endangered primate in relation to both environmental and anthropogenic factors.We selected the Sanje mangabey (Cercocebus sanjei, Mittermeier et al., 2006) for this study, an IUCN-endangered, predominantly ground-dwelling, and frugivorous primate first described in 1979 and occurring only in two separated forests in the Udzungwa Mountains of south-central Tanzania (Ehardt et al., 2005;McCabe et al., 2019;Rovero et al., 2006).The area is of outstanding importance for biodiversity conservation (Burgess et al., 2007;Rovero, Menegon, et al., 2014) and a high-priority area for primate conservation in Tanzania (Davenport et al., 2013).
Despite several ecological studies, robust knowledge on population abundance and habitat association remains preliminary, with studies rarely targeting the whole species range with consistent methodology.Earlier estimations of population abundance were based on crude extrapolations of home-range data from a few groups (e.g., Ehardt, 2001;Ehardt et al., 2005;Rovero et al., 2009).More recently, Paddock et al. (2020) used acoustic detections to estimate density using a distance sampling approach, hence taking into account imperfect detection.This method enabled group detection using the distinctive mangabey whoop-gobble long-call, which is given by males in each group in the mornings, allowing for remote sensing without the need for visual observations of this elusive species.However, effort was limited to a single survey across a spatially and numerically limited set of listening posts, and the analysis did not consider spatial covariates.Similarly, habitat associations have been studied with camera traps, but only for parts of the species range (Martin et al., 2015;Oberosler et al., 2020a;Rovero, Martin, et al., 2014) and considered only a few environmental and anthropogenic variables.Studies showed that camera traps are very efficient at detecting this ground-dwelling species (Hegerl et al., 2017;Martin et al., 2015;Oberosler et al., 2020a), which is instead very rarely sighted by human observations (Barelli et al., 2013;Rovero et al., 2012).
Here, we aimed to provide more accurate quantification of abundance and density of this elusive and iconic species by integrating camera-trapping and acoustic data to overcome previous limitations, particularly the difficulties associated with collecting systematic and spatially comprehensive density data from acoustic detections.Thus, we used camera-trap data collected consistently in both forests, and hence across most of the species' range, to calibrate spatially explicit occupancy estimates to density from acoustic surveys, and used the regression to derive density estimates across the species' range.Moreover, by integrating visual counts of group size, we also predicted individual densities and total population abundance.Lastly, we also aimed to determine spatial patterns of occurrence in relation to both habitat and potential anthropogenic disturbance to better assess habitat associations and vulnerability of this endangered and iconic primate.

Study area
The Udzungwa Mountains in south-central Tanzania (7 40 0 -8 40 0 S and 35 10 0 -36 50 0 E; Figure 1) represent the largest block within the Eastern Arc Mountains and the richest in biodiversity (Rovero, Menegon, et al., 2014).Part of the area is protected as Udzungwa Mountain National Park (1990 km 2 ), established in 1992, whereas the remaining part is preserved either as nature reserve or forest reserve (Rovero et al., 2009).The two forests where the Sanje mangabeys occur are Mwanihana (MW, 167 km 2 ) within the Udzungwa Mountain National Park, and Uzungwa Scarp Nature Reserve (USNR, 372 km 2 ), which is approx.150 km to the southwest (Figure 1).These areas are both east-facing escarpment slopes ranging in elevation between 290 and 2100 m above sea level (asl), with a gradient of vegetation cover from lowland deciduous to montane evergreen and upper montane bamboo forests.Lowland zones in both areas are dominated by regenerating and secondary forest due to past logging and degradation, whereas the higher elevations are mainly characterized by undisturbed, closed-canopy forest.The two forests share similar ecological characteristics (Appendix S1: Table S1), but different management regimes: MW has effective law enforcement being in a well-protected national park, while USNR is less efficiently protected (Rovero et al., 2012) and more degraded (Oberosler et al., 2020a;Rovero, Menegon, et al., 2014).In fact, the nature reserve is a forest island surrounded by small villages and suffers from logging, habitat destruction, and illegal hunting, whereas MW has minor anthropogenic disturbance from the villages located east of the park boundary.The listening posts were located randomly inside cells of a systematic grid, with a minimum distance between locations of 2 km (Figure 1), which is in line with the known home-range size of the target species (Ehardt et al., 2005).The distance to each mangabey vocalization heard was estimated by three independent observers and the distance and direction mapped to estimate the position of the group.Estimates included in the analyses were truncated at 1 km as it was unlikely calls were accurately detectable over this distance.Data collection was carried out in the morning, between 07:00 and 09:00, when the Sanje mangabey is known to call at the highest frequencies (Ehardt et al., 2005), and every survey was conducted once (further details in Paddock et al., 2020).

Data collection: Acoustic and camera-trapping surveys
Camera-trapping surveys were conducted in both forests in the same period as the acoustic survey (i.e., July-November 2017).Data in MW were collected through the Tropical Ecology Assessment and Monitoring (TEAM) Network project (Rovero & Ahumada, 2017), which follows a standardized and systematic protocol for monitoring medium-to-large terrestrial mammals, and the same protocol was extended in USNR.Thus, for each forest, we placed three consecutive arrays of 20 camera traps (Reconyx HC500, Reconyx Inc., Holmen, WI, USA) for a total of 60 camera stations per forest (Figure 1).Cameras were located based on a systematic grid of 2-km 2 cell size and deployed in the field for a minimum of 30 consecutive days.Inside each cell, camera traps were placed on the closest wildlife trail to the centroid of the cell and attached to trees at approximately 50 cm from to the ground, facing the presumed trail at approximately 2-3 m distance.The total sampling effort yielded 1839 camera days for MW and 1792 for USNR (Appendix S1: Table S1) from 117 camera-trap sites overall, as three were damaged, with a mean effort per camera of 31 days for both study areas.

Data collection: Covariates
We derived six covariates that could affect the species spatial distribution both in terms of density and occupancy and used them as proxies of environmental and anthropogenic factors.Thus, we considered: (1) habitat cover, extracted from the GLC2000 dataset (European Commission, 2004) and based on Food and Agriculture Organization of the United Nations habitat classification (i.e., broadleaf evergreen forest and deciduous forest), (2) elevation, extracted from a 90-m digital elevation model raster, (3) ground slope, (4) distance from the camera location to the closest river, and (5) distance to the closest forest border.Values were extracted using the built-in tools in the open-source software QGIS (QGIS Development Team, 2019).Based on Pearson's correlation coefficient (threshold = 0.5), these variables were not collinear.

Group density estimates
We reanalyzed data from Paddock et al. (2020) with the aim of building density estimation models with account of anthropogenic and environmental covariates.We used the package "Distance" (Miller et al., 2019) in R (R Core Team, 2018) to model the probability of detection and estimate Sanje mangabey groups' density from distance data.We built and compared several model combinations using different detection functions and covariates, to find the best fitting detection function.The best model was selected using Akaike information criteria (AIC) where the best model was considered at ΔAIC = 0. We did not use adjustment terms as covariates were involved (Miller et al., 2019).Given the limited sample size (N = 28 listening posts), we only used elevation, ground slope, and distance to the closest reserve border as density covariates (see next section on assumptions of covariate relationship with density and occupancy).We averaged covariate value across a 1-km buffer around each listening post since Sanje mangabeys are generally heard within such area.Group densities were calculated for each study forest and then extrapolated considering the total available forest area (MW = 150.59km 2 and USNR = 314.48km 2 ; Marshall et al., 2010).

Occupancy modeling
Sanje mangabey occupancy (ψ) accounting for imperfect detection (p) was estimated using single-species singleseason occupancy modeling (MacKenzie et al., 2002).Occupancy is the proportion of sites where the target species is expected to occur when the likelihood to record it is smaller than 1, namely when the nondetection of the species does not imply its absence from the site (MacKenzie et al., 2002).We built a single model for both forests.The input data consisted of a matrix of sampled sites i by the sampling occasions j, filled with species detections (1) and nondetections (0).This consisted in a matrix of 117 sites and 137 sampling occasions.We scored as NA a site that was not sampled during a specific occasion.We first built a null model with constant detection and occupancy, and then selected the best supported model for the detection probability by fitting distance to the closest forest border and habitat type.In fact, previous studies assessed that the former could influence the detection probability of a range of species in either direction (Greco & Rovero, 2021;Rovero, Martin, et al., 2014), likely triggering increasing shyness or boldness.Moreover, detectability could also be affected by the density of forest floor vegetation in different habitat types, leading to changes in the efficiency of camera at detecting animals (Greco et al., 2021).We did not use camera-trap effort as potential covariate for the detection probability because our effort was homogeneous across cameras.We then used the best supported detection model (Appendix S1: Table S2) to model occupancy.Occupancy probability was modeled as a function of (1) habitat type, (2) elevation, (3) ground slope, (4) distance to the closest reserve border, and (5) distance to the closest river.We expected occupancy to increase away from the forests' borders, where anthropogenic disturbance may be higher, and that would increase in association with an optimal environment for the Sanje mangabeys, such as closer to watercourses, in association with broadleaves evergreen forests, and at mid-elevation.We also expected averaged site-specific occupancy to be higher in the national park compared with the nature reserve (Oberosler et al., 2020a;Paddock et al., 2020).Finally, we ranked all the occupancy models with different combinations of variables against the null model.Models were run using the package "unmarked" (Fiske & Chandler, 2011) in R. Model selection and ranking was based on the AIC score, with ΔAIC < 2 representing best supported models.If multiple models resulted top ranked, we followed Dormann et al.'s (2018) procedure to averaging predictions: we computed predictions on link scale for each best supported model, back-transformed them to response scale, and averaged predictions.Using the best models, we predicted occupancy across the forest areas by using a grid of 100 × 100 m pixels populated with site covariate values (Rovero & Spitale, 2016).This was done by using the "predict" function and plotted the result on a map using the "levelplot" function in the R package "lattice" (Sarkar, 2008).

Density to occupancy calibration
Given that listening post locations might not perfectly overlap with just one pixel of predicted occupancy, we averaged predicted occupancy values from the nine pixels centered on each listening post.We then regressed density values from survey sites on these predicted occupancy values to determine the relationships.We used generalized linear models with a gamma distribution and log-link function for each population within the two forests to test the relationship between the two metrics, with occupancy values as the independent variable and density values as the dependent variable (N = 13 sites for MW and N = 9 sites for USNR).Contrary to MW, USNR had sites with no acoustic records, and we assumed that these reflected nondetections rather than absence.Hence, we did not include them in the analyses.We chose a gamma distribution as the data are continuous and positive (Bolker, 2008) and because we hypothesized that the relationship might not be linear.This is because as occupancy increases along the 0-1 gradient, density may increase at a higher rate given it is not bounded.Based on the calibration results, we derived density (±SE) estimates from the occupancy values at the 100 × 100 m resolution, using the function "predict" in the package "stats" (Chambers & Hastie, 1992).We averaged density values to predict group density in each forest and multiplied them by the forest area to predict group abundance.Lastly, we derived individual density and abundance by using an average count of 39.2 individuals per group in MW (N = 5) and 31.7 in USNR (N = 3; Paddock et al., 2020).
Acoustic density estimations were best fitted by a half-normal detection function and with the distance to the closest reserve border as the variable retained (goodness of fit = 0.69; Appendix S1: Table S2).Occupancy values at the listening posts successfully predicted density estimates from acoustic surveys for both MW (β = 1.90 ± 0.92 SE; p = 0.03; R 2 = 0.32) and USNR (β = 2.49 ± 0.93 SE; p = 0.03; R 2 = 0.36; Figure 3).The spatially explicit densities derived from such relationships (Figure 4) were significantly higher in MW than in USNR (Welch two-sample t test: t = 16.65;df = 51,988; p value <0.001), and resulted in an estimated mean ± SE density of 0.26 ± 0.05 groups/km 2 and 0.24 ± 0.06, respectively.These estimates translated into 39.66 ± 8.28 mangabey groups across MW and 77.94 ± 18.01 across USNR.Predicted population size using available mean group counts resulted in 1555 ± 324.58 and 2471 ± 570.88 individual monkeys in MW and USNR, respectively, corresponding to an estimated individual density of 10.32 ± 2.16 individuals/km 2 in MW and 7.86 ± 1.82 in USNR (Table 3).

DISCUSSION
By integrating density estimation from acoustic surveys with spatially explicit occupancy predictions, we derived density with account for spatial heterogeneity across the entire range of the endangered Sanje mangabey, hence refining the population estimates by Paddock et al. (2020), which were based on acoustic surveys at a limited sample of sites.Occupancy modeling also allowed us to determine habitat associations across the species' range and assess potential vulnerability of this nonhuman primate species to anthropogenic disturbance.
The use of occupancy as a cost-effective population assessment tool when abundance estimation is complex has long been recognized (e.g., MacKenzie et al., 2006;Noon et al., 2012;Wilson & Schmidt, 2015) and rests on an assumed positive relationship between occupancy density (MacKenzie & Nichols, 2004).Such relationship is inherently complex because occupancy asymptotes at 1.In our study, the design of 2-km 2 grid cell size used for camera-trapping broadly matched the documented home range of Sanje mangabey (Ehardt et al., 2005;Mwamende, 2009).Hence, given that the grid cell size was close to the target species' home range, occupancy and abundance should broadly correspond (Linden et al., 2017), which likely underlines the occupancy-density relationship we found.Indeed, several studies highlighted how the spacing of detectors in relation to the home-range and movement pattern of the target species influenced the relationship (Rogan et al., 2019;Tempel & Gutierréz, 2013), with appropriate designs being those where such spacing matches the target species' home-range size, so to assume independence of detections (Clare et al., 2015;Tempel & Gutierréz, 2013;Wilson & Schmidt, 2015).In addition, the use of camera-trapping to estimate occupancy in continuous habitats is problematic due to the inherent violation of the closure assumption, given that camera traps are point detectors and animals regularly move in and out of sampled "sites" (Neilson et al., 2018).However, this issue makes estimated occupancy from camera detections a biased approximation of true occupancy, especially when low-density species move fast across large areas, hence saturating occupancy when density may be low (Lewis et al., 2015;Neilson et al., 2018;Parsons et al., 2017).
Findings from studies that investigated both forests found that the northern forest appears well protected compared with the southern one (Oberosler et al., 2020b), with the population in the nature reserve forest that has suffered decades of poor protection and, as a result, it is highly threatened and with lower density (Paddock et al., 2020).Our predictions of average group density per forest are broadly in line with those presented by Paddock et al. (2020) in that density was significantly higher in MW than in USNR; however, we estimated relatively higher density in USNR than from the acoustic distance sampling data alone (i.e., 0.22 vs. 0.15 groups/km 2 and T A B L E 1 Results of the occupancy models (top five and null model) for the Sanje mangabey (Cercocebus sanjei) in two forests, Mwanihana and Uzungwa Scarp Nature Reserve in Tanzania, detected by camera traps in 2017.

Occupancy models
No 69.6 vs. 45.9 groups for ours and the earlier study, respectively).This divergence may be caused by the spatially limited sample of listening posts used by Paddock et al. (2020) in comparison to the diffuse grid of camera traps we deployed, as well as the single acoustic survey event conducted for that study that may have caused a lower overall detection of the target species.The higher group density we estimated in USNR compared with previous studies translated into a higher overall abundance of approximately 4.026 ± 875.61 individuals in both forests (vs.3167 ± 436.62 estimated by Paddock et al., 2020), 61% of which are estimated to occur in USNR given its much larger area than MW.We caution that the conversion   from groups to individuals' abundance may be biased by the relatively small sample of group counts available (N = 8; Paddock et al., 2020) and needs therefore to be considered with this limitation in mind for future improvements.However, population estimates with methods that account for imperfect detection (Paddock et al., 2020) and with data from across the species' range (this study) were lacking (see review in Paddock et al., 2020); hence, our study provides for an important and comprehensive contribution to refine the species' conservation status and as a baseline for future studies.
The lower group density in USNR compared with MW mirrors results from earlier occupancy analyses from camera-trapping data (Hegerl et al., 2017;Oberosler et al., 2020aOberosler et al., , 2020b) ) and is likely explained by the elevated anthropogenic encroachment that has affected USNR for decades, mainly in the form of illegal hunting, tree cutting, and charcoal making (Oberosler et al., 2020a;Rovero et al., 2012).Such disturbance is seemingly more impactful on the strictly arboreal and folivorous colobine monkeys that underwent severe declines in USNR, being the target of selective hunting techniques (Rovero et al., 2012(Rovero et al., , 2015)), while the predominantly ground-dwelling and frugivorous mangabeys appear not to be specifically targeted by hunters and better adapted to exploit more heterogeneous and lightly disturbed habitats (Ehardt et al., 2005;Lloyd, 2017;Rovero, Martin, et al., 2014).Yet, their lower density in USNR than MW (with a seemingly lower group size in the former forest; Paddock et al., 2020) and the results of occupancy analyses suggest that the mangabey is vulnerable to human disturbance.In fact, the decreasing occupancy in proximity to protected area borders coincides in the insulated USNR with proximity to densely human-populated areas, which surround the entire forest block.Indeed, the grid of camera traps distributed across most of the species' range allowed us to study spatial variation in occurrence, and hence abundance.Elevation and gross habitat type emerged as important predictors consistently for both populations.Sanje mangabeys displayed a general preference for sub-montane forest at mid-elevation, with higher occurrence roughly between 500 and 1000 m asl.These areas of presumed optimal habitat for the species generally consist of evergreen forest with greater mean basal area than elsewhere (Cavada et al., 2016), indicating predominance of closed-canopy and old-growth forest.Nonetheless, our results indicated that deciduous forest zones were also used, likely reflecting the great niche flexibility of the Sanje mangabey, which is known to exploit mosaic forests, thanks to its diet predominantly based on seeds and invertebrates (Ehardt et al., 2005;Lloyd, 2017;McCabe et al., 2019;Rovero et al., 2012).
In conclusion, the approach we tested, based on calibrating occupancy to independent density estimations, appears valuable for spatially explicit density estimations of poorly detected species such as the Sanje mangabey, provided that the sampling design is appropriate to relate the two inherently different metrics.Our analysis provided for the first comprehensive assessments of the abundance and vulnerability of this primate, identifying priority areas for conservation, especially the forest areas at lower elevation and close to human settlements.While the relatively higher abundance of Sanje mangabeys in the less protected USNR than earlier reported gives cause for optimism, this forest continues to suffer from illegal human encroachment, calling for the need of increased protection and mitigation of anthropogenic threats.

F
I G U R E 1 Map of the study area in Tanzania (central panel and inset), with the two forest blocks (Mwanihana forest and Uzungwa Scarp Nature Reserve) where the Sanje mangabey (Cercocebus sanjei) occurs (left and right panels).In these two forests, black and white dots represent camera-trap locations (white are those where mangabeys were detected), whereas crosses represent the listening posts where acoustic surveys were made (yellow are those where mangabeys were detected).

F
I G U R E 2 Bivariate predictions of estimated occupancy probability with SE from the averaging of the predictions of the three best supported models for the Sanje mangabey (Cercocebus sanjei) in relation to: (a) elevation, (b) forest type, and (c) distance to the closest forest border.Additionally, (d) shows the box plot for site-specific occupancy values divided by forest, where central bold lines depict medians, boxes show interquartile range with lower and upper quartile values, and whiskers indicate minimum and maximum values.T A B L E 2 Parameter estimates from the averaging of the best supported models of Sanje mangabey (Cercocebus sanjei): occupancy (ψ) and detection (p) probability estimates across the two study areas (Mwanihana forest and Uzungwa Scarp Nature Reserve) in Tanzania.

F
I G U R E 3 Relationship between occupancy and acoustic density estimates derived in Mwanihana (left, N = 13 listening posts) and Uzungwa Scarp forests (right, N = 9), Tanzania, where the Sanje mangabey (Cercocebus sanjei) was studied.The graph represents the result of a generalized linear model with gamma distribution and log-link function.Shaded area represents SE.F I G U R E 4 Maps of the predicted group density of the Sanje mangabey (Cercocebus sanjei) in Mwanihana forest (left) and Uzungwa Scarp forest (right), Tanzania.
Francesco Rovero conceived the study.Ilaria Greco and Christina Lynette Paddock analyzed the data, with the contribution of Francesco Rovero.Christina Lynette Paddock, Gr ainne Michelle McCabe, Claudia Barelli, Steven Shinyambala, and Arafat S. Mtui conducted fieldwork.Ilaria Greco and Francesco Rovero wrote the manuscript and all authors contributed critically to the drafts.
Model ranking is based on the Akaike information criterion (AIC), with models having ΔAIC < 2.00 considered best supported (appearing in boldface). Note: T A B L E 3 Estimated group density, number of groups, population size, and individual density of the Sanje mangabey (Cercocebus sanjei) in Mwanihana forest (MW) and Uzungwa Scarp Nature Reserve (USNR), Tanzania.