Introduction

At a global scale, increasing deforestation rates and other anthropogenic drivers of landscape transformation cause habitat loss and fragmentation (Laurance et al. 2014). This imposes immense pressure on wildlife populations, especially for highly specialized species (Frank and Amarasekare 1998) or those with large home ranges and dispersal abilities, such as medium- or large-bodied felids (Sunquist and Sunquist 2019). Species adaptability to these changes, and populations’ survival, often depends on how easily individuals can move through the landscape in search of food, better quality habitat and mates (Zeller et al. 2012). Small and isolated wildlife populations, without the possibility to disperse and increase their gene pool, are those most threatened by extinction risk (Smith 1993; Wikramanayake et al. 2011; Goossens et al. 2016). Therefore, it is crucial to identify and protect highly utilized habitats, as well as dispersal corridors linking population strongholds. Only well-connected landscapes with large enough core habitat can facilitate frequent gene flow, leading to genetically diverse and stable populations (Manel et al. 2003; Storfer et al. 2007; Thapa et al. 2018). Understanding and mapping population connectivity and species dispersal routes are further important in the current context of globally increasing human-wildlife conflict (Cushman et al. 2018). Habitat fragmentation and increased poaching pressure may push individuals to the periphery of protected areas or even to move out into human-dominated lands, elevating the likelihood of conflict (Nyhus and Tilson 2004; Macdonald et al. 2012).

When assessing local or regional landscape connectivity, either for human-conflict mitigation or conservation planning, it is crucial to base assessments on rigorous empirical data and reliable modelling methods. This is particularly important for Southeast Asia, which has one of the highest rates of biodiversity loss worldwide, primarily caused by forest conversion to smallholder farmland and large industrial monoculture plantations (Miettinen et al. 2011; Gaveau et al. 2016).

Due to challenging landscapes characteristics and thick canopy cover, movement data for target species are limited (Grassman et al. 2005; Mohamad et al. 2015; Hearn et al. 2019). Therefore, we need robust data for creating landscape resistance layers for defining corridors and core areas. Mateo-Sánchez et al. (2015) and Keeley et al. (2016) found that, in the absence of data on movement behaviour or population genetic structure, which is the case for our three study species in Sumatra, spatial information on habitat use can be a useful surrogate for deriving a resistance to movement layer. Mateo-Sánchez et al. (2015) and Keeley et al. (2016) showed that a negative exponential function best describes the relationship between habitat suitability and resistance values. Unlike many other connectivity modelling approaches, both cumulative resistant kernels and factorial least-cost paths, as implemented in UNICOR (Landguth et al. 2012), account for species dispersal abilities, which is crucial for accurately predicting landscape-scale connectivity patterns (Cushman et al. 2013a, b). Dispersal abilities are not well known for our focal species; however, some studies have defined relationships between maximum species dispersal distances and home range size (Bowman et al. 2002; Whitmee and Orme 2013).

Current understanding of wild cat ecology in the Kerinci Seblat (KS) landscape, and in Sumatra in general, is mostly limited to the critically endangered Sumatran tiger (Panthera tigris sumatrae). Several studies inform on the population status of mesopredators in Sumatra - clouded leopard (Neofelis diardi diardi) (Haidir et al. 2018, 2020), Asiatic golden cat (Catopuma temminckii) and marbled cat (Pardofelis marmorata) (McCarthy 2013; Pusparini et al. 2014; Sunarto et al. 2015). However, a recent study by Struebig et al. (2018), although focused on the Sumatran tiger, suggested that small-medium sized felids moving through the Kerinci Seblat landscape may encounter more interactions with humans than do tigers, again highlighting the knowledge gap on species dispersal patterns and meta-population connectivity. Therefore, applying a multi species approach focused especially on mesopredator species with diverse body sizes, ranges and habitat requirements (Schuette et al. 2013; Lesmeister et al. 2015; Moreira-Arce et al. 2016) can provide protected area managers with crucial information to prioritize management scenarios, leading to better assessments of particular interventions (Sauer et al. 2013).

Within Southeast Asia, Indonesia is reported to have lost six million hectares of primary lowland forest from 2000 to 2012, equivalent to 470,000 ha per year (Margono et al. 2014). The most recent report on national deforestation by the Indonesian Ministry of Environment and Forestry recorded a loss of 223,000 ha in 2018 (KLHK 2019), with Sumatra accounting for 25% (59,000 ha) of this loss. In Sumatra, the Kerinci Seblat (KS) landscape is a biodiversity stronghold. Within the landscape, an intact area covering 1.39 million hectares, the Kerinci Seblat National Park (KSNP) is one of the largest protected areas in Southeast Asia, but its elongated shape makes it susceptible to deforestation pressures. These include complete forest clearance resulting from smallholder land conversion to agriculture, which is precipitated by the creation of logging roads and presence of a large-scale road network that increases access to remoter forest areas (Linkie et al. 2006; Gaveau et al. 2009; Margono et al. 2014). Thus, empirically-based detection of important core habitats and ecological corridors for felid species of conservation priority in tropical landscapes is urgent (Linkie et al. 2006).

