Using occupancy-based camera-trap surveys to assess the Critically Endangered primate Macaca nigra across its range in North Sulawesi, Indonesia

Abstract

Primates are one of the most threatened groups of mammals. Understanding their patterns of population occurrence and abundance, especially in response to threats, is critical for informing conservation action. The crested black macaque Macaca nigra is the only Critically Endangered species of Sulawesi's seven endemic macaques. Little is known about its distribution or its response to deforestation and hunting. We conducted a camera-trap survey across the entire species range using an occupancy-based analytical approach to (1) establish the first range-wide baseline of occurrence, (2) investigate how environmental and anthropogenic factors influence occurrence, (3) identify priority conservation subpopulations, and (4) test the efficacy of the sampling and analytical protocol for temporal monitoring of M. nigra using occupancy as the state variable. From 9,753 camera-trap days, M. nigra was detected on 473 days at 77 of the 111 camera locations. Species occupancy was 0.66 and highest inside protected areas and closed canopy forest. We identified eight distinct subpopulations, based on distribution and forest fragment size. To inform future monitoring, we used a power analysis to determine if our effort would allow us to detect inter-annual occupancy declines of 10%, and found that 90 camera locations surveyed for 3 months (8,100 camera days) across three consecutive seasons is the effort required to detect such change with 80% certainty. Our study underscores the importance of well-managed protected areas and intact forests for the long-term survival of the crested black macaque, and tests the effectiveness of camera traps to monitor primates at the landscape scale.

Introduction

Globally, an estimated 60% of primate species are at risk of extinction (i.e. categorized as threatened on the IUCN Red List), with 75% of all primate species being in a state of decline (Estrada et al., 2017). The threats to primate species are broad and wide-ranging, with many acting in synergy to hasten population declines, but hunting for trade and habitat loss from logging and agricultural expansion are the primary drivers (Cowlishaw & Dunbar, 2000; Estrada et al., 2017). To understand the effects of threats on primate populations, bespoke monitoring programmes are critical.

Monitoring of primate populations provides forewarning of declines and, by recognizing causal factors, a good monitoring programme can identify mitigating solutions. However, monitoring primates using abundance metrics can be resource intensive and time consuming at large spatial scales (Walsh & White, 1999; Cavada et al., 2017), especially if the focal species are hunted and therefore wary of humans, or live in remote and inaccessible locations (Fashing & Cords, 2000). To address this, several recent studies have recommended using occupancy as the state variable, rather than abundance, as a more cost-effective option (Sanje mangaby Cercocebus sanjei, Rovero et al., 2014; indri Indri indri, Keane et al., 2012; Alaotra reed lemur Hapalemur alaotrensis, Guillera-Arroita et al., 2010).

Occupancy modelling corrects for detectability issues by using presence/absence (detection/non-detection) data to estimate the proportion of sites occupied by the species of interest (Mackenzie et al., 2002), and as such is considerably less labour intensive and less costly to implement than count-based metrics (Royle & Nichols, 2003; Joseph et al., 2006). To this end, camera traps produce data that are suited to occupancy analysis and therefore present a viable survey method, especially for elusive primates living in landscapes with difficult terrain. Here, we use an occupancy based camera-trap approach to conduct a baseline assessment for the crested black macaque Macaca nigra, a hard-to-detect and mostly terrestrial species of forest-dwelling primate endemic to the northernmost peninsula of the Indonesian island of Sulawesi.

The crested black macaque exemplifies the plight of many forest-dwelling primates. Resource extraction and land requirements for agriculture in Sulawesi have resulted in the loss of most of its lowland forest habitat (Margono et al., 2014), creating a fragmented habitat mosaic that is probably restricting gene flow between subpopulations (Evans et al., 2003). Added to this are retaliatory killings in response to foraging in crops, and hunting for bushmeat consumption, which in turn makes the species shy and secretive in large parts of its range.

