Driven to the edge: Species distribution modeling of a Clawed Salamander (Hynobiidae: Onychodactylus koreanus) predicts range shifts and drastic decrease of suitable habitats in response to climate change

Abstract Climate change is one of the major threats to global amphibian diversity, and consequently, the species distribution is expected to shift considerably in the future. Therefore, predicting such shifts is important to guide conservation and management plans. Here, we used eight independent environmental variables and four representative concentration pathways (RCPs) to model the current and future habitat suitability of the Korean clawed salamander (Onychodactylus koreanus) and then defined the dispersal limits of the species using cost distance analysis. The current habitat suitability model generated using the maximum entropy algorithm was highly consistent with the known distribution of the species and had good predictive performance. Projections onto years 2050 and 2070 predicted a drastic decrease of habitat suitability across all RCPs, with up to 90.1% decrease of suitable area and 98.0% decrease of optimal area predicted from binary presence grids. The models also predicted a northeastward shift of habitat suitability toward high‐elevation areas and a persistence of suitability along the central ridge of the Baekdudaegan Range. This area is likely to become a climatic refugium for the species in the future, and it should be considered as an area of conservation priority. Therefore, we urge further ecological studies and population monitoring to be conducted across the range of O. koreanus. The vulnerability to rapid climate change is also shared by other congeneric species, and assessing the impacts of climate change on these other species is needed to better conserve this unique lineage of salamanders.

Among vertebrates, amphibians are especially prone to the detrimental effects of climate change. Their ecology is closely tied to requirements for humidity and temperature for cutaneous respiration (Carey & Alexander, 2003), and many species undergo biphasic life histories with an aquatic stage at least during their larval phase (Duellman & Trueb, 1994). These aquatic ecological traits and related habitat requirements can be easily disturbed by rapid climate change (Bickford et al., 2010). With 40% of all amphibian species already threatened (IUCN, 2020), the effects of climate change are likely to exacerbate the problems caused by other known threats including habitat destruction, infectious diseases, overexploitation, and pollution (Beebee & Griffiths, 2005;Gallant et al., 2007;Lips, 1998;Lips et al., 2006;Rowly et al., 2016;Sodhi et al., 2008;Wake, 2012).
Therefore, predicting range shifts as a consequence of climate change is important to better guide management, policy-making, and conservation decisions, as well as to determine important areas for conservation.
Correlative species distribution models (SDMs) have become the main method for conducting such predictions ). Generated by correlating species' distribution and a set of environmental predictors, SDMs provide a representation of continuous habitat suitability across the landscape Guisan & Zimmermann, 2000). The SDMs can then be projected onto different space to locate suitable habitats outside the species' known distributions (Pyron et al., 2008;Raffini et al., 2020;Srivastava et al., 2021), as well as to different time periods to predict range shifts under climate change (Borzée, Andersen, et al., 2019;Hu et al., 2020;Mothes et al., 2020). The development of user-friendly interfaces and statistical packages such as the software MaxEnt have greatly enhanced the ability to implement such projections , opening the doors for this method to be applied to other species with presumed vulnerability to climate change.
Here, we applied correlative species distribution modeling to estimate current and future habitat suitability of the Korean clawed salamander (Onychodactylus koreanus), a species of hynobiid salamander endemic to the Korean Peninsula (Poyarkov et al., 2012;Lee & Park, 2016). This species is primarily found in the mountainous regions of the Korean Peninsula, mainly along the mountains branching from the Baekdudaegan Range, in close association with clear low-temperature streams and forests ( Figure 1). Due to their characteristics and ecology, members of the genus Onychodactylus are highly likely to become threatened by climate change (Suk et al., 2018), similarly to other Korean salamanders (Borzée, Andersen, et al., 2019). First, unlike other hynobiids, the members of this genus are lungless, making them more stringently bound to moist and cool environments (Kuzmin, 1995;Poyarkov et al., 2012).
Second, unlike most plethodontids, a diverse family of lungless salamanders, Onychodactylus goes through aquatic reproduction and larval phase (Park, 2005). This requires unpolluted and lowtemperature mountainous or subterranean water bodies with high levels of dissolved oxygen (Hong, 2017;Jung, 2020;Kuzmin, 1995;  (e) Lee & Park, 2016), environments that are prone to degradation by climate change. Previous studies on plethodontids have predicted range shifts and the loss of diversity hotspots due to climate change (Borzée, Andersen, et al., 2019;Milanovich et al., 2010;Parra-Olea et al., 2005), and similar patterns are expected for Onychodactylus.
However, while the detrimental effects of climate change on the genus have been speculated (Suk et al., 2018), no study has applied SDMs to investigate potential range shifts and future habitat suitability in response to climate change to our knowledge.
The effects of climate change on the distribution of O. koreanus will also have conservation implications for the future.
In  Figure 2). While the occurrence dataset is naturally correlated temporally due to ecological characteristics of the species (e.g., cryptic amphibian with increased activity and detection during rainy months), the risk of temporal correlation leading to biased occurrence sampling is unlikely in our case. This is because our target species is a strict habitat specialist that occupies narrow environments along mountain streams and forests immediately adjacent to them, such that the observation of larvae in a given area is a good indicator of the presence of adults in the same area, and vice versa (Hong, 2017;Lee & Park, 2016). Therefore, accounting for spatial sampling bias was more important in our case than considering temporal bias, and we accounted for spatial sampling bias by spatial occurrence thinning and modified background selection (see Model preparation section below). and 2020, strictly enforcing temporal match between occurrence and climate data would mean that only a small subset of occurrence points could be used for modeling. As this can lead to biased model estimates, we consider the use of the full occurrence dataset to be more beneficial than enforcing temporal match between data, while acknowledging the presence of the temporal mismatch.

