The study site for all climatic predictions was restricted to latitudes of 20°30’ – 22°30'N and longitudes of 106° – 108°30'E (Fig. 1), with a background area of 43,400 km², covering all known distributions of Q. acanthophora in northeast Vietnam and its potential occurrence in southern China (Dubois and Ohler 2009; Hecht et al. 2013; Nguyen et al. 2025).

Two coordinates of Q. acanthophora from the Global Biodiversity Information Facility (GBIF) were downloaded (GBIF.org 2026), using the “rgbif” package (Chamberlain et al. 2024) in R v4.6.0 (R Core Team 2026). Additionally, interviews with local people and forest rangers identified 32 streams likely occupied by Q. acanthophora, and 42 coordinates from 16 streams were subsequently verified by our direct observations during field surveys (Pham et al. in press). Only one occurrence locality of Q. acanthophora was randomly selected within each square kilometer, using the “dismo” (v1.3-8) package in R (Aiello-Lammens et al. 2015; Hijmans et al. 2022). Such spatial filtering can improve the quality of prediction models by decreasing geographical bias and autocorrelation effects (Radosavljevic and Anderson 2014; Kiedrzynski et al. 2017). Finally, 20 representative occurrence records of Q. acanthophora were employed for further analyses (Fig. 1).

To predict the current potential distribution of Q. acanthophora, we obtained 19 bioclimatic variables from the Worldclim 2.1 at 30 arc seconds resolution (Fick and Hijmans 2017). For future projections, we used three global circulation models (GCMs) from the sixth Coupled Model Intercomparison Project (CMIP6): Meteorological Research Institute Earth System Model (MRI-ESM2-0) developed by the Meteorological Research Institute in Japan (Eyring et al. 2016; Yukimoto et al. 2019), Interdisciplinary Research on Climate version 6 (MIROC6) developed by the Atmosphere and Ocean Research Institute, the National Institute for Environmental Studies, and the Japan Agency for Marine-Earth Science and Technology, Japan (Tatebe et al. 2019), and MR) developed by the Beijing Climate Center, China Meteorological Administration (Wu et al. 2019). These climate models are able to capture regional monsoon dynamics, to simulate ocean-atmosphere interactions, and to predict tropical cyclone behavior. Three climate change scenarios of the Shared Socioeconomic Pathways (SSPs) were downloaded to assess the different level of climate change impacts over 20-year periods (i.e. 2050s: 2041–2060; 2070s: 2061–2080; 2090s: 2081–2100): SSP2-4.5 (“middle of the road”: temperatures are projected to rise by 2.7 °C by the end of the century), SSP3-7.0 (“regional rivalry”: to rise by around 3.2 °C to 4.1 °C) and SSP5-8.5 (“fossil-fuelled development”: to rise by around 4.7 °C to 5.1 °C) (O’Neill et al. 2016; Riahi et al. 2017). In this study, we only presented one projected map for each future climate scenario, by averaging the raster cell values of suitability across the three global circulation models (GCMs).

To select climatic conditions for predictions, we first employed 19 variables to import into the MaxEnt software v. 3.4.4 (Phillips et al. 2006), and applied the default settings to identify the variables that contribute most to the model (Suppl. material 1: table S2). Next, a pairwise Pearson correlation analysis of selected variables was conducted to mitigate the influence of potential multicollinearity, using the “usdm” (v2.1-7) package in R to identify and eliminate highly correlated environmental variables with coefficient values larger than |0.75| (Dormann et al. 2013; Naimi 2015). Our final dataset comprised six variables, including Mean Diurnal Range (Bio2), Isothermality (Bio3), Temperature Seasonality (Bio4), Max Temperature of Warmest Month (Bio5), Mean Temperature of Coldest Quarter (Bio11), Precipitation of Driest Quarter (Bio17), which were finally selected for subsequent analyses (Suppl. material 1: table S3).

Occurrence data of Q. acanthophora and six selected climate variables were imported into the MAXENT software to predict the potential distribution at present and to assess its shifts in all future scenarios.

For the species within a restricted distribution, defining the study background to include environmentally dissimilar areas beyond the accessible region of Q. acanthophora may artificially inflate the AUC (Area Under the Curve) and specificity values by increasing the contrast between presence and background locations. Thus, a buffer layer (18,492 km²) with a radius of 40 km around occurrence points with dissolved boundaries was created by using the R packages “dismo” and “raster” (v3.6-30). The buffer layer is expected to encompass an extent of occurrence (EOO) of Q. acanthophora within the study area (Fig. 2A) (IUCN SSC Amphibian Specialist Group 2017; Nguyen et al. 2024; Pham et al. in press). We randomly sampled 10,000 background points within the buffer area using the “randomPoints” function of the “dismo” package (Hijmans et al. 2022), and this approach can reduce the potential for artificial inflation of AUC values and provides a more realistic evaluation of model discrimination. Next, we applied the Jackknife method to generate a default model. We then used this default model to test model hyperparameters with six feature combinations (L, Q, H, LQ, QH and LQH; where L = linear, Q = quadratic, H = hinge) and regularization values ranging from 0.5 to 5 with a 0.5 increment, using the “gridSearch” function of the “SDMtune” package (v1.3.1) in R. A total of 60 candidate models were evaluated using the AIC test (Akaike Information Criterion), and the optimal model was selected with the lowest AICc value (i.e., delta. AICc = 0) (Warren and Seifert 2011). Based on these filtering criteria, the final SDM for Q. acanthophora was generated with the feature class of LQH and the regularization value of 1.5.

