Extraordinary diversity among allopatric species in the genus Goniurosaurus (Squamata: Eublepharidae): understanding niche evolution and the need of conservation measures

Given the high degrees of adaptation to specific microhabitats and restricted-range endemism, Goniurosaurus (Tiger geckos) serves as a unique model to study the complex evolution in lizards. Using phylogenetic analyses, we estimated the first divergence date of Goniurosaurus to the Eocene (~ 45.3 mya). The diversification within four monophyletic species groups began in the mid-Miocene between ~ 13.4 and 7.7 mya and continued to at least the early Pleistocene (~ 2 mya). Their ancestor was predicted to originate somewhere in contiguous continental Eastern Asia, whereas the current regions in which each monophyletic Goniurosaurus species group radiated are respectively their own ancestral regions. Together with factors of altitudinal gradient and climate conditions, we reconstructed relevant niche models of Goniurosaurus including ancestral reconstructions. Consequently, low elevations were predicted to be the most probable ancestral state for Goniurosaurus and all its groups as well. Both climatic niche conservatism and divergence have shaped the extraordinary species richness of allopatric Chinese and Vietnamese tiger geckos. In terms of endangerment, Goniurosaurus has been considered one of the most susceptible lizard groups under severe human impacts, especially climate change. The assessments of their niche evolution can provide a science-based pre-signal of vulnerability, thereby improving the efficacy of conservation measures to safeguard species of Goniurosaurus in the future. Accordingly, almost all closely related species of Goniurosaurus in China and Vietnam were identified with a high rate of niche conservatism, which should be included in conservation priorities under potential impacts of climate change.