Numerous presence-only studies have explored the distribution and status of M. nigra (Mackinnon & Mackinnon, 1980; Sugardjito et al., 1989; Lee, 1997; Rosenbaum et al., 1998; Melfi et al., 2007; Palacios et al., 2011; Kyes et al., 2013) and have estimated the population size to be 4,000–6,000 individuals (Riley, 2010). Sharp declines of up to 80% over 40 years in M. nigra populations have indicated the seriousness of the threats facing the species (Supriatna & Andayani, 2008). Macaca nigra is now categorized as Critically Endangered on the IUCN Red List (Supriatna & Andayani, 2008) and is one of the top 25 most threatened primates (Schwitzer et al., 2017). Although previous work has highlighted the vulnerability of M. nigra, the limited survey efforts yielded incomplete knowledge on distribution, status and responses to anthropogenic threats at a regional scale. In this study, we use camera-trap data to (1) investigate the influence of environmental and anthropogenic factors on M. nigra occupancy across its native range (2) provide an empirical baseline of its status, using occupancy as the state variable of interest, and (3) evaluate the efficacy of our protocol to monitor primates at the landscape scale.

Study area

Macaca nigra is native to North Sulawesi province, where its range extends from the northern tip and some of the small surrounding islands (Lembeh, Manadotua and Talise) to the Onggak-Dumoga River and the south-east landscape of Bogani Nani Wartabone (Johnson et al., 2019). Within this area, we used knowledge of the habitat preferences of M. nigra (O'Brien & Kinnaird, 1997; Rosenbaum et al., 1998; Palacios et al., 2011) and recent landcover maps (KLHK, 2015) to determine the extent of potential habitat that would be ecologically capable of supporting a macaque population and use this to define our study area (Fig. 1). These included sites characterized by forest cover and/or scrub habitats. North Sulawesi has altitudes of 0–1,995 m and has a wet season during October–May and a dry season during June–September.

Fig. 1 Camera-trap placements across North Sulawesi province in relation to protected and unprotected areas and the main habitat types. The line represents the boundary between Macaca nigra and Macaca nigrescens (Johnson et al., 2019).

Methods

Data collection: camera trapping

We first divided the potential habitat into sampling units (termed sites) of 2 × 2 km (n = 796 sites). This size was chosen as it is close to the largest recorded M. nigra home range of 4.06 km² (O'Brien Kinnaird, 1997), but larger than the purported average of 2.16 km² (Riley, 2010). A sample of 115 of these sites was then selected to be surveyed using camera traps (64 Reconyx HC500, Holmen, USA; 50 Cuddeback Black Flash, DePere, USA; one Bushnell Natureview HD Live View, Overland Park, USA). This was done in two phases. Firstly, during January 2016–October 2017, we surveyed 67 sites. This included all 28 sites across the Tangkoko Nature Reserve, and 39 sites selected with random systematic sampling within the larger region of the southern Bogani Nani Wartabone National Park landscape (Fig. 1). Secondly, during March–July 2018, we randomly selected 48 sites from across patches of potentially suitable habitat between these protected areas. Although implemented in different years because of the number of cameras available, we conducted the actual survey phases within a relatively short period of time, minimizing the risk of violating the assumption that the population was closed.

Because of logistical constraints, camera traps were not uniformly deployed within each site. Rather, we deployed cameras in the part of the site most accessible, provided it was > 50 m from the site's border. Our grid-based approach ensured sufficient spacing between camera traps. We positioned cameras along wildlife trails, attaching them to tree trunks at a height of c. 50 cm, with the sensor aimed parallel to the ground facing the monitoring area. On average each camera was deployed for a minimum of 3 months. Four locations, however, produced no data as the cameras were either stolen or malfunctioned; we therefore present data collected from 111 camera traps.

Data collection: covariates

Occupancy of sites by M. nigra was expected to vary across North Sulawesi according to habitat and anthropogenic factors (Rosenbaum et al., 1998). We therefore derived a spatial dataset of covariates that could potentially explain any heterogeneity in M. nigra occupancy. These were presence of forest cover (Yes/No), elevation, slope, normalized difference vegetation index (NDVI), distance to roads and settlements, protected area (Yes/No), level of anthropogenic disturbance in the landscape (measured by the Human Footprint Index; Venter et al., 2016), and distance from the edge of continuous forest (Table 1). All covariates were extracted for the exact camera location, except for NDVI, slope and elevation. As a measure of landscape roughness, and therefore accessibility, we calculated the mean of slope and elevation across each site. We calculated NDVI using red and near infrared spectral bands in Landsat 8 imagery (Hansen et al., 2013), averaged across each site. All spatial datasets for covariates were calculated and derived in ArcGIS 10.2 (Esri, Redlands, USA).

Table 1 Definition and predicted effect of covariates used to model variation in detectability and occupancy of Macaca nigra across sites.

Data analyses: occupancy modelling

We discretized our sampling data, following Rovero & Spitale (2016), into 24 consecutive 5-day sampling occasions (of 1, 5 and 10 days, a 5-day duration best facilitated model convergence). For each site and each occasion, a 1 indicated detection (a photograph) and 0 indicated non-detection (no photograph) of M. nigra.

To assess the influence of ecological and anthropogenic factors on occupancy, we fitted single-season occupancy models to the data using the package unmarked (Fiske & Chandler, 2011) in R 3.5.1 (R Development Core Team, 2010). Single-season occupancy models estimate both the probability that a site is occupied (ψ) and the probability that the species is detected if it is present (p) using a maximum likelihood approach (Mackenzie et al., 2002). All numerical covariates were first standardized into z-scores and assessed for collinearity using Pearson's rank correlation. As no covariates were found to exhibit collinearity > 0.6, no covariates were excluded in any of the candidate models (Supplementary Table 1).

We next examined the effect of covariates (Table 1) on our parameter of interest, ψ. However, we expected that a number of covariates that influence occupancy of M. nigra may also affect the abundance of M. nigra and therefore detectability (Royle & Nichols, 2003). Additionally, there are other factors that could also influence the detectability. We therefore began our model selection process by explicitly accounting for p. Following recommendations of MacKenzie et al. (2006), we identified a suitable covariate structure for p whilst holding the covariate structure for ψ constant. Covariate structure for p was assessed using the Akaike information criterion corrected for small sample sizes (AICc). We then fixed the covariate structure for p with the covariate structure that had the lowest AICc (Burnham & Anderson, 2002) and proceeded to assess the role of our covariates on ψ, ranking the competing models according to their AICc values. If multiple models were within 2 ΔAICc points, they were considered equally supported (Burnham & Anderson, 2002) and estimates were instead derived by averaging coefficients from across these multiple models using the MuMIn package (Barton, 2012) in R. We used a parametric bootstrap to check the adequacy of our model fit using the χ ² approach on the most saturated model (Burnham & Anderson, 2002).

Once the best (or averaged) occupancy model was identified, we used the corresponding occupancy probability to estimate the proportion of the 796 sites likely to be occupied by M. nigra. We then mapped predicted occupancy probability of M. nigra across its native range using the approach suggested by Rovero & Spitale (2016). This approach involved fitting our best model to a data frame containing a grid of covariate values at the scale of 1 × 1 km across the range of M. nigra. Maps of estimated occupancy were created in ArcMap 10.2 (Esri, Redlands, USA).

Finally, we used predicted occupancy to identify regions that may be important to the species based on the presence of characteristics associated with a high occupancy probability. We classified these as spatially distinct areas with > 50 km² of continuous forest cover and predicted occupancy > 0.7. We defined spatially distinct as being separated by > 100 m of habitat not capable of supporting the species. Although separation of forest patches as we define it here does not necessarily preclude movement of M. nigra individuals between populations, it enables easy demarcation of distinct forest blocks that may be exposed to different threats.

Data analyses: optimal study design and power analysis

We evaluated the optimal effort for an occupancy-based camera-trap monitoring protocol for M. nigra by first estimating the number of sites (s) and occasions (K, 5-day sampling periods) required to achieve a desired level of precision for the ψ and p estimators. These were calculated using estimates from the best model and the corresponding asymptotic variance of the occupancy indicator (Equation 1; Mackenzie & Royle, 2005):

(1)$$TS = s \times K = \displaystyle{{K\times {\rm \psi} } \over {var\lpar {\hat{\rm \psi}} \rpar }}\left[{\lpar {1- {\rm \psi} } \rpar + \displaystyle{{\lpar 1-p^\ast \rpar } \over {\,p^\ast {-}\, K \times p \times {\lpar {1-p} \rpar }^{K-1}}}} \right]$$

where TS is total effort, s is the number of sites, K is the number of occasions, ψ the probability of occupancy, p the probability of detection provided the site is occupied, and p* the probability that the species is detected at least once after K occasions. As the number of occasions can be increased free of additional cost in a camera-trap survey (at least until a camera needs to be serviced, in this case after c. 3 months, or K = 18) we defined the optimum monitoring effort for M. nigra as one that achieves the target standard error of 0.05 in our estimates whilst minimizing s. We did this by resolving the Equation for s and assuming a duration of 3 months for the surveys.

As smaller absolute declines in occupancy become more detectable as the number of seasons of monitoring increases (Beaudrot et al., 2018), we calculated the number of seasons needed to detect a 10% annual occupancy decline if the optimal monitoring effort (calculated by Equation 1) is chosen as the long-term monitoring protocol (survey effort = s × K). To do this, we simulated data for an annual 10% decline across a time series. We used the estimates from our best model to determine the initial input parameters for occupancy and detection and set colonization to zero, so as to simulate a decline. We then took the generated data and fitted it to a dynamic occupancy model without covariates (Mackenzie et al., 2003). These were then assessed to determine the number of seasons required to detect an annual decline of 10% with 80% confidence. Our methods follow those detailed by Beaudrot et al. (2018) and all simulations and analyses were done in R, using code available on Github (Ahumada, 2017).

Results

Detection and occupancy

From 9,749 camera-trap days, M. nigra was detected on 473 separate days in 71 of 111 sites, yielding a naïve occupancy estimate of 0.64. These 71 camera locations confirmed the presence of the species in 12 spatially distinct forest fragments, distributed across North Sulawesi, indicating a greater distribution than previously known.

Investigating covariate influence on detection probability, whilst holding occupancy constant, revealed p as a positive function of being inside both a protected area and forest, and a negative function of the Human Footprint Index (Supplementary Table 2, Supplementary Fig. 1). This model gave a detection probability of 0.28 ± SE 0.025 and was subsequently used to explore which combination of covariates best explained M. nigra occupancy.

From the different combinations of covariates, the most parsimonious occupancy model also included protected area and forest (Table 2). M. nigra occupancy was higher in forests (β coefficient 1.14 ± SE 0.5) and inside protected areas (0.69 ± SE 0.13). However, five models were ranked < 2 ΔAICc units and thus considered equally supported (Table 3, Fig. 2). We therefore estimated detectability and occupancy by averaging across these models, resulting in p = 0.66 and ψ = 0.28 (Table 2). By taking this estimate for occupancy and extrapolating it to all sites across the landscape (n = 796 or 3148 km²), we estimated the area of occupancy of M. nigra to be 2,101 km². Finally, our goodness-of-fit test (χ ² = 1869, P = 0.65) indicated no evidence for a lack of fit for the most saturated model, suggesting that the more parsimonious models provide an adequate description of the data. According to our criteria of predicted occupancy > 0.7 and continuous forest cover > 50 km², we identified eight distinct regions that probably contain important subpopulations (Fig. 3, Supplementary Table 3, Supplementary Fig. 2).

Fig. 2 The functional relationships of predicted occupancy, with standard error, of M. nigra, with covariates included, in the best-fit models (Table 3): (a) protective area status, (b) forest cover (bush/scrub or closed canopy forest), (c) Euclidean distance from forest edge (negative values for cameras in non-forested habitat), and (d) Normalized Difference Vegetation Index (NDVI).

Fig. 3 Predicted occupancy across the mainland native range of M. nigra. Eight potentially important landscapes, defined as spatially distinct continuous forest blocks > 50 km² with predicted occupancy > 0.7, are identified by numbers (Supplementary Fig. 2, Supplementary Table 3).

Table 2 Model selection results for covariate effects in determining occupancy probability of M. nigra across its native range in North Sulawesi. The top ranked models are shown as those with ΔAIC c < 3 followed by the null and average model. All models depicted (ΔAICc < 2) were included in the averaged model.

¹ PA, protected area; HFI, Human Footprint Index.

² Number of parameters.

³ Akaike information criterion corrected for small sample size.

⁴ Relative difference in AICc values compared to top-ranked model.

⁵ AICc model weight.

⁶ Estimate of occupancy ± SE.

⁷ Cumulative AICc model weight.

⁸ Estimate of detection probability ± SE.

Table 3 Summary of conditional model averaged parameters based on the best-supported models identified in Table 2. Estimates of the β coefficient are reported for standardized covariates (scaled to mean = 0 and unit variance of 2). See Table 2 for the models with covariates for both ψ and p.

¹ PA, protected area; HFI, Human Footprint Index.

* Significant.

Optimal study design and power analysis

Considering the estimates from our averaged model and assuming constant probabilities, an optimal survey design would have 123 camera-trap sites, each with 8 repeats (1.5 months when each repeat constitutes a duration of 5 days) to achieve a target standard error of 0.05 (Fig. 4a). However, for camera traps further repeats are easily gained, without incurring further cost, by simply delaying the retrieval of the cameras. Therefore, it is possible to reduce the number of sites needed by increasing the number of repeats, without additional cost, and thereby reducing overall survey costs. Under this option, if the number of repeats is increased to 18 (90 days, or c. 3 months), the number of camera-trap sites could be reduced to 90 (Fig. 4b). As the cheapest means of surveying M. nigra with suitable precision, this is the protocol that we advocate. Simulating the power provided by this protocol, if applied across multiple seasons, 3 years would be sufficient to detect a 10% decline in occupancy with 80% certainty.

Fig. 4 Total survey effort required (sites × repeats) (a), and total number of unique camera sites required (b) to achieve a given precision in the occupancy estimate as a function of the number of repeats required per camera. Repeats are represented by months; a single repetition has a duration of 5 days and 1 month is therefore six repeats. Curves were created based on the averaged model (p = 0.28 and ψ = 0.66) and on the asymptotic properties of the maximum likelihood estimates. The curves correspond to different standard errors in the occupancy estimate.

Discussion

Our study demonstrates an extension in the known range of M. nigra compared to previous data. We show that camera-trap data coupled with occupancy analysis can provide a robust baseline assessment of a predominantly ground-dwelling forest primate of conservation concern. We estimated an occupancy of 0.66 (2,101 km²) for M. nigra, which represents the first baseline for future regional monitoring of the species. Estimated occupancy was only slightly higher than the naïve occupancy (0.64), which is probably an artefact of a long survey duration, which resulted in a high cumulative detection probability, thereby reducing its influence on occupancy. Detection probability, however, varied across sites, and was higher in closed canopy forests and in areas with a lower Human Footprint Index. We attribute this to a preference by M. nigra for forests and a sensitivity to human disturbance, which cause a reduced abundance of the species in non-forested and disturbed habitats, with less abundant species typically being less detectable (Royle & Nichols, 2003). However, this could also be a consequence of M. nigra being more elusive and cautious in more disturbed areas.

Although occupancy was significantly higher in closed canopy forest than in bush/scrub habitats, we were not able to test for effects of forest type (pristine vs degraded/logged) because of the limited resolution of available satellite imagery. We used NDVI as a remotely sensed proxy for vegetation structure and green biomass (Myneni et al., 1995) but, although retained in the best models, it was not significant. This is probably the result of the overriding importance of forest presence over the quality of that forest in determining occupancy. The dependence of M. nigra on forest is nevertheless apparent, highlighting the species’ susceptibility to forest conversion to other land uses.

Along with a preference for forest cover, estimated occupancy was higher within protected areas, which is probably a result of both the continuous extent of optimal habitat and lower anthropogenic disturbance. Deforestation observed within protected areas during 2001–2017 was much lower than outside (26.4% outside protected area vs 6.2% loss within protected area; Hansen et al., 2013), a factor that probably implies reduced hunting and disturbance. However, protected areas currently cover only 20% of the area estimated to be occupied by the species. Considering that M. nigra, like many other primates (Estrada et al., 2017), is vulnerable to anthropogenic activities, the species remains at risk across 80% of its range.

Our modelling approach predicted eight spatially distinct regions that are likely to support important subpopulations of M. nigra (Fig. 3, Supplementary Table 3, Supplementary Fig. 2). Five of these landscapes represent the first scientific record of the species' presence (Supplementary Table 3). In combination with evidence of range expansion (Johnson et al., 2019), this finding is of conservation relevance and contradicts speculations that Tangkoko contains the last viable population (Supriatna & Andayani, 2008; Palacios et al., 2011; Kyes et al., 2013; Engelhardt et al., 2017). Although Tangkoko does appear to hold an important and viable population (Engelhardt et al., 2017) it comprises just 3% of the total area occupied by the species (Supplementary Table 3). Therefore, although we acknowledge that population density and viability of subpopulations will vary across the species' range, the significance of these additional areas for the species cannot be ignored.

Our discovery of previously undocumented populations, in addition to the increased range size of 220 km² recently reported (Johnson et al., 2019), results in an area of occupancy of 2,101 km² and an extent of occurrence of 7,810 km². The Critically Endangered status of M. nigra was a result of suspected severe population declines. As we only provide a baseline, using new survey methodology, we can neither support nor refute this. However, with a greater range than previously thought, our findings may influence the threatened status of the species and therefore have implications for any future review of the species' Red List status. With this in mind, we caution that the species remains highly threatened, particularly outside protected areas, and emphasize the need for standardized monitoring that follows the protocol we recommend. This will identify trends in the population and lead to a clearer picture of the species' status. In the meantime, we stress the need for increased conservation management in the landscapes identified as important (Supplementary Table 3).

In addition to contributing valuable information on the distribution and status of the species, our study facilitates the design of a targeted monitoring protocol, which we define as 90 sites surveyed with camera traps for 90 days (18 × 5-day occasions) per season. Testing the power of this design under simulation, we found it would be possible to detect an annual change in occupancy of 10% (a decline we consider sufficiently severe to initiate a conservation intervention), with 80% certainty if implemented over 3 seasons. By demonstrating the robustness of our suggested protocol in meeting the objectives of such a monitoring programme, researchers can be confident in investing resources in its long-term implementation (Yoccoz et al., 2001; Legg & Nagy, 2006).

Meaningful monitoring studies, especially at large spatial scales, are considerable investments of limited conservation resources. Procurement of the cameras required for a study such as this constitutes a high cost, which must be justified given the availability of other presence/absence survey methods such as line transect surveys. We were able to use data collected during this study to conduct a basic feasibility comparison between the two methods and found that the cost efficiency of camera traps was far greater than that of a more labour intensive transect method, even after considering the cost of cameras (Supplementary Table 4). As such, although initial costs are high, camera traps may provide a cheaper alternative for monitoring certain species.

In conclusion, our study provides a robust baseline for M. nigra occupancy and, through simulation, indicates how the monitoring of a generally elusive and Critically Endangered primate through camera trapping could be feasible at a landscape level. We emphasize the necessity of such a monitoring approach for understanding population responses to ecological and anthropogenic factors and hence informing conservation efforts. The baseline we report forms the basis on which temporal trends in M. nigra occupancy could be detected, and the study design we suggest provides a robust means of detecting potential declines within a suitable timeframe. In demonstrating the potential of camera-trap data analysed in an occupancy framework to monitor M. nigra, we add to a growing body of literature that promotes occupancy as a feasible alternative to abundance-based metrics for monitoring wild primates (Guillera-Arroita et al., 2010; Keane et al., 2012; Gerber et al., 2014; Rovero et al., 2015).