Current and Future Habitat Suitability of the Endangered Java Slipper Orchid (Paphiopedilum javanicum) in Sundaland, Indonesia

2.1 Study area

This study was conducted across the Sundaland region of Indonesia, focusing on the islands of Java, Bali, and Lombok. The study area extends approximately from 105° to 116° E longitude and 6° to 9° S latitude, encompassing a wide range of climatic, topographic, and ecological conditions. Elevation ranges from coastal lowlands to montane environments exceeding 3,500 m above sea level, particularly in the volcanic mountain system of Java and Bali.

The climate of the region is tropical monsoonal, characterized by distinct wet and dry seasons, with mean annual temperatures ranging from approximately 20–30℃ and annual precipitation varying from 1,500 to over 4,000 mm depending on elevation and exposure. Vegetation types include lowland tropical rainforests, montane forests, subalpine vegetation, and fragmented secondary forests, reflecting strong environmental gradients across the islands.

The study area encompasses numerous protected areas, including national parks, nature reserves, and wildlife sanctuaries, that serve as critical refugia for biodiversity conservation in Sundaland.

2.2 Study species

Paphiopedilum javanicum is a terrestrial plant. Leaves up to 30 × 5 cm, tessellate; dark green spots often coalesce. Inflorescences consisting of a single flower held on a ± 40 cm long peduncle, which is usually longer in the wild than when the plant is cultivated. Flowers up to 10 cm broad across the petals; dorsal sepal ovate, acute, green with white margins, longitudinally striped dark brown. Petals long-oblong held nearly horizontal but sometimes drooping, greenish, finely spotted reddish brown in the lower two-thirds, white or pink in the apical third, lip dull brownish green [3].

2.3 Species occurrence data

Species occurrence records of Paphiopedilum javanicum were compiled from herbarium collections of Bali Botanic Garden and Bogor Botanic Garden and supplemented with records from the Global Biodiversity Information Facility (GBIF). All records were carefully cleaned prior to modelling. Duplicate records and points with identical geographic coordinates were removed. To minimize the effect of spatial sampling bias and spatial autocorrelation, spatial thinning was applied using a minimum distance threshold of 5 km between occurrence points. Records collected before 1970 were excluded to ensure consistency with contemporary climate data. After the data cleaning process, a total of 48 high-quality occurrence points were retained and used in the final modelling.

2.4 Environmental variables

