Frontiers of Biogeography Modeling the environmental refugia of the endangered Lichtenfelder’s Tiger Gecko (Goniurosaurus lichtenfelderi) towards implementation of transboundary conservation

Climate change has potential effects on global biodiversity by shifting the optimal distribution of terrestrial organisms, particularly species with narrow distributions. Goniurosaurus lichtenfelderi, a forest-dwelling lizard, is found on both the island and mainland sites of northern Vietnam and southern China. The species is categorized as Vulnerable in the IUCN Red List and was recently listed in CITES Appendix II and the Vietnam Government’s Decree 06 in 2019 due to severe anthropogenic impacts on its populations. In this study, we employ Maxent species distribution modeling with climatic and vegetation cover data to identify the potential distribution of G. lichtenfelderi. We also used this approach to assess future climate impacts on the potential distribution under different climate change scenarios. Our model predicts that the potential distribution of G. lichtenfelderi will shrink significantly under future scenarios and even vanish in the entire study area under novel environmental conditions of the BCC-CSM 1-1 – RCP 8.5 scenario by the 2070s. Overall, the current potential distribution is expected to shift towards higher latitudes within the next decades. The forecasted maps provide useful guidelines to implement conservation strategies to mitigate synergistic impacts from climate change and other negative anthropogenic activities. In the context of the potentially severe impacts, the border areas between China and Vietnam, Yen Tu Mountain Range, Bai Tu Long National Park, and their surroundings should be considered core refugia for the species, where conservation measures need to be prioritized in the future.


