https://doi.org/10.19184/geosi.v8i1.31862 Research Article Species Distribution Modelling Using Bioclimatic Variables on Endangered Endemic Species (Bubalus depressicornis and Bubalus quarlesi) Septianto Aldiansyah 1 * , Khalil Abdul Wahid 2 1 Department of Geography, Faculty of Mathematics and Natural Sciences, Universitas Indonesia, Depok, West Java, 16424, Indonesia 2 Department of Geography, Faculty of Social Science, Universitas Negeri Malang, Jl. Semarang 5, Malang, East Java, 65145, Indonesia *Corresponding author, Email address : septiantoaldiansyah863@gmail.com INTRODUCTION The populations of Lowland Anoa (Bubalus depressicornis) and Mountain Anoa (Bubalus quarlesi) in Sulawesi have declined by less than 2,500 adults with a population loss rate of up to 20% in the last 14 -18 years (Burton et al., 2016a; Burton et al., 2016b; Arini et al., 2020) caused by poaching (O'Brien & Kinnaird, 1996; Burton et al., 2005), and deforestation (Mustari, 2019; Arini et al., 2021), thus making anoa move to specially protected areas (Burton et al., 2005; Steinmetz et al., 2014; Mustari, 2019). The IUCN categorizes the status of these two anoas as Endangered Species (EN). This status indicates a high level of threat and has the potential to become extinct in the future if ABSTRACT Sulawesi Island is an island located in the Wallacea area. Most of the fauna on the island of Sulawesi is a transitional fauna from Australia and Asia. This study aims to model the potential distribution of the species Bubalus depressicornis and Bubalus quarlesi using famous models in the present and in the future as a result of climate change phenomena throughout the island of Sulawesi and beyond their natural habitat. The parameters used are bioclimatic variables and in-situ presence data. The method used is Maximum Entropy by comparing the GLM, SVM, and RF algorithms. The model is evaluated with reference to the values of AUC, COR, TSS, Deviance, and observation data. The RF model is quite good in modeling the distribution of B. depressicornis and B. quarlesi species with AUC values of 0.92 and 1, COR values of 0.59 and 0.84, TSS values of 0.87 and 1, and Deviance values of 0.37 and 0.08, respectively, while the results of data observations show values of 80% and 84%. B. depressicornis was most affected by bio14=0.665, while B. quarlesi was most affected by bio2=0.525, which means that this endemic species is suitable to live in a tropical climate with a warm and wet climate throughout the year, where the difference in temperature at night and during the day is very large. In the future, B. depressicornis and B. quarlesi are estimated to be compatible in an area of 143,281.78 km2 (81%) and 136,892.89 km2 (77%) of the Sulawesi. Keywords: Species Distribution Model; Bubalus depressicornis; Bubalus quarlesi; Bioclimatic; Climate Change ARTICLE INFO Received : 21 June 2022 Revised : 3 February 2023 Accepted : 4 March 2023 Published : 18 April 2023 . Geosfera Indonesia p-ISSN 2598-9723, e-ISSN 2614-8528 available online at : https://jurnal.unej.ac.id/index.php/GEOSI Vol. 8 No. 1, April 2023, 1-18 © 2023 by Geosfera Indonesia and Department of Geography Education, University of Jember. This is an open access article under the CC-BY-SA license (https://creativecommons.org/licenses/by-sa/4.0/) 1 https://doi.org/10.19184/geosi.v8i1.31862 https://orcid.org/0000-0002-5432-8322 https://orcid.org/0000-0002-4732-9079 https://jurnal.unej.ac.id/index.php/GEOSI https://creativecommons.org/licenses/by-sa/4.0/ 2 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 conservation action is not taken immediately. Anoa is included on CITES Appendix I, namely: one type of animal that is strictly prohibited to be hunted or traded in whole or in part. The price of anoa meat can reach IDR 22,500/kg (USD $1.57) (Burton et al., 2005). Anoa is classified as a rare and protected wild animal based on the Decree of the Minister of Agriculture of the Republic of Indonesia No: 421/KPTS/UM/8/1970, Minister of Agriculture Number: 90/KPTS/2/1972, Law Number 5 of 1990 and Government Regulation Number 7 of 1999. Anoa occupies the top position for five criteria, including endemicity, population status, habitat conditions, threats, and species management status. Species Distribution Modeling (SDM) is the most popular species distribution analysis method worldwide for predicting the probability of a species based on environmental, bioclimatic, and topographical variables (Byeon et al., 2018; Jung et al., 2019; Yoon & Lee, 2021). This model has been identified as one of the best global predictive methods when estimating the similarity of climatic conditions to the presence of known species (Phillips et al., 2006). The performance and quality of the SDM model will be better when including bioclimatic variables (partly or completely) in the prediction modeling (Truong et al. 2017; Ahmed et al., 2020; Ahmadi et al., 2020). The SDM analysis using the Bioclimatic variable has been able and proven to be able to model estimates of the distribution of living things in the future (Morales-Barbero & Alvarez, 2018; Jung et al., 2019). This is supported by more than 75% of research that has adopted SDM related to terrestrial ecosystems and always uses part or all of the Bioclimatic variables provided by WorldClim (Booth, 2018; Dyderski et al., 2018). Zhang et al. (2023) modeled the global distribution potential of Gentiana rhodantha using bioclimatic data in China which found that the degree of habitat fragmentation would increase, but the suitable area would decrease and its distribution would shift northeastward in the future. Hosni et al. (2022) found a significant correlation between SDM and the ecology of the pest Galleria mellonella using only bioclimatic data. Although SDM exhibits a wide distribution, there is an economic impact of the pest on Honeybee's population and global distribution. Bandara et al. (2022) also found SDM's ability to model the distribution of the species Kerivoula malpasi and Kerivoula picta in Sri Lanka using bioclimatic data and found the distribution of K. malpasi to be fragmented and likely to be spatially adjacent to K. picta in the future. Several previous anoa studies focused on areas with relatively narrow areas (e.g., Irawati & Arini, 2012; Wardah et al., 2012; Ranuntu & Mallombasang, 2015; Priyono et al., 2020). As a result, these studies are limited to research areas, exploration outside natural habitats, and requires time and money. This limitation is what this research tries to overcome which makes the study different from previous studies. This research is intended to explore the use of climate and environmental data as well as the in-situ distribution of species obtained from open and public data sources using the SDM for a wider scale. This model is quite effective when used to determine the current distribution and possible future potential distributions if it only uses climate data (Elith et al., 2011; Booth, 2018; Ahmed et al., 2020; Amiri et al., 2020; Andersen et al., 2022). This study aims to: 1) model the potential distribution of Bubalus depressicornis and Bubalus quarlesi not only in their natural habitat but also outside their natural habitat. Ex-situ conservation will be possible after the discovery of new areas outside the original habitat. 2) comparing several famous regression models to get an overview of the distribution of species as well as testing a regression model that is suitable for modeling the distribution of endemic species. 3) modeling the prediction of species distribution as a result of future climate change phenomena to see how much influence climate change has on the distribution of this Anoa species in its natural habitat. METHODS Study Area This study focuses on the entire area of Sulawesi Island and surrounding islands which is located between 01.7°N-05.8°S and 112.7°E-125.3°E. This area is the eleventh largest island 3 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 in the world with a projected area of 177,518.41 km2 by WGS 1984, located on the east side of Kalimantan, west of the Maluku Islands, south of the Lesser Sunda Islands and north of Mindanao (Philippines). The study area in this research is presented in the following Figure 1. Figure 1. Location of (a) Modelled area, in (b) Sulawesi, (c) Indonesia. Model Selection Maximum Entropy (MaxEnt) is the method used in this research (Phillips et al., 2006) which has proven successful in predicting the distribution of various fauna species using only species presence data (Guo & Liu, 2010; Elith et al., 2011). There are three algorithm models adopted, namely: Generalized Linear Model (GLM), Support Vector Machine (SVM), and Random Forest (RF). Each algorithm can transform nonlinear data between abnormal distributions and derived variables from independent variables (McCullagh & Nelder, 1989; Shabani et al., 2016), redundant data reading, errors during modeling, dependence on distributional assumptions (Ravand & Baghaei, 2016; Hamidi et al., 2018; Schmidt & Finan, 2018) or with large amounts of data (Howard et al., 2014), which can reduce over-fitting when predicting distributions (Ren et al., 2015; Hill et al., 2017). Species Distribution Model Data on the presence of species and environmental parameters that reflect the relevance of habitat variables will be explored using MaxEnt (Phillips et al., 2006; Elith et al., 2011). MaxEnt method runs on RStudio Software (RStudio Team, 2020) by utilizing available packages. The data used are 19 WorldClim climate data for an average of 1970-2000 with a resolution of 1 km (Worldclim, 2020) and the Climate Model Intercomparison Project version 5 (CMIP5) which is the average prediction of climate conditions for 2080-2100 (Taylor, 2009) with an accuracy of up to 0.5 km and is quite good at simulating geographic and temperature distributions (Kamruzzaman et al., 4 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 2021). Representative Concentration Pathways (RCP) data on CMIP5 used is 8.5 (high scenario) to model the distribution of species in 2092 (Taylor, 2009). WorldClim and CMIP5 spatial resolution data are reduced to 2.5 km according to the standard factor approach (Wilby et al., 2004). The model is built from worldclim data called "bio" which is called from the RStudio software. Data management is done by excluding climate data that has a high correlation value (vifcor = 0.7) and then knowing the relative importance of these variables. Prediction of future climatic conditions using climatological data models with the same variables and with a spatial resolution of up to 2.5 km (Wilby et al., 2004) for the next 70 years (Taylor, 2009). There were 296 observations for B. depressicornis and 15 observations for B. quarlesi recorded in database of The Global Biodiversity Information Facility (GBIF). However, only 15 data for B. depressicornis and 4 for B. quarlesi were used after ignoring duplicate data, coordinate errors, and data that did not have coordinates. It is assumed that each presence of data represents 1 species. The MaxEnt model is built on predictors called the “biome”. Replication (cross-validation) was applied 100 times. The resulting model is named "m". Distribution modeling is carried out based on "biom" and "m" data. The species distribution model has then been processed into ArcMap 10.4.1 software by dividing the distribution model class into 2 classes using interval equal method to see how far the suitability of climatic conditions for these species is. Model Evaluation The model of each algorithm was evaluated and compared by looking at the value of Receiver Operating Characteristics - Area Under Curve (ROC-AUC) with cross-validation (Shabani et al., 2018), Correlation (COR), True Skill Statistics (TSS) (Fourcade et al., 2018) and Deviance (Agresti, 2018). The AUC value represents the model validation based on the classification of positive and negative values for variables and presence data. A good AUC value if value is above 0.7 (Shabani et al., 2018). The evaluation was also carried out on the distribution model that had been classified using the results of observations from Mustari (2019). In addition, predictions of the distribution of B. depressicornis (Burton et al., 2016a), and B. quarlesi (Burton et al., 2016) from the IUCN Red List were also used. The research flow chart can be seen in Figure 2. 5 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Figure 2. Research flow chart RESULTS AND DISCUSSION Endemic Species Distribution Model This prediction model only displays scatter predictions based on presence data in pixels only. In Figures 3a and 4a, the GLM model does not have a low level of sensitivity to the location of the presence of low data on all species, this can be seen from the almost even distribution in all regions and generalized areas that are not species habitats. In contrast to GLM, the SVM model is slightly sensitive to non-habitat areas (low distribution = red) in B. depressicornis (Figure 3b) and B. quarlesi (Figure 4b) when viewed from the point of distribution of species presence. In Figure 3c, the RF model shows the level of sensitivity to topographic variations and selects areas that show the presence of B. depressicornis (points of presence) as the distribution area of this species, such as Rawa Aopa Watumohai National Park (RAWNP) (41) conservation area in Southeast Sulawesi Province. The RF model also highlights the conservation area of the Bogani Nani Wartabone National Park (BNWNP) (5). Areas of distribution of B. quarlesi such as the conservation area of the Gandang Dewata National Park (TNGD) (32) in Central Sulawesi are also 6 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 highlighted (Figure 4c). Both of these areas are B. depressicornis habitats in Sulawesi. In addition, this model also highlights the Lore Lindu National Park (TNLL) (25) conservation area. The BNWNP (5) was also highlighted, but with a moderate level of distribution. Figure 3. Distribution Model of Bubalus depressicornis: (a) GLM Model; (b) SVM model; (c) RF models. Figure 4. Distribution Model of Bubalus quarlesi: (a) GLM Model; (b) SVM model; (c) RF models. Determining the type of climate that is suitable for both species must go through a test of relative importance derived from variables. Based on the correlation test on the selection of existing variables. Seven variables have the potential to affect the distribution of B. depressicornis, namely Bio2: diurnal range, Bio8: temperature of wettest quarter, Bio9: temperature of driest quarter, Bio14: precipitation of driest month, Bio15: precipitation seasonality, Bio18: precipitation of warmest quarter, and Bio19: precipitation of coldest quarter. Likewise with B. quarlesi, but without the bio14 variable. The variables above represent trends in average annual temperature and annual rainfall over a period of years; seasonal range of temperature and annual rainfall; and the temperature of the coldest and warmest months, and the rainfall of the wet and dry months which are extreme environmental factors/limiting the movement of species. Based on the correlation, it shows that the distribution of B. depressicornis is influenced by bio14=0.665 and also influenced by bio9=0.036 and bio19=0.036. The training correlation shows that the distribution of B. quarlesi is most influenced by bio2=0.525 with the least influence by bio15=0.175 and bio19=0.1755. The contribution of each parameter is presented in Figure 5. 7 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 The results of the relative importance of the variables indicate that these two species are suitable for living in a tropical climate. The climate is usually warm and wet all year round and the rains are heavy and continuous. One day in an equatorial climate can be very similar to the next day. However, the gap between the hottest and coldest temperatures in 24 hours is larger than the change in the average temperature throughout the year (McKnight, 2000). Anoa has adapted perfectly to the geology, physics, climate, and landscape of the island of Sulawesi (Mustari, 2019). The typical vegetation that inhabits most of the forests on Sulawesi Island is very suitable for the life of anoa (Broto, 2015) and cannot be separated (Mustari, 2019), because it can reduce the intensity of sunlight (Aldiansyah, 2022) and maintain a macro-climate that is the habitat of anoa. Figure 5. Relative Bioclimatic Variable Importance of Bubalus depressicornis and Bubalus quarlesi. Besides being influenced by the structure and geological conditions in the Pleistocene Period (Nugraha & Hall, 2018), the endemicity of this species is influenced by climate, soil, the shape of the earth's surface, and factors from other living things (Whitten & Henderson, 2012). The endemicity of mammals in Sulawesi requires ideal climate and weather conditions, especially air temperature and rainfall so that they can support the survival of species in nature (Whitten & Henderson, 2012). A study conducted by Kasim (2002) and Jahidin (2003) explained that anoa will faint and even die if left under the scorching sun without any protection. Model Accuracy Based on data analysis, it shows that the distribution of SVM and RF models is the best, although in some cases both models have advantages and disadvantages. In B. depressicornis species, SVM and RF showed strong AUC values (Table 1). When referring to the accuracy results, SVM is the best, but based on observations from several studies, SVM makes mistakes. For example, Buton Island, which is the original habitat of B. depressicornis (Martin et al., 2012; Mustari, 2019), such as those found in Lambusango Wildlife Sanctuary (WS) and North Buton Nature Reserve (NR) (Mustari, 2019). In the SVM model, these two habitats are classified as areas with a low distribution of B. depressicornis, even though the area is the original habitat of this species (Mahmud, 2009). Although RF is not better than SVM, RF is quite sensitive in classifying the distribution of this species when viewed from the appearance of the topographical structure which is the habitat requirement of this species (Mustari, 2003; Mustari, 2019; Aldiansyah, 2022). The accuracy values of each model based on the AUC, COR, TSS, and Deviance values are presented in Table 1. 8 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Table 1. Model Accuracy Evaluation Species Model GLM SVM RF AUC Bubalus depressicornis 0.88 0.89 0.92 Bubalus quarlesi 0.99 1 1 COR Bubalus depressicornis 0.5 0.86 0.59 Bubalus quarlesi 0.85 1 0.84 TSS Bubalus depressicornis 0.8 0.83 0.87 Bubalus quarlesi 0.99 1 1 Deviance Bubalus depressicornis 3.41 0.26 0.37 Bubalus quarlesi 1.03 0.4 0.08 The distribution of the B. depressicornis species on Sulawesi Island can still develop, given the many factors that cause this species to thrive outside its natural habitat. In addition to climatic and topographical factors, other factors that can affect the distribution of B. depressicornis species are the abundance of feed types, and the availability of habitat components needed by animals, such as vegetation types and water availability (Ranuntu & Mallombasang, 2015). evenly distributed throughout the anoa's original habitat which can minimize the possibility of competition between other mammals for food (Broto, 2015), as well as species breeding activities outside their natural habitat (ex-situ) carried out in captive breeding centers or conservation forest areas (Mayasari et al., 2018). Therefore, any area that has the above factors has the potential to become a distribution area for the B. depressicornis species on the island of Sulawesi. The distribution modeling of B. quarlesi using the GLM, SVM and RF models also showed the same results as the B. depressicornis species. Some built-up areas in West Sulawesi Province and Agricultural/Plantation Areas in South Sulawesi Province were classified as having high B. quarlesi distribution. Although it has a difference in the COR and Deviance values of the RF model to the SVM model. The RF model gives quite realistic model results when compared to the SVM model. Overall, the model used displays very good results, but it can be affected by many factors such as an error such as pixel size, number of samples, and area of study. The author suspects that there are some limitations in this study, such as the pixel resolution that uses 2.5 km so that it is possible to generalize the surrounding pixels if the unit is changed to 1 km or more in detail. The distribution of the two species is very minimal and uneven. The author is also limited to three species distribution model algorithms. The selection of variables on several factors that influence their distribution in nature is not taken into account such as human activities which are a real threat to this species (O'Brien & Kinnaird, 1996; Burton et al., 2005; Mustari, 2019; Arini et al., 2020; Aldiansyah, 2022). Research conducted by Collins & McIntyre (2015) on thirty studies on modeling the distribution of species in the world that the GLM model is only used in 43% and 20% of RF. Other models that are mostly applied such as Artificial Neural Network (ANN), Boosted Regression Trees (BRT), BIOMOD, Classified Tree Analysis (CTA), Flexible Discriminant Analysis (FDA), General Additive Model (GAM), Generalized Boosted Model (GBM), Multivariate Adaptive Regression Splines (MARS), Mixture Discriminant Analysis (MDA), and Surface Range Envelopes (SRE). Different results will be obtained when using different distribution models (Shabani et al., 2016) and different results will be obtained if using the same model based on the species and study area studied. In future types of research, it is hoped that the model will be tested more diverse and the data more evenly distributed across the region. 9 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 New Distribution Model Both of our distribution models are based only on 15 (B. depressicornis) and 4 (B. quarlesi) presence points. The local distribution of the two anoa species is difficult to identify because these species occur in forest patches at different elevations or in a sympatric/parapattric manner (Burton et al., 2005). Based on the skull and bone records of B. quarlesi and B. depressicornis, as well as morphological descriptions, it is only found in the Central Sulawesi region, northern of Buton Island (Burton et al., 2005), the northern peninsula, and along the southeast peninsula (Groves, 1969). The B. depressicornis species is currently distributed on the northern peninsula of Sulawesi Island, as far east as BNWNP, in the central region, throughout the eastern and southeastern parts of the island, but this species is no longer found on the southwestern peninsula. (Burton et al., 2005). The new distribution modeled by bioclimatic variables as well as in-situ species distribution data obtained from GBIF gives promising results. Figures 6a and 6c present a distribution map of species grouped with at least 50% probability. This area includes B. depressicornis covering an area of 105,417.19 km2 (59%) and B. quarlesi covering an area of 90,139.42 km2 (51%) of the total area of Sulawesi (Fern Green color shows the distribution area). About 80% of the areas matched the presence of B. depressicornis (Figure 8a), whereas 84% of the areas matched the presence of B. quarlesi (Figure 8b). Area match values based on AUC were 0.92 for B. depressicornis and 1 for B. quarlesi (Table 1). This value shows strong results. Some conservation areas that were and are still natural habitats for B. depressicornis and B. quarlesi are classified correctly and incorrectly based on the distribution and important habitats of anoa such as conservation areas, and protected forests in Sulawesi, presented in Figure 8. Validation The distribution models of B. depressicornis (Burton et al., 2016a), and B. quarlesi (Burton et al., 2016b) in Figures 6b & 6d, have similar patterns and regions. The B. depressicornis and B. quarlesi are estimated to inhabit or fit in at 97,800.7 km2 and 70,967.1 km2, respectively. This area is smaller than the estimated area of this study. The predicted distribution model also has corridors, B. depressicornis also shows connectivity from the side of the corridor between West Sulawesi and South Sulawesi, and B. quarlesi between West Sulawesi and Central Sulawesi which did not exist in previous studies (Figures 6a & 6c). As for B. depressicornis which has a corridor that connects from South Sulawesi around Bone Bay to Southeast Sulawesi, from the data released by IUCN, this corridor is not predicted but can be identified in the new distribution model (Figure 6a & 6b). The IUCN in 2016 also recorded that there is a possibility of B. depressicornis experiencing extinction in Southeast Sulawesi and this data has similarities with the new distribution model shown in Figure 7. However, this statement is not entirely true, because in June 2022 two B. depressicornis were seen which were identified based on visible physical characteristics, namely the Anoa parent looks solid black while the cubs are yellowish brown in the vicinity of the nickel mining area of Konawe Regency (Yunus, 2022). The new distribution model covers several important distribution areas and habitats (Mustari, 2019). Areas confirmed to have B. depressicornis, such as: Mount Tinombala (17), Takolekaju Mountains (24), Nantu-Boliyohuto WS (7), RAWNP (41), Tanjung Peropa WS (43), Tanjung Amolengo WS (44), Tanjung Batikolo WS (45), North Buton NR (49) (Figure 8a). The areas confirmed to have B. quarlesi include: Mayoa-Mampi-Seko (25), Mount Pompangeo (26), Faruhumpenai NR (31), Forest in the Toraja Mountains (37), and GDNP (32) (Figure 8b). Other areas confirmed to have a mixture of B. depressicornis and B. quarlesi both from DNA fossils and the presence of species such as the Panua NR (6), LLNP (11), Kambuno Katena Protected Forest (38), Latimojong Mountains (35), Tanjung Polewali WS (42), Mangolo Nature Park (46), Grand Forest Park Tanjung Nipa Nipa (47) (Figures 8a & 8b). 10 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Figure 6. Anoa Distribution Prediction: (a) A new model of Bubalus depressicornis; (b)Prediction of Bubalus depressicornis from IUCN (2016); (c) New model Bubalus quarlesi; (d) Prediction of Bubalus quarlesi from IUCN (2016). Figure 7. Predicted distribution of Bubalus depressicornis and probability of extinction: (a) Sulawesi; (b) Southeast Sulawesi. 11 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Figure 8. The prediction of the distribution of Anoa overlaps with the results of observations for the species identified by Mustari (2019) in Sulawesi: (a) Bubalus depressicornis; (b) Bubalus quarlesi. Future Species Distribution Scenarios Distribution scenarios for B. depressicornis and B. quarlesi are still viable in most areas of Sulawesi if viewed from future climatic conditions. It is estimated that B. depressicornis and B. quarlesi still fit in an area of 143,281.78 km2 (81%) and 136,892.89 km2 (77%) of Sulawesi's area. However, both of these species experienced a reduction in the distribution area. The area, including Buton Island and the surrounding islands, is an area that is predicted to be unsuitable as a habitat for B. depressicornis from a climatic aspect. Other areas include the Pati-Pati WS (13), RAWNP (41), Tanjung Peropa WS (43), Tanjung Amolengo Wildlife Refuge (44), Tanjung Batikolo WS (45), Tanjung Polewali WS (42), North Buton NR (49), and Lambusango WS (50) also cannot become a habitat for B. depressicornis due to climate change (Figure 9a). The Lambuyan- Pangimanan WS (14) is not suitable as a habitat for B. quarlesi (Figure 9b). According to Leclerc et al. (2020), Sulawesi is potentially very vulnerable to climate change. In addition to being a threat to island ecosystems (Leclerc et al., 2018), climate change will have an impact on populations in nature (Leclerc et al., 2020). Sulawesi has never been identified by climate change for the mammal category (Pacifici et al., 2018) or other taxonomic groups (Foden et al., 2013). Anoa tends to be picky in habitat and food as a result of their inability to tolerate even small climate changes (Leclerc et al., 2020). If mammals that reproduce quickly are more at risk of extinction (González-Suárez et al., 2013), anoa that reproduces slowly will be more vulnerable to extinction because they are known to only give birth to 1 child/birth (Mustari, 2019; Mustari, 2020). This became serious after this study was conducted. This study can provide an overview of the first actions for conservation efforts in their natural habitat. 12 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Figure 9. Distribution Prediction of Anoa 2092: (a) Bubalus depressicornis; (b) Bubalus quarlesi. CONCLUSION Three models can visualize the distribution of species based on presence data and bioclimatic variables only. The best model based on the comparison of statistical tests and observational data shows that the RF model is quite good for modeling species distribution. Variables that have a high correlation with climate are bio2, bio8, bio9, bio14, bio15, bio18, and bio19. The results of the relative importance show that B. depressicornis and B. quarlesi are only suitable for living in a tropical climate with a warm and wet climate throughout the year, where the difference in temperature at night and during the day is very large. There is a corridor connecting the B. quarlesi area, which is between West Sulawesi and South Sulawesi, and B. quarlesi between West Sulawesi and Central Sulawesi. Ex-Situ conservation actions can be carried out on these two species because most of Sulawesi's areas have suitable habitats from a climate perspective. However, as climate change occurs, some areas that are native habitats are becoming more vulnerable. In the future, B. depressicornis is estimated to be compatible in an area of 143,281.78 km2 or equivalent to 81% of the total area, while B. quarlesi is estimated to be suitable for an area of 136,892.89 km2 or equivalent to 77% of total area of Sulawesi. The habitat areas that are thought to be unsuitable include the Pati-Pati WS, RAWNP, Tanjung Peropa WS, Tanjung Peropa Wildlife Reserve (WR), Tanjung Batikolo WS, Tanjung Polewali WS, Preserve North Buton NR, Lambusango WS on B. depressicornis, and Lambuyan-Pangimanan WS on B. quarlesi. ACKNOWLEDGMENTS The authors extend their appreciation to Dean of Faculty of Mathematics and Natural Sciences, Universitas Indonesia. 13 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 DECLARATIONS Conflict of Interest The authors declare that in the research and preparation of this article, there are no conflict of interests related to certain organizations, institutions, and individuals or groups. Ethical Approval On behalf of all authors, the corresponding author states that the paper satisfies Ethical Standards conditions, no human participants, or animals are involved in the research. Informed Consent On behalf of all authors, the corresponding author states that no human participants are involved in the research and, therefore, informed consent is not required by them. DATA AVAILABILITY Data used to support the findings of this study are available from the corresponding author upon request. REFERENCES Agresti, A. (2018). An introduction to categorical data analysis. John Wiley & Sons. Ahmadi, K., Alavi, S. J., Amiri, G. Z., Hosseini, S. M., Serra-Diaz, J. M., & Svenning, J. C. (2020). The potential impact of future climate on the distribution of European yew (Taxus baccata L.) in the Hyrcanian Forest region (Iran). International Journal of Biometeorology, 64(9), 1451– 1462. https://doi.org/10.1007/s00484-020-01922-z. Ahmed, N., Atzberger, C., & Zewdie, W. (2020). Integration of remote sensing and bioclimatic data for prediction of invasive species distribution in data-poor regions: a review on challenges and opportunities. Environmental Systems Research, 9(1), 1-18. https://doi.org/10.1186/s40068-020-00195-0. Aldiansyah, A. (2022). Kajian Spasial Konservasi Anoa (Bubalus Quarlesi Ouwens 1910) Pada Skala Lanskap Di Taman Nasional Gandang Dewata, Sulawesi Barat. Tesis. Universitas Indonesia. Amiri, M., Tarkesh, M., Jafari, R., & Jetschke, G. (2020). Bioclimatic variables from precipitation and temperature records vs. remote sensing-based bioclimatic variables: Which side can perform better in species distribution modeling?. Ecological Informatics, 57, 101060. https://doi.org/10.1016/j.ecoinf.2020.101060. Andersen, D., Litvinchuk, S. N., Jang, H. J., Jiang, J., Koo, K. S., Maslova, I., ... & Borzée, A. (2022). Incorporation of latitude-adjusted bioclimatic variables increases accuracy in species distribution models. Ecological Modelling, 469, 109986. https://doi.org/10.1016/j.ecolmodel.2022.109986. Arini, D.I.D., Christita, M., Sheherazade, N., Mayasari, A., Suryaningsih, R., & Simamora, A. T.A.J. (2020). A review of anoa conservation efforts in Sulawesi, Indonesia. IOP Conference Series: Earth and Environmental Science, 533(1), 012003. https://doi.org/10.1088/1755- 1315/533/1/012003. https://doi.org/10.1016/j.ecoinf.2020.101060 https://doi.org/10.1016/j.ecolmodel.2022.109986 https://doi.org/10.1088/1755-1315/533/1/012003 https://doi.org/10.1088/1755-1315/533/1/012003 14 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Bandara, A. M. J., Madurapperuma, B. D., Edirisinghe, G., Gabadage, D., Botejue, M., & Surasinghe, T. D. (2022). Bioclimatic Envelopes for Two Bat Species from a Tropical Island: Insights on Current and Future Distribution from Ecological Niche Modeling. Diversity, 14(7), 506. https://doi.org/10.3390/d14070506. Booth, T.H. (2018). Why understanding the pioneering and continuing contributions of BIOCLIM to species distribution modelling is important. Austral Ecology, 43(8), 852- 860.https://doi.org/10.1111/aec.12628. Broto, B.W. (2015). Struktur dan komposisi vegetasi habitat anoa (Bubalus spp.) di Hutan Lindung Pegunungan Mekongga, Kolaka, Sulawesi Tenggara. Seminar Nasional Masyarakat Biodiversitas Indonesia, 1 (3), https://doi.org/10.13057/psnmbi/m010339. Burton, J.A., Hedges, S., & Mustari, A.H. (2005). The taxonomic status, distribution and conservation of the lowland anoa Bubalus depressicornis and mountain anoa Bubalus quarlesi. Mammal Review, 35(1), 25-50. Burton, J., Wheeler, P., & Mustari, A.H. (2016a). Bubalus depressicornis. The IUCN Red List of Threatened Species, 2016 (e.T3126A46364222), 14. https://doi.org/http://dx.doi.org/ 10.2305/IUCN.UK.2016-2.RLTS.T3126A46364222.en. Burton, J., Wheeler, P., & Mustari, A.H. (2016b). Bubalus quarlesi. The IUCN Red List of Threatened Species 2016: e.T3128A46364433. The IUCN Red List of Threatened Species 2016, 8235, 14. https://doi.org/10.2305/IUCN.UK.2016-2.RLTS.T3128A46364433.en. Byeon, D. hyeon, Jung, S., & Lee, W. H. (2018). Review of CLIMEX and MaxEnt for studying species distribution in South Korea. Journal of Asia-Pacific Biodiversity, 11(3), 325–333. https://doi.org/10.1016/j.japb.2018.06.002. Collins, S.D., & McIntyre, N.E. (2015). Modeling the distribution of odonates: a review. Freshwater Science, 34(3), 1144-1158. https://doi.org/10.1086/682688. Dyderski, M. K., Paź, S., Frelich, L. E., & Jagodziński, A. M. (2018). How much does climate change threaten European forest tree species distributions?. Global Change Biology, 24(3), 1150- 1163. https://doi.org/10.1111/gcb.13925. Elith, J., Phillips, S. J., Hastie, T., Dudík, M., Chee, Y. E., & Yates, C. J. (2011). A Statistical Explanation of MaxEnt for Ecologists. Diversity and Distributions, 17(1), 43-57. https://doi.org/10.1111/j.1472-4642.2010.00725.x. Foden, W.B., Butchart, S.H., Stuart, S. N., Vié, J.C., Akçakaya, H.R., Angulo, A., ... & Mace, G.M. (2013). Identifying the world's most climate change vulnerable species: a systematic trait- based assessment of all birds, amphibians and corals. PloS One, 8(6), e65427. https://doi.org/10.1371/journal.pone.0065427. Fourcade, Y., Besnard, A. G., & Secondi, J. (2018). Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics. Global Ecology and Biogeography, 27(2), 245-256. https://doi.org/10.1111/geb.12684. https://doi.org/10.1111/aec.12628 https://doi.org/http:/dx.doi.org/10.2305/IUCN.UK.2016-2.RLTS.T3126A46364222.en https://doi.org/http:/dx.doi.org/10.2305/IUCN.UK.2016-2.RLTS.T3126A46364222.en https://doi.org/http:/dx.doi.org/%2010.2305/IUCN.UK.2016-2.RLTS.T3126A46364222.en https://doi.org/10.2305/IUCN.UK.2016-2.RLTS.T3128A46364433.en https://doi.org/10.1111/j.1472-4642.2010.00725.x https://doi.org/10.1111/j.1472-4642.2010.00725.x https://doi.org/10.1371/journal.pone.0065427 https://doi.org/10.1371/journal.pone.0065427 https://doi.org/10.1111/geb.12684 15 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Froese, G.Z., Contasti, A. L., Mustari, A.H., & Brodie, J.F. (2015). Disturbance impacts on large rain- forest vertebrates differ with edge type and regional context in Sulawesi, Indonesia. Journal of Tropical Ecology, 31(6), 509-517. DOI: https://doi.org/10.1017/S0266467415000450. González-Suárez, M., Gómez, A., & Revilla, E. (2013). Which intrinsic traits predict vulnerability to extinction depends on the actual threatening processes. Ecosphere, 4(6), 1-16. https://doi.org/10.1890/ES12-00380.1. Groves, C.P. (1969). Systematics of the anoa (Mammalia, Bovidae). Beaufortia, 17, 1-12. Guo, Q., & Liu, Y. (2010). ModEco: an integrated software package for ecological niche modeling. Ecography, 33(4), 637-642. https://doi.org/10.1111/j.1600-0587.2010.06416.x. Hamidi, O., Tapak, L., Abbasi, H., & Maryanaji, Z. (2018). Application of random forest time series, support vector regression and multivariate adaptive regression splines models in prediction of snowfall (a case study of Alvand in the middle Zagros, Iran). Theoretical and Applied Climatology, 134(3), 769-776. https://doi.org/10.1007/s00704-017-2300-9. Hill, L., Hector, A., Hemery, G., Smart, S., Tanadini, M., & Brown, N. (2017). Abundance distributions for tree species in Great Britain: A two‐stage approach to modeling abundance using species distribution modeling and random forest. Ecology and Evolution, 7(4), 1043-1056.. https://doi.org/10.1002/ece3.2661. Hosni, E. M., Al-Khalaf, A. A., Nasser, M. G., Abou-Shaara, H. F., & Radwan, M. H. (2022). Modeling the potential global distribution of honeybee pest, Galleria mellonella under changing climate. Insects, 13(5), 484. https://doi.org/10.3390/insects13050484. Howard, C., P.A. Stephens, J.W. Pearce‐Higgins, R.D. Gregory & S.G. Willis (2014). Improving species distribution models: the value of data on abundance. Methods in Ecology and Evolution 5(6): 506–513. https://doi.org/10.1111/2041-210X.12184. Irawati, D., & Arini, D. (2012). Anoa dan Habitatnya di Sulawesi Utara. Balai Penelitian Kehutanan Manado. www.Bpk-Manado.Litbang.Dephut.Go.Id. [19 Jan 2022]. Jahidin. 2003. Populasi dan Perilaku Anoa Pegunungan (Bubalus (Anoa) quarlesi Ouwens) di Taman Nasional Lore Lindu. Tesis. Institut Pertanian Bogor. Bogor. Jung, J. M., Byeon, D. hyeon, Jung, S., & Lee, W. H. (2019). Effect of climate change on the potential distribution of the common cutworm (Spodoptera litura) in South Korea. Entomological Research, 49(12), 519–528. https://doi.org/10.1111/1748-5967.12398. Kamruzzaman, M., Shahid, S., Islam, A.R.M., Hwang, S., Cho, J., Zaman, M., ... & Hossain, M. (2021). Comparison of CMIP6 and CMIP5 model performance in simulating historical precipitation and temperature in Bangladesh: a preliminary study. Theoretical and Applied Climatology, 145(3), 1385-1406. https://doi.org/10.1007/s00704-021-03691-0. Kasim, K. 2002. Potensi anoa (Bubalus depressicornis dan Bubalus quarlesi) sebagai alternatif satwa budidaya dalam mengatasi kepunahannya. Disertasi. Program Pascasarjana Institut Pertanian Bogor, Bogor. Leclerc, C., Courchamp, F., & Bellard, C. (2018). Insular threat associations within taxa worldwide. Scientific reports, 8(1), 1-8. DOI: https://doi.org/10.1038/s41598-018-24733-0. https://doi.org/10.1017/S0266467415000450 https://doi.org/10.1890/ES12-00380.1 https://doi.org/10.1890/ES12-00380.1 https://doi.org/10.1111/j.1600-0587.2010.06416.x https://doi.org/10.1002/ece3.2661 https://doi.org/10.1111/2041-210X.12184 http://www.bpk-manado.litbang.dephut.go.id/ https://doi.org/10.1038/s41598-018-24733-0 16 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Leclerc, C., Courchamp, F., & Bellard, C. (2020). Future climate change vulnerability of endemic island mammals. Nature communications, 11(1), 1-9. DOI: https://doi.org/10.1038/s41467-020- 18740-x. Mahmud, W. (2009). Kelimpahan Anoa Dataran Rendah (Bubalus depressicornis) dan Faktor yang Mempengaruhinya di Hutan Lambusango Pulau Buton, Sulawesi Tenggara. Skripsi. Universitas Gadjah Mada. Martin, T.E., Kelly, D.J., Keogh, N. T., Heriyadi, D., Singer, H.A., & Blackburn, G.A. (2012). The avifauna of Lambusango Forest Reserve, Buton Island, south-east Sulawesi, with additional sightings from southern Buton. Forktail, 28, 107-112. Mayasari, A., Arini, D.I.D., Suryawan, A., Suryaningsih, R., Vernia, R.E., Abinawanto, A., Suryanda, A., & Bowolaksono, A. (2018). Sexual behaviour of lowland anoa (Bubalus depressicornis) in the captivity of anoa breeding centre (ABC) Manado. MATEC Web of Conferences, 197(620), 1–5. https://doi.org/10.1051/matecconf/201819706007. McCullagh, P. & J.A. Nelder (1989). Generalized Linear Models. Chapman and Hall, London, 511pp. McKnight, T.L. (2000). Climate zones and types. Physical geography: a landscape appreciation. Morales-Barbero, J., & Vega-Álvarez, J. (2019). Input matters matter: Bioclimatic consistency to map more reliable species distribution models. Methods in Ecology and Evolution, 10(2), 212– 224. https://doi.org/10.1111/2041-210X.13124. Mustari, A.H. (2003). Ecology and Conservation of Lowland Anoa, Bubalus depressicornis, in Sulawesi, Indonesia. University of New England. Mustari, A.H. (2019). Ekologi, Perilaku, dan Konservasi Anoa. PT Penerbit IPB Press. Nugraha, A. M. S., & Hall, R. (2018). Late cenozoic palaeogeography of Sulawesi, Indonesia. Palaeogeography, Palaeoclimatology, Palaeoecology, 490, 191-209. https://doi.org/10.1016/j.palaeo.2017.10.033. O'Brien, T. G., & Kinnaird, M. F. (1996). Changing populations of birds and mammals in North Sulawesi. Oryx, 30(2), 150-156. https://doi.org/10.1017/S0030605300021530. Pacifici, M., Visconti, P., & Rondinini, C. (2018). A framework for the identification of hotspots of climate change risk for mammals. Global Change Biology, 24(4), 1626-1636. https://doi.org/10.1111/gcb.13942. Phillips, S.J., Anderson, R.P. & Schapire, R.E. (2006) Maximum entropy modeling of species geographic distributions. Ecological Modelling 190(3): 231–259. https://doi.org/10.1016/j.ecolmodel.2005.03.026. Priyono, D.S., Solihin, D.D., Farajallah, A., & Purwantara, B. (2020). The first complete mitochondrial genome sequence of the endangered mountain anoa (Bubalus quarlesi) (Artiodactyla: Bovidae) and phylogenetic analysis. Journal of Asia-Pacific Biodiversity, 13(2), 123–133. https://doi.org/10.1016/J.JAPB.2020.01.006. https://doi.org/10.1038/s41467-020-18740-x https://doi.org/10.1038/s41467-020-18740-x https://doi.org/10.1051/matecconf/201819706007 https://doi.org/10.1111/2041-210X.13124 https://doi.org/10.1016/j.palaeo.2017.10.033 https://doi.org/10.1016/j.palaeo.2017.10.033 https://doi.org/10.1017/S0030605300021530 https://doi.org/10.1111/gcb.13942 https://doi.org/10.1111/gcb.13942 https://doi.org/10.1016/j.ecolmodel.2005.03.026 https://doi.org/10.1016/j.ecolmodel.2005.03.026 https://doi.org/10.1016/J.JAPB.2020.01.006 17 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Ranuntu, R., & Mallombasang, S. (2015). Studi Populasi Dan Habitat Anoa (Bubalus sp) di Kawasan Hutan Lindung Desa Sangginora Kabupaten Poso. Mitra Sains, 3(2), 81-94. Retrieved from http://mrtg.untad.ac.id/index.php/MitraSains/article/view/129. Ren, S., Cao, X., Wei, Y., & Sun, J. (2015). Global refinement of random forest. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (pp. 723-730). Ravand, H., & Baghaei, P. (2016). Partial least squares structural equation modeling with R. Practical Assessment, Research, and Evaluation, 21(1), 11. https://doi.org/10.7275/d2fa-qv48. RStudio Team (2020). RStudio: Integrated Development for R. RStudio, PBC, Boston, MA URL http://www.rstudio.com/. Schmidt, A. F., & Finan, C. (2018). Linear regression and the normality assumption. Journal of Clinical Epidemiology, 98, 146-151. https://doi.org/10.1016/j.jclinepi.2017.12.006. Shabani, F., Kumar, L., & Ahmadi, M. (2018). Assessing accuracy methods of species distribution models: AUC, specificity, sensitivity and the true skill statistic. Global Journal of Human Social Science, 18(1), 6-18. Shabani, F., Kumar, L., & Ahmadi, M. (2016). A comparison of absolute performance of different correlative and mechanistic species distribution models in an independent area. Ecology and Evolution, 6(16), 5973-5986. https://doi.org/10.1002/ece3.2332. Steinmetz, R., Srirattanaporn, S., Mor‐Tip, J., & Seuaturien, N. (2014). Can community outreach alleviate poaching pressure and recover wildlife in South‐East Asian protected areas?. Journal of Applied Ecology, 51(6), 1469-1478. https://doi.org/10.1111/1365-2664.12239. Taylor, K.E. (2009). A summary of the CMIP5 experiment design. Retrieved from http://cmip-pcmdi. llnl. gov/cmip5/docs/Taylor_CMIP5_design. pdf. Truong, T. T. A., Hardy, G. E. S. J., & Andrew, M. E. (2017). Contemporary remotely sensed data products refine invasive plants risk mapping in data poor regions. Frontiers in Plant Science, 8(May). https://doi.org/10.3389/fpls.2017.00770. Wardah, W., Labiro, E., Massiri, S. D., Sustri, S., & Mursidin, M. (2012). Vegetasi kunci habitat Anoa di Cagar Alam Pangi Binangga, Sulawesi Tengah. Jurnal Penelitian Kehutanan Wallacea, 1(1), 1-12. http://dx.doi.org/10.18330/jwallacea.2012.vol1iss1pp1-12. Whitten, T., & Henderson, G. S. (2012). Ecology of Sulawesi. Tuttle Publishing. Wilby, R.L., Charles, S.P., Zorita, E., Timbal, B., Whetton, P., & Mearns, L.O. (2004). Guidelines for use of climate scenarios developed from statistical downscaling methods. Supporting material of the Intergovernmental Panel on Climate Change, available from the DDC of IPCC TGCIA, 27. Worldclim. (2020). Bioclimatic Variable. Retrieved from https://www.worldclim.org/data/bioclim.html. [25 Jan 2022]. http://mrtg.untad.ac.id/index.php/MitraSains/article/view/129 https://doi.org/10.1016/j.jclinepi.2017.12.006 https://doi.org/10.1002/ece3.2332 https://doi.org/10.1111/1365-2664.12239 http://dx.doi.org/10.18330/jwallacea.2012.vol1iss1pp1-12 https://www.worldclim.org/data/bioclim.html 18 Septianto Aldiansyah & Khalid Abdul Wahid / Geosfera Indonesia 8 (1), 2023, 1-18 Yoon, S., & Lee, W. H. (2021). Methodological analysis of bioclimatic variable selection in species distribution modeling with application to agricultural pests (Metcalfa pruinosa and Spodoptera litura). Computers and Electronics in Agriculture, 190, 106430. https://doi.org/10.1016/j.compag.2021.106430. Yunus, S.R. (2022). Habitat Terganggu, Anoa Berkeliaran di Lokasi Perusahaan Tambang di Konawe. Kompas ID Online News June 7th 2022. https://www.kompas.id/baca/nusantara/2022/06/07/habitat-terganggu-anoa-muncul-di- lokasi-perusahaan-tambang-di-konawe. [20 June 2022]. Zhang, H., Sun, X., Zhang, G., Zhang, X., Miao, Y., Zhang, M., ... & Huang, L. (2023). Potential Global Distribution of the Habitat of Endangered Gentiana rhodantha Franch: Predictions Based on MaxEnt Ecological Niche Modeling. Sustainability, 15(1), 631. https://doi.org/10.3390/su15010631. https://doi.org/10.1016/j.compag.2021.106430 https://www.kompas.id/baca/nusantara/2022/06/07/habitat-terganggu-anoa-muncul-di-lokasi-perusahaan-tambang-di-konawe https://www.kompas.id/baca/nusantara/2022/06/07/habitat-terganggu-anoa-muncul-di-lokasi-perusahaan-tambang-di-konawe METHODS Study Area Model Selection Species Distribution Model Model Evaluation Model Accuracy New Distribution Model Validation Future Species Distribution Scenarios CONCLUSION