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., Reference Estrada, Garber, Rylands, Roos, Fernandez-Duque, Di Fiore and Li2017). 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, Reference Cowlishaw and Dunbar2000; Estrada et al., Reference Estrada, Garber, Rylands, Roos, Fernandez-Duque, Di Fiore and Li2017). 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, Reference Walsh and White1999; Cavada et al., Reference Cavada, Ciolli, Barelli and Rovero2017), especially if the focal species are hunted and therefore wary of humans, or live in remote and inaccessible locations (Fashing & Cords, Reference Fashing and Cords2000). 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., Reference Rovero, Martin, Rosa, Ahumada and Spitale2014; indri Indri indri, Keane et al., Reference Keane, Hobinjatovo, Razafimanahaka, Jenkins and Jones2012; Alaotra reed lemur Hapalemur alaotrensis, Guillera-Arroita et al., Reference Guillera-Arroita, Lahoz-Monfort, Milner-Gulland, Young and Nicholson2010).
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., Reference Mackenzie, Nichols, Lachman, Droege, Royle and Langtimm2002), and as such is considerably less labour intensive and less costly to implement than count-based metrics (Royle & Nichols, Reference Royle and Nichols2003; Joseph et al., Reference Joseph, Field, Wilcox and Possingham2006). 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., Reference Margono, Potapov, Turubanova, Stolle and Hansen2014), creating a fragmented habitat mosaic that is probably restricting gene flow between subpopulations (Evans et al., Reference Evans, Supriatna, Andayani and Melnick2003). 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, Reference Mackinnon and Mackinnon1980; Sugardjito et al., Reference Sugardjito, Southwick, Supriatna, Kohlhaas, Baker, Erwin and Lerche1989; Lee, Reference Lee1997; Rosenbaum et al., Reference Rosenbaum, O'brien, Kinnaird and Supriatna1998; Melfi et al., Reference Melfi, Tasirin, Jonas, Yosep, Sabintoe, Houssaye and Kambey2007; Palacios et al., Reference Palacios, Engelhardt, Agil, Hodges, Bogia and Waltert2011; Kyes et al., Reference Kyes, Iskandar, Onibala, Paputungan, Laatung and Huettmann2013) and have estimated the population size to be 4,000–6,000 individuals (Riley, Reference Riley2010). 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, Reference Supriatna and Andayani2008). Macaca nigra is now categorized as Critically Endangered on the IUCN Red List (Supriatna & Andayani, Reference Supriatna and Andayani2008) and is one of the top 25 most threatened primates (Schwitzer et al., Reference Schwitzer, Mittermeier, Rylands, Chiozz, Williamson, Macfie and Cotton2017). 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., Reference Johnson, Hilser, Andayani, Hunowu, Linkie, Patandung and Bowkett2019). Within this area, we used knowledge of the habitat preferences of M. nigra (O'Brien & Kinnaird, Reference O'brien and Kinnaird1997; Rosenbaum et al., Reference Rosenbaum, O'brien, Kinnaird and Supriatna1998; Palacios et al., Reference Palacios, Engelhardt, Agil, Hodges, Bogia and Waltert2011) 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.
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 km2 (O'Brien Kinnaird, Reference O'brien and Kinnaird1997), but larger than the purported average of 2.16 km2 (Riley, Reference Riley2010). 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., Reference Rosenbaum, O'brien, Kinnaird and Supriatna1998). 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., Reference Venter, Sanderson, Magrach, Allan, Beher, Jones and Watson2016), 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., Reference Hansen, Potapov, Moore, Hancher, Turubanova, Tyukavina and Townshend2013), averaged across each site. All spatial datasets for covariates were calculated and derived in ArcGIS 10.2 (Esri, Redlands, USA).
Data analyses: occupancy modelling
We discretized our sampling data, following Rovero & Spitale (Reference Rovero, Spitale, Rovero and Zimmermann2016), 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, Reference Fiske and Chandler2011) 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., Reference Mackenzie, Nichols, Lachman, Droege, Royle and Langtimm2002). 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, Reference Royle and Nichols2003). 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. (Reference Mackenzie, Nichols, Royle, Pollock, Bailey and Hines2006), 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, Reference Burnham and Anderson2002) 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, Reference Burnham and Anderson2002) and estimates were instead derived by averaging coefficients from across these multiple models using the MuMIn package (Barton, Reference Barton2012) in R. We used a parametric bootstrap to check the adequacy of our model fit using the χ 2 approach on the most saturated model (Burnham & Anderson, Reference Burnham and Anderson2002).
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 (Reference Rovero, Spitale, Rovero and Zimmermann2016). 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 km2 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, Reference Mackenzie and Royle2005):
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., Reference Beaudrot, Ahumada, O'Brien and Jansen2018), 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., Reference Mackenzie, Nichols, Hines, Knutson and Franklin2003). 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. (Reference Beaudrot, Ahumada, O'Brien and Jansen2018) and all simulations and analyses were done in R, using code available on Github (Ahumada, Reference Ahumada2017).
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 km2), we estimated the area of occupancy of M. nigra to be 2,101 km2. Finally, our goodness-of-fit test (χ 2 = 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 km2, we identified eight distinct regions that probably contain important subpopulations (Fig. 3, Supplementary Table 3, Supplementary Fig. 2).
1 PA, protected area; HFI, Human Footprint Index.
2 Number of parameters.
3 Akaike information criterion corrected for small sample size.
4 Relative difference in AICc values compared to top-ranked model.
5 AICc model weight.
6 Estimate of occupancy ± SE.
7 Cumulative AICc model weight.
8 Estimate of detection probability ± SE.
1 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.
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 km2) 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, Reference Royle and Nichols2003). 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., Reference Myneni, Hall, Sellers and Marshak1995) 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., Reference Hansen, Potapov, Moore, Hancher, Turubanova, Tyukavina and Townshend2013), 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., Reference Estrada, Garber, Rylands, Roos, Fernandez-Duque, Di Fiore and Li2017), 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., Reference Johnson, Hilser, Andayani, Hunowu, Linkie, Patandung and Bowkett2019), this finding is of conservation relevance and contradicts speculations that Tangkoko contains the last viable population (Supriatna & Andayani, Reference Supriatna and Andayani2008; Palacios et al., Reference Palacios, Engelhardt, Agil, Hodges, Bogia and Waltert2011; Kyes et al., Reference Kyes, Iskandar, Onibala, Paputungan, Laatung and Huettmann2013; Engelhardt et al., Reference Engelhardt, Muniz, Perwitasari-Farajallah and Widdig2017). Although Tangkoko does appear to hold an important and viable population (Engelhardt et al., Reference Engelhardt, Muniz, Perwitasari-Farajallah and Widdig2017) 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 km2 recently reported (Johnson et al., Reference Johnson, Hilser, Andayani, Hunowu, Linkie, Patandung and Bowkett2019), results in an area of occupancy of 2,101 km2 and an extent of occurrence of 7,810 km2. 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., Reference Yoccoz, Nichols and Boulinier2001; Legg & Nagy, Reference Legg and Nagy2006).
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., Reference Guillera-Arroita, Lahoz-Monfort, Milner-Gulland, Young and Nicholson2010; Keane et al., Reference Keane, Hobinjatovo, Razafimanahaka, Jenkins and Jones2012; Gerber et al., Reference Gerber, Williams and Bailey2014; Rovero et al., Reference Rovero, Mtui, Kitegile, Jacob, Araldi and Tenan2015).
Acknowledgements
We thank the Ministry of Science and Technology for permission to conduct this research in Indonesia; Bogani Nani Wartabone National Park management authority for research permission; BKSDA North Sulawesi, in particular Pak Yakub Ambagau, for support of our research; our field assistants Nicodemus Malir, Hanisa Silayar, Fandy Pongantung, Adha Yakseb, Omega Judistira, Idham Pohi, Muhamad Ihwan, Sukrisno and Luciano Gawina; Graeme Gillespie for providing advice during the conception of this project; and two anonymous reviewers for their useful comments. Funding was provided by the Global Environmental Facility, as administered through the project Enhancing the Protected Area System of Sulawesi, Fondation Segré and Wildlife Reserves Singapore Conservation Fund.
Author contributions
Study design and fieldwork: CJ, RR, LN, HH, AP, WP, IH; data analysis: CJ, FR; writing: CJ, AB, ML, HH, FR.
Conflicts of interest
None.
Ethical standards
This research abided by the Oryx guidelines on ethical standards.