The environmental variables used in the analysis were derived from WorldClim 2.1 (https://www.worldclim.org/) at a spatial resolution of 30 arcseconds (~1 km). The dataset included 19 bioclimatic variables (BIO1-BIO19), which represent annual trends, seasonality, and extreme or limiting environmental factors. Additionally, elevation data from WorldClim 2.1 were incorporated to account for topographic influences on species distribution.

Prior to model construction, multicollinearity among environmental variables was assessed using Pearson correlation analysis. When pairs of variables showed high correlation (|r| > 0.7), only the variable considered more ecologically relevant to the species was retained in the model to reduce redundancy and overfitting. This procedure ensured that the final set of predictor variables provided independent and robust contributions to model performance.

To project future habitat suitability, we utilized downscaled and bias-corrected climate projections from WorldClim 2.1. Future climate scenarios were based on the UKESM1-0-LL model under the Shared Socioeconomic Pathway (SSP) 5-8.5. The selected time periods included (1) 2041-2060 (UKESM1-0-LL SSP5-8.5, 30 arcseconds), (2) 2061-2080 (UKESM1-0-LL SSP5-8.5, 30 arcseconds), and (3) 2081-2100 (UKESM1-0-LL SSP5-8.5, 30 arcseconds).

2.5 Data analysis

All analyses were performed using Ecocommons (https//www.ecocommons.org.au/), a web-based ecological modelling platform that integrates the Biodiversity and Climate Change Virtual Laboratory (BCCVL). The modelling was conducted using the Maximum Entropy (MaxEnt) algorithm, which is widely used for species distribution modelling due to its ability to predict suitable habitats based on presence-only data. In this study, 60 percent of the species presence points were randomly selected for modelling, and the remaining 40 percent were used to evaluate the model. The jackknifing method was used to evaluate the contribution of the various bioclimatic variables to the model. The response curves demonstrating the relationships between the probability of plant presence and the environmental variables were evaluated [16]. Model outputs were exported from EcoCommons and mapped in ArcGIS 10.4. The ROC curve evaluates model performance by plotting sensitivity (true positive rate) against 1-specificity (false positive rate). Higher AUC values indicate better discrimination; AUC values greater than 0.7 are commonly interpreted as acceptable, greater than 0.8 as good, and greater than 0.9 as excellent.

The distribution analysis was divided into five classes based on the probability of the suitable area ranging from 0 to 1. The classification that is used here is not suitable (0-0.2) and currently not suitable (0.2-0.4), low suitable (0.4-0.6), high suitable (0.6-0.8), and very high suitable (0.8-1) [17]. The current and predicted analyses were calculated for each island based on the appropriate classification and presented as percentages of the total area. Thus, the changes were calculated using Excel to provide the changes in areas that increase or decrease between the current and future prediction [17].

Our findings revealed that the AUC value for the training data is 0.963 in both the current and future scenarios. This result is lower compared with research [14], which has an AUC training value of 0.973. Overall, the high AUC values across all time periods suggest that the model effectively differentiates suitable habitats for P. javanicum, supporting its use for forecasting future distribution trends under climate change scenarios. The AUC value will be higher for species with a narrow distribution range [11].

The environmental variables influencing the distribution of species are derived from the jackknife test. The higher the contribution, the greater the variable’s influence on the final prediction [16]. The determining factor in habitat modelling is species-specific. In selecting variables used to build a Maxent model, it is necessary to consider ecological relevance and data suitability [14]. Our findings revealed that temperature and elevation appear to have a strong influence on the habitat suitability of P. javanicum in Java, Bali, and Lombok. The result is similar to the findings of a previous study conducted on the same species in Mount Lawu, which highlighted the importance of temperature-related variables [14].

The dominance of Bio5 is closely linked to the physiological characteristics of Paphiopedilum orchids. Many species in this genus are adapted to shaded, moist, and thermally stable forest understories. Their leaves are relatively thin and not structurally adapted to withstand prolonged heat or intense radiation. In addition, successful germination and early development depend on symbiotic fungi, which are themselves sensitive to soil temperature and moisture fluctuations. High maximum temperatures during the warmest month may therefore limit not only adult plant survival but also recruitment processes. These physiological constraints explain why P. javanicum is rarely found in lowland or seasonally hot environments and instead concentrates in montane habitats with moderated thermal extremes [18]. Meanwhile, P. javanicum was also found in Mount Lawu with a temperature of 21.8-29.7℃ and an elevation of 1,231 – 1,825 m a.s.l. [14].

Elevation also contributed substantially to the model, reinforcing the role of montane environments as thermal refuges. Elevation gradients integrate multiple environmental factors, including lower temperatures, higher humidity, frequent cloud cover, and reduced thermal amplitude. These conditions create microclimates that buffer organisms from regional warming. Studies on other Paphiopedilum species similarly report strong associations with mid- to high-elevation forest zones where temperature extremes are limited and canopy cover maintains stable humidity [19].

The current potential distribution pattern shows that highly suitable habitats are concentrated in mountainous regions of West Java, Central Java, East Java, Bali, and Lombok. This is consistent with field observations indicating that P. javanicum typically occurs between approximately 1,200 and 2,500 m a.s.l. [14]. Although research [3] reported occurrences as low as 750 m, most populations are associated with cooler montane forest. The preference for such elevations reflects the species’ sensitivity to heat stress and its dependence on shaded, moist habitats.

Importantly, several mountain complexes were identified as remaining suitable across all time periods: Mount Rinjani, the mountain areas around Ulun Danu Beratan, Mount Raung, Mount Argopuro, Bromo Tengger Semeru National Park, Mount Arjuna, Mount Kawi, Mount Liman, Mount Lawu, Mount Merbabu National Park, Mount Merapi National Park, Mount Sumbing, Mount Sindoro, Mount Prau, Mount Slamet, Mount Ceremai National Park, Mount Halimun Salak National Park, Mount Gede-Pangrango National Park, Mount Cikuray, and Mount Karang. These areas can be considered long-term climate refugia for several reasons.

First, these mountain systems possess high elevation ranges, allowing upward shifts of suitable thermal zones as regional temperatures increase. Species such as P. javanicum can persist by moving to higher altitudes where maximum temperatures remain within physiological tolerance. Second, many of these areas have complex topography, producing diverse microhabitats with varying slope aspects, canopy cover, and moisture regimes. Such heterogeneity increases the probability that localized cool and humid pockets will persist even under broader climatic warming. Third, most of these mountains are located within protected areas or forest reserves, where land-use change is relatively limited compared to the surrounding lowlands. Reduced anthropogenic disturbance enhances their role as stable refuges for temperature-sensitive species.

The persistence of suitability in these montane landscapes highlights their importance as conservation priorities. Protecting forest integrity in these areas will not only conserve existing populations but also provide future habitat as climate envelopes shift. Conservation strategies should therefore focus on preventing habitat degradation, controlling illegal collection, and maintaining canopy cover and watershed functions that regulate local microclimates. In addition, these refugial zones may serve as priority sites for population reinforcement or assisted migration if low-elevation populations decline.

Species distribution modeling (SDM) represent a valuable tool to develop predictions to understand the future distribution of plant and evaluate appropriateness of the habitat for a particular species [16]. Our findings also revealed that under future climate scenarios, habitat suitability shows a shifting trend. By 2040, the distribution of highly suitable areas remains relatively stable, but shows minor expansion in some regions of West Java and Bali. However, by 2060 and 2080, a noticeable decline in suitable habitats is observed in parts of Central Java and eastern Java, indicating the potential impact of climate change on the species' distribution.

The contrasting trends in habitat suitability among regions reflect underlying differences in topography and regional climate sensitivity. West Java consistently retains the highest proportion of suitable habitat across all time periods, which can be attributed to its extensive mountainous terrain and complex topography, which create persistent cool, moist microclimates. These conditions are less vulnerable to surface warming and provide climatic refugia for montane species under future climate change.

In contrast, Central Java and East Java exhibit more pronounced declines in habitat suitability, particularly in later projections (2060–2080). These regions are dominated by lower elevations and broader lowland areas, where increasing temperatures and changes in precipitation patterns under the SSP5-8.5 high-emission scenario intensify thermal stress and reduce habitat suitability for temperature-sensitive species such as P. javanicum. This spatial divergence highlights the buffering effect of topographic complexity in West Java and underscores the importance of high-elevation landscapes as long-term refugia for conservation planning.

The percentage area of P. javanicum suitable habitat shows that West Java remains the most critical region for P. javanicum conservation. West Java has many mountains and hills, and the topography influences the environmental conditions in the region. Several studies reported the occurrence of P. javanicum in Gede-Pangrango National Park [20], Cakrabuana Mountain, Majalengka [21], and Jambu Mountain, Sumedang [22].

The identification of priority sites for species conservation is a crucial concern in plant conservation management. The application of SDM helps us to identify with a high predictive accuracy. Thus, we predict the potential distribution of Paphiopedilum javanicum based on bioclimatic variables, which allow us to highlight areas that could be transformed. The establishment of protected areas and/or the extension of the existing ones is surely the first step towards effective protection for these species [23]. The model indicates that the current total suitability area of P. javanicum habitat in the protected area is 301,872.51 Ha (3.99%) of the total area, while the projected area shows an increase in the year 2040 and 2060, but then a decrease in the year 2080. The suitable habitat is a high-altitude area that has suitable environmental conditions to support the P. javanicum growth. Policymakers should consider the influence of anthropogenic activities in the conservation management design of P. javanicum, because this plant has economic value, even though it grows in mountain areas and not inside protected areas, but many people take it from nature, so the number is limited.

In this study, the model primarily relies on bioclimatic variables and elevation, while other potentially important factors, such as soil properties, vegetation structure, canopy cover, and biotic interactions, were not explicitly incorporated due to data limitations at the regional scale. Nevertheless, climate remains the dominant driver of large-scale species distributions, particularly under future climate change scenarios. The use of spatially filtered occurrence data and variable selection procedures helped minimize model uncertainty. Future studies integrating edaphic variables, land cover, and species interactions would further improve the ecological realism of the predictions.