The RSPL, situated in Bohol Province within the Central Visayas Region of the Philippines, spans approximately 10,453 hectares and represents the largest remaining contiguous natural habitat on the island (Fig. 1). The site is characterised by a heterogeneous landscape of forested limestone hills, open grasslands and karst formations, with land cover dominated by open forests and interspersed with brush, cultivated areas and patches of built-up land (Birding2Asia 2014; Aureo et al. 2020). Elevation within the RSPL ranges from 0 to 870 metres above sea level, with higher altitudes concentrated in the central and eastern sections, corresponding to steeper and more rugged terrain (Aureo and Decena 2023). These upland areas form part of Bohol’s critical recharge zones, supplying water to numerous springs and rivers. The region experiences a Type IV climate, based on the Modified Coronas Classification (MCC) system, with rainfall relatively evenly distributed year-round and an annual precipitation average of around 3,756 mm (Lantican 2001; BirdLife International 2014). The predominant soil type is the Batuan–Faraon complex, which, while suitable for agriculture in flatter areas, is subject to fertility constraints and erosion risks on steeper slopes. Approximately 62% of the terrain consists of slopes between 0–18%, allowing for various land uses, whereas steeper zones require targeted management to maintain ecological stability and water quality. The RSPL thus plays a central role in sustaining both biodiversity and ecosystem services within Bohol’s complex topography.

The study implemented an integrative spatial modelling framework to identify suitable habitats and HCVAs for threatened species in RSPL (Fig. 2). The methodology followed a structured pipeline consisting of four stages: data input, analytical processing, output generation and area identification. Environmental preparation involved standardising bioclimatic and topographic data layers, with climatic variables interpolated and terrain variables, such as the Terrain Ruggidness Index, derived from a DEM. Species occurrence records underwent spatial thinning to mitigate sampling bias before habitat suitability was modelled using MaxEnt. Model performance and data integrity were rigorously assessed using the VIF, Jackknife test, AUC and TSS and sensitivity analysis (MTP vs. maxSSS). Through composite richness mapping and threshold sensitivity analysis, individual habitat suitability maps were converted into binary layers and aggregated to produce a species richness map. Finally, HCVAs were delineated by applying statistical thresholds to the richness distribution, providing a repeatable, quantitative basis for supporting site-based conservation decisions.