Introduction
A high proportion of terrestrial organisms from around the world will face the risk of extinction due to severe effects of climate change (Huggett 2004, Thomas et al. 2004, Parmesan 2006, Huey et al. 2009, Nogués-Bravo et al. 2010, Sandel et al. 2011, Monastersky 2014, Pimm et al. 2014, Urban 2015. Wiens (2016) documented the local extirpation of almost half of 976 investigated species during the 20 th century, especially in tropical regions. The IUCN listed more than 300 terrestrial species from Pacific Island nations as threatened by climate change (Taylor and Kumar 2016). Monastersky (2014) indicated that climate change accounts for approximately 7% of declines in animal populations and is expected to become more severe over time. Urban (2015) predicted that climate change presents an extinction risk to 9% of species in Asia and suggested amphibians (13%) and reptiles (9%) are the most vulnerable vertebrate groups. Narrowly-distributed species and habitat specialists with unique ecological niches are considered to be particularly vulnerable since they are less capable of responding to impacts of climate change compared to wide-ranging species (Huggett 2004, Sandel et al. 2011, Pimm et al. 2014, Markle and Kozak 2018. As ectotherms, reptiles are generally susceptible to environmental changes. Some remote forest-dwelling species show a limited dispersal capability, and as such climate change is expected to affect those severely (Araújo et al. 2006, López-Alcaide and Macip-Ríos 2011, Fitzgerald et al. 2018, Powers and Jetz 2019, Vicente Liz et al. 2019). However, actual impacts of global climate change on reptiles have not been well understood because the distribution of most reptiles fall within remote and restricted areas, as shown for many species in Vietnam (Guisan and Hofer 2003, Segurado and Araújo 2004, Sterling et al. 2006, Urban 2015, Ngo et al. 2016, van Schingen et al. 2016a, Le et al. 2017).
Lichtenfelder`s Tiger Gecko, Goniurosaurus lichtenfelderi (MOCQUARD, 1897), is one of 24 species of the eublepharid genus Goniurosaurus, characterized by a high level of local endemism and narrow distribution (Orlov et al. 2008, Liang et al. 2018, Ngo et al. 2019a, Qi et al. 2020a, b, Zhu et al. 2020a, b, Uetz et al. 2021). This habitat specialist was first described from an island located in the Gulf of Tonkin in northern Vietnam (Orlov et al. 2008). Subsequently, the tiger gecko has been recorded from several provinces across the mainland of northeastern Vietnam and was very recently discovered in Chongzuo, southern China (Grismer 2000, Orlov et al. 2008, Nguyen et al. 2009, Gawor et al. 2016, Zhu et al. 2020b). Due to the restricted distribution, G. lichtenfelderi is potentially vulnerable to anthropogenic threats such as habitat loss and climate change, resembling other tiger geckos (Yang and Chan 2015, Ngo et al. 2016, Ngo et al. 2019b. Consequently, G. lichtenfelderi is assessed as Vulnerable (VU) in the IUCN Red List for Threatened Species (Nguyen 2018). In addition, as locally over-harvested for the international pet trade, the species has been included in CITES (the Convention on International Trade in Endangered Species of Wild Fauna and Flora) Appendix II and the Vietnam Government's Decree No. 06/2019/ND-CP (Group IIB) in 2019 (Ngo et al. 2019b). However, appropriate conservation measures have not yet been undertaken to safeguard in situ populations of this highly threatened species.
In this study, we predicted the potential distribution of G. lichtenfelderi using distribution records and environmental variables (including climate and vegetation cover data). The maximum entropy algorithm (MAXENT) was employed to predict the species' current potential distribution and evaluate alterations of suitable habitat for the species under different scenarios of climate change. We hypothesize that the contemporary suitable area of G. lichtenfelderi will contract under any scenario of climate change. Furthermore, we propose potential core refugia of G. lichtenfelderi based on projected suitable habitats, to improve the efficacy of conservation measures.

Study area
The study area (within 20-24 0 N and 104-109 0 E) with elevations ranging from 1 m to 2,139 m a.s.l. (Fig. 1), encompassing the entire known distributions of four captured Goniurosaurus species in northern Vietnam, was selected based on previous studies, direct observations, and interviews with local community members (Grismer et al. 1999, Vu et al. 2006, Orlov et al. 2008, Ziegler et al. 2008, Orlov et al. 2020, Zhu et al. 2020b. To obtain occurrence localities of the target species, we conducted field surveys in Hai Duong, Bac Giang and Quang Ninh provinces (both in the mainland and on offshore islands), northern Vietnam. Two new distribution records of G. lichtenfelderi from Ba Vi National Park (NP), Ha Noi Capital, northern Vietnam (Orlov et al. 2020) and Chongzuo, China (Zhu et al. 2020b) were also included ( Fig. 1). Based on previous descriptions of its natural habitat, G. lichtenfelderi is only found along small streams situated within evergreen forests on granitic formations (Grismer 2000, Orlov et al. 2008. Macro-environmental characteristics of the study region are presented in detail in Supplementary  Table S1.

Data collection
Occurrence data of G. lichtenfelderi were compiled from literature and direct observations in the field from 2014 to 2020. Only data with precise location information were included in the analysis apart from one pseudo-presence from Chongzuo, China, for which exact information was not available (Zhu et al. 2020b). Using pseudo-data likely affects the accuracy of SDMs by potentially introducing a sampling bias. Rather than being randomly selected, the pseudo-presence point was intentionally placed in a forest with dense vegetation in Chongzuo spotted based on satellite imagery, in proximity to the recorded distribution of Vietnamese populations in the border area. This pseudo-point is highly expected to be within the distribution of − or at least in close proximity to − the natural population of G. lichtenfelderi in China.
Coordinates of captured individuals were recorded with a GPS Garmin 64. The data will be shared upon reasonable request. A total of 190 occurrence records in the WGS84 projection of G. lichtenfelderi were initially selected. However, several records were removed by using the spatial filtering in packages "dismo" and "sp" in R v 3.1.2 (R Core Team 2018). Only one occurrence locality of G. lichtenfelderi was randomly selected within each 1 km square. Such spatial filtering can improve the quality of prediction models by decreasing geographical bias, autocorrelation effects, and uncertainty (Veloz 2009, Radosavljevic and Anderson 2014, Kiedrzyński et al. 2017. For macro-climatic data, we initially extracted 19 bioclimatic variables in the study area. The current (version 1.4) and future climatic conditions from the Coupled Model Intercomparison Project Phase 5 (CMIP5) were obtained from Worldclim (https://www.worldclim.org/) (Hijmans et al. 2005). To predict the future potential distribution of the target species, we selected future climatic data from two general circulation models (GCM) including Community Climate System Model version 4 (CCSM4) (Gent et al. 2011) and Beijing Climate Center -Climate System Model 1-1 (BCC-CSM 1-1) (Wu et al. 2014) by 2050s (average 2041 -2060) and by 2070s (average 2061 -2080). Two climate change scenarios of representative concentration pathways (RCPs): RCP 4.5 and RCP 8.5, representing intermediate and the most severe levels of accumulation of greenhouse gas concentrations on the future climate, were used for each model, respectively (Moss et al. 2010, Van Vuuren et al. 2011. For vegetation cover data, we extracted the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) within the study area using Moderate Resolution Imaging Spectroradiometer (MODIS) sensor data in recent five years from 2015 to 2019 (https://earthexplorer. usgs.gov/). The mean, maximum, minimum, median, range, standard deviations (STD), and mean, maximum, minimum of warmest and coldest quarters for the enhanced vegetation index (EVI) and normalized difference vegetation index (NDVI) were generated using Quantum GIS software (QGIS Version 3.12.0, Development Team. 2020. Available at http://qgis. osgeo.org [Downloaded 25 March 2020]). These environmental variables are listed in Supplementary  Table S1.
To remove the influence of multicollinearity on the modeling process, a pairwise Pearson correlation analysis of 36 environmental variables was conducted using the package "raster". Highly correlated variables with Pearson's correlation coefficients (r) larger than │0.7│ were removed. Consequently, 28 variables were excluded and only 11 remaining variables were included in the final analyses (Supplementary Table  S2). Six out of these selected variables were replaced in accordance with each future climate model, while five remaining variables of vegetation cover were used for the vegetation model. The climatic and vegetation cover variables were used separately to model the potential distributions in order to account for both broad-scale climatic factors delimiting the species potential distribution and fine-scale habitat availability within this potential distribution. This procedure allowed us to disentangle the potential impact of climate change from those of land-use changes, which are largely unknown in the near future.

Species distribution modeling
Occurrence data of G. lichtenfelderi and six selected climate variables under current conditions and future scenarios, and five remaining ones under current vegetation cover conditions were separately imported into the Maxent software (Phillips et al. 2006) v. 3.4.1 (https://biodiversityinformatics.amnh. org/opensource/maxent/ [accessed March 2020]) to predict suitable habitats. For a rare species with a restricted distribution, the selection of a study background area that includes very environmentally distant areas from the presences of G. lichtenfelderi can inflate the AUC scores (Lobo et al. 2008) (see more below). Thus, we generated a buffer background region for model building by circling an area with a radius of 75 km around occurrence points with dissolved boundaries using the R packages "dismo" and "raster". The selected buffer layer is expected to encompass all potential refugia of G. lichtenfelderi within the study area ( Fig. 1) and improves the accuracy of discriminate use of AUC in predicted models. We did project models to areas beyond the selected buffer site within the study area to assess alternative potential refugia in the context of future climate change.
As only a small number of species presence locations were available, using the default automatic configuration may not be always appropriate and may lead to model overfitting (Shcheglovitova andAnderson 2013, Muscarella et al. 2014). In order to minimize model overfitting, we used the ENMeval package in R and employed the jackknife method to determine the optimal model parameter configuration by varying the "feature" and "regularization multiplier" parameters for both the climatic and vegetation models (e.g. Supplementary Fig. S1) (Warren and Seifert 2011, Muscarella et al. 2014, Radosavljevic and Anderson 2014. The optimal model parameter configuration was identified based on evaluation criteria including the AUC test statistic, omission rate at two different thresholds (OR-MTP, OR-10) and delta.AIC c (Warren and Seifert 2011, Shcheglovitova and Anderson 2013, Muscarella et al. 2014, Radosavljevic and Anderson 2014, Morales et al. 2017). In particular, AUC is an effective threshold-independent index to distinguish presence from absence or background (0.5 indicates that the performance of the model is no better than random, while values closer to 1.0 indicate better model performance) (Fielding and Bell 1997, Phillips et al. 2006, Peterson et al. 2011). In addition, the OR-MTP (minimum training presence omission) and OR-10 (10% training omission) rates are thresholddependent metrics to evaluate model overfitting when their values are greater than the expectation of zero and 10%, respectively (Fielding and Bell 1997, Peterson 2011, Radosavljevic and Anderson 2014. The AIC (Akaike Information Criterion) corrected for small sample sizes reflects both model goodness-of-fit and complexity, and the model with the lowest AIC c value (i.e. delta.AIC c = 0) is considered the best model (Warren and Seifert 2011).
The minimum training presence threshold was selected in Maxent to convert continuous predicted values to binary maps of suitable and unsuitable areas (Pearson et al. 2007, Phillips et al. 2006. To train the data, we ran 100 bootstrap replications by sampling with replacement from the presence points (Phillips et al. 2006). Default values were applied for the maximum number of interactions (500), the maximum number of background points (10,000), and the prevalence of the species (0.5) (Phillips and Dudík 2008). Lastly, we selected the cloglog output format, in which values ranging from 0 to 1 indicate the probability of suitable environmental conditions for the species. In particular, values closer to 1 suggest a greater probability of species occurrence. The contribution of each variable to the current climatic and vegetation models was determined by measuring the permutation importance (Phillips and Dudik 2008).
To evaluate the accuracy of the climatic and vegetation maps, the multivariate environmental similarity surfaces (MESS) analysis was developed in recent versions of Maxent (Elith et al. 2010). MESS analysis compares the environmental similarity of variables and identifies the similarity of any given pixel to a reference set of pixels of chosen predictor variables within the whole study site. It is used to determine areas with nonanalog (novel) environments where the model presents extrapolations compared to highly similar (interpolated) areas (Fitzpatrick and Hargrove 2009, Elith et al. 2010, 2011. In this study, similarity/ novelty was classified into four levels: "< -10" (high extrapolation), "-10 -0" (low extrapolation), "0 -10" (low interpolation) and "> 10" (high interpolation). We also performed the MESS analysis to create maps for all future climate scenarios and to evaluate the alteration of novel climatic conditions from the reference conditions under the current climatic model within the whole study site.
Maxent was first run using default parameter settings to find potential new populations of G. lichtenfelderi using six selected climate variables and only 25 representative occurrence points collected from 2014 to 2019. The model was computed by selecting "auto features" and the regularization multiplier of "1.00". Afterwards, a short survey was conducted in Quang Ninh Province in July 2020 by interviewing local community members to confirm new populations of G. lichtenfelderi based on results from the initial prediction model. New coordinates were combined with previous locations to improve the model. Because low omission rates are desirable when searching for new populations (Peterson 2006, Lobo et al. 2008, we also identified optimal configurations by varying the "feature" and "regularization multiplier" parameters in the ENMeval package to improve the predictability of 25-occurrence model.

Core distribution areas
The core refugia for G. lichtenfelderi were identified within highly suitable areas in terms of climate, which have values above 10% training presence cloglog threshold (high occurrence probability). This is a stricter criterion for converting to a binary map with smaller suitable habitats (Radosavljevic and Anderson 2014). We also identified buffer refugia with the least-suitable environmental conditions with values above the minimum presence cloglog threshold (medium occurrence probability) (Phillips et al. 2006, Radosavljevic and. They were all combined with suitable areas in the vegetation model. To identify priority areas and assess the effectiveness of the protected areas to safeguard the species, we collected all shapefiles of nature reserves and national parks within the occurrence region of G. lichtenfelderi from the website of https://www.protectedplanet. net/ (accessed June 2019) (Fig. 1). The projected environmental refugia for the target species were afterwards overlaid with the protected area layer.

Results
Goniurosaurus lichtenfelderi was recorded in Yen Tu mountain range from southern Quang Ninh Province, through Bac Giang to Hai Duong provinces, and three non-contiguous areas in Ha Noi City, two islands in Quang Ninh Province, northern Vietnam and Chongzou, southern China. The species was documented in four protected areas, namely Bai Tu Long National Park (NP) and Dong Son -Ky Thuong Nature Reserve (NR) in Quang Ninh Province, Tay Yen Tu NR in Bac Giang Province, and Ba Vi NP in Ha Noi City (Fig. 1).
Using 25 occurrence records (including the pseudopresence from China) and current climatic variables, the model using default parameter settings predicted a new potential distribution of G. lichtenfelderi in northeastern border areas between China and Vietnam, which was beyond its initially known range (Fig. S2A). Based on the prediction results, we conducted extensive interviews with local community members in the northern districts of Quang Ninh Province. A total of 11 new records of G. lichtenfelderi were documented; none of these records were located within a protected area. We produced two additional, optimized models of LHQ -1.5 and H -1.5 configurations using the 25 initial records, selected following evaluation metrics (Fig. S1A, Table S3, Supplementary data). Predicted maps of the optimal models highlighted suitable areas more apparently in northern Quang Ninh Province, compared to the prediction of the default model (Fig. S2).
Using a total of 36 representative occurrence records as an input, we finally selected the optimal current climatic model of H -1.5 and the vegetation model of LQ -1.0 following our evaluation criteria (Fig. S1B, C). Specifically, both Delta.AIC c values of the models were zero; average.test.AUCs reached high values (0.94 and 0.87, respectively); and the values of OR-MTP and OR-10 were lower than 0.04 and 0.11, respectively (Table S3). According to the results of MESS analysis, we measured high similarity (interpolated) habitat covering 10,042 km 2 (approximately 75%) in the current climate model and 7,513 km 2 (approximately 56%) in the vegetation model within the selected buffer, with very small areas of extrapolation (high novelty) in both of these models (Fig. S4A, B).
Four important climate variables (Mean diurnal temperature range (Bio2), Mean temperature of driest quarter (Bio 9), Precipitation of warmest quarter (Bio 18) and Precipitation seasonality (Bio 15)) accounted for 96% of permutation importance for the current climatic model (Fig. S3A). For the vegetation model, three variables of NDVI of STD, NDVI of Mean Coldest Quarter and NDVI of Minimum Warmest Quarter accounted for 85.9% of permutation importance (Fig. S3B).
Under the current climate conditions, the potential distribution of G. lichtenfelderi mainly covers the sites of occurrence and their surroundings. Moreover, we recorded potentially suitable habitat in coastal areas and other islands in Quang Ninh Province, northern Vietnam, and a few additional sites in the border area of Fangchenggang, southern China ( Fig. 2A, B). In particular, climatically suitable habitat consists of a total area of 16,133 km 2 , of which approximately 13,408 km 2 (83%) are within the buffered selection (Fig. 4). Under the vegetation conditions, suitable habitat only consists of a total of 7,629 km 2 (57%) within the climatically suitable area in the selected buffer (Fig. 2C, D).
With respect to future projections for the 2050s and 2070s under the scenarios of RCP 4.5 and RCP 8.5, each predicted significantly smaller suitable areas compared to that determined by the current climate model, and the range contraction reached the highest level under RCP 8.5 by the 2070s (Fig. 3, 4). In fact, the climatically suitable areas under future scenarios fluctuated considerably depending on RCPs and periods examined (Fig. 3, 4). The largest climatically suitable area under the CCSM4 -RCP 4.5 scenario by the 2050s covers 6,209 km 2 in the selected buffer (46% suitable area under the current climate conditions), while under the BCC-CSM 1-1 -RCP 8.5 scenario by the 2070s, the potential distribution of G. lichtenfelderi was not recorded under novel climatic conditions (0 km 2 in the selected buffer) (Fig. 3D, 4). Using a MESS analysis to compare the current climate model with the future ones, we found that high similarity (interpolated) habitats within the potential climatic distribution of G. lichtenfelderi were partly replaced by high novelty (extrapolated) habitats under the future scenarios (Fig.  S4A, S5). Especially, under the RCP 8.5 scenarios by the 2070s, novel climatic conditions were documented in most of the study area (Fig. S5D, H). In general, the potential distribution of G. lichtenfelderi under future scenarios tends to shift towards higher latitudes in border areas between China and Vietnam (Figs 2, 3).
We also consider suitable areas as potential refugia for the long-term persistence of G. lichtenfelderi in Frontiers of Biogeography 2022, 14.1, e51167 © the authors, CC-BY 4.0 license 6 the context of current environmental conditions and impacts of future climate change. In particular, the green shapes, which represent highly suitable habitats, mostly cover areas in Bai Tu Long NP, Yen Tu Mountain Range, and border areas in both China and Vietnam, and the yellow areas with lower levels of suitability cover surrounding sites (Fig. 5). A few scattered green patches are also recorded on some islands and coastal areas from Quang Ninh Province (Fig. 5).

Modeling
To predict the potential distribution of G. lichtenfelderi, we used Maxent models with known and newly collected occurrence data and environmental variables (climate and vegetation). The current models of climate and vegetation performed very well based on model performance and evaluation metrics. Further, model tuning was performed to ensure model selection balance low AIC c and omission rates (Warren and Seifert 2011, Shcheglovitova and Anderson 2013, Muscarella et al. 2014, Radosavljevic and Anderson 2014, Morales et al. 2017. Following the MESS analysis, the majority of the potential distribution of G. lichtenfelderi is in areas of high similarity (interpolation), while areas of high novelty (extrapolation) are beyond the potential distribution of the current models (e.g. Supplementary Fig. S4A, B). These results confirmed that the optimal models performed well, are not overfit, and are meaningful to use for future projections.
Eleven focal environmental factors (climate and vegetation) were selected to predict the potential distribution of G. lichtenfelderi. Being ectothermic, the selected climate variables may influence the tiger geckos' spatial distribution, physiology, behavior and reproductive biology (López-Alcaide and Macip-Ríos 2011, Ngo et al. 2019a, Vicente et al. 2019. Furthermore, all G. lichtenfelderi specimens were recorded along small streams and none was found outside evergreen forests (Ngo et al. pers. obs). Thus, a high level of humidity and dense vegetation coverage might associate with the natural history of G. lichtenfelderi. Ngo et al. (2019a) already confirmed the importance of these environmental factors for the microhabitat requirement of another congener, G. catbaensis. There is a significant gap in the knowledge about habitat requirements of G. lichtenfelderi and more studies yet need to be undertaken to understand its role in the ecosystem. Thus, the specific biological relationships and mechanisms between the selected environmental variables (macro-climate and vegetation cover) and the species remain unclear.

Finding new populations
Beyond the known distribution of G. lichtenfelderi, the initial and two optimal models predicted the potential habitat in a few small patches of northern  We recorded sympatric occurrences of these two lizards by interviews with local community members in at least five locations. The discovery of new occurrence locations of G. lichtenfelderi from the predicted maps could provide further justification for the use of the geographically close pseudo-presence points, as we included in this study from Chongzuo, China. In our case, with so few presence points to work from, selecting a carefully vetted pseudo-presence point did not detract from model performance and likely significantly improved model results to predict the population of G. lichtenfelderi in China. The updated model using 36 occurrence points predicted other potential distribution sites of G. lichtenfelderi in Fangchenggang, southern China, a few islands and coastal areas in Quang Ninh Province, and a few small areas in Bac Ninh and Vinh Phuc provinces, northeastern Vietnam ( Fig. 2A, B). Over the last two decades, intensive surveys have never documented the species in Tam Dao National Park, Vinh Phuc Province (Nguyen et al. 2009). However, additional overlooked populations of G. lichtenfelderi will likely be discovered from one or more of the remaining unstudied potential areas.

Potential refugia
To date, only a few studies used the Maxent approach to predict the potential distribution of lizard species in Vietnam (e.g. Goniurosaurus catbaensis and Shinisaurus crocodilurus) under contemporary and future conditions (van Schingen et al. 2016a, b, Le et al. 2017, Ngo et al. 2019a). The potential distribution of these species was predicted to be diminished due to the impact of climate change. The projected distribution range of G. lichtenfelderi is consistent with the trend of contraction under climate change. Novel climatic conditions will arise within the current potential distribution of G. lichtenfelderi, even the currently realized niche of G. lichtenfelderi will be totally replaced following the BCC-CSM 1-1 -RCP 8.5 model projection by the 2070s (Figs 3, 4). Overall, the potential distribution of G. lichtenfelderi tends to shift towards higher latitudes, following a similar  trend in many terrestrial species (Walther et al. 2002, Deutsch et al. 2008, Lenoir and Svenning 2013. Under contemporary climatic conditions, the potential distribution of G. lichtenfelderi mostly covers occurrence areas and their surroundings ( Fig. 2A). The current map shows that the potential range of G. lichtenfelderi does not overlap with the occurrence extent of each remaining Goniurosaurus congeners in Vietnam, except for G. catbaensis from the Gulf of Tonkin, northern Vietnam (Fig. 1, Fig. 2A). However, the potentially sympatric patches only occur in a few small karst islands in Ha Long Bay, Quang Ninh Province and Cat Ba NP, Hai Phong City, where the target species has never been recorded, as it has only been found in granite-dwelling forests (Ziegler et al. 2008, Gawor et al. 2016, Ngo et al. 2019a. It is likely that G. lichtenfelderi has a distinct climate niche from other congeners of the G. luii group in Vietnam. A macroclimatic niche divergence may explain the evolutionary split between the granite-adapted G. lichtenfelderi species and the karst-dwelling G. luii group (Orlov et al. 2008, Ziegler et al. 2008), but these hypotheses should be further investigated in the future.
Since all Goniurosaurus species are likely narrowlydistributed habitat specialists and adapted to unique ecological niches, we only identified stable refugia for G. lichtenfelderi under the impacts of climate change, rather than identifying colonization areas (Orlov et al. 2008, Ngo et al. 2019a. A potential refugium, which is highly suitable in the context of current environmental conditions and can be relatively stable under climate change scenarios, can potentially increase landscape connectivity to facilitate genetic exchange between subpopulations (Littlefield et al. 2017). Thus, we suggest that key regions safeguard wild populations of G. lichtenfelderi include the contiguous Yen Tu Mountain Range (Refugium-1), Bai Tu Long NP (Refugium-2), northern parts of Quang Ninh Province,, and border areas in Fangchenggang and Chongzuo, China (Refugium-4) (Fig. 5). Within these refugia, three protected areas (Bai Tu Long NP, Tay Yen Tu NR and Dong Son -Ky Thuong NR) contain a significant proportion of the area of habitat suitability (approximately 365 km 2 ).
Thus, we highly recommend that these three protected areas and the green regions in Figure 5 be prioritized as core areas. Their surrounding regions bounded with yellow patches can function as buffer zones; for example, Hai Duong Province and the area between Dong Son -Ky Thuong NR and Tay Yen Tu NR can serve as corridors connecting the core localities (Fig. 5). Van Schingen et al. (2016a) also highlighted the importance of key forest areas within Yen Tu Mountain Range and the border region in Quang Ninh Province to maintain genetic exchange and the long-term persistence of another threatened lizard, Shinisaurus crocodilurus. We assume that other sympatric species will also be affected by the contraction of suitable habitats due to the impact of climate change and they should be safeguarded in the same manner.
We highlight the importance of the two borderregion refugia in Quang Ninh Province, northern Vietnam (Refugium-3), and Fangchenggang and Chongzuo, southernmost China (Refugium-4), in relation to in situ conservation actions. In particular, under the potential threat of climate change, the currently suitable refugia of G. lichtenfelderi are likely to shrink and move towards higher latitudes in the border areas between China and Vietnam within the next decades. Moreover, the two refugia shape a contiguously suitable and unfragmented habitat for the species between the two countries (Fig. 5). Therefore, research and management collaborations between authorities of China and Vietnam are essential to protect the remaining wild populations of G. lichtenfelderi.
With respect to the primary region in Bai Tu Long NP, surrounding islands and coastal areas in Quang Ninh Province, increasing sea levels and catastrophic events caused by the global climate change may be a potential threat to the island-populations of G. lichtenfelderi. As Ngo et al. (2019a) documented in a case study of G. catbaensis, a sub-population seemed to be wiped out from a former occurrence site on Cat Ba Island, Hai Phong City after an extreme flood event. In fact, the target species only occurs in restricted areas along small streams on granitic rocks, thus, microhabitat degradation is expected to be exacerbated in the summer by heavy rains (Ngo et al. pers. obs).

Anthropogenic impacts and conservation
Previous observations have already shown the degradation and fragmentation of habitats in the Yen Tu Mountain Range caused by direct anthropogenic threats, namely deforestation as a consequence of coal mining, timber logging and tourism development, and waste deposition (van Schingen et al. 2014, 2015. Habitat fragmentation likely generates potential barriers that prevent genetic exchange through dispersal between G. lichtenfelderi populations. These anthropogenic impacts will make the species even more vulnerable to the impacts of climate change. In fact, no G. lichtenfelderi individual was observed along stream sections in close proximity to the residential and agricultural regions (Ngo et al. pers. obs). We assume that the actual habitat loss could be worse than the predicted loss under climate change scenarios.
In addition, all narrowly distributed tiger geckos are extremely vulnerable to over-exploitation (Yang andChan 2015, Ngo et al. 2019b). A large number of wild animals of G. lichtenfelderi were over-harvested to meet high demand from the local and international pet markets. Particularly, G. lichtenfelderi was the most common tiger gecko species in the international pet trade from 1999 to 2018 comprising 43.6% (7,281 individuals) of all imported individuals into the United States (US) (Ngo et al. 2019b).
Synergistic effects of habitat degradation and overexploitation will certainly lead to an exacerbation of climate change impacts in the future. The species may even go extinct before climatic effects become more apparent. Thus, mapping optimal refugia and identifying key regions for G. lichtenfelderi will help to Frontiers of Biogeography 2022, 14.1, e51167 © the authors, CC-BY 4.0 license 10 enhance the effectiveness of stringent conservation strategies to mitigate the impacts of climate change and other human disturbances. Furthermore, we propose to establish a species and habitat conservation area for the highly threatened species G. lichtenfelderi and Shinisaurus crocodilurus in unprotected areas in the border area between China and Vietnam.

Conclusion
In this study, we predicted the potential distribution of G. lichtenfelderi by using the maximum entropy algorithm in the Maxent software under current and future environmental scenarios across its range in northern Vietnam and southern China. This resulted in successful field surveys documenting new occurrence records for the species. Contemporary potential distributions were mainly recorded in the occurrence areas of G. lichtenfelderi and their surroundings. Suitable refugia for G. lichtenfelderi are expected to decline significantly and shift towards higher latitudes under climate change scenarios by the 2050s and 2070s.
Border areas between China and Vietnam, Yen Tu Mountain Range and granitic habitats on islands in Quang Ninh Province, including three protected areas, are considered core refugia for the implementation of conservation measures. Unprotected areas in the border region between China and Vietnam covering a high proportion of suitable habitat for both endangered species of G. lichtenfelderi and S. crocodilurus are proposed to be set aside as a species and habitat conservation area in the future.
The results of the present study are expected to form a baseline for formulating and implementing appropriate in situ conservation measures to safeguard wild populations of G. lichtenfelderi and facilitate research management collaborations between authorities from China and Vietnam to mitigate impacts of climate change and other human disturbances.

Author Contributions
HNN, HQN, TQP and LRG conducted field surveys. HNN and DR computed the species distribution models. HNN, TQN, LRG, DR and TZ wrote the manuscript.

Supplementary Material
The following materials are available as part of the online article at https://escholarship.org/uc/fb Figure S1. Four evaluation criteria resulting from three models made across a range of feature-class combinations and regularization multipliers. Figure S2. Predicted suitable habitat for Goniurosaurus lichtenfelderi under current climate conditions using 25 occurrence points with a set of different configurations. Figure S3. Permutation importance in the current prediction models. Figure S4. Multi-Environment Similarity Surface (MESS) maps of novel habitat for current models. Figure S5. Multi-Environment Similarity Surface (MESS) maps of novel habitat for the circulation models of climate change. Table S1. Current climate and vegetation cover conditions within the study site. Table S2. Multicollinearity test results showing Pearson correlation coefficients of eleven selected environmental variables used in species distribution modeling. Table S3. Evaluation metrics of Maxent ENMs generated by ENMeval for the optimal climatic and vegetation models.