This study surveyed the main forest habitat types in the Kerinci Seblat landscape aiming to: (1) spatially predict habitat use patterns of clouded leopard, Asiatic golden cat and marbled cat using single-species occupancy models; (2) model, map and quantify high-density movement and dispersal corridors for the three focal species; and, (3) model deforestation risk, a main driver of habitat loss and fragmentation in KS landscape, to determine its potential effects on population connectivity of the studied felids.

Methods

Study area

The study area encompasses 16,000 km2 Kerinci Seblat landscape, which stretches across the west-central Sumatran section of the Bukit Barisan mountain range that runs the length of the island (Fig. 1). The landscape consists of 15 districts that are predominantly covered by Kerinci Seblat National Park (KSNP) and Batanghari Protection Forest, and other land-use types: dry agricultures, rubber, coffee, palm oil, and cocoa plantations, and paddy field across four provinces: Jambi, West Sumatra, South Sumatra and Bengkulu. The Kerinci Seblat landscape is listed as a National Strategic Area under Indonesian Government Act No. 26/2008 on National Spatial Planning because of its high environmental and biodiversity values (MoPWH 2017). The three national parks of Kerinci Seblat, Bukit Barisan Selatan and Gunung Leuser form the Tropical Rainforest Heritage Site, a natural UNESCO World Heritage Site (IUCN) (IUCN 2004).

Fig. 1
figure 1

Study areas in Kerinci Seblat landscape, black dots indicate camera trap locations with study area names and district capitals

Full size image

Camera trap surveys

The camera trap surveys sampled seven study sites inside and adjacent to KSNP, spanning from the north to the southern-most extent of the KS landscape: Kambang (KM), Bungo (BG), Muara Hemat (MH), Sipurak (SP), Renah Kayu Embun (RKE), Ipuh (IP) and Karang Panggung (KP) - inside and adjacent to KSNP (Fig. 1 and Table S1). The four surveys in BG, SP, RKE and IP aimed at repeating camera trapping previously undertaken in 2004 and 2010 (see Linkie et al. 2008; and Wong et al. 2013). These sites cover the main forest and land-use types of the landscape. However, camera trap deployments in MH (n = 143 single camera placements), KM (130) and KP (106) aimed at sampling the forest-farmland interface, using a strip-shaped camera trap polygon (15–18 km long, 3–5 km wide and spanning 27–32 km2). The distance between camera trap stations ranged from 0.4 to 0.7 km.

Surveys in BG, SP, RKE and IP covered the same areas studied by Haidir et al. (2018). A total of 292 camera trap stations, using a combination of Cuddeback Ambush IR (Non Typical Inc., WI, USA) and Panthera IV camera trap units (Panthera Foundation, NY, USA), were set with gaps ranging from 0.8 to 1.4 km, covering a roughly circular area of 60–70 km2. At each trap station, paired cameras were set except for RKE where only ~ 75% were in pairs. Each camera was placed on a pole/tree next to forest trails at a height of 40–60 cm above the ground, c.a. 2–2.5 m from the centre of the trails. No bait or lure was applied at trap sites. Two field teams of five to six personnel, checked the units fortnightly to clean the cameras and replace memory cards and batteries. In total, camera polygons covered 18 villages and represented a mosaic of forest-farmland-forest or forest-farmland-settlement. All surveys took place between June 2014 and April 2016 (Table S1).

Species habitat use

In order to estimate species occupancy (habitat use), a single species single season occupancy model was used, based on four main assumptions: (i) occupancy state is closed, where species occupancy and detection probability at all sites (camera trap locations) remained constant over a survey period but may change between surveys; (ii) sites and replicates are spatially and temporally independent, where detecting species of interest at a site is independent of detecting the species at other sites or during other time intervals, (iii) site and survey covariates that influence occupancy are quantified and incorporated in the model calculation, and (iv) factors that influence detection probability are explained through incorporating site covariates and survey covariates within the analyses (MacKenzie et al. 2002), although in this study the (iv) assumption is considered constant across sites (Linkie et al. 2007).

Our large-scale spatial study did not consider certain finer-scale temporal covariates. Camera traps were active for 24 h/day over several months and all data were used, so we did not consider either daily activity budget or date as a covariate. It is possible that our study animals adjusted their daily activity based on weather, which for our study area would be rainy or dry, as the temperature is fairly constant throughout the year being an equatorial rainforest. However, reviewing the scientific literature for similar studies, we decided to follow an occupancy modelling approach that is typically used studies covering large spatial scales (Brodie et al. 2015a; Espinosa et al. 2018; Penjor et al. 2019). Therefore, during model development, we tested eight landscape covariates that were considered likely to influence the spatial behaviour of the marbled cat, golden cat and clouded leopard (McCarthy 2013; Haidir et al. 2018). We included elevation (elev) and slope (slope) using data obtained at 30 m resolution from the Shuttle Radar Topographic Mission (SRTM) (Rabus et al. 2003). We obtained NDVI (Normalized Difference Vegetation Index) (vegcov) data using Global Forest Change data version 1.6 for year 2018. Combination of cloud free Landsat 8 OLI composite images over the year 2018 was used to generate NDVI as a ratio between the red and near infrared values. This dataset was first published in 2013 (version 1.0, see Hansen et al. 2013) and then updated each year (currently version 1.7 which provides data from 2000 to 2019). During the study period, there might have been some changes in vegetation cover, but these would have been minor. This dataset was first published in 2013 (version 1.0, see Hansen et al. 2013) and then updated each year (currently version 1.7 which provides data from 2000 to 2019). Euclidean distance to forest edge (fordist) was calculated based on official forest cover data from BAPLAN (Indonesia Ministry of Environment and Forestry’s Planning and Mapping Centre), data obtained from year 2014. Euclidean distances to villages (vildist), distance to major roads (national and provincial roads; roadist) and distance to rivers (rivdist) were calculated based on spatial layers from BAKOSURTANAL (Indonesia Land Survey and Mapping Agency) for the year 2018. All layers were projected to UTM 47 Mercator Southern Hemisphere Projection and re-sampled to 250 m resolution following Macdonald et al. (2018b) study on clouded leopard in Borneo.