Three types of replication procedures in MaxEnt with bootstrap (100 replicates), cross-validation and subsampling were tested (Phillips et al. 2006). The optimal procedure was selected based on its performance across multiple evaluation metrics of AUC, omission rate, sensitivity, and specificity. Consequently, the replication method of bootstrap was selected for all predictions (Suppl. material 1: table S1). Lastly, we selected the cloglog output format, in which values ranging from 0 to 1 represent cloglog-transformed suitability. The AUC score was used to evaluate the predictive performance, ranging from 0 to 1 (perfect performance). The contribution of each variable to the current climatic model was determined by measuring the permutation importance and contribution percentage (Phillips and Dudík 2008). The minimum training presence cloglog threshold and 10 percentile training presence cloglog threshold were both selected in Maxent to convert continuous predicted values to maps of suitable and highly suitable areas, respectively, which were separated to unsuitable areas (Phillips et al. 2006; Peterson 2011).

To assess the predictive capabilities, multivariate environmental similarity surface (MESS) analyses were further used to quantify areas characterized by potential extrapolation errors (Elith et al. 2010). MESS analysis compares the environmental similarity of any given grid cell projection area to the environmental background, which was used to train the models. We also presented one projected MESS map only for each future climate scenario, by averaging the raster cell values across the three global circulation models (GCMs). In particular, similarity/novelty was classified into three levels: “< 0” (Extrapolation), “=0” (Margin) and “>0” (Interpolation).

The priority conservation areas for Q. acanthophora, which were identified within suitable habitats in terms of climates, were overlaid with the available layer of tropical evergreen rain forests within restricted latitudes of 21°–22°N and longitudes of 106°–108°E (JAXA EORC 2025). For assessing the conservation priorities, information on protected areas was included as shapefiles of nature reserves and national parks in Vietnam, which were further overlaid with highly suitable habitats of Q. acanthophora at present, using the higher threshold of 10% training presence cloglog. All maps were generated using Quantum GIS software (QGIS version 4.0.1), showing the study site, layers of suitable habitats, tropical evergreen rain forests and protected areas.

Quasipaa acanthophora only inhabits evergreen forests of northern Vietnam within the tropical monsoon zone, characterized by a seasonal variation in climatic conditions (Dubois and Ohler 2009; Hecht et al. 2013; Nguyen et al. 2025). The MaxEnt model indicated that temperature-related climatic variables, particularly the Max Temperature of Warmest Month (Bio5) during the rainy season in northern Vietnam and Temperature Seasonality (Bio4), play an important role in predicting the potential distribution of Q. acanthophora, accounting for 90% of the contribution. Amphibians, including Q. acanthophora, in tropical forests are highly sensitive to temperature changes and seasonal variability that strongly influence their physiological performance and distribution patterns (Buckley and Jetz 2007; Sunday et al. 2012; Scheffers et al. 2016). Pham et al. (in press) demonstrated that species’ temperature is significantly positively correlated with environmental conditions, including ambient temperatures with a range from 18.6 °C – 28.0 °C across four surveyed months (i.e., March, May, June and November). Humidity along the streams where the species occurs was relatively high, ranging from 68% to 99% (mean ± SE: 86.0 ± 1.2%) (Pham et al. in press). These variables likely reflect the species’ dependence on mild climatic conditions associated with streams in the mountain forests.

The climatic model predicted suitable habitats, covering the natural distribution of Q. acanthophora across the mountain forests of Bac Ninh, Quang Ninh and Lang Son provinces, northeast Vietnam. Based on the predicted distribution of optimal environmental conditions, forest regions in Cao Bang and Thai Nguyen provinces (Vietnam), as well as Guangxi Province (China) should be prioritized for future field surveys to identify new populations of Q. acanthophora. Based on previous predictions for the potential distribution of Goniurosaurus lichtenfelderi, a tiger gecko species occurring in stream-associated habitats in syntopy with Q. acanthophora, extending potentially to the northeastern Vietnam – southwestern China border region, Ngo et al. (2022) subsequently conducted field surveys and confirmed its presence in the predicted areas.