Environmental and species occurrence data were sourced from authoritative local, national and international repositories to support the spatial analyses conducted in this study (Suppl. material 1). Climatic variables, including temperature and precipitation data (Suppl. material 2), were obtained from the Philippine Atmospheric, Geophysical and Astronomical Services Administration (https://www.pagasa.dost.gov.ph), which provides long‑term climate data and derived bioclimatic layers for the Philippines (Department of Science and Technology - Philippine Atmospheric, Geophysical and Astronomical Services Administration 2023). Species occurrence records were compiled from the Global Biodiversity Information Facility (GBIF 2023) available at https://www.gbif.org/, supplemented by local data from the Biodiversity Assessment and Monitoring System (BAMS) reports of the Department of Environment and Natural Resources (Department of Environment and Natural Resources 2019), as well as field survey data collected by the Central Visayas Agriculture, Aquatic and Natural Resources Research and Development Consortium (CVAARRDEC). These species observations were vetted for coordinate accuracy and taxonomic consistency prior to inclusion in the analyses. In addition, topographic variables, such as elevation, slope, aspect and terrain ruggedness (Suppl. material 3), were derived from the ALOS Digital Elevation Model (https://www.eorc.jaxa.jp/ALOS/), which provides a 12.5‑metre resolution terrain dataset ideal for fine‑scale ecological and spatial analyses (JAXA 2019). Moreover, the RSPL boundary shapefile was obtained from World Database on Protected Areas available at https://www.protectedplanet.net/en and was verified using spatial data from the DENR. Together, these multi‑source datasets (Suppl. material 4) enabled a robust assessment of species richness, habitat distribution and topographic variability across the protected area.

In this study, spatial interpolation of 19 bioclimatic variables was performed using IDW to generate continuous proxy environmental gradients aligned with the spatial resolution and extent of the reference DEM. The input data consisted of point‑based climatic observations (e.g. Bio_1–Bio_19) stored in CSV format, each containing geographic coordinates and associated climate values. The IDW method was implemented using the gstat package in R (Pebesma 2004), with a power parameter of p = 2, based on the following formula:

Z(x0) = [Σ(Z(xi)/d(x0, xi)^p)] / [Σ(1/d(x0, xi)^p)],

where Z(x0) is the estimated value at an unsampled location x0, Z(xi) is the observed value at location xi, d(x0, xi) is the Euclidean distance between x0 and xi and p is the inverse distance weighting exponent.

A value of p = 2 was selected as it is a common choice that ensures closer points have a stronger influence on the interpolated value. Since only three sampling points were available, the interpolation was constrained to the three nearest neighbours (nmax = 3), minimising extrapolation artefacts and overfitting. The resulting raster was then projected and resampled to match the spatial resolution, extent and coordinate reference system (CRS) of the DEM. IDW was selected over geostatistical methods such as kriging due to the extremely limited number of input points, which are insufficient to reliably fit a semi-variogram model (Li and Heap 2014). IDW is a deterministic method that performs robustly under data‑sparse conditions and does not assume statistical properties of spatial autocorrelation, making it a practical and reproducible approach for interpolating environmental surfaces in data‑limited contexts.

To derive key topographic variables for SDM, a digital elevation model was processed using terrain analysis methods implemented in the terra package in R. The variables generated included elevation, slope, aspect and the topographic ruggedness index (TRI), all of which can shape microclimates and species distribution patterns (Wilson et al. 2007). These topographic predictors capture terrain complexity and environmental gradients, which are critical for explaining species distributions, especially in mountainous or heterogeneous landscapes (Körner 2007; Elith and Leathwick 2009). Elevation was extracted directly from the DEM and represents the height above sea level in metres. Slope was calculated as the maximum rate of elevation change between a cell and its neighbours and is expressed in degrees. The slope (S) is calculated as:

S = atan(sqrt((dz/dx)² + (dz/dy)²)),

where dz/dx and dz/dy are the partial derivatives of elevation in the x and y directions, respectively.

Aspect, representing the compass direction that a slope faces, is calculated as:

A = atan²(dz/dy, -dz/dx).

The Topographic Ruggedness Index (TRI) is calculated as:

TRI = (1/n) Σ|z0 - zi|,

where z0 is the elevation of the central cell, zi are the elevations of surrounding cells and n is the number of surrounding cells.

A MaxEnt modelling approach (version 3.4.4) was employed to predict the potential habitat suitability of the target species using presence-only occurrence data and a suite of environmental predictors (Phillips et al. 2006a; Phillips et al. 2006b; Phillips et al. 2009). Occurrence records were initially cleaned by removing duplicate entries and erroneous coordinates. To mitigate the effects of spatial autocorrelation and sampling bias, we applied spatial thinning and "bias file", resulting in a final dataset of 83 spatially independent occurrence records (Pearson et al. 2006).

To ensure reproducibility and address the complexities of species-environment relationships, the model was parameterised using Linear, Quadratic, Product, Threshold and Hinge (LQPTH) features (Suppl. material 5). We utilised the logistic output format and enabled clamping and extrapolation to account for environmental values outside the training range. Model performance and variable importance were evaluated using response curves and jackknife tests, respectively. To improve calibration, presence samples were added to the background and all output grids were generated in .asc format.

To minimise model over-fitting and improve the interpretability of environmental responses, we implemented a variable reduction process to refine the initial set of 23 bioclimatic and landscape predictors. We subjected the initial variables to a VIF analysis, iteratively excluding those with a VIF > 5 to eliminate extreme redundancy (Dormann et al. 2012). This screening ensured that the final models were calibrated using a robust subset of five independent predictors — Annual Mean Temperature, Annual Precipitation, Elevation, Slope and Terrain Ruggedness Index (TRI) — all of which maintained VIF values < 4.0, well below the conservative threshold of 5.0 (Naimi et al. 2013). The relative influence of these predictors was further evaluated using Percent Contribution and Jackknife tests of variable importance, where a variable was considered a "primary driver" if it yielded the highest regularised training gain in isolation or if its omission resulted in the most substantial decrease in gain.

To represent the available environment, 10,000 background points were randomly generated across the RSPL. Model performance was evaluated using a dual-metric approach: the Area Under the Receiver Operating Characteristic Curve (AUC) to assess general discrimination and the True Skill Statistic (TSS) to evaluate accuracy independently of prevalence. To ensure reliability across varying sample sizes, we utilised an adaptive validation framework: 5-fold cross-validation for species with n < 10 and a Jackknife (Leave-One-Out) procedure for data-limited taxa (n < 10).

Final models were trained using the fully thinned dataset to produce continuous habitat suitability maps. To transform continuous probability surfaces into binary presence-absence maps, we applied the maxSSS threshold. While MTP was initially considered to minimise omission, a sensitivity analysis revealed that MTP inflated predicted presence by up to 74% in certain taxa by encoding spatial noise into the binary maps (Suppl. material 6). The maxSSS threshold was selected as it optimises the trade-off between sensitivity and specificity and is widely recognised as one of the most robust thresholds for conservation planning (Liu et al. 2013). All modelling procedures were conducted using the dismo, terra, raster and usdm packages (Naimi et al. 2013; Hijmans et al. 2017; Hijmans 2020; R Core Team 2023).

The species richness data were analysed and visualised to assess the distribution and determine appropriate conservation thresholds (Table 1). We aggregated taxa to identify 'multi-taxa hotspots' that represent high-integrity ecosystems beneficial to the broader community, rather than managing for a single flagship species. The species richness values for each grid cell were compiled into a numeric vector and used to generate both the cumulative distribution and the frequency distribution graphs. The cumulative distribution was created using the ecdf() function, which calculates the empirical cumulative distribution of species richness across the landscape. Vertical lines were added to mark the statistical thresholds (Mean ± 1 SD, CRFD inflection point, 80^th, 90^th and 95^th percentiles) using the abline() function in the ggplot2 package, with distinct line styles and colours for clarity. The frequency distribution was generated as a histogram using hist() and complemented with a density plot overlay (lines(density(...))) to highlight the distribution’s shape and modality. The thresholds were similarly annotated with vertical lines.

To identify HCVAs, five statistical thresholds were assessed: the Mean ± 1 SD, the inflection point of the Cumulative Richness Frequency Distribution (CRFD), the 80^th, 90^th and the 95^th percentiles. These thresholds were selected to capture a range of richness intensities across the landscape, from common areas (Mean ± 1 SD) to rare hotspots (90^th and 95^th percentiles). The CRFD inflection point was computed to locate the richness value at which the trend shifts from common to rare (Kukkala and Moilanen 2012). In selecting an appropriate threshold, both statistical rigour and ecological relevance were considered, especially within the context of an existing protected area (Hansen et al. 2020). Each HCVA map was treated as a binary surface (1 = HCVA, 0 = Non‑HCVA) and species occurrences were extracted to assess whether they fell within HCVA areas. The sensitivity of each HCVA map, defined as the proportion of observed species occurrences correctly identified within HCVA areas (True Positive Rate), was calculated using the formula:

Sensitivity = Number of occurrences within HCVA / Total occurrences x 100.

To estimate the uncertainty of these sensitivity values, 95% confidence intervals were calculated using the normal approximation for binomial proportions (Agresti and Coull 1998), an approach well established for assessing spatial prediction performance (Fielding and Bell 2002). This approach provides a robust and interpretable measure of map effectiveness for conservation planning, aligning with best practices in species distribution modelling (Fawcett 2006). Recent studies also emphasise the use of threshold‑based metrics, such as sensitivity, specificity and AUC, in selecting spatial prediction maps for conservation prioritisation (Buckland ST et al. 2022; Bracken FSA et al. 2022; Thomas et al. 2024). By focusing on sensitivity and its associated uncertainty, this method allows for a direct and transparent assessment of how well HCVA thresholds capture actual species occurrences, making it a valuable tool for guiding site‑level conservation and management.

The MaxEnt models for the 13 species yielded high predictive scores, with a mean AUC of 0.905 and a mean TSS of 0.717 (Table 2). Performance remained statistically robust across the ensemble, with AUC values ranging from 0.860 (E. mindanense) to 0.950 (G. crinigera). The high TSS values (mean 0.717) further confirm that the models maintained high accuracy independent of prevalence, effectively distinguishing suitable habitats across the RSPL. However, the robust statistical performance of the models must be interpreted through the lens of the data constraints inherent to high-diversity, data-poor tropical landscapes. A primary limitation is the reliance on bioclimatic variables interpolated from only three local sampling points. These layers should be viewed as proxy gradients reflecting broad landscape transitions rather than high-fidelity climate surfaces.

The multi-collinearity screening successfully reduced the environmental suite to five independent predictors, all maintaining VIF values < 4.0 (Table 2). Across the 13 species, Annual Mean Temperature (29.2%) and Elevation (28.2%) emerged as the most frequent primary drivers of habitat suitability (Fig. 3). However, the percentage contribution values revealed high niche specificity amongst taxonomic groups. For instance, temperature was the dominant constraint for restricted endemics, contributing 62.3% to the model for G. crinigera. Annual Precipitation was the critical driver for the amphibian P. guentheri (60.2%), reflecting the physiological dependency of karst-dwelling amphibians on moisture regimes. Elevation was the primary predictor for dipterocarp trees, notably P. steerii (54.9%) and D. grandiflorus (45.2%). Furthermore, terrain complexity — represented by TRI and Slope — contributed significantly to the models for M. fascicularis (33.1%) and M. schreibersii (35.6%), respectively.

Species distribution models (SDMs) for thirteen threatened species (Fig. 4) revealed distinct spatial patterns of habitat suitability across the RSPL. Each species exhibited a unique spatial signature, reflecting ecological specificity and differential responses to environmental gradients within the landscape. Highly restricted and fragmented suitability patterns were identified for Gallicolumba crinigera, Dipterocarpus grandiflorus and Shorea guiso. These species are primarily concentrated in the western and northern sectors of the RSPL, where suitable habitat is confined to intact forest interiors. Similarly, localised hotspots of high suitability were identified for Hydrosaurus pustulatus, Macaca fascicularis and Todiramphus winchelli, suggesting a heightened vulnerability to site-specific disturbances within their narrow ranges. Meanwhile, intermediate spatial patterns were observed for Pitta steerii, Hopea acuminata and Miniopterus schreibersii, with habitat suitability concentrated primarily in the southern and central corridors of the protected area. In contrast, a broader and more continuous suitability gradient was predicted for Platymantis guentheri, Edolisoma mindanense and Voacanga megacarpa. The expansive distributions of these three species across the central and western portions of the landscape indicate higher ecological plasticity compared to the more restricted endemics. These results highlight a stark spatial heterogeneity across the RSPL, where specialised "suitability islands" for forest-dependent birds and dipterocarps contrast with the broader landscape-level suitability predicted for certain herpetofauna and flora.

The species richness underscores the spatial heterogeneity of habitats within the RSPL, highlighting areas of high ecological significance (Fig. 5). The northern and north-western sections of RSPL, which register the highest richness classes (10–12 species), align with findings from other tropical forest ecosystems that identify intact, mature forest areas as hotspots of species diversity (Myers et al. 2000; Jarnevich et al. 2015). In contrast, southern and eastern portions of RSPL, classified within lower richness classes (0–3 species), likely reflect degraded or fragmented habitats that offer fewer niches and lower resource availability for threatened species.

The distribution of species richness across the study area is characterised by a multimodal frequency pattern, with a notable concentration of lower richness values (Fig. 6, left). The cumulative distribution indicates that the landscape is predominantly composed of areas supporting 2–5 species, after which the cumulative proportion increases more gradually. A critical inflection point is observed at the CFRD threshold of six species, representing approximately the 70^th percentile of the landscape. Beyond this point, the availability of high-richness areas decreases significantly, confirming that locations with high species counts are spatially restricted and relatively rare. The frequency distribution (Fig. 6, right) highlights a primary peak in pixel density at a richness count of 3, followed by a secondary peak at 5, suggesting a landscape dominated by moderate biodiversity. However, the distribution exhibits a sharp decline following the 90^th percentile (~ 9 species) and 95^th percentile (10 species) thresholds. These statistical markers delineate the most ecologically significant hotspots, which represent less than 10% of the total study area.

While the 90^th and 95^th percentiles effectively pinpoint rare hotspots, these thresholds capture only the very top fraction of species-rich areas, making them suitable for highly-targeted conservation priorities (Fig. 7). However, in a context where the area already has a formal protected status, it is vital to consider thresholds that balance ecological significance with spatial coverage, ensuring the protection of both rare hotspots and areas of moderate richness that may serve as vital habitats or corridors (Hansen et al. 2020). The CRFD threshold yielded the highest sensitivity (90%), followed by standard deviation and top 20% (both 83%). The top 5% threshold performed the least effectively (29%), highlighting the trade-off between spatial precision and ecological inclusivity (Fig. 8). The CRFD inflection point (~ 6 species) emerges as an ideal threshold for ths purpose, as it marks the point at which species richness transitions from common to rare within the landscape. Similar inflection points analyses have been used to justify thresholds for identifying HCVAs, as they capture ecologically meaningful richness patterns, while still representing a more inclusive range of habitats (Kukkala and Moilanen 2012). Thus, using the CRFD inflection point allows for a balanced approach — one that captures significant hotspots without excluding moderately rich habitats, aligning well with conservation best practices and ensuring long‑term resilience of biodiversity within the protected area.