We used photographic evidence from camera trap surveys, sorted into two-week sampling occasions adopting the approach of previous studies by ( Linkie et al. 2007; Haidir et al. 2018), these detection data were then converted into detection matrices for four species: clouded leopard, golden cat, marbled cat and leopard cat. Detection matrices were developed through ‘camtrapR’ package in R (Niedballa et al. 2016). However, due to low detection of the leopard cat (< 5 photographs), that species was excluded from the analysis.

For the three focal species: clouded leopard, golden cat and marbled cat, we applied a single-species, single-season occupancy approach, psi (\(\hat{\psi}\)), which incorporates a function of detection probability (Mackenzie 2006; Kéry et al. 2013). All variables were extracted at the camera trap station location and, before the analysis, these were inspected for collinearity by calculating Pearson’s correlation (Dormann et al. 2013). From the pair of highly correlated variables (|r|>0.7; Booth et al. 1994), we excluded the one with higher AICc value (Akaike’s Information Criterion corrected for small sample sizes; Burnham et al. 2002) in a univariate model.

We assembled and tested a set of 10 candidate and biologically realistic models per species (Table S2). Expecting a parabolic relationship between occupancy and elevation, this covariate was tested in its quadratic term, while for the others, non-correlated covariates, we used normal (non-quadratic) relationships (Penjor et al. 2019). The 10 candidate models were then compared using the AICc and all models with ΔAICc ≤ 4 were averaged using model weights (Burnham et al. 2002) in the ‘wiqid’ package (Meredith 2018). We used ΔAICc ≤ 4 so that we could consider the influence of a greater number of models and their covariates for the final species ’model averaging’(see Richards 2005; and Symonds and Moussalli 2011). To generate predictive maps of habitat use for each species, we applied the averaged model coefficients to the raster representing each of the final model covariates (Rhodes et al. 2009; Banner and Higgs 2017).

To test the performance of the predicted habitat use models based on occupancy following Gould et al. (2019), higher quality data are required, which would substantially increase the cost of the surveys. Due to lack of sufficient data, instead we validated the final models using a Bayesian five-fold cross-validation with 80% of the data used for training and 20% used for validation (Petracca et al. 2018; Penjor et al. 2019). We decided to use Bayesian fivefold, over 10-fold which has high computational demands. To assess model accuracy, we calculated the proportion of expected detection ŷi to observed detection yi. (proportion of true positive and negative observations), sensitivity (proportion of true positive observations correctly identified) and specificity (proportion of true negatives correctly identified).

Landscape connectivity models

Species occupancy modelling generated predicted probabilities for habitat use that reflect the likelihood of each cell being used by a focal species (Gould et al. 2019; Penjor et al. 2019). To model landscape connectivity for each species, we transformed the habitat use layers into ‘resistance to movement’ values. Mateo-Sánchez et al. (2015) and Zeller et al. (2012) found that, in the absence of data on movement behaviour or population genetic structure, which is the case for our three species in Sumatra, spatial information on habitat use can be a useful surrogate for deriving a resistance to movement layer. Mateo-Sánchez et al. (2015) and Keeley et al. (2016) showed that a negative exponential function best describes the relationship between habitat suitability and resistance values. Therefore, we transformed habitat use (H) into resistance (R) ranging from 1 (low landscape resistance to movement) to 100 (high resistance to movement) using the following equation:

$$R = (\exp \left(-1 \times \frac{H}{C}\right)\times 100+Y$$
(1)

Equation 1 transformation of habitat use layer (predictive habitat use) into resistance layer, where R is the resistance (cost) value, H is habitat use (occupancy, ψ), C is the factor used to determine the shape of the curve (0.15), and Y is the value used to convert minimum resistance to 1 (Y for clouded leopard = + 0.67, golden cat = + 0.42 and marbled cat − 0.07).

Landscape features such as settlements, water bodies and roads, likely constituting a barrier to the movement and dispersal of the three felid species, were not included in the habitat use model. Therefore, we ‘burnt in’ these features to the final resistance layers (reclassifying the pixels). For this, we considered large water bodies, settlements and major roads to be impermeable to felid movement and, therefore, assigned them a resistance of 100. Field observations have shown that these species can occasionally cross rivers (Haidir 2016) and, so, major rivers were given a medium resistance value of 50, unless the base resistance value was higher, in which case the higher value was retained.

Landscape connectivity models, which predict patterns of population connectivity, require animal source locations, reflecting the actual distribution and density of the target population (Cushman et al. 2018). The number of source points corresponding to the number of individuals in our study area was informed by the most recent published literature on density estimates of the three focal species in Southeast Asian countries. In summary, we generated a total of 206 source locations for clouded leopard (Haidir et al. 2020), 260 for golden cat and 346 for the marbled cat (Hearn et al. 2016; Rustam et al. 2016; Naing et al. 2017). We then distributed these sets of source points in the landscape, retaining 2 km minimum distance between the points, probabilistically to the habitat use of each species. This was done in three steps: (i) we first created a uniform random raster with values between 0 and 1 with identical extent and cell size to the area of interest; (ii) the random raster was then subtracted by the predicted occupancy layer, so that the positive values of the resulting layers reflect areas with higher than the random probability of species occurrence; and lastly; and, (iii) from the positive values of the resulting layer we randomly selected a set of cells with minimum gap between the pixels of 3 km for clouded leopard and 2 km for both golden cat and marbled cat. The number of the selected cells representing individuals source locations followed the estimated number of individuals (see Fig. S1 for illustration).

To map and quantify core movement areas and dispersal corridors for each of the focal species we applied cumulative resistant kernel (Compton et al. 2007) and factorial least-cost path approaches to the source points and resistance surface layers described above (Cushman et al. 2009). For both methods, we used the UNICOR program (Landguth et al. 2012). The cumulative resistant kernel identifies the main pattern of synoptic connectivity and core habitat areas (Kaszta et al. 2018) by predicting the total movement density across the landscape. This is calculated by summing all individual least-cost kernels from all dispersal source points (Compton et al. 2007). The factorial least-cost path approach complements the cumulative resistant kernels by defining narrow linkages in the landscape where the movement pattern is constrained. Therefore, the least-cost paths set up with larger dispersal thresholds can be applied to define dispersal corridors. The final network of paths is computed by summing the least-cost paths between all possible pairs of points (Cushman et al. 2013b). Therefore, combinations of the factorial least-cost path and cumulative resistant kernel has the advantage over other methods, such as circuit theory, by explicitly including dispersal thresholds that enable the method to realistically reflect the limited distances that can be traversed by real organisms and to explore scale dependent relationships with varying dispersal ability (Cushman and Landguth 2012; Cushman et al. 2013b; Kaszta et al. 2020).

Empirical dispersal abilities of the three felids are unknown. Therefore, to estimate the dispersal thresholds for connectivity analyses we applied Bowman equation where dispersal is calculated as 40 × (home range size0.5) in highly suitable habitat (Bowman et al. 2002). With the average home range of clouded leopard varying from 23 to 45 km2, the golden cat from 33 to 48 km2 and marbled cat from 2 to 29 km2 (Grassman 2004), such estimated dispersal thresholds vary between 57 and 277 km. Based on these estimates, and for the purpose of this study, we chose a lower threshold for cumulative resistant kernels of 125,000 cost units (CU) and an upper threshold of 250,000, reflecting distances of 125 km and 250 km respectively, in a uniform landscape of resistance equal to 1 (Hearn et al. 2019; Kaszta et al. 2020). To model long-distance movement and dispersal, beyond the extent of locally connected populations, we defined the threshold to 1,250,000 CU, which is five times higher than the upper threshold defined for the cumulative resistance kernels (Kaszta et al. 2019; Kaszta et al. 2020). This allows prioritization of linkage areas beyond the dispersal ability of most individuals, which are important for long-term connectivity of meta-populations (Elliot et al. 2014; Cushman et al. 2018; Kaszta et al. in press). To better reflect local gradients of movement density, the final three least-cost paths density layers (further referred to as LCP) were smoothed by calculating the focal mean of a 5 km radius.

To compare the landscape connectivity of the three felids we calculated Pearson’s correlation and the averaged absolute difference between the cumulative resistance kernels layers (across both dispersal thresholds for sensitivity analysis) and between the three least-cost paths layers.

To identify the areas of high conservation priority defined by zones of high core areas overlap between the three felids, we followed three different approaches. First, we used a threshold to define core areas of high-density movement for all three felid species (Cushman et al. 2018; Kaszta et al. 2019; Kaszta et al. 2020). For that, we selected the value of the 55th percentile of the cumulative resistant kernel raster with the lowest maximum value across all three species (i.e. the golden cat at a dispersal threshold of 250,000 CU and threshold value of 1.64). We reclassified the cumulative resistant kernel rasters into binary layers (‘0’/‘1’), with areas equal or above the threshold were reclassified as 1, and below that value as 0. We then summed up all the binary layers (‘0’/’1’) representing core areas of each felid to delineate the inter-species core areas overlap. Second, to avoid the arbitrary selection of a threshold for defining core areas, we normalized the values of cumulative resistant kernel surfaces across all species by rescaling the original cumulative resistance layers to continuous values ranging between 0 and 1. We then summed up the three rescaled layers. The third approach of calculating core areas overlap across species was through to multiplying the three rescaled cumulative resistance kernels layers from the previous steps.

We also identified the most important corridors predicted to be jointly used by all three species. To maintain only significant connections, we used the 20th percentile value of the species with the lowest maximum LCP (i.e. 0.24 for marbled cat) (Cushman et al. 2018; Kaszta et al. 2020) and converted the LCP layers into binary layers with 1 assigned to the values equal or above the threshold. These binary LCPs were then summed up and overlapped with a layer of the protected area network to calculate the proportion of the most important corridors that is legally protected.

To assess how much of the core areas and corridors are contained within the protected area network in the study area, we calculated for each species the proportion of the resistant kernels for each species being officially protected. Lastly, we further calculated for all three species the proportion of the joined connectivity for all three species (sum and multiplication of the rescaled resistant kernels layers) that lie within the protected area network.

Deforestation risk

To estimate rates of deforestation and model forest loss, we performed a time-series forest cover analysis using ArcGIS 10.1. The forest cover maps were obtained from BAPLAN (Indonesia Ministry of Environment and Forestry’s Planning and Mapping Centre) between the years 2000 and 2017. We calculated the extent of forest cover in 2000 as the landscape baseline and mapped the amount of forest cleared for the years 2003, 2006, 2009, 2011, 2014 and 2017, to determine the rates and locations of loss.

To model forest loss risk, we created two layers that depict the historical record of deforestation in the study area. The first layer was the deforestation data (binary, deforested = 1; forest = 0) for 2003–2014 that was used to train the deforestation model. The second layer contained information on deforestation (also binary) occurring from 2014 to 2017 and was used to validate the predicted deforestation model. To sample the landscape, we created 10,000 random points with a minimum distance of 3 km between each point to avoid spatial autocorrelation. We used Moran's I test in 'spdep' (Bivand and Wong 2018) package in R to test for spatial autocorrelation in the best model’s residuals’. To maintain the proportion between forested and deforested regions in the study area (Kaszta et al. 2017), a total of 89 points were selected from deforested areas and another 717 points from forested areas (from 2003 to 2014). For consistency, we selected 30 points from deforested areas in 2017 and 688 points from forested areas, in 2017, to validate the predictive model. All calculations were performed in R (R Core Team 2017).

To investigate the landscape factors driving deforestation, and to predict future forest loss and fragmentation patterns, we developed 10 candidate generalized linear models (GLMs) with a binary response variable (1: deforested, 0: forest). We used the same GIS covariate layers as for the small felid species’ habitat use analysis, testing for collinearity prior to model inclusion. We, therefore, used the following non-correlated variables: elevation, slope, distance from roads, and distance from villages (see Table 2).

Model parsimony was assessed using AIC corrected for small sample size (AICc). Beta coefficient values of the final model with the lowest AICc were used to predict future deforestation risk. Based on the deforestation model coefficients we generated a predictive map of probability of deforestation for the period between 2015 and 2026 by using a percentage tree cover layer updated by the deforestation from 2003 to 2014.

To assess the predictive strength of the models, we calculated the area under the ROC curve (AUC), sensitivity and specificity by comparing the observed values (deforested areas in 2017) with predicted deforestation values (generated from the best deforestation model using data from 2003 to 2014).

To assess which core habitats of the three focal felids might be most affected by future deforestation risk, we summed the probability of deforestation within core habitats of each species and in the areas of high core areas and corridors overlap amongst the three cats. We then compared these values to the deforestation probability predicted for the entire study area. The same procedure was applied to calculate the percentage of deforestation predicted to occur in the protected area network.

The whole processes and workflow that indicate steps taken in this study starting from field data collection to final results are visually shown in Fig. 2.

Fig. 2
figure 2

Workflow of the whole process in defining core areas and corridors from occupancy approach

Full size image

Results

Focal species habitat use

From the 55,856 combined camera trap nights recorded, there were 211 clouded leopard detections, 137 golden cat detections and 50 marbled cat detections. The model-averaged (\(\hat{\psi}\)) estimate with the highest density interval (HDI) was 0.18 (0.12–0.27) for clouded leopard, 0.29 (0.19–0.42) for golden cat and 0.21 (0.06–0.36) for marbled cat with top candidate models of respective species presented in Table 1. The southwest quadrant (Ipuh) had the highest predicted habitat use for all species at \(\hat{\psi }\) > 0.50, followed by BG, SP, RKE and MH at \(\hat{\psi }\) ~ 0.40, and BG and KM (\(\hat{\psi }\) ~ 0.30), with Batanghari (north-eastern) at \(\hat{\psi }\) ~ 0.57.

Table 1 Standardized coefficient (β) for predicted habitat use of clouded leopard, golden cat, and marbled cat based on the top ranked and averaged occupancy models
Full size table

Clouded leopard habitat-use was found to be explained by two top models within ∆AICc ≤ 4 (Table S2) that included elevation, slope, distance to villages, vegetation cover, tree cover and distance to forest edge (Table 1). Our results showed that the species responded negatively to higher elevation and positively to the increased tree cover and vegetation cover. Higher distances to both forest edge and villages were also found to have a strong positive impact on clouded leopard habitat use.

Golden cat habitat use was explained by seven plausible models (∆AICc ≤ 4, Table S2) and the averaged model showed a significant negative relationship with elevation (Table 1). This species was found to prefer areas closer to the forest edge and with gentler slopes.

Marbled cat habitat use was explained by six plausible models (Table S2) and the averaged model (∆AICc ≤ 4, Table 1) revealed a negative species response to higher elevation and a positive response to slope, distance to villages and forest edge, vegetation cover and tree cover.

Model validation indicated that model accuracy for each species was above 0.75 (95% CRI 0.73–0.91) and the AUC above 0.66. The specificity values varied from 0.86–0.94 and sensitivity from 0.10 to 0.28 (Table S3).

Populations connectivity: core areas and dispersal corridors

The overlap for all three species is shown in Fig. 3, while the cumulative resistant kernels layers for the tested dispersal threshold of 125 kCU and 250 kCU, as well as LCP (corridors) layers for clouded leopard, golden cat and marbled cat, are presented in Fig. S1.

Fig. 3
figure 3

Predicted species occupancy (top), resistance layer (middle), and core areas and corridors (bottom). First row indicates species-wise occupancy values from low (blue) to high (red) for clouded leopard (a), golden cat (b) and marbled cat (c). Second row is landscape resistance, same respective species, with darker colour indicating lower resistance. The third row is predicted core areas and corridors from low density (blue) to high (red), both core areas and corridors are based on a lower dispersal threshold of 125 k cost units for the three species

Full size image

The total value of connectivity (sum of kernel density pixels) for both dispersal thresholds was higher for golden cat and double that of clouded leopard, with the marbled cat having the lowest total predicted density of movement (Table 2). The mean value of kernel density was also highest for the golden cat (2.8 and 1.8 for dispersal thresholds 125 kCU and 250 kCU, respectively) and lowest for the marbled cat (1.4 and 0.7; Table S4).

Table 2 Standardized coefficients (β) for the final deforestation model with the lowest AICc
Full size table

The sensitivity analysis between the two tested dispersal thresholds showed a high correlation between cumulative resistant kernels (r \(\ge\) 0.94) and the averaged absolute difference was low (AAD < 0.75), indicating that the resistant kernels correspond to each other well in the two thresholds (Table S5).

The individual population connectivity patterns of the three species showed that generally clouded leopard and golden cat had substantially larger core areas in comparison to marbled cat, despite the dispersal threshold (Fig. S1). A binary core areas map, based on 55th percentile threshold applied over the resistant kernels at lower dispersal threshold of 125 k CU, showed that clouded leopard had one predominant core area (4000 km2) and several smaller core habitats (< 1000 km2; Fig. 3). Golden cat core areas were substantially larger with two main cores (8500 km2 and 3400 km2) and several smaller core habitats (< 1500 km2). Marbled cat, with the smallest and most fragmented core areas amongst the three cats, was predicted to have one larger core area (1500 km2), with several much smaller cores (50–300 km2) all distributed far (> 20 km) from each other. LCPs of the three felids revealed that all these core areas are potentially still connected by long-distance dispersal corridors, but some of those linkages are weak, especially for the marbled cat (Fig. 3 and Fig. S1). Clouded leopard corridors were the strongest of those predicted amongst the three species.

Patterns of connectivity and species overlap

The AAD between kernel density surfaces (Table S5) was generally small across all species and both dispersal thresholds (AAD < 2). However, the largest difference between kernel densities reflect by high AAD and low correlation was found between kernel densities of marbled cat and golden cat (AAD = 1.92 and 1.33, r = 0.29 and 0.37 for dispersal thresholds 125 kCU and 250 kCU, respectively; Table S5). The lowest difference in connectivity patterns, as shown by the lowest AAD and highest correlation between kernel densities, was found between the clouded leopard and marbled cat (AAD = 0.69 and 0.45, r = 0.72 and 0.76 for dispersal thresholds 125 kCU and 250 kCU, respectively). Dispersal corridors (LCPs) showed much lower correlation (< 0.03) and higher AAD (6.6–17.7) than kernel density layers (Table S6). However, similarly to kernel density surfaces the highest correlation and smallest AAD was reported between the clouded leopard and marbled cat (AAD = 6.64, r = 0.03) and the highest difference was found between golden cat and marbled cat (AAD = 17.75, r = 0.03).

The analysis of spatial overlap between the core areas defined by the threshold of 55th percentile revealed that 7–11% (approximately 1250–1650 km2, depending on the tested dispersal thresholds) of the total core areas were potentially utilized by all three species (Fig. 3). These habitats were mostly (88%) located inside the protected area network. Core habitat overlap of two species only was much larger and represents 18–22% of all core areas. However, only 30–44% (depending on the dispersal threshold) of this habitat was inside the protected area network. The majority of the core areas (71% for both dispersal thresholds) were inhabited by only one species and approximately 21% of it was located in the protected area network (Fig. 3). The core areas defined by the binary overlap of all three species covered 4924 km2, with the majority (54%) of all three species core areas encompassed inside Kerinci Seblat National Park. The remaining 46% of the three species overlapped located outside KSNP was found to lie within Batanghari Protection Forest (north-eastern of the landscape; 23%), along the border in the western KSNP 15%, in south-western KSNP (further from KSNP forest, > 10 km) 6%, and the remaining (2%) were scattered in forested areas in the north western of the landscape.

The sum and multiplication of rescaled resistant kernel layers of the three felids indicated similar areas of the highest importance for joint conservation of the three felids, when comparing to the results of binary core areas overlap (Fig. 3). The proportion of the total value of the summed-up kernel density layers within the protected area network was 41–42%, depending on the dispersal threshold. The proportion of the multiplied kernel density layers, indicating the potential presence of all three species, being located within the protected area network, was 87–94%.

Overlap in dispersal corridors (LCPs) revealed that only 10% of predicted corridors were used in common by all three felid species, and 68% of all corridors, whether shared or not, lay within the protected area network (Fig. 3) and 46% of these two-species corridors were within the protected area network.. The single species corridors represented the vast majority of the corridors network (63%), with 32% of these single-species corridors lay within the protected area network (Fig. 3).

Deforestation risk

The most parsimonious deforestation model showed high accuracy in its predictive power (AUC = 0.80; Table S7) and model validation indicated that the model had a sensitivity of 0.11 and specificity of 0.97 (using threshold of 0.5). The final deforestation model was not affected by spatial autocorrelation (Moran’s I = 0.17, p > 0.06).

The top model, which ranked much higher than the second top model (ΔAICc = 8.1; Table S7), included four of the five tested variables: elevation, slope, distance from roads and distance from villages (Table 1). The model predicted a lower probability of deforestation in areas that were higher elevation, of steeper slope and greater distance to roads and villages (Fig. 4).

The deforestation risk model predicted that 1052 km2 of the landscape was at risk of clearance with 32% of the total deforestation probability predicted to occur within protected areas (Fig. 5).

Fig. 4
figure 4

Core areas and dispersal corridors overlap between the three felids for dispersal threshold of 125 k cost units (a), and prediction map of future probability of deforestation (b)

Full size image
Fig. 5
figure 5

Multi-species connectivity. Core areas overlap for the three species for dispersal threshold 125 k CU (panels a, b and c) and 250 k (panels d, e and f). Overlap of the binary core areas defined by a resistant kernels threshold value (panels a and d), sum of the 0–1 rescaled resistant kernels for the three species (panels b and e), multiplication of the 0–1 rescaled resistant kernels (panels c and f)

Full size image

The layer of future deforestation risk predicted that approximately 4% of the total deforestation will occur within both the clouded leopard and golden cat core habitats and, 1.6–2.8% within marbled cat core areas (calculations were based on the binary core areas and for two dispersal thresholds; Table S8). Furthermore, 2.7–3.7% of the total predicted deforestation will affect the joint core areas and 4.5% was predicted to occur within corridors jointly utilized by the three species. Additionally, 5–6% of the deforestation was predicted to affect core areas and 5% to occur in corridors jointly utilized by two felids (Fig. 5; Table S8).

Discussion

This study presents a quantitative framework to assess landscape connectivity for populations of small- and medium-size felid species, and the effects of future habitat loss, in Sumatra. However, the methodology proposed here can also be applied to other areas and other species, and can ultimately assist efficient management of practical, on-the-ground and multi-species conservation efforts. Based on an extensive camera trap surveys we identified habitat use and landscape connectivity for populations of clouded leopard, golden cat and marbled cat. We identified the most important core habitats and dispersal corridors for each of them, as well as key habitats for joint conservation of these felids. Using six intervals (data from year 2000 as baseline: 2003, 2006, 2009, 2012, 2015, and 2017) to sample the deforestation throughout the 17 years’ overall (2000–2017) period of deforestation data, we have also modelled the probability of future deforestation risk as a major threat to these species in the landscape. Hence, we detected key habitats, the loss of which might, in the longer-term, affect the stability not only of populations of the three felids but also of other species ecologically linked to them.

Focal species habitat use

All three felids showed a significant and non-linear relationship with increasing elevation, avoiding lowland areas but also disfavouring high altitude areas. Similar findings were reported by McCarthy et al. (2015) on golden cat, and Sunarto et al. (2015) for all three species in Sumatra, and Hearn et al. (2018) for clouded leopard and marbled cat in Borneo.

The results of our occupancy models support several assessments of Southeast Asian small- and medium-sized wild cats that indicate the sensitivity of these cats to habitat disturbances (Rayan and Mohamad 2009; Sunarto et al. 2015; Brodie et al. 2015b; Granados et al. 2016). In line with findings of previous studies (Brodie et al. 2015a; Tan et al. 2017; Macdonald et al. 2018a; Penjor et al. 2018; Hearn et al. 2019; Haidir et al. 2020), our analysis revealed that predicted habitat use for all three felids is confined to densely forested areas and further from human disturbances. Similarly, golden cat and marbled cat, although found to be less sensitive to forest cover changes, and not specifically confined to old-growth forest, were highly associated with areas that have increased tree/canopy cover, including agroforestry land, forest plantations and selectively logged forest (Wearn et al. 2013; McCarthy et al. 2015).

Population connectivity

The connectivity analysis revealed that the golden cat population in our landscape was the most connected of the three felid species, with large and strong core areas representing meta-population strongholds, linked by a diffused network of dispersal corridors (Fig. 4). The largest golden cat core habitat stretches across the western KSNP, with more than half of it being located outside of the park boundary. The second core area covers Batanghari Protection Forest. These two main core areas are connected through a network of corridors that are partially located in the northern KSNP.

The clouded leopard population, although predicted to be much less connected than that of the golden cat, was also predicted to have two well-defined and relatively robust core population strongholds (Fig. 3). The largest one was located almost entirely inside KSNP and the other in the Batanghari Protection Forest. They are linked by a network of dispersal corridors in northern KSNP, with the strongest (highest LCP values) corridor overlapping one of the golden cat’s dispersal corridors, highlighting its high importance for conservation.

Clouded leopard landscape connectivity patterns were spatially and quantitatively most similar to those of the marbled cat. The main marbled cat core area, located in central KSNP, although much smaller than that of the clouded leopard, overlapped with the largest core area of the latter (Fig. 3). Our analysis predicted the marbled cat population to be the least connected amongst the three species, with small patchy core habitats, which makes this species the most susceptible to landscape disturbance.

Conservation effectiveness and future habitat loss

Our results showed that large parts of predicted core areas are inside state forests but some are not legally protected by any of the protection categories defined by the IUCN, from 12% of the core habitats jointly utilized by the three felids, up to 70% used by only two species, and 79% used exclusively by one of the species. This is similar for dispersal corridors, where 47% of corridors are jointly suitable for the three species, 54% of corridors jointly suitable for two species, and 68% of corridors used by only one of the species, are located outside of the protected area network. The analysis of overlap for core areas and corridors identified areas of high conservation importance where core habitats of all three species overlap—one of them is central KSNP and the second is Batanghari Protection Forest.

Batanghari Protection Forest and its surrounding area also have a resident population of Sumatran tigers (Dinata 2008). It is partly categorized under Indonesian law as ‘limited production forest’ and watershed protection forest. However, the protection status of this area is not recognized by the IUCN classification scheme and, in practice, this status and limited conservation attention have resulted in forest clearance and the development of a local road network (Sloan et al. 2019), which impacts the predicted core habitats and corridors. The area between KSNP and Batanghari was identified in our analysis as a network of key corridors used by all three felid species and connecting the most important core habitats. We recommend, for future infrastructure development in this area, the implementation of strategic environmental assessments be conducted that take into account the importance of this key area for population connectivity of the three species.

Our deforestation risk model indicated that areas of lower elevation and gentler slope, located close to roads and villages were most susceptible to deforestation. This trend is similar to findings from deforestation modelling in Borneo (Cushman et al. 2017). Most importantly, our model showed that, over the 10 years covered by our predictions, 335 km2 (equivalent to ~ 2% of the overall landscape) of high deforestation risk in the 16,000 km2 Kerinci Seblat landscape is predicted to take place within protected areas, although the areas of the highest probability of deforestation are mostly located outside the protected areas network. Furthermore, up to 10% of the predicted deforestation may occur within core habitats and up to 5% of the deforestation is likely to be in the dispersal corridors utilized by the three or two felid species, which are the areas of the highest conservation priority.

Scope and limitations

The validation of the occupancy models showed high accuracy and high AUC and specificity values. However, the sensitivity of these models was low, indicating the good ability of our models to predict non-occurrence of the three species but relatively low capacity to predict their presence. This is most probably due to the low detectability rate of these three rare felids during our surveys, with marbled cat probability of detection being the lowest. The low detectability might constitute a potential limitation to the habitat use predictions made by our models, and further surveys with increased trapping effort and improved sampling design should be carried out in the Kerinci Seblat landscape.

Management implications

The current status of the Kerinci Seblat landscape as one of 73 National Strategic Areas of importance to regional sustainability should promote the integration of a landscape-based development approach for wildlife conservation (MoPWH 2017). The role of remnant forests adjacent to KSNP in facilitating population connectivity shows the need for such an approach in maintaining viable populations of small- and medium-sized felid species. Our multi-species models offer suggestions on further designation of core areas and both physical and functional corridors within the national park and other protected area management units. The results and maps produced in this study provide information for national park managers in identifying priority areas for conservation, i.e. forest patrols and other support to areas that are likely used for felid movement. The KSNP authority, as the main government institution, is the epicentre of future collaborations with local government agencies and stakeholders in managing the landscape.

Although, a relatively large portion of the predicted small-and medium-sized cat dispersal corridors and some of the most important core habitats, are located in the protected area network, this conservation status does not guarantee prevention of deforestation risks. The majority of the remaining forest is currently within the protected area network (67%) and an additional ~ 30% is within state-owned forest areas, therefore relatively little additional forest remains outside of the protected forest network. Consequently, our deforestation risk model predicts that one third of the predicted deforestation may occur within the boundaries of the protected area network. Poor et al. (2019) found that the role of protected areas in effectively curbing deforestation should be strengthened, although protected areas in Sumatra already have the advantage of practical protection afforded by its location at rather high elevation typified by rugged slopes that discourage deforestation through the expense of access.

With the current pace of infrastructure development, managers of conservation areas face two substantial challenges in (1) mitigating the impacts of accelerated deforestation rates in central Sumatra; and, (2) ensuring national policies and conservation programs are well-communicated and integrated into provincial and district government plans (Poor et al. 2019; Sloan et al. 2019). In the meantime, national government restrictions on new road developments and/or expansion of existing roads that transect national park forests, e.g. those listed in the UNESCO Tropical Rainforest Heritage of Sumatra sites (GoRI 2019), should be well-communicated to the local authorities and stakeholders.

We propose that both the results and methodological framework of this study serve as a guideline for future development planning to achieve a smart and scientifically-guided compromise between development and wildlife conservation goals (Rayan and Linkie 2015). The methodology applied here, and our results, may be transferable to identify connectivity corridors for other sensitive species of concern and/or other areas where major infrastructure developments are taking place (MoPWH 2017), e.g. where the trans-Sumatra highway and major road upgrade projects are being planned or have already begun (CMoEA 2011).