We examined 97 preserved specimens of Goniurosaurus, including 25 specimens of G. bawanglingensis, four of G. sinensis (= G. kwanghua; = G. hainanensis sensu stricto in results), 38 of G. hainanensis from the east side of the island (= G. cf. lichtenfelderi in results) and 13 of G. hainanensis (= G. hainanensis sensu stricto in results) from the west side of the island, and 17 of G. zhoui. Due to the poor preservation of some specimens, only 23 of G. bawanglingensis, 24 of G. hainanensis from the east side of the island, 13 of western clade G. hainanensis, and 18 of G. zhoui were included in the subsequent statistical analyses. The specimens examined in this study are deposited in ECNUEast China Normal University; HNNU Hainan Normal University; LSUHCLa Sierra University Herpetological Collection; MCZMuseum of Comparative Zoology, Harvard University; MVZMuseum of Vertebrate Zoology, University of California, Berkeley; SYSSun Yat-sen University.

The morphological characters examined followed Grismer et al. (2002), Zhou et al. (2018), Zhou et al. (2019), and Zhu et al. (2020b). External measurements were taken with digital calipers (Neiko 01407A Stainless Steel 6-Inch Digital Caliper, USA) to the nearest 0.1 mm on the right side of each individual by Shuo Qi. Abbreviations of morphological characters are as follows: SVL snout-vent length, from tip of snout to vent; AG axilla to groin length, from posterior edge of forelimb insertion to anterior edge of hind limb insertion; ID internarial distance, distance between nares; HL head length, from the tip of snout to posterior edge of occiput; HW maximum head width; IO interorbital distance, distance between posteriormost points of eyes; AD diameter of auditory meatus; SL snout to eye distance, measured from tip of snout to anteriormost point of eye; ED greatest diameter of eye; FLL forelimb length, from axilla to the tip of the fourth finger; HLL hind limb length, from groin to the tip of the fourth toe. SPL supralabials; IFL infralabials; N nasal scales surrounding nare; IN internasals; PostIN granular scales bordering the internasals;PM postmentals; GP gular scales bordering postmentals; CIL eyelid fringe scales or ciliaria; PO preorbital scales, number of scales in a line from the posterior margin of external naris to the anterior margin of the bony orbit; GST granular scales surrounding dorsal tubercles; PTL paravertebral tubercles between limb insertions; DTR longitudinal dorsal tubercle rows at midbody; MB scales around midbody; PP precloacal pores; PAT postcloacal tubercles. Bilateral scale counts are given as left/right.

All statistical analyses were conducted using R version 4.4.2 (R Core Team 2024). A multiple factor analysis (MFA) using the R package FactorMineR (Husson et al. 2017) and visualized using the Factoextra package (Kassambara and Mundt 2017) was used to compare differences and similarities in morphospace of individuals from the G. lichtenfelderi group. The MFA used a concatenated data set comprised of 12 size-corrected morphometric characters (SVL, AG, ID, HL, HW, IO, AD, SL, ED, EE, FLL, HLL) and 15 meristic characters (SPL, IFL, N, IN, PostN, PM, GP, CIL, PO, GST, PTL, DTR, MB, PP, PAT). To remove potential effects of allometry in the morphometric characters (see Chan and Grismer 2022), measurements were size-corrected using the following equation: X_adj = log(X)-β[log(SVL)-log(SVL_mean)], where X_adj = adjusted value; X = measured value; β = unstandardized regression coefficient for each population; and SVL_mean = overall average SVL of all populations (Thorpe 1975, 1983; Turan 1999; Lleonart et al. 2000). MFA is a global, unsupervised, multivariate analysis that incorporates qualitative and quantitative data (Pagès 2015) simultaneously, making it possible to include different data types in a nearly total morphological evidence environment. All data types are standardized, preventing one data type from overleveraging the output. The MFA dataset is in Table S1.