| Model preparation
In addition to the bioclimatic variables, we included two topographic variables (slope and elevation), three forest cover variables (needleleaf forest cover, broadleaf forest cover, and mixed/other trees cover), two variables representing human activities (urban landscape cover and human footprint), and a raster representing the distance to the nearest water bodies, as these variables are highly relevant to the ecology of O. koreanus (Kuzmin, 1995;Poyarkov et al., 2012;Lee & Park, 2016 (Hijmans, 2019) and rgdal (Bivand et al., 2021).
Maximum entropy (implemented in the software MaxEnt) is a robust method of species distribution modeling with presence-only data . However, one of the major issues with this type of data is sampling bias (Kramer-Schadt et al., 2013;Phillips et al., 2009). As a presence-only algorithm, MaxEnt requires a set of background points for model construction (Merow et al., 2013).
By default, 10,000 background points are randomly selected across the entire extent of the study, but this does not take sampling bias into account (Phillips et al., 2009). As our occurrence data had both regional (for GBIF and NES data) and national (between R. Korea and D.P.R. Korea) biases (Syfert et al., 2013), we employed target group background sampling (Merow et al., 2013;Phillips et al., 2009 Shin, Min, et al., 2021).
After spatially thinning this dataset using a 1-km distance parameter in the R package spThin (Aiello-Lammens et al., 2015), 4,523 unique data points remained (https://doi.org/10.17632/ yc3nw 4d9f2.2; Shin, Min, et al., 2021). We converted these data points into a density raster reflecting sampling intensity across the Korean Peninsula and used this raster for bias correction in subsequent MaxEnt modeling.
To reduce the clustering of occurrence points, we spatially thinned our occurrence data so that each data point is at least one kilometer apart, thereby putting only one occurrence point within each 30-arc-second (~1 km) grid cell of environmental layers. This procedure was conducted in R using the package spThin (Aiello-Lammens et al., 2015) and reduced the number of data points from 922 to 187 occurrence data used to calibrate the models (https:// doi.org/10.17632/ yc3nw 4d9f2.2; Shin, Min, et al., 2021). After reducing spatial clustering, the occurrence points still encompassed the totality of the known range of O. koreanus. Also, considering the known capability of MaxEnt algorithm to predict habitat suitability even with small sample sizes (Pearson et al., 2007), our spatially thinned occurrence dataset is sufficiently large to generate models with good predictive ability.
We used a two-step approach to select environmental variables to calibrate our models. First, we conducted five bootstrap replicated runs in MaxEnt version 3.4.1 (https://biodi versi tyinf ormat ics. amnh.org/open_sourc e/maxen t/) with the spatially thinned occurrence dataset, 27 environmental layers, auto features, and a default regularization multiplier (= 1). For background point selection, we used the density raster to sample 10,000 target group background points. From this initial run, we retained nine variables with contribution scores greater than 3%. To reduce multicollinearity among the retained variables, we then conducted a Spearman's correlation test in the R package ntbox (Osorio-Olvera et al., 2020) and removed highly correlated variables if | coefficient | > 0.7. This resulted in the following set of eight environmental variables used to estimate the current habitat suitability of O. koreanus: bio1 (annual mean temperature), bio3 (isothermality), bio4 (temperature seasonality), bio13 (precipitation of wettest month), bio14 (precipitation of driest month), slope, needleleaf forest cover, and mixed/other trees cover.
As O. koreanus is a lungless ectotherm living in and around mountain streams, its distribution is expected to be strongly governed by precipitation, humidity, temperature, and vegetation cover. Therefore, we considered the selection of the eight environmental variables to be appropriate. While human-related variables were eventually excluded from the model calibration due to the initial filtering criterion, the effects of human activities can still be captured by other environmental variables (e.g., low forest cover in highly populated area).
Nevertheless, we included human-related variables to calculate cost rasters for the estimation of dispersal limits as human activities are still important to determine the species' distribution and dispersal limits (see Estimation of dispersal limits section below).
To optimize model parameters for calibration, we tested a combination of six MaxEnt feature classes (Linear, Linear Quadratic, Hinge, Linear Quadratic Hinge, Linear Quadratic Hinge Product, and Linear Quadratic Hinge Product Threshold), and 16 regularization multipliers (ranging from 0.5 to 8 with a 0.5 increment) under two cross-validation schemes: 10-fold cross-validation and spatial block cross-validation. While the former method is commonly used to generate and test SDMs, we also tested the latter method as its use has been recommended when model transfer to different environmental conditions is the goal (Muscarella et al., 2014;Kass et al., 2021). We used 10,000 target group background points to generate candidate models under each cross-validation method. To select optimal model parameters from the sets of candidate models, we applied and compared two selection criteria.
First, we selected the model with a combination of feature classes and regularization multiplier that had the lowest AICc value, describing model fit and complexity (Warren & Seifert, 2011;Muscarella et al., 2014). Second, we selected the model that had the lowest 10% omission rate as well as the highest test AUC (sequential selection method; Kass et al., 2021). For both crossvalidation methods, the optimal parameters based on AICc were the LQHPT features combined with a regularization multiplier of 1.5. On the other hand, based on the sequential selection method, the optimal parameters under the spatial block cross-validation was the L feature with a regularization multiplier of 4.5, while the optimal parameters under 10-fold cross-validation were the H feature combined with a regularization multiplier of 3.5. To select the final optimal model parameters, we projected the habitat suitability predicted from these three parameters across the landscape. As the northern distribution limit of O. koreanus in D.P.R.
Korea remains uncertain due to the presence of closely related species in the country (Poyarkov et al., 2012;Borzée, Litvinchuk, Ri, et al., 2021), selecting model parameters that predicted habitat suitability in D.P.R. Korea with minimum overprediction was key to generate accurate current habitat suitability model. According to these criteria, we selected the model parameters with the lowest AICc (LQHPT 1.5) to calibrate the final current habitat suitability model. We conducted the model tuning procedure using R packages dismo  and ENMeval (Muscarella et al., 2014;Kass et al., 2021).

| Current habitat suitability
We used MaxEnt for our modeling because of the algorithm's ability to explicitly account for spatial bias in the occurrence dataset (via a density raster; Shin, Messenger, et al., 2021). We calibrated the current habitat suitability model with 10,000 target group background points and 10-fold cross-validation (Thapa et al., 2018). We applied all MaxEnt feature classes (LQHPT) with a regularization multiplier of 1.5, as determined by the ENMeval run. We selected cloglog as the output format (Borzée, Andersen, et al., 2019;Di Pasquale et al., 2020;McDonald et al., 2021;Phillips et al., 2017), and we took the average of the 10 replicated models for the final model evaluation. We also used a jackknife analysis to assess the contribution of each variable to the model.

| Future projections
For future projections we used the Community Climate System been widely applied to predict future habitat suitability of terrestrial vertebrates (Hu et al., 2020;Morovati et al., 2020;Sutton et al., 2015) and was also effective in predicting future habitat suitability of While bioclimatic variables were projected onto future conditions, we maintained topographical and land-cover variables constant throughout all projections due to the lack of future scenarios corresponding to climate change for these variables. Although vegetation cover is likely to shift with climate change (Choi et al., 2011), thereby introducing a bias in our analyses, adding these variables was still favorable to better characterize the ecological niche of O. koreanus compared to using bioclimatic variables only.

| Estimation of dispersal limits
One of the most significant limitations of SDMs is the inability of the modeling algorithms to explicitly account for the dispersal ability of organisms, and this can lead to unrealistic estimates of the current and future habitat suitability (Guisan & Thuiller, 2005;Ofori et al., 2017;Zhang, Mammola, et al., 2020). Although the dispersal ability of O. koreanus is unknown, it is likely to be limited as the species is a lungless ectotherm and as its distribution is tightly linked to forest environments. Thus, we tried to overcome this limitation by estimating the cost of dispersal through various landscape features (Sahlean et al., 2014).
To do so, we first generated a set cost rasters, the values of which represent the difficulty (cost) of moving through the land-  (Kuzmin, 1995;Park, 2005;Song & Lee, 2009;Poyarkov et al., 2012;Hong, 2017;Lee & Park, 2016;Jung, 2020; and our knowledge of the species' biology. We also relied on the literature for specific information, such as the known upper elevation limit for the species (Song & Lee, 2009

| Current habitat suitability
Based on evaluation metrics, our model is acceptable as an ad-

| Dispersal limits
Our estimation of dispersal limits under the current climate condition was largely matching with the current known distribution of the species. Also, the transition of habitat suitability in central D.P.R.
Korea predicted by the current habitat suitability model was generally matching with the northern dispersal limits estimated from the cost distance analysis.
Regarding future projections, all RCPs across both time frames where the habitat is environmentally suitable but unreachable. It is also notable that while future dispersal limits fluctuated slightly between time periods and scenarios, no great expansion or contraction of dispersal limits was predicted (Figure 4).

| Area of suitable and optimal habitats
The total area of current suitable habitat estimated within the dispersal limits was 32,593.1 km 2 and that of the current optimal habitat was 10,729.4 km 2 . As with continuous habitat suitability, the binary estimation of suitable and optimal habitats was highlighted by a drastic contraction. The models predicted sharp decreases of both suitable and optimal habitats between the present and 2050, with further decreases predicted for 2070 ( Figure 5)

| D ISCUSS I ON
In agreement with our hypothesis, the models predicted drastic Increased annual mean temperature and increased fluctuations in precipitation are likely driving the loss of suitable habitat and range shifts predicted in our models. Such environmental changes triggered by rapid climate change are likely to cause habitat degradation by altering a suite of important habitat requirements for the species. For example, the species is highly dependent on lowtemperature mountain forests receiving abundant precipitation (Poyarkov et al., 2012;Lee & Park, 2016). However, future climate change scenarios for R. Korea indicate an increased severity of droughts (Kim et al., 2013), as well as an increase in extreme precipitation events (Sung et al., 2018). Therefore, increasing fluctuations in temperature and precipitation are likely to create habitat conditions unsuitable for the species by causing desiccation and flooding of habitats and affecting the survival of both larvae and adults. In addition, the breeding environment of the species is strongly associated with unpolluted, low-temperature mountain streams with high levels of dissolved oxygen (Hong, 2017;Jung, 2020;Lee & Park, 2016;Park, 2005), and therefore, larvae of the species are highly dependent on these requirements (Hong, 2017;Jung, 2020;Kuzmin, 1995). The amount of dissolved oxygen decreases with increasing temperature in the aquatic environments (Ficklin et al., 2013;Matear et al., 2000), and such decreases in dissolved oxygen can negatively impact amphibian communities (Bickford et al., 2010). Therefore, warming climate can result in degradation of the mountain stream environments by reducing the amount of dissolved oxygen. The combination of these factors is particularly negative to the larvae considering the extended aquatic larval stage of the species (2-3 years; Lee & Park, 2016). Furthermore, climate change can also cause climatic disturbances and subsequent degradation of subterranean environments used by the species for reproduction (Mammola et al., 2018(Mammola et al., , 2019Park, 2005;Zhang, Mammola, et al., 2020).

F I G U R E 5
In addition, complex interactions between climate change and other factors not accounted for in our models are likely to have negative impacts on the population dynamics of the species. For example, while our models assumed constancy of land-cover variables through time due to the lack of applicable future scenarios, forest environments will most likely be affected by climate change (Choi et al., 2011;Joshi et al., 2012). In R. Korea, a drastic contraction of the forest volume has been predicted under future climate conditions (Kwak et al., 2012). In addition, cool temperate forests that are the principal habitats of O. koreanus have been predicted to be replaced by subtropical vegetation by 2100 due to climate change (Choi et al., 2011). The specific impacts of replacement and shift of vegetation on population dynamics are difficult to assess at present, but it will most likely lead to serious habitat degradation and population decline due to the loss of key habitats (Cushman, 2006).
Finer scale environmental variables that are important in determining habitat suitability, such as soil pH, are also likely to change with climate warming. Acidic soils negatively affect amphibian communities through disturbances in osmoregulation, growth, and respiration (Sugalski & Claussen, 1997;Wyman & Jancola, 1992), and the habitat requirements of O. koreanus reflect avoidance of environments with low soil pH (Hong, 2017;Jung, 2020). Soil acidification can be exacerbated by altered precipitation patterns caused by climate change (Rengel, 2011), resulting in decreased habitat suitability. While microhabitat features and plasticity can provide buffering effects to climate change (Huey et al., 2018;Scheffers et al., 2014), habitat degradation at multiple levels is still likely to put severe pressure on the species in the future.
Other anthropogenic stressors are additional threats to the species. Widespread threats such as deforestation, pollution, road construction, and conversion of landscapes to agriculture put a strong pressure on amphibian populations by hampering dispersal, reproduction, foraging, and other behaviors associated with species survival (Beebee & Griffiths, 2005;Lee & Miller-Rushing, 2014;Shin, Jeong, et al., 2020). Such development of the habitat is ongoing and causing widespread environmental degradation across the Korean Peninsula (Lee & Miller-Rushing, 2014;Park et al., 2014).
This is particularly more problematic for R. Korea, where the rate of development is high (Borzée, 2021;Lee & Miller-Rushing, 2014) (Iskandar & Erdelen, 2006). Therefore, we urge additional ecological studies and population monitoring on the species to be conducted and, ideally, the establishment of long-term population monitoring programs. Such programs can be implemented in the form of regular surveys conducted by researchers and can also benefit from the participation of citizen scientists (Borzée et al., 2018;Borzée, Baek, et al., 2019;Deutsch et al., 2017).
Our results are also applicable to other species of Onychodactylus.
While we used O. koreanus as an example, the potential negative impacts of climate change on future habitat suitability are not unique to the species. All ten known species of Onychodactylus are lungless and share largely similar environmental requirements for survival (Poyarkov et al., 2012;. Therefore, all of these species are likely to face similar challenges under rapid climate change (Pyron, 2018). Onychodactylus represents an ancient lineage of salamanders (Frost et al., 2006;Pyron & Wiens, 2011), and assessing the potential impacts of climate change and other anthropogenic pressures on the genus is required to better understand and conserve the unique evolutionary diversity of this lineage.

ACK N OWLED G M ENTS
We thank the associate editor Dr. Luciano Bosso and two anonymous reviewers for helpful comments that greatly improved the initial draft of this study. We also thank Il-Kook Park for technical advice on GIS visualizations of MaxEnt outputs and also for helpful discussions on the concepts of ecological niche modeling.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this study are available from the Mendeley Data repository (https://doi.org/10.17632/ yc3nw 4d9f2.2; Shin, Min, et al., 2021) and the appendices accompanying this manuscript.

A PPE N D I X A
The original and reclassified values of landscape variables used to create cost rasters. We reclassified the original raster values based on our knowledge of the species. The "thresholds" specified for reclassification of habitat suitability output models are based on maximum training sensitivity plus specificity thresholds computed in MaxEnt runs. The cost rasters resulting from these reclassified rasters were used as inputs for cost distance analyses to estimate dispersal limits of Onychodactylus koreanus under each climate change scenario (see text for methods).