Distribution and habitat assessments of the Slender racer, Orientocoluber spinalis, for the registration of nationally endangered species in the Republic of Korea

Conservation assessments are essential for preserving biodiversity. However, many reptile species have not been evaluated owing to data deficiencies. The Slender racer (Orientocoluber spinalis) is threatened in four out of six inhabiting countries. However, despite its apparent rarity and data deficiency, the International Union for Conservation of Nature (IUCN) has classified it as a Least Concern. In this study, we combined field surveys, habitat analysis, and ecological niche models (ENMs) to identify the critical habitat characteristics of O. spinalis, evaluate its distribution status in the Republic of Korea, and register it as a nationally endangered species. Across the country, we found a few small populations on the mainland but large populations on the islands. Orientocoluber spinalis is mainly found in low-altitude ecotone habitats between grasslands and forests. Based on previous genetic and climatic studies, we propose designating it as an endangered species to conserve this species in protected areas such as national parks, and its non-isolated mainland populations can be preserved as source populations.


Results
Data description. Initially, we collected 246 individual location data points, comprising 154 from databases and personal records and 92 from our direct field surveys (Table 1). After applying the criteria, we used 117 data points from databases and personal data and six from field collection data. Among the 246 individual location data points, 172 (69.9%) were located on the islands and 74 (30.1%) on the mainland. During the field surveys, we observed snakes mainly on Oeyeon Island (n = 51), Udo Island (n = 25), Ui Island (n = 8), and in Woraksan National Park (n = 6; Fig. 1). We filtered overlapping data points within a 1 km radius to obtain 123 data points (70 islands and 53 mainlands) for downstream analyses. Many island points were excluded from the analysis because of the high density of island populations relative to the mainland populations.
Habitat characteristics. The average values of the four environmental variables were as follows: altitude of 110.1 ± 159.4 m, 7.2 ± 6.1° slope, 13.0 ± 1.2 °C annual mean temperature, and 1,368.4 ± 337.6 mm annual precipitation. Habitat of O. spinalis (r = 100 m) was primarily composed of forest (37.7%; 1.13 ± 0.99 ha), followed by Table 1. Location data by 17 administrative districts and the habitable area in the ensemble model. To reduce bias, only one data was selected among the overlapped data within a radius of 1 km and used to analyze the habitat characteristics and for building the ecological niche models habitable area. The percentage indicates the habitable area relative to each partial area.  Distribution characteristics. Among the three base models, random forest (RF) model showed the highest accuracy and transferability with the highest area under the receiver operating characteristic curve (AUC) and true skill statistic (TSS) values, followed by boosted regression trees (BRT, also known as the generalized boosting method/GBM) and maximum entropy (MaxEnt; Table 3). The ensemble model had the most accurate predictability with AUC of 0.966 and TSS of 0.805 compared to the three other ENMs. Further, the ensemble model indicated that altitude had the greatest contribution (34.6%), followed by the annual mean temperature (21.3%), distance to forest (20.0%), slope (13.3%), distance to grassland (5.5%), and annual precipitation (5.4%).
The ensemble model predicted that O. spinalis had only 11,217 km 2 as a suitable habitat area, which represented 13.5% of the total evaluated area of 83,049 km 2 ( Table 1). The habitable area for the mainland area was 8462 km 2 (10.8%) out of 78,422 km 2 , while the island area had a habitable area of 2755 km 2 (59.5%) out of 4627 km 2 . Habitable areas were mainly located in the western and southern coastal regions, comprising great agricultural plains and islands in the southwestern sea, while the areas were limited to the northern and eastern mainland regions (Fig. 3). In addition, habitable areas were found in several large cities, such as Gwangju, Busan, and Ulsan. www.nature.com/scientificreports/  Figure 2. Eight land cover types of Orientocoluber spinalis habitat. The areas of each land cover type within a 100 m radius of O. spinalis location data (a) and distribution density map (b) which shows the distance from the location data to each major land cover type. This graphs were generated using the package ggplot2 (https:// ggplo t2. tidyv erse. org) in R v.4.0.5 software package (R Core Team, https:// www.r-proje ct. org/). www.nature.com/scientificreports/

Discussion
The combined application of field surveys and high-resolution land cover analysis enabled us to characterize the key habitat components of O. spinalis and develop empirical high-resolution ENMs. The ensemble model showed the highest AUC and TSS values, providing reliable results. This approach allowed for more precise identification of the distribution pattern of O. spinalis and a realistic threat assessment. The predominant land cover types surrounding the location point of O. spinalis were forest, grassland, and dry cropland. Among these, the closest land cover type to the location point was grassland, followed by forest and bare land. These results suggest that O. spinalis primarily inhabited grassland, cropland, or bare land adjoining the forest edge. Based on our field surveys and a previous radiotelemetry study 38 , O. spinalis was predominantly observed in grasslands adjacent to forests. The average altitude of the individual location data that made the greatest contribution to ENMs was 110 m a.s.l. Based on the topographical results and habitat analyses, we inferred that the preferred habitat of O. spinalis consisted of grassland, cropland, or bare land adjoining the edges of relatively flat lowland forests. Such habitat usage characteristics are often reported in O. spinalis 38 and various other grassland snake species, including the Eastern racer (Coluber constrictor), Grass snake (Natrix natrix), and the Common garter snake (Thamnophis sirtalis) [39][40][41][42] .
Despite the ENM, indicating large suitable habitats for O. spinalis in the southwestern and southeastern parts of the Korean Peninsula, particularly in agricultural plains like the Naju and Gimhae Plains, our field observations revealed low occurrences of O. spinalis in these areas. We speculate that O. spinalis inhabits lowland adjoining forests for two possible reasons. First, anthropogenic activities can fundamentally change suitable habitats into hostile ones. Irreplaceable reptile habitats are closely related to anthropogenic pressure 43 . Most land in the plains has experienced arable land rearrangements, such that the land's topography has greatly changed. Following paddy cultivation, the basking and feeding areas were highly disrupted. Such disturbances could also negatively influence populations of small snakes or lizards, which are the major components of O. spinalis diet 33,44 . Second, habitable areas for O. spinalis are conveniently accessible for anthropogenic activities. Since mountainous forests account for 64.5% of the total land area in the Republic of Korea 45 , anthropogenic development is largely concentrated in the lowlands. Therefore, the grassland area in the Republic of Korea decreased by approximately 50% in 2020 compared to 1995 46 . These disturbances would have led to further habitat fragmentation and destruction. Apart from agricultural plains, we documented several individuals of O. spinalis from protected mainland areas, such as Woraksan and Byeonsanbando National Parks and Seonunsan Provincial Park. Despite the lack of major components, including large grasslands or dry croplands, offered by key habitats, these protected areas provided suitable basking and foraging sites. In addition, even though suitable areas were found in several large www.nature.com/scientificreports/ urban areas, we rarely found O. spinalis at the outermost edge of the cities, indicating that cities were unsuitable for O. spinalis. Contrary to the limited observations in mainland areas, we recorded numerous individuals of O. spinalis on the islands. We suppose several reasons for this observation. Firstly, many islands with small mountainous forests, adjacent grasslands, or dry croplands are highly habitable for O. spinalis, as shown in our models. These environments provided sufficient resources, such as basking and foraging areas for O. spinalis 38 . For example, the Japanese keelback (Hebius vibakari) and Tsushima smooth skink (Scincella vandenburghi), which are known prey, are enough for islands 44 . Secondly, anthropogenic activities in islands are limited owing to the low human population density compared to mainland areas. In addition, dry cropland, which can be used as an alternative for O. spinalis to grasslands, is used as arable land 40,47 , instead of rice paddies, which are unsuitable in small islands. Third, suitable habitats are usually limited to small islands, thereby increasing the observation chances of O. spinalis.
The priority of conservation efforts of either island or mainland populations of O. spinalis is a key question. We identified habitable areas on the islands, agricultural plains, and a few forested areas and observed many, very few, and few individuals in the order. However, despite the abundance of island populations in the southwestern sea, island populations pose several challenges in terms of long-term conservation. First, island populations are vulnerable to demographic variability 48,49 . They are isolated from each other by the sea barrier, limiting gene flow among them and often having very low genetic diversity 50,51 . In our previous study, 27 individuals on Oeyeon Island had only one haplotype for the mitochondrial Cytb and ND4 genes 52 . Second, the island population could be highly vulnerable to the expected climate change, which could move the entire O. spinalis distribution to the north 30 . Since island populations cannot shift their latitudinal distribution range, progressive climate change could accelerate the growth of the entire island population in the Republic of Korea in the near future. Although the mainland, Woraksan, and Byeonsanbando National Parks, have more than five populations within 500 km 2 , the population density is very low. In our previous study, mainland populations, such as the Woraksan population, had relatively high genetic diversity, resulting in gene flow between subpopulations 52 . Although anthropogenic disturbances such as forest roads limit certain gene flows in mountainous mainland populations, the situation is better than that among island populations. Considering our current and previous results, we propose that for the long-term survival of O. spinalis in the Republic of Korea, more attention should be paid to mainland mountainous populations located in protected areas than island populations.
The mainland Republic of Korea has very few populations of O. spinalis, with stable populations reported only from Woraksan and Byeonsanbando National Parks and Seonunsan Provincial Park. Despite high population densities on islands, O. spinalis populations are isolated from each other and have low genetic diversity, with low priority for conservation purposes. Currently, the Korean O. spinalis population is threatened as to be designated as Endangered under the IUCN assessment criteria 53 , based on our current results on key habitats and distribution rates and previous results on low genetic diversity and large climate change impacts 30,52 . Therefore, we propose designating O. spinalis as an endangered species in the Republic of Korea for its long-term conservation. Future conservation efforts should be directed primarily toward valuable mainland mountainous populations, including reducing anthropogenic disturbances in the lower ranges of montane forests. Furthermore, as this is the first distribution and habitat assessment on O. spinalis worldwide, we expect that this approach can be effectively applied not only to the Republic of Korea, but also to other Asian countries such as Russia, Mongolia, and Kazakhstan where O. spinalis is also threatened.

Methods
Location data collection. We extracted and compiled location data for O. spinalis from various Korean research institutions, such as the National Institute of Ecology (n = 52) and the Korea National Park Service (n = 31); public databases, such as Inaturalist (https:// www. inatu ralist. org; n = 14) and Naturing (https:// www. natur ing. net; n = 3); and from several Korean reptile researchers' records (n = 40). We then applied the following four criteria to the compiled dataset to filter out unreliable or inaccurate data points: (1) data points recorded prior to 2000, (2) inaccurate coordinates, including regional centroids (e.g., downtown) or coordinates recorded on the sea, (3) no detailed coordinate information to the fifth decimal place, and 4) data without the name of the recorder and collection date. Additionally, we screened 106 data points collected through field surveys conducted between 2020 and 2022. During the field surveys, we considered each live individual, shed skin, and carcass as a unique data point. We arranged the obtained location data based on administrative districts comprising eight metropolitan cities and nine provinces in the Republic of Korea (Table 1, Fig. 1).
Since spatial autocorrelation may occur in adjacent location data 54 , especially in the field survey sites, we removed overlapping data points within a 1 km radius and finally selected 123 location data for the analysis 30,55,56 . To select location data, we used QGIS v.3.4.7 57 .
Habitat characteristics analysis. We considered 12 variables to characterize the habitat of O. spinalis: eight land cover types, two topographic variables, and two climatic variables. For land-cover types, we downloaded seven raster layers, each representing the cover of urban areas, cropland, forest, grassland, wetland, bare land, and freshwater bodies, at a 1 m × 1 m grid resolution (Korea Ministry of Environment; https:// egis. me. go. kr/; Fig. 1). Considering the distinctly different ecological roles of dry cropland and rice paddy as animal habitats, we further subdivided the cropland into rice paddy and dry cropland, resulting in a final set of eight land cover types; land cover data near the Northern Limit Line (NLL; near the border between South and North Korea) were not provided due to military security reasons (Fig. 1).
To identify the important land cover types for O. spinalis habitat, we calculated the area of eight land cover types within a 100 m radius of each location and the distances between each location and the eight major land www.nature.com/scientificreports/ cover types (Table 2). We considered that the area (approximately 3.14 ha) within a 100 m radius was sufficient to represent the habitat of O. spinalis, as the species' home range is approximately 2.11 ha 38 . We defined the shortest distance from the location point for each land-cover type as the distance to that land-cover type. In addition, we visualized the distribution density of snakes according to distance from each land cover type using the package ggplot2 58 in R v. 4.0.5 software package 59 .
To extract topographic variables (altitude and slope), we downloaded numerical topographic data at a 1 m × 1 m grid resolution (National Geographic Information Institution; http:// map. ngii. go. kr/). Based on the contour lines of the topographic map, we constructed a digital elevation model and extracted the altitude and slope.
For climatic variables, we selected annual mean temperature (AMT) and annual precipitation (APP), which are the most commonly used variables in distribution models implemented for reptiles 60 . We produced these climatic variables by applying the inverse distance-weighted interpolation method with average climatic data collected from 1991 to 2020 from 95 meteorological stations provided by the Korea Meteorological Administration (https:// data. kma. go. kr/). We used QGIS v.3.4.7 for the generation of environmental variables 57 .
Ecological niche modeling. Environmental variables. To identify and evaluate suitable habitats for O. spinalis in the Republic of Korea, we built ENMs using six environmental variables: distance to forest (DTF), distance to grassland (DTG), altitude, slope, AMT, and APP. Although many different environmental variables can be considered to generate ENMs, the specificity of presence data increases with increasing variables, leading to a high false negative rate 61 . Therefore, we only included the six variables with low multicollinearity (| r |< 0.6) 62 . We used DTF and DTG as indicators of the major habitats of the eight land-cover types. These two types constituted most of the habitat areas and were closest to the location data. A previous study also suggested that forests and grasslands are major factors in O. spinalis habitat 38 . Since area variables are difficult to apply in ENMs, we used only distance variables. Snakes generally prefer grasslands along forest edges and avoid deep forest interiors 39,42 ; therefore, we treated the DTF as having a higher positive value further away from the forest and a higher negative value deeper in the forest interior. However, we treated DTG as having only positive and zero values because grasslands in the Republic of Korea do not have a large continuous area and are less negatively affected by depth than forests. For the remaining variables, we selected two geographic (altitude and slope) and two climatic (AMT and APP) variables because they are the most basic variables for ENMs applied to animals 60 . Owing to the heavy computational operating load with high-resolution data of 1 m × 1 m covering the Republic of Korea, we resampled the variables to a relatively lower grid resolution of 120 m × 120 m, which is the highest spatial resolution that can be handled by the biomod2 package 63 .
Selecting ENMs. We generated the final ENM for O. spinalis within the ensemble-modeling framework. The ensemble model is a multimodal approach that addresses the shortcomings of individual ENM algorithms, reinforces their strengths 63,64 , and cross-validation 60 . Ensemble methods have been frequently used to determine habitat suitability for various snake species 65,66 . The proper selection of individual base ENMs is critical for effective habitat prediction using the ensemble model. We selected random forest (RF), boosted regression trees (BRT), and maximum entropy (MaxEnt) as base models and relatively recently developed machine-learning algorithms for interpolation with high reliability [67][68][69][70] . Random forest, an approach that combines several randomized decision tree predictors of independent samples 71,72 , has been shown to have high classification accuracy and high stability against small perturbations of data 71,73 . The BRT combines many simple regression decision trees using boosting techniques, is robust to nonlinear relationships, and can process interaction effects between predictors 74,75 . MaxEnt estimates the target distribution by calculating the convergence value (maximum entropy) based on its experimental average, which enables it to incorporate interactions between predictors and calculate the optimal probability distribution 76 . ENMs were created using the package biomod2 63 in R version 4.0.5 59 .
Model evaluation. After building the ENMs using three different algorithms, we combined these models into a final ensemble model weighting the true skill statistic (TSS) 77 cutoff value of 0.6. We generated individual ENMs in 15 bootstrap replicates and 5,000 iterations, using 123 presence and 1,000 pseudo-absence data. We randomly selected 25% of the location data for the model evaluation. We used the area under the receiver operating characteristic curve (AUC) and the TSS to evaluate the predictive performance of the ENMs. The AUC suggests the optimal correlation between sensitivity and 1-specificity and is typically used to evaluate calculations based on machine learning methods 78 . The TSS was calculated as sensitivity + specificity -1, presenting a simple and intuitive measure of ENM performance 77 . The AUC value between 0.5 and 0.7 indicates moderate predictive performance, between 0.7 and 0.9 indicates high model performance, and AUC greater than 0.9 indicates excellent model performance 79 . Meanwhile, the TSS value between 0.4 and 0.6 indicates moderate predictive performance, between 0.6 and 0.7 indicates high model performance, and TSS greater than 0.7 indicates excellent model performance 79 . We used the value maximizing the sum of sensitivity and specificity (MaxSSS) 80 to threshold the continuous model prediction output into a binary presence-absence map. The MaxSSS threshold draws a reliable interpretation because maximizing the sum of sensitivity and specificity is equivalent to minimizing the sum of false negatives and false positives 81 . We then overlaid the distribution of protected areas in the Republic of Korea (national parks, provincial parks, and county parks; Korea Database on Protected Areas; http:// www. kdpa. kr/) on the ENM to visually examine the overlap between the protected areas, location data, and suitable habitats predicted by the ENMs (Fig. 3). www.nature.com/scientificreports/

Data availability
The datasets generated and analyzed in the current study are available from the corresponding author upon reasonable request. www.nature.com/scientificreports/