A PERMANOVA analysis from the vegan package 2.5–3 in R (Oksanen et al. 2020) was used to determine if the centroid locations and group clusters of each species/population from the MFA were statistically different from one another (Skalski et al. 2018) based on the load scores of dimensions 1–5. Using load scores as opposed to raw data, which are normally used, allows for the incorporation of the categorical characters which cannot be run untransformed in a PERMANOVA. All load scores for the PCAs, however, were used. The analysis calculates a Euclidean (dis)similarity matrix using 50,000 permutations. A pairwise post hoc test calculates the differences between the populations, generating a Bonferroni-adjusted p value and a pseudo-F ratio (F statistic). A p < 0.05 is considered significant and larger F statistics indicate more pronounced group separation. A rejection of the null hypothesis (i.e., centroid positions and the spread of the data points [i.e., clusters] are no different from random) signifies a statistically significant difference between species/populations.

One-way analysis of variance (ANOVA) was performed on a dataset coded for species to examine statistically significant mean differences (p < 0.05) among characters using car R package. Character means showing significant differences were subjected to a Tukey HSD test to determine which pairs of species differed significantly for those specific characters. Violin plots and the inserted/independent boxplots were generated to visualize the range, frequency, mean, 50% quartile, and degree of differences between the dependent variables for datasets with statistically significant mean differences using ggplot2 R package.

Mitochondrial and nuclear DNA. Seventy-eight newly collected samples of Goniurosaurus geckos were used in this study, including 53 from the G. lichtenfelderi group, 18 from the G. luii group, and 7 from the G. yingdeensis group. All Goniurosaurus specimens were collected during field surveys from 2014 to 2023. Prior to February 5^th, 2021—the date on which the newly revised List of National Key Protected Wild Animals in China (National Forestry and Grassland Administration and Ministry of Agriculture and Rural Affairs of the People’s Republic of China 2021) came into effect, specimens were retained as vouchers, and tissues for DNA extraction were mainly taken from muscle or tail tips. After this date, specimens were no longer retained but were released at the site of capture following non-lethal sampling, with either oral swabs or tail tip tissue collected for molecular analyses. This protocol was implemented to ensure full compliance with updated national wildlife protection regulations, which prohibit the long-term removal of protected species from the wild. DNA was extracted from each tissue sample using a standard extraction kit (Tiangen Biotech, Beijing, China). Two mitochondrial genomic fragments with 495 base pairs (bp) of 16S ribosomal RNA (16S) and 396 bp of cytochrome b (cyt b) and two nuclear genes with 369 base pairs (bp) of C-mos and 990 bp of Rag1 were amplified in this study. Primers used in this study are listed in Table 1. PCR amplifications were performed in a 20 μL reaction volume with the following cycling conditions: an initial denaturing step at 94°C for 5 min, 35 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 30 s, extension at 72°C for 1 min, and a final extension step of 72°C for 10 min (Liang et al. 2018). PCR products were purified with spin columns and then sequenced with forward primers using BigDye Terminator Cycle Sequencing Kit as per the guidelines on an ABI Prism 3730 automated DNA sequencer by Wuhan Tianyi Huiyuan Bioscience & Technology Inc.

For the Vietnamese samples, DNA was extracted and amplified at Hanoi University of Science (HUS), Hanoi, Vietnam, following the protocol described in Ngo et al. (2021).