In addition, islands of Bai Tu Long National Park, Ha Long Bay (Quang Ninh Province), and Cat Ba National Park (Hai Phong City) are suggested as the potential distribution of Q. acanthophora. However, herpetofaunal surveys have not recorded the species in these insular habitats (Gawor et al. 2015; Ngo pers. obs.). Geographic isolation (i.e., mainland and islands) and particular restriction to stream microhabitats likely limit the species’ dispersal and colonization to these islands (Heaney 2007; Whittaker et al. 2008; Lomolino et al. 2010; Pham et al. in press).

Future climate change scenarios project a significant contraction in suitable habitats of Q. acanthophora, restricted to a very small region along the China–Vietnam border. Under the most extreme scenario with an increase of the global temperature, approximately 4.7 °C –5.1 °C, during the period 2060–2100, climatically suitable habitats for the species are projected to shrink severely. The pattern of contraction under climate change was noted for all tiger geckos (Goniurosaurus spp.) in northeastern Vietnam as well (Le et al. 2017; Ngo et al. 2021a, 2022, 2025). These projections underscore the potentially negative impacts of climate change not only on Q. acanthophora and Goniurosaurus species, but also on other threatened herpetofauna (including Tylototriton spp. and Shinisaurus crocodilurus) occurring in the study site (Bernardes et al. 2013; Nguyen et al. 2022). A previous study on Asian warty newts (Paramesotriton spp.) in northern Vietnam has shown that climate change can produce contrasting responses among genetic lineages, with the spatial distribution of lineage-specific populations either contracting or expanding, and reduced gene flow increasing vulnerability in certain regions (Tran et al. 2025). Thus, integrating population genetic data with species distribution models should therefore be considered in future studies for the target species of Q. acanthophora.

A total of 20 occurrence records were used to predict the potential distribution of Q. acanthophora and its projected changes under future climate scenarios. However, the number of occurrence records was very limited because the species inhabits a highly restricted distribution (Nguyen et al. 2024; Pham et al. in press), which may reduce the robustness of the model predictions. To improve the accuracy of models, we restricted the background points to a smaller buffer surrounding the occurrence records rather than using the entire study area. This approach reduces the artificial inflation of AUC and specificity values that can arise from the broad background extent and provides a more realistic assessment of model discrimination. We also averaged projections across three global circulation models (GCMs) to reduce uncertainty associated with a single climate model only. Additionally, three MaxEnt replication procedures—bootstrap (100 replicates), cross-validation, and subsampling were compared, and the best-performing procedure was selected for the final model. Finally, we selected the optimal feature class and regularization multiplier through model tuning rather than relying on the default MaxEnt settings.

Like many ectothermic species with a restricted distribution, Q. acanthophora is predicted to be particularly vulnerable to climate change (Luedtke et al. 2023; Hayden and Blom 2024; Li et al. 2024; Ngo et al. 2025). Extinction risk to the species may take place, even more rapidly under direct human impacts, such as overharvesting for food consumption and forest degradation (e.g., conversion to agriculture, mining, tourism development) (Pham et al. in press), and such impacts have been recorded in other endemic species in northeast Vietnam as well (e.g., Tylototriton vietnamensis, Goniurosaurus spp., Shinisaurus crocodilurus) (Bernardes et al. 2013; van Schingen et al. 2015, 2016; Ngo et al. 2021b, 2022, 2025). In fact, the current climate change under these impacts has increased the frequency, severity, and damage of natural disasters (e.g., storms, floods) on these species in northeast Vietnam (OCHA 2024; Kieu and Tran 2025; World Bank 2025; Pham et al. in press).

Most amphibian species, particularly narrowly distributed endemics (including Q. acanthophora), are unlikely to cope with climate change by shifting their distribution to newly suitable stream habitats within evergreen forests under future climate scenarios. This limitation is primarily due to narrow distribution ranges, specialized habitats and limited dispersal ability (Menin et al. 2011; Sunday et al. 2011). Furthermore, future models have not incorporated land-use change, yet continued human activities are likely to cause further declines in the natural forest in Vietnam. As a result, even areas projected to remain climatically suitable may no longer provide optimal forest conditions for the target species in the future. The preservation of forest habitats is therefore critical, as they maintain stable microclimatic conditions and buffer Q. acanthophora against broader macroclimatic changes projected under future climate change scenarios. As noted by Ngo et al. (2025), forest protection is also suggested as a conservation measure for the persistence of threatened herpetofauna in northeastern Vietnam, particularly under the combined pressures of future climate change and ongoing anthropogenic impacts.

Identifying core habitats is a crucial step to enhance the effectiveness of conservation measures to protect wild populations and natural habitats of the target species – Q. acanthophora. Based on the current prediction of potential distribution combined with layers of tropical evergreen rain forests and protected areas, we propose two areas, including Tay Yen Tu Nature Reserve (Bac Ninh) and Dong Son–Ky Thuong Nature Reserve (Quang Ninh), serving as the priority conservation areas of Q. acanthophora (Fig. 6).