Introduction
It has been acknowledged that ecological conditions are the major factors that determine the distribution of species as well as being the main driver of speciation (Darwin and Wallace 1858;MacArthur 1984;Orr and Smith 1998;Schluter 2009;Glor and Warren 2011). In particular, the largest geographic range in which a species can sustain a stable population size, is frequently limited by its physiological tolerances to ecological factors defined by its fundamental niche (Grinnell 1917;Kearney and Porter 2004). Together with biotic interactions, dispersal barriers and chance events, species actually occupy only fractions of their fundamental niche, termed realized niche that can explain the high diversity of complex range-restricted species, especially in reptile groups distributed in the tropics (Hutchinson 1957;Sterling et al. 2006;Heaney 2007;Mittelbach et al. 2007;Soberón 2007;Böhm et al. 2013;Stein et al. 2014;Steinbauer et al. 2016). In endangerment concerns, such reptiles are at high risk of extinction, following Criterion B and D2 in the IUCN Red List relative to restricted geographic range (Böhm et al. 2013;IUCN 2021). Against the background of increasing anthropogenic impacts, their survival fate is severely imperiled (Thuiller et al. 2005;Ohlemüller et al. 2008;Böhm et al. 2013;Cahill et al. 2013;Trew and Maclean 2021;Cox et al. 2022). Consequently, ecological alternations (e.g., climate change) have been globally acknowledged to be a looming factor of extinction, rather than being a driver of speciation in the past (Thuiller et al. 2005;Schluter 2009;Lavergne et al. 2013;Bellard et al. 2014;Smith et al. 2019).
Given the potential for niche evolution, the level of adaptability in a species to different environmental conditions can be estimated to assess its chances of survival. Accordingly, natural populations could be subjected to ecology-based divergent selection across the landscape under different abiotic factors (e.g. climate) as non-physical barriers favoring reproductive isolation and consequent speciation in new habitats by accumulating genetic difference (niche divergence). Consequently, the evolved lineages were detected with distinct evolutionary responses from their ancestor and a high level of physiological tolerances to niche alternations (Ahmadzadeh et al. 2013;Krehnwinkel et al. 2015;Rato et al. 2015;Ahmadi et al. 2021). On the other hand, closely related taxa tend to occupy more similar niches than distant phylogenetic congeners and retain ancestral niches over time (niche conservatism), when their speciation was mainly related to the allopatric distribution patterns by variance. Such species may be more susceptible to climate change due to constraints in genetic variation caused by previous niche conservatism limiting the likelihood of adaptation to changing conditions. Therefore, the looming threat containing novel conditions may exceed their narrower physiological tolerance, consequently leading to extirpation (Wiens and Graham 2005;Hadly et al. 2009;Wiens et al. 2010;Peterson 2011;Lavergne et al. 2013;Pyron et al. 2014;Ahmadi et al. 2021). Understanding the potential for niche evolution is particularly important in evaluating the status of endangered species and thereby proposing science-based measures of conservation (Richardson and Whittaker 2010;Wiens et al. 2010;Moritz and Potter 2013). Recently, species distribution models (SDM) have been employed to identify patterns and mechanisms of niche evolution, enabling researchers to test causal relationships between ecological traits and diversification processes (Evans et al. 2009;Hof et al. 2010;Salamin et al. 2010;Broennimann et al. 2012;Lavergne et al. 2013;Rato et al. 2015;Ahmadi et al. 2021).
Due to the high degrees of adaptation to specific microhabitats and restricted-range endemism, tiger geckos belonging to the genus Goniurosaurus Barbour, 1908 can serve as a unique model to study the complexity of evolution in lizards (Honda et al. 2014; 1 3 Liang et al. 2018;Ngo et al. 2021b;Grismer et al. 2021). Recently, systematic issues of Goniurosaurus were resolved and ongoing discoveries of new species have been updated based on integrative taxonomic approaches using morphology and phylogenetics (Liang et al. 2018;Grismer et al. 2021;Ngo et al. 2021b;Zhu et al. 2022). In particular, 25 tiger geckos have been assigned to one of four monophyletic species groups: the kuroiwae group with six exclusively insular species from the Ryukyu Archipelago, Japan; the lichtenfelderi group with five species and the luii group with nine species (whereof G. kadoorieorum is recently considered a synonym of G. luii) distributed throughout insular and mainland sites in China and Vietnam; and the yingdeensis group with five continental species in southern China ( Fig. 1) (Nguyen et al. 2009;Nguyen 2011;Liang et al. 2018;Qi et al. 2020;Zhu et al. 2020aZhu et al. , 2020bZhu et al. , 2021Grismer et al. 2021;Ngo et al. 2021bNgo et al. , 2022a. Although most of their distributions are adjacent to each other, none of them have been recorded to be in sympatry (Orlov et al. 2008;Yang and Chan 2015;Zhu et al. 2020aZhu et al. , b, 2021Ngo et al. 2021b). Oceanic archipelagos and disjunct karst ecosystems could serve as physical barriers and/or provide unique microclimate zones that might impede dispersal and hence gene flow among the Goniurosaurus populations (Clements et al. 2006;Sterling et al. 2006;Heaney 2007). Research presented by Liang et al. (2018) indeed underlined the central role of Hainan Island, which formulated species diversification within the lichtenfelderi group. Using an analysis of ancestral habitat reconstruction, Grismer et al. (2021) shed light on the evolution of habitat preference in Goniurosaurus. Accordingly, karstic ecosystems were highly supported as the most probable ancestral condition for Goniurosaurus as well as for three of the four species groups (e.g. kuroiwae, luii and yingdeensis), whereas the lichtenfelderi group most likely evolved in non-karstic habitats.
Goniurosaurus is considered one of the most threatened genera under ongoing severe human impacts (e.g. unsustainable exploitation and habitat loss) (Ngo et al. 2019). Furthermore, as highly endemic ectotherms with limited dispersal capacity, potential factors-including small population size, stochastic events and climate change can also push Goniurosaurus species to the brink of extinction (Ngo et al. 2019(Ngo et al. , 2021a(Ngo et al. , b, 2022a. To safeguard their populations, legal regulations listing all Tiger gecko species in CITES (Convention on International Trade in Endangered Species of Wild Fauna and Flora) appendices have recently come into force (Ngo et al. 2019;CITES Notification No. 2020/068), whereof 17 species have been considered threatened in the IUCN Red List (IUCN 2021). In spite of a wealth of studies on systematics, as well as assessments of demography and ecology (Liang et al. 2018;Ngo et al. 2019Ngo et al. , 2021aNgo et al. , b, 2022aGrismer et al. 2021), understanding of the potential of niche evolution in Goniurosaurus remains to be investigated in order to identify priorities in species conservation planning.
Following evolutionary hypotheses, cladogenesis in Goniurosaurus species from Eastern Asia was assumed to be strongly related to past orogenic processes that facilitated the geographic separation (e.g., disjunct limestone habitats, islands, rivers) which consequently led to reproductive isolation (Honda et al. 2014;Liang et al. 2018;Ngo et al. 2021bNgo et al. , 2022b. However, the simple vicariance might play only an indirect role and cannot fully explain the mechanism of speciation within Goniurosaurus. Combining a dated phylogenetic tree of all Goniurosaurus species together with their distribution data, altitudinal gradient and sets of climate conditions, we herein test evolutionary models that may explain their extraordinary richness in allopatry and high levels of endemism. In terms of species conservation planning, we also intend to identify closely related species of Goniurosaurus preserving climatic conditions over time (niche conservatism), that allows us to include them in priorities of conservation protection due to the particular vulnerability to climate change in the future.
We constructed Maximum Likelihood (ML), Bayesian Inference (BI), and Bayesian Evolutionary Analysis by Sampling Trees (BEAST) phylogenetic trees using a concatenated data set composed of 3070 base pairs (bp) of the mitochondrial genes, 16 s (633 bp) and cytb (1075 bp), and the nuclear genes, CMOS (472 bp) and RAG1 (890 bp), from 103 specimens of 24 species of Goniurosaurus with varying degrees of sequence coverage across the samples. Concatenation followed the comparison of separate gene trees to confirm there were no major discordances. One species, Eublepharis macularius, served as the outgroup (Grismer 1988;Jonniaux and Kumazawa 2008) to root all trees. See Grismer et al. (2021) for all sequence data and GenBank accession numbers.
A Maximum likelihood (ML) analysis partitioned by gene was implemented using the IQ-TREE webserver (Nguyen et al. 2015;Trifinopoulos et al. 2016) preceded by the selection of substitution models using TIM2 + F + I + G4 for 16 s and cytb and HKY + F for CMOS and RAG1. To avoid over parameterization, protein coding genes were not partitioned by codon. One-thousand bootstrap pseudo-replicates via the ultrafast bootstrap (UFB: Hoang et al. 2018) approximation algorithm were employed, and nodes having UFB values of 95 and above were considered strongly supported (Minh et al. 2013). We considered nodes with values of 90-94 as well-supported. A Bayesian inference (BI) analysis was carried out in MrBayes 3.2.3. (Ronquist et al. 2012) on XSEDE using the CIPRES Science Gateway (Cyberinfrastructure for Phylogenetic Research; Miller et al. 2010) employing default priors and models of evolution that most closely approximated those selected by the BIC and used in the ML analysis. Two independent Markov chain Monte Carlo (MCMC) analyses for each data set were performed-each with four chains, three hot and one cold. The MCMC simulations ran for 100 million generations. Trees were sampled every 10,000 generations, and the first 10% of the trees from each run from each data set were discarded as burn-in. The parameter files from both runs were checked in Tracer v1.6 (Rambaut et al. 2014) to ensure that convergence and stationarity of all parameters had effective sample sizes (ESS) well-above above 200. Postburn-in sampled trees from each respective run were combined and 50% majority-rule consensus trees were constructed. Nodes with Bayesian posterior probabilities (BPP) of 0.95 and above were considered highly supported. We considered nodes with values of 0.90-0.94 as well-supported.
An input file was constructed in Bayesian Evolutionary Analysis Utility (BEAUti) v. 2.4.6 using a relaxed lognormal clock with unlinked site models, linked trees and clock models, and a Yule prior and run in BEAST on CIPRES (Cyberinfrastructure for Phylogenetic Research; Miller et al. 2010). bModelTest was used to numerically integrate over the uncertainty of substitution models of each gene while simultaneously estimating phylogeny using MCMC analyses. MCMC chains were run for 100,000,000 generations and logged every 10,000 generations. The BEAST log file was visualized in Tracer v. 1.6.0 (Rambaut et al. 2014) to ensure effective sample sizes (ESS) were well-above 200 for all parameters. A Maximum clade credibility tree using mean heights at the nodes was generated using TreeAnnotator v.1.8.0 (Rambaut and Drummond 2014) with a burn-in of 1000 trees (10%). Nodes with BPPs of 0.95 and above were considered strongly supported. We considered nodes with values of 0.90-0.94 as well-supported.

Ancestral range estimation
The resulting BEAST tree was used to estimate the ancestral range at each node using the R package BioGeoBEARS (Matzke 2014). We followed the dating scheme of Jonniaux and Kumazawa (2008) in order to compare it with the BEAST and treePL (Smith and O'Meara 2012) analyses of divergence times of Liang et al. (2018). BioGeoBE-ARS allows for both probabilistic inferences of ancestral geographic ranges and statistical comparisons of range expansion from different models in a likelihood framework employing the Akaike Information Criterion (AIC) to allow the data to choose the best fitting model. Available models in BioGeoBEARS include a likelihood version of the parsimony-based Dispersal Vicariance Analysis DIVA ("DIVALIKE") (Ronquist 1997), the likelihood-based Dispersal-Extinction Cladogenesis (DEC) model of the LAGRANGE program (Ree and Smith 2008), and the Bayesian-based BayArea ("BAYAREALIKE") (Landis et al. 2013). Additionally, each model incorporates founder-effect speciation (+ J) which is purported to be particularly important when reconstructing biogeographic scenarios of insular lineages (Matzke 2014). For a geography input file, we classified Goniurosaurus species into distinct areas of the four monophyletic groups' current geography following the division of Liang et al. (2018), including (1) Guangdong Province, China (C) with four species of yingdeensis group (Qi et al. 2020); (2) Hainan Island, China (H) with four species of lichtenfelderi group (except for G. lichtenfelderi); (3) Ryukyu Archipelago, Japan (R) with six species of kuroiwae group; and (4) Vietnam and China (V) from Red River fault in northern Vietnam, northward to southern Guizhou Province, China, with nine species of luii group and G. lichtenfelderi. Each species occurs in only a single region and as such was allowed to assign only a single selected area in the analysis (Liang et al. 2018;Grismer et al. 2021;Ngo et al. 2021b).

Occurrence records
A total of 223 occurrence records of 16 Goniurosaurus species in China and Vietnam were collected from our field surveys during the two last decades by ourselves and others ( Fig. 1). In order to improve the quality of predictions based on our models, which may be affected by geographic bias as results of duplicates representing pseudo replicates, we randomly selected only one record within each 1 km square by using spatial filtering functions in R v 3.1.2 (RStudio Team 2018; Aiello-Lammens et al. 2015). Consequently, 160 representative records of three Goniurosaurus species groups (including lichtenfelderi-44 records, luii-79 records and yingdeensis-38 records, see details in Supplementary information Table S1) were used for further analyses.

Ancestral elevational range reconstruction
Based on significant changes of vegetation as Bain and Hurley (2011) suggested for the study area, the elevational ranges were divided into three different levels, namely low (below 300 m), medium (310-800 m) and high (above 800 m). The coding of elevation level for each Goniurosaurus species, which was determined from the literature and field observations of the authors, was extracted from coordinate data of occurrences using GPS equipment ( Supplementary Information Fig. S1). The classification 1 3 was mapped onto the BEAST tree of Goniurosaurus using a stochastic character mapping (SCM) implemented in the R package phytools (Revell 2012) in order to derive probability estimates of the ancestral states at each node. A transition rate matrix was identified that best fit the data by comparing the corrected Akaike Information Criterion (AICc) in the R package ape (Paradis and Schliep 2018). Three transition rate models were considered: a 2-parameter model having different rates for every transition type (the ARD model); a single-parameter model with equal forward and reverse rates between states (the symmetrical rates SYM model); and a single rate parameter model that assumes equal rates among all transitions (ER). Lastly, an MCMC approach was used to sample the most probable 1000 trait histories from the posterior using the make. simmap() command and then summarized them using the summary() command.

Climatic niche comparisons
In terms of environmental data, 19 bioclimatic variables for current conditions with a spatial resolution of 30 arc-second were obtained from Worldclim (http:// www. world clim. org/, version 2.0) (Fick and Hijmans 2017). To mitigate the influence of potential multicollinearity of predictors on the models, we computed a pairwise Pearson correlations to identify and eliminate highly correlated climatic variables with r coefficient values larger than │0.7│ (Dormann et al. 2013). Our final data set comprised seven climatic variables including Mean Diurnal Range (BIO2), Isothermality (BIO3), Max Temperature of Warmest Month (BIO5), Mean Temperature of Wettest Quarter (BIO8), Precipitation of Wettest Month (BIO13), of Driest Month (BIO14) and Warmest Quarter (BIO18), which were selected for subsequent analyses.
To compare the macro-climatic niche among Goniurosaurus groups and sister taxa of each group, all pairwise tests were computed by using the package "ecospat" in R software (Cola et al. 2017). The niche equivalency and two-way similarity tests, presented in detail by Warren et al. (2008) and Broennimann et al. (2012), were calculated between pairs of Goniurosaurus groups. The degree of overlap in climatic-niche space was further evaluated using Schoener's D as recommended by Rödder and Engler (2011), which varies from 0 (no overlap) to 1 (complete overlap) (Schoener 1970;Warren et al. 2008). However, the niche overlap indices among species-pairs cannot be reliably estimated due to the limited availability of records with less than 10 , which is the case in most Tiger geckos. Following a jackknifing approach proposed by  to overcome the limitation, all available Goniurosaurus records of each monophyletic group were pooled into a dataset and then repeated "leave-one-species-out" to calculate Schoener's D values. We subsequently used the value of 1-D as the relative degree of overlap between the omitted one and all remaining species.
Each macro-climatic niche was described in an orthogonal space combining the first two principal components in the PCA-env approach (Broennimann et al. 2012). The overlap level of given macro-climatic niche spaces was evaluated by direct observations on the same environmental background. The proportion of each principal component (PC) explaining climatic variance was calculated by the proportion of its eigenvalue to the sum. To define the available macro-climatic space, we constructed a minimum convex polygon (MCP) covering a total of selected Goniurosaurus records using the package "adehabitatHR" in R software (Calenge 2006).

Species distribution modelling
For species for which only one or two occurrence records were available (namely G. araneus, G. gezhi, G. liboensis, G. kwangsiensis, G. kadoorieorum, G. gollum and G. zhoui), we estimated their potential distributions by applying a circular buffer with a radius of 10 km around occurrence points with dissolved boundaries in R. In terms of SDM, the potential distributions of remaining tiger geckos were predicted by using the ensemble of small models (ESMs) combining six modelling techniques: artificial neural network (ANN), classification tree analysis (CTA), generalized linear models (GLM), generalized additive models (GAM), generalized boosting regression models (GBM), and maximum entropy modelling (MAXENT.Phillips). All models were calibrated using the "biomod2" package (Breiner et al. 2015;Thuilleret et al. 2016;Cola et al. 2017). Based on the seven selected climatic variables, we computed a principal component analysis (PCA) in R to derive an orthogonal niche space. Only principal components (PCs) with eigenvalues > 1 were subsequently considered, which were subsequently employed in the SDM (Rato et al. 2015). Although the potential distribution of each species was projected within the entire background, all models were only trained with presence-only data and 10,000 random points selected within the MCP area. Subsets of 70% training data and 30% test data were used to model and evaluate the SDMs using the two performance indices AUC and TSS (Fielding and Bell 1997;Allouche et al. 2006). Values closer to 1.0 indicate better model performance (Breiner et al. 2015).

History of niche occupancy
Following an approach proposed by Evans et al. (2009), the evolutionary history of climatic niche occupancy of Goniurosaurus species in China and Vietnam was reconstructed throughout predicted niche occupancy (PNO) profiles in the "phyloclim" package. The suitable probability derived from ESMs of each Goniurosaurus species was re-scaled to the sum of 1, and corresponding PC scores binned into 1000 grid cells together with designed ones of occurrence-limited species in order to obtain the PNO profile for each PC. From these PNO profiles, we resampled 1000 times using the generalized least squares method with an assumption of Brownian motion evolution to estimate the climatic spectrum on the BEAST dated tree of Goniurosaurus species in China and Vietnam (Rato et al. 2015;Ahmadzadeh et al. 2016;Ahmadi et al. 2021).

Ancestral climatic state estimation
We tested for phylogenetic signals in climatic niche axes based on climatic PCs and our dated phylogenetic tree, using Blomberg's K (Blomberg et al. 2003) and Moran's I (Münkemüller et al. 2012) in the "phylosig" function of the package "phytools" in R. A phylogenetic signal was identified with a significance level of p < 0.05, indicating that Blomberg's K and Moran's I are significantly different from zero. In particular, values of K range from zero to infinity, where K < 1 indicates weak phylogenetic signal, and K > 1 strong signal, implying that sister taxa are, respectively, more or less divergent in environmental traits. In terms of Moran's I, the value reaches to + 1 indicating a perfect cluster of 1 3 similar traits, whereas the value of − 1 indicates the highest autocorrelation of dissimilar traits.
We illustrated the ancestral state of climatic PCs to assess the pattern of climatic evolution of Tiger geckos by using the "contMap" function of the "phytools" package in R software (Revell 2012(Revell , 2013. For climatic input data, we averaged the climatic spectrum of each climatic PC from the PNO profiles for each Tiger gecko. The best model to describe trait evolution was ranked based on the AICs tests for four commonly proposed models: Brownian motion (BM), Ornstein-Uhlenbeck (OU), Early Burst (EB) and Lambda, using the "fitContinuous" function of "geiger" package (Pennell et al. 2014). We computed the maximum likelihood (ML) to predict the climatic state following the equation from Felsenstein (1985).

Phylogenetic relationships and divergence dating
The time-calibrated BEAST analysis recovered a phylogeny with well-supported nodes (BPP ≥ 90) throughout the tree (Fig. 2). The divergence dates are similar to those calculated Fig. 2 The ancestral geographical range of Goniurosaurus species reconstructed under BioGeoBEARS using DEC + J model on the dated BEAST tree. Colors indicate the current biogeographic unit of each species. Pie charts with colors at the nodes indicate the estimated probability of each area inherited from the ancestor after the immediate cladogenetic process by Liang et al. (2018) and fall well within the range of their calculated highest posterior density (HPD).
Based on our dated phylogeny, the first divergence of Goniurosaurus was estimated to happen during the Eocene at approximately ~ 45.3 mya and the radiation continued across East Asia up until the Pliocene. Diversification within four monophyletic species groups began in the mid-Miocene between ~ 13.4 and 7.7 mya (million years ago) and continued to at least the early Pleistocene ~ 2 mya (Fig. 2).

Ancestral reconstruction of geographical ranges and elevation
The BioGeoBEARS model indicated that the DEC + J model recovered the best fit to the data and most likely to infer the correct ancestral range at each node being that it had the lowest AIC-wt score (Supplementary Information Table S2). The range analyses estimated that the respective ancestral regions for each group are the same regions in which each group currently inhabits. Liang et al. (2018) considered this to be a recent dispersal event. However, the geographical origin of Goniurosaurus is not clearly specified. The analysis indicated that all the deep nodes outside those of the species groups, recovered an essentially equal probability that Goniurosaurus originated anywhere across the combined areas. Given the young age and relatively recent geographic isolation of the Ryukyu Archipelago and Hainan Island from the late Mioence (Sterling et al. 2006;Honda et al. 2014;Liang et al. 2018), as for now, we assume that the Goniurosaurus evolved in continental East Asia. With additional fieldwork and the potential discovery of new species in the north of Guangdong Province, the region of origin could be further refined.

History of altitudinal transitions
All Goniurosaurus species have been recorded in areas ranging from 40 to 1000 m elv ( Supplementary Information Fig. S1). The AIC scores for the three transition models were ARD = 47.0085, SYM = 43.60873, and ER = 40.07807. The SCM using the best ER model predicted that the "low" level is the most probable ancestral state at the crown tree and all group nodes, except for the "high" level at the common node of three lineages (namely G. kwanghua, G. hainanensis, G. lichtenfelderi) (Fig. 3).

Niche overlap comparison
Climatic niche spaces of three Goniurosaurus groups in China and Vietnam were illustrated in the PCA-env analysis and its first two dimensions accounted for 65.5% of overall variance in the climatic data set (whereof, PC1: 38.6%-explained by Mean Diurnal Range (BIO2), Precipitation of Driest Month (BIO14) and Warmest Quarter (BIO18); and PC2: 26.9%-related by Isothermality (BIO3), Max Temperature of Warmest Month (BIO5), Mean Temperature of Wettest Quarter (BIO8), Precipitation of Wettest Month (BIO13)) (Fig. 4). The analysis revealed no overlaps among their climatic spaces, especially the niche space of the yingdeensis group which is entirely separated from the remaining groups (Fig. 4). The niche overlap between group-pairs in terms of Schoener's D only ranges from nearly 0 (yingdeensis and two remaining groups) to very limited overlap with D = 0.15 (lichtenfelderi and luii groups). Furthermore, the climatic spaces occupied by the three 1 3 groups were neither significantly similar (p > 0.05) nor equivalent (p = 1.00) in all grouppair comparisons ( Supplementary Information Fig. S2, Table S3).
The Jackknifing approaches omitting records of each species from the pooled data of each group documented low "Schoener's D"-high "1-D" values in some cases, including G. lichtenfelderi (1-D = 0.998) omitted from the lichtenfelderi group, G. catbaensis (0.844) and G. liboensis (0.418) from the luii group, and G. varius (0.61) from the yingdeensis group (Table 1). Our results of the PCA-env analysis further noted in these cases that the niche space of each group shrunk significantly ( Supplementary Information Fig. S3).

History of climatic evolution
The first four climatic PCs from the PCA analysis for SDM account for 90% of overall climatic variance, whereof PC1: 39.3% was mainly associated with Mean Temperature of Wettest Quarter (BIO8), PC2: 23.9% with Precipitation of Driest Month (BIO14), PC3: 18.1% with Max Temperature of Warmest Month (BIO5) and PC4: 7.8% with Isothermality (BIO3). These PCs variables were employed together with geographic records for the SDMs. The indices of TSS and AUC indicated that the fitted EP-ESMs of nine selected Tiger geckos performed well ( Supplementary Information Fig. S4). Using  Supplementary Information Fig. S5), we reconstructed on the BEAST dated tree the ancestral niche occupancy profiles of Chinese-Vietnamese Goniurosaurus species, which presented only the niche evolution pattern of divergence in given Goniurosaurus groups and intermixed with the conservatism pattern among their lineages ( Supplementary Information Fig. S5).
Phylogenetic signals in climatic variables were only detected in PC1 and PC4 (p-values ≤ 0.05), but these signals were relatively low according to Blomberg's K values, being < 1 and Moran's I values only reached 0.16 (Supplementary Information  Table S4). The AICc test suggested the Brownian Motion (BM) model was the most likely scenario (Supplementary Information Table S4). Consequently, the ancestral state estimation also recovered a mix of both climatic patterns (Fig. 6). In particular, the pattern of divergence was documented among three Goniurosaurus groups in PC4 and PC2, and between the yingdeensis group and two remaining groups in PC1 and PC3. Given among their lineages, almost Goniurosaurus species followed a conservative Fig. 4 Climatic niche space occupied by Goniurosaurus groups according to the PCA-env analysis including the lichtenfelderi group (granite-stream adapted group)-luii group and yingdeensis group (karst adapted group), and contribution of climatic variables presented along the first two axes (Correlation circle: PC1 and PC2). The solid (100%) and dashed contour (50%) lines illustrate the available macro-climate space pattern, except for G. lichtenfelderi of the lichtenfelderi group in four PC variables, and G. catbaensis in the first three PCs, G. kwangsiensis in PC2 and PC3 and G. liboensis in PC1 of luii group (Figs. 5, 6).

Evolution
Given estimations of the ancestral range in this study and Liang et al. (2018), the origin of the ancestral range of Goniurosaurus was predicted to be somewhere in continental Eastern Asia, including the Ryukyu Archipelago and Hainan Island prior to the Eocene, as opposed to a specified location. On the other hand, our analyses revealed that the current regions in which each monophyletic Goniurosaurus group radiated, are respectively their own ancestral region of origin. Dispersals of populations of insular groups (namely kuroiwae and lichtenfelderi) were therefore considered related to prior tectonic events following the late Paleocene, such as the collision of India and Eurasia and the clockwise rotation of the Philippine Sea plate (Seno et al. 1993;Ota 1998;Teng and Lin 2004;Sterling et al. 2006;Liang et al. 2018). Subsequent seaways served as a dispersal barrier that prevented gene flow among populations. Therefore, vicariant evolution resulted in the formation of  (Grismer et al. 2021). We used ecological niche models to elucidate the role of the Grinnellian niche (e.g., climate conditions) in the diversification of Goniurosaurus in China and Vietnam. Accordingly, the niche divergence in all PCs based on climatic variables was strongly correlated with the formation of the yindeensis group by the middle Oligocene (Figs. 4,5,6). On the other hand, the dominant climate of the tropical monsoon toward southeastern China and northern Vietnam was caused by the uplift of Tibetan Plateau as a result of India's collision with Eurasia as well (Sterling et al. 2006;Miao et al. 2013). The subsequent envelopment of monsoon tropical conditions since the Oligocene, might have also driven the cladogenesis of the two remaining groups (e.g., lichtenfelderi and luii). In particular, we documented the very low overlap between their niche spaces in the PCA-env analysis (Fig. 4) and the discordance in their climatic spectrums (e.g., PC2 and PC4) in the PNO estimation (Figs. 5,6). Furthermore,  Grismer et al. (2021) estimated that karst ecosystems are the most probable ancestral habitat for most Goniurosaurus groups, except for the lichtenfelderi group which mostly inhabits nonkarst environments (Grismer et al. 2021).
Hainan Island is strongly supported as the ancestral range of the lichtenfelderi group (Liang et al. 2018). A mixture of karst and granite formations on Hainan Island that may have promoted the phylogenetic split between G. bawanglingensis and the most recent common ancestor of G. zhoui and G. kwanghua, following the evolution of habitat preference (Liang et al. 2018;Grismer et al. 2021). The emergence of other species within the lichtenfelderi group might generally be associated with climatic cycles during the Pliocene and early Pleistocene (Herzschuh et al. 2016). Accordingly, low-to-high transitions along the elevation gradient toward highland refugia on Hainan Island occurred during the interglacial mid-Pliocene, and thereby tiger geckos (the common ancestor of G. kwanghua, G. hainanensis and G. lichtenfelderi) avoided much warmer climatic conditions. The latest speciation event between G. hainanensis and G. lichtenfelderi in non-karst habitats during the early Pleistocene might relate to altitudinal partitioning as well, but vice versa. In particular, the Hainan populations of Goniurosaurus dispersed overwater or overland to occupy low elevations in the Vietnamese mainland during sequential glacial periods of sea level decline (Li et al. 2010;Herzschuh et al. 2016;Liang et al. 2018). Following a pattern of vicariance, these populations were gradually isolated by increasing the sea levels, and eventually evolved into G. lichtenfelderi (~ 2 mya) at both island and mainland sites in Vietnam.
Besides the ocean serving as barrier, the latest speciation event in the lichtenfelderi group is also supported by the non-physical barrier of climatic niche limits. We documented the pattern of niche divergence in all selected climatic conditions between G. lichtenfelderi and Hainan populations of Goniurosaurus (Figs. 5, 6). The role of Grinnellian niche in cladogenetic processes was further assessed in other groups of Goniurosaurus. In general, the low phylogenetic signals within Goniurosaurus were strongly influenced by various climate spectra in the luii group. In particular, three Goniurosaurus species of the luii group followed the pattern of climatic divergence, including G. catbaensis in most climatic conditions (except for Isothermality of PC4), G. kwangsiensis in Precipitation of Driest Month (PC2) and Max Temperature of Warmest Month (PC3), and G. liboensis in Mean Temperature of Wettest Quarter (PC1) (Figs. 5, 6). On the other hand, the remaining species in the luii group and all species in the yingdeensis group adapted to their convergent gradients of all climatic conditions, representing a high degree of niche conservatism (Figs. 5,6). Accordingly, the radiation of their common ancestors in Vietnam and China mainland might be enveloped within homogenous niche spaces. The subsequent speciation events among these populations of Goniurosaurus might be related to potential barriers created by disjunct karst ecosystems, and/ or large rivers or canyons offering various micro-niches (Clements et al. 2006;Chen et al. 2014;Qi et al. 2020;Grismer et al. 2021;Ngo et al. 2021bNgo et al. , 2022bZhu et al. 2021). For example, the Zuojiang and Ky Cung rivers might establish an eastern boundary between G. luii and G. araneus in China and G. huuliensis in Vietnam, respectively (Ngo et al. 2022b). Given a case of uncertain taxonomy in Goniurosaurus, the high level of niche conservatism agreed well with the morphological and genetic similarities between G. kadoorieorum and G. luii, confirming the junior synonym status as proposed by Grismer et al. (2021) and Ngo et al. (2021b).

Conservation
Due to high degrees of endemism and small population sizes, all species of Goniurosaurus are particularly vulnerable to human impacts (Ngo et al. 2019(Ngo et al. , 2022bGrismer et al. 2021). In this study, our results suggest that the closely related species within the yingdeensis group, the lichtenfelderi group on Hainan Island and some in the luii group, show patterns of niche conservatism and hence appear to be more susceptible to climate change (Wiens and Graham. 2005;Hadly et al. 2009;Wiens et al. 2010;Peterson 2011;Lavergne et al. 2013;Pyron et al. 2014;Ahmadi et al. 2021). Their probability of survival in the future has been even reduced decisively by the loss of genetic diversity that likely occurs in most Goniurosaurus species due to population declines (e.g., overexploitation) and habitat destruction (e.g., forest fragmentation, quarrying of limestone karst and infrastructure) (Frankham et al. 2002;Ngo et al. 2019;2022b, c). To avoid the extinction pressure of climate change, a natural response of species are range shifts towards optimal environments (Sinervo et al. 2010). However, the range-restricted tiger geckos are unable to migrate larger distances to climatic refugia, which were forecasted to be greatest towards higher latitudes in proximity to temperate regions, due to either habitat specialization or lacking forest corridors (Williams et al. 2007;Ngo et al. 2021aNgo et al. , 2022a. The dense forest canopy, wide elevation gradient and complex karst formations likely construct unique micro-refugia containing favorable micro-climate conditions which could buffer species against suboptimal climates (Clements et al. 2006;Sterling et al. 2006;Ohlemüller et al. 2008;Suggitt et al. 2018). However, this last survival chance is currently getting narrower due to the intensive loss of micro-refugia as the result of deforestation and quarrying in karst formations (Clements et al. 2006;Queiroz et al. 2013;CEPF 2020;Ngo et al. 2022b, c). Therefore, higher conservation priorities and stringent protection should be immediately implemented to prevent these niche-conservative species of Goniurosaurus from entering an extinction vortex. However, this modelling-based assessment does not necessarily indicate that the remaining Goniurosaurus species can overcome environmental changes to adapt. In particular, both G. lichtenfelderi and G. catbaensis showing the pattern of niche divergence while G. huuliensis following the pattern of niche conservatism, all three of them are potentially impacted by climate change (Le et al. 2017;Ngo et al. 2021aNgo et al. , 2022a. Although approaches to identify the pattern of niche evolution of Goniurosaurus are necessary for conservation priorities, the species' endangered level should be further assessed based on comprehensive biological knowledge of population status, ecological characteristics, and human impacts. The combination of biological background data is expected to significantly improve the efficacy of conservation plans and actions to safeguard the future existence of the species of Goniurosaurus.

3
QTJP01.02/23-25). Field surveys were partially funded by Cologne Zoo, the Rufford Foundation (Project: 30597-1). Cologne Zoo is a partner of the World Association of Zoos and Aquariums (WAZA): Conservation Project 07011, 07012 (Herpetodiversity Research, Amphibian and Reptilian Breeding and Rescue Stations). The research of Hai Ngoc Ngo in Germany is funded by the German Academic Exchange Service (DAAD).

Data availability
The data that support the findings of this study are available from the corresponding author [Hai Ngoc Ngo], upon reasonable request.

Declarations
Conflict of interest Authors declare no conflict of interest.
Informed consent All the authors give the consent for publication of the article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.