A total of 561 sequences were used in the phylogenetic analyses, including 543 sequences of Goniurosaurus species and 18 outgroup sequences (see Table S2 for detail). Among them, 280 newly sequenced data have been deposited in GenBase (https://ngdc.cncb.ac.cn/genbase), and all other sequences were downloaded from GenBank (https://www.ncbi.nlm.nih.gov/genbank).

ddRAD library preparation and sequencing. To further evaluate putative species boundaries and clades recovered in the Sanger sequencing analyses in a more genomically comprehensive context, a small subset of samples was selected for double-digest restriction site-associated DNA sequencing (ddRAD-seq). Genomic DNA was extracted from 76 individuals of the G. lichtenfelderi group, as well as an outgroup sample of G. yingdeensis, using a universal DNA extraction kit (GenStar, Beijing, China), following the manufacturer’s instructions. DNA quality and concentration were assessed by agarose gel electrophoresis and a Qubit fluorometer (Thermo Fisher Scientific, USA). Library construction and sequencing were conducted by Guangzhou Jierui Bioscience & Technology Inc. Briefly, 200 ng of genomic DNA from each sample was digested with EcoRI and PstI (New England Biolabs, USA) and ligated to adapters containing unique barcodes and compatible overhangs using T4 DNA ligase (New England Biolabs). Barcoded samples were pooled in equal volumes, and fragments in the 350–550 bp range were isolated through agarose gel electrophoresis and purified (Omega Bio-tek, Nocross, USA). The pooled library was PCR-amplified and sequenced on an Illumina NovaSeq platform using 150 bp paired-end reads.

ddRAD-seq bioinformatics and SNP calling. Raw reads were demultiplexed and quality-filtered with the process_radtags module in Stacks v2.6 (Catchen et al. 2013; Rochette et al. 2019), discarding low-quality reads and adapter sequences. Reads were truncated to 135 bp and had to be at least 140 bp long to be retained. De novo assembly was performed using ustacks (minimum coverage m = 2, maximum mismatches per locus M = 10), and loci were merged across individuals with cstacks (number of mismatches allowed = 10). Each sample was then aligned against the catalog using sstacks. The tsv2bam and gstacks modules carried out genotype calling under a maximum-likelihood framework, after which populations filtered out low-quality and potentially linked SNPs. Following these steps, a total of 5301 high-quality SNPs, spanning 5416 RAD loci were retained for downstream population structure analyses and phylogenetic inference of the G. lichtenfelderi group.

Sequence-based phylogenetic inference. DNA sequences were aligned using the MAFFT algorithm (Katoh and Standley 2013) and trimmed with gaps partially deleted in Gblocks (Talavera and Castresana 2007). We generated up to 2242 nucleotides of aligned sequence data including partial sequences for two mitochondrial genes (16S, 487 bp; cyt b, 396 bp) and two nuclear genes (C-mos, 369 bp; Rag1, 990 bp). For Bayesian phylogenetic inference in MrBayes 3.2.4 (Ronquist et al. 2012), PartitionFinder v2.1 (Lanfear et al. 2017) was used to determine the best-fit partitioning scheme and nucleotide substitution models. Two independent MCMC runs were conducted with 10,000,000 generations each, sampling every 1000 generations and discarding the first 25% as burn-in, achieving a potential scale reduction factor (PSRF) < 0.005. Maximum likelihood analysis was performed in IQ-TREE v2.0 (Minh et al. 2020) using the BIC criterion implemented in ModelFinder to select the best substitution model (Kalyaanamoorthy et al. 2017). A bootstrap consensus tree was inferred from 1000 replicates. Node support was evaluated using Bayesian posterior probabilities (BPP), ultrafast bootstrap (UFB) values, and SH-aLRT (approximate likelihood ratio test) values. Nodes with BPP ≥ 0.95, UFB ≥ 95, and SH-aLRT ≥ 80 were considered strongly supported.

To further examine genetic divergence, uncorrected pairwise sequence distances for the cyt b gene were calculated in MEGA 6 (Tamura et al. 2013) using the Kimura 2-parameter (K2P) model. Mean genetic distances were computed both within and between species. Haplotype networks were constructed using PopART v1.7 (Leigh and Bryant 2015) with the Median Joining Network algorithm. All phylogenetic trees were visualized and edited in TVBOT v2.6.1 (Xie et al. 2023).

SNP-based Structure and ML Tree Inference. We used a model-based clustering approach implemented in STRUCTURE v2.3.4 (Pritchard et al. 2000). Analyses were conducted under an admixture model with correlated allele frequencies. The SNP dataset, formatted for STRUCTURE, was analyzed across a range of K values to determine the most likely number of genetic clusters. For each K, we performed multiple independent runs with a burn-in of 100,000 iterations followed by 1,000,000 MCMC iterations. Results were averaged across replicates, and the optimal K was identified using the ΔK method of Evanno et al. (2005), based on the rate of change in log-likelihood values [LnP(D)]. Individuals were then assigned to clusters according to their estimated ancestry proportions. We also verified the consistency of clustering patterns using alternative methods, although STRUCTURE served as the primary tool for inference. These results revealed clear genetic subdivisions and admixture patterns within the G. lichtenfelderi group.

To reconstruct phylogenetic relationships, we reconstructed a maximum likelihood (ML) tree based on the SNP dataset. A concatenated alignment including one SNP per RAD locus was generated in PHYLIP format using Stacks output. Tree inference was performed in IQ-TREE v2.0, which selected the best-fit substitution model via ModelFinder. Branch support was assessed with 1000 ultrafast bootstrap replicates. The resulting tree provided a well-resolved phylogeny of the G. lichtenfelderi group. An outgroup from a congeneric species outside the group was used for rooting, allowing interpretation of evolutionary direction. All analyses were conducted with default parameters unless specified otherwise, and the final tree was visualized using TVBOT v2.6.1.

Dating analyses. We estimated divergence times using PAML v4.9j (Yang 2007) based on a concatenated dataset of four genetic markers. A relaxed molecular clock model (clock = 2) was applied to account for rate variation among lineages. The F84 model was used to model nucleotide substitution rates, as determined by prior model testing. The approximate likelihood method (usedata = 2) was employed to improve computational efficiency for large datasets. Among-site rate heterogeneity was modeled using a gamma distribution with a shape parameter (alpha = 0.5) and five discrete categories (ncatG = 5).

Priors were set following standard practice: the mean substitution rate prior (rgene_gamma) was set with α = 2 and β = 20 (corresponding to 0.1 substitutions/site/100 million years ago, Mya), and the rate variance across branches prior (sigma2_gamma) was defined as α = 1 and β = 10. MCMC chains were run for 6,000,000 iterations, sampling every 10 iterations, with a burn-in of 1,000,000, resulting in 500,000 samples used for posterior estimation. Convergence was verified using standard diagnostics.

As there are no reliable fossil records for Goniurosaurus, we applied secondary calibrations based on divergence estimates within Eublepharidae (Agarwal et al. 2022). Additionally, three fossil calibrations from Gekkota (Agarwal et al. 2020) were used: Burmese amber fossils for crown Gekkota (offset = 99 Mya); Pygopus hortulanus for the stem of Pygopus (offset = 23 Mya); the divergence between Phelsuma inexpectata and P. ornata (uniform prior: 0.05–5 Mya). We also included the earliest fossil record of Tarentola (offset = 11.6 Mya) as a genus-level constraint. Final divergence estimates are reported as medians with 95% highest posterior density (HPD) intervals in millions of years ago (Mya). Sample information for dating analyses see Table S3.

Phylogenetic inference from mitochondrial and nuclear DNA. In Figure 1A, both Bayesian Inference (BI) and Maximum Likelihood (ML) phylogenetic analyses were conducted based on a concatenated dataset comprising two mitochondrial (16S, cyt b) and two nuclear (C-mos, Rag1) gene fragments. These combined analyses consistently recovered the same interspecific relationships within the G. lichtenfelderi group, providing robust support for species delimitation. Both methods strongly supported the monophyly of G. bawanglingensis and G. zhoui, with high nodal support values (BPP = 0.99/1.00; UFB = 99/100; SH-aLRT = 97.8/100), thereby validating their status as distinct evolutionary lineages. In contrast, G. hainanensis is recovered as polyphyletic, forming two distinct clades—eastern and western. Goniurosaurus sinensis was strongly supported but was embedded within the western clade of G. hainanensis (BPP = 1.00; UFB = 98; SH-aLRT = 94.6). The sister relationship between G. lichtenfelderi and the eastern clade of G. hainanensis received weak support (BPP = 0.61; UFB = 79; SH-aLRT = 68.9), indicating that this relationship should be interpreted cautiously. Although G. sinensis is strongly supported it is embedded within the western clade of G. hainanensis (BPP = 100), and together they form a monophyletic lineage with G. lichtenfelderi and the eastern clade of G. hainanensis. Based on mitochondrial cyt b genetic distances (Table 2), G. bawanglingensis and G. zhoui exhibit substantial genetic divergence from members of the G. lichtenfelderi species group. The mean interspecific genetic distances were 17.52–19.81% for G. bawanglingensis and 14.24–16.25% for G. zhoui, respectively, supporting their recognition as distinct evolutionary lineages. The mean genetic distance between G. lichtenfelderi and G. sinensis was 7.20%, while G. lichtenfelderi showed a distance of 7.69% to the western clade of G. hainanensis and 6.13% to the eastern clade. In contrast, the western clade of G. hainanensis was genetically very close to G. sinensis, with a mean distance of only 2.45%, whereas its distance to the eastern clade of G. hainanensis was 7.69%.

Haplotype network analysis based on mitochondrial sequence (cyt b) revealed a clear pattern of genetic differentiation within the G. lichtenfelderi group (Fig. 1B). Goniurosaurus bawanglingensis and G. zhoui each formed distinct and compact haplotype clusters, consistent with their monophyly recovered in the phylogenetic trees, thereby supporting their validity as independent species. In contrast, the haplotypes of G. hainanensis, G. sinensis, and G. lichtenfelderi showed complex relationships. The western clade of G. hainanensis shared highly similar and closely connected haplotypes with G. sinensis, without clear separation, indicating minimal genetic differentiation between them. The eastern clade of G. hainanensis exhibited haplotypes that were very similar to those of G. lichtenfelderi; however, no intermixing was observed, with G. lichtenfelderi maintaining a distinct and well-separated cluster.

Taken together, the phylogenetic results suggest that the taxon currently referred to as G. hainanensis on Hainan Island may comprise two deeply divergent evolutionary lineages, but this inference requires further validation with higher-resolution genomic data.

Phylogenetic inference from SNP data. The maximum likelihood (ML) tree inferred from concatenated SNP data (one SNP per locus) recovered five distinct clades, each corresponding to putatively independent lineages (Fig. 2A). The outgroup sample (G. yingdeensis) was included in the ddRAD-seq dataset for rooting purposes but is not displayed in the phylogenetic trees presented in this study. The phylogeny exhibited strong support at nearly all nodes (black dots: bootstrap >90%), and clearly resolved the relationships among geographically isolated populations. Individuals from the same locality or morphologically defined species formed cohesive genetic clusters, such as G. bawanglingensis, G. zhoui, G. lichtenfelderi, and the eastern (G. cf. lichtenfelderi) and western clades (G. hainanensis sensu stricto) of G. hainanensis. Compared to traditional sequence-based markers, the SNP data provided higher resolution, revealing subtle genetic divergences even between geographically proximate populations. The SNP phylogeny also supports the paraphyly of G. hainanensis on Hainan Island, consistent with results from mitochondrial and nuclear genes.

Model-based clustering implemented in STRUCTURE, evaluated across K = 2–10, iden­ti­fied K = 3 as the optimal number of genetic clusters according to the ΔK method (Fig. 2B). Under this optimal clustering scheme (K = 3), the individuals of the G. lichtenfelderi group were assigned to three genetic clusters that only partially correspond to phylogeny. The first cluster (red) includes both the eastern and western clades of G. hainanensis along with continental populations of G. lichtenfelderi, reflecting genetic admixture or shared ancestry signals at the population level, rather than distinct population structuring. The second cluster (yellow) corresponds primarily to G. bawanglingensis, and the third cluster (blue) to G. zhoui. However, the admixture plot reveals signs of recent gene flow between G. bawanglingensis and G. zhoui, particularly in populations from locality Exianling (EX). Individuals from this site exhibit varying proportions of ancestry from both genetic clusters, indicating potential hybridization or recent introgression. The uneven admixture proportions among individuals from the same site further suggest that this genetic mixing may be ongoing and not the result of ancient shared ancestry.

Divergence dating. The MCMCtree analysis (Fig. 3) based on concatenated mitochondrial (16S, cyt b) and nuclear (C-mos, Rag1) genes. Effective sample sizes (ESS) for all parameters exceeded 200, indicating good mixing and convergence of the analysis.

Divergence time estimates place the most recent common ancestor (MRCA) of Eublepharidae at 82.9 Mya (95% HPD: 93–73 Mya). The divergence between Goniurosaurus and the clade comprising Eublepharis and two African genera occurred at 59.4 Mya (95% HPD: 67–53 Mya). Within Goniurosaurus, G. kuroiwae group diverged around 43.1 Ma (95% HPD: 53–34 Mya), followed by G. yingdeensis group at 34.2 Mya (95% HPD: 42–26 Mya), and finally, the divergence between the G. lichtenfelderi group and G. luii group occurred at 28.7 Mya (95% HPD: 36–22 Mya). Within the G. lichtenfelderi group, G. bawanglingensis is the earliest diverging lineage, with an estimated divergence time of 14.7 Mya (95% HPD: 20–10 Mya); followed by G. zhoui at 10.9 Mya (95% HPD: 15–7 Ma); Subsequently, two major clades emerged: one comprising G. lichtenfelderi and the eastern clade of G. hainanensis, and the other comprising the western clade of G. hainanensis and G. sinensis, with a divergence time of 5.3 Mya (95% HPD: 8–3 Mya). The divergence between G. lichtenfelderi and the eastern G. hainanensis clade occurred at 3.8 Mya (95% HPD: 6–2 Mya), while the western G. hainanensis clade diverged from G. sinensis much more recently, at 1.1 Mya (95% HPD: 2–0.3 Mya).

Statistics for morphometric data. The MFA and subsequent PERMANOVA recovered statistically different morpho-spatial differences among all species pairs except the eastern and western clades of G. hainanensis (Fig. 4; Table 3). Dimension 1 in the MFA accounted for 18.3% of the variation in the data set and loaded most heavily for the meristic data and dimension 2 accounted for an additional 16.9% of the variation. The PERMANOVA analysis indicated that each species’ centroid position is statistically different from others in the morphological space, except that between the eastern and western clades (adjusted p value = 1). However, the ANOVA analyses still indicated that the eastern and western clades differ significantly in head length (HL, p = 0.003), head width (HW, p = 0.019), hind limb length (HLL, p = 0.022), Internarial distance (ID, p = 0.002), and number of postmental scales (PM, p = 0.006) (Table 4). This can also be seen in the violin plots with inserted boxplots for morphometric data and the boxplots with inserted scatter plot for meristic data (Figs 5, 6). These results suggested that the eastern and western clades of G. hainanensis are morphologically different from each other in some characters despite the relatively smaller divergences compared to the G. bawanglingensis and G. zhoui, that is consistent with the conclusion of molecular analyses.

The results presented here provide a robust genetic framework for understanding species boundaries within the G. lichtenfelderi group. Both sequence- and SNP-based trees robustly support the distinctiveness of G. bawanglingensis and G. zhoui as independent evolutionary lineages. However, G. hainanensis is consistently recovered as polyphyletic, comprising two deeply divergent clades on Hainan Island. The western clade, which includes samples from the type locality Wuzhishan, clusters closely with G. sinensis in both sequence- and SNP-based phylogenies, with minimal genetic distance and recent divergence time estimates. Given this strong genetic continuity and the historical context of the type locality, we propose to treat G. sinensis as a junior synonym of G. hainanensis.

In contrast, the eastern clade of G. hainanensis forms a distinct lineage that is genetically closer to continental populations of G. lichtenfelderi than to its western counterpart. Moderate levels of genetic divergence exist between the eastern clade and G. lichtenfelderi, comparable to those observed between other recognized species within the genus (e.g., G. chengzheng vs. G. gezhi; G. luii vs. G. huuliensis). However, due to the lack of comprehensive morphological data for continental G. lichtenfelderi collected by the same observer and the absence of high-coverage genomic evidence to further evaluate the depth and consistency of divergence, we restrict the distribution of G. lichtenfelderi sensu stricto to the Asian mainland and tentatively treat the Hainan Island lineage as G. cf. lichtenfelderi. Future studies incorporating extensive morphological datasets and genome-wide data will be essential to determine whether this lineage warrants formal taxonomic recognition.