Genomic DNA was extracted from muscle tissue samples, using a DNA extraction kit from Tiangen Biotech (Beijing) Co., Ltd. Primers used for 16S were r16S-5L (5’- GGTMMYGCCTGCCCAGTG -3’) and 16Sbr-H (5’- CCGGTCTGAACTCAGATCACGT-3’) (Palumbi et al. 1991), for cyt b the primers were L14731 (5’- TGGTCTGAAAAACCATTGTTG-3’) (Honda et al. 2014) and H15149m (5’- GCMCCTCAGAAKGATATTTGYCCTCA-3’) (Chambers and MacAvoy 1999), for CMOS the primers were FU-F (5’- TTTGGTTCKGTCTACAAGGCTAC -3’) and FU-R (5’- AGGGAACATCCAAAGTCTCCAAT -3’) (Gamble et al., 2008) , and for RAG1 the primers were R13 (5’- TCTGAATGGAAATTCAAGCTGTT -3’) and R18 (5’- GATGCTGCCTCGGTCGGCCACCTTT -3’) (Groth and Barrowclough 1999). The PCR procedure was performed with an initial denaturation at 94 °C for 5 min, 35 cycles of 94 °C for 30 s, 55 °C for 30 s and 72 °C for 1 min, followed by a final extension at 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 Shanghai Majorbio Bio-pharm Technology Co., Ltd.

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, 16S (633 bp) and cyt b (1075 bp), and the nuclear genes, CMOS (472 bp) and RAG1 (890), from 103 specimens of 23 species of Goniurosaurus with varying degrees of sequence coverage across the samples (Table 1). Concatenation followed the comparison of separate gene trees to confirm there were no major discordances. One species, Eublepharis macularius, served as an outgroup (Grismer 1988; Jonniaux and Kumazawa 2008) to root the trees. Sequence data and GenBank accession numbers are listed in Table 1.

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 16S and cyt b and HKY+F for CMOS and RAG1. To avoid over parameterization, protein coding genes were not partitioned by codon. One-thousand bootstrap pseudoreplicates 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. Post-burn-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 (Huelsenbeck et al. 2001; Wilcox et al. 2002). 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 Markov chain Monte Carlo (MCMC). 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 (Huelsenbeck et al. 2001; Wilcox et al. 2002). We considered nodes with values of 0.90–0.94 as well-supported.

The BEAST tree was converted to newick format and pruned using the drop.tip () command (Paradis and Schliep 2018) in the R package ape [v.3.4.3] to include only the earliest diverged individual of each species. Habitat preference (karst or non-karst; see below) was mapped onto the tree using 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) values in the R package ape (Paradis and Schilep 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.

The coding of habitat preference for each species was determined from the literature and field observations of the authors (Table 2). A species’ habitat preference was coded as “karst” if it had a strong association with karstic habitats. Many such species may range into forested areas or even other rock types (granite or volcanic) if nearby but their densities are much greater in karstic areas. A species was coded as “non-karst” if detected only in forested areas or forested areas with other rock types (e.g. granite). These species never show any strong preference for karstic microhabitats even if such habitats exist within their range.

An ecomorphological analysis was conducted using four of the five recorded species from Vietnam (Grismer et al. 1999; Vu et al. 2006; Orlov et al. 2008; Ziegler et al. 2008; Nguyen et al. 2009; Nguyen 2011; Orlov et al. 2020) for which there existed a substantially large morphometric data set (Ngo et al. 2021): the karst-adapted G. catbaensis Ziegler, Nguyen, Schmitz, Stenke, and Rösler, 2008, G. huuliensis Orlov, Ryabov, Nguyen, Nguyen, and Ho, 2008, and G. luii Grismer, Viets, and Boyle, 1999 of the luii group and the granite stream-adapted G. lichtenfelderii (Mocquard, 1897) of the lichtenfelderi group (Figs 1, 2). A total of 486 live individuals and 54 museum specimens of four species were examined for morphological data, comprising 194 individuals of G. catbaensis (21 juveniles, 93 females, and 80 males), 80 individuals of G. huuliensis (two juveniles, 46 females, and 32 males), and 88 individuals of G. luii (11 juveniles, 43 females, and 34 males) of the luii species group and 178 individuals of G. lichtenfelderii (14 juveniles, 72 females, and 92 males) of the lichtenfelderi group.

Measurements were taken with dial calipers to the nearest 0.1 mm on the right side of each individual. Abbreviations are as follows: snout-vent length (SVL), from tip of snout to vent; axilla to groin length (AG), from posterior edge of forelimb insertion to anterior edge of hind limb insertion; maximum body width (BW), greatest width of torso, taken at level of midbody; maximum body height (BH), from dorsal surface of body to belly; internarial distance (ID), distance between nares; head length (HL), from the tip of snout to posterior edge of occiput; maximum head width (HW); cheek height (CH), from posterior edge of labial to top of head at parietal region; interorbital distance (IO), distance between posteriormost points of eyes; diameter of auditory meatus (AD); snout to eye distance (SL), measured from tip of snout to anteriormost point of eye; diameter of eye (ED), greatest diameter of eye; eye to ear distance (EE), from posterior margin of eye to posterior margin of ear; forelimb length (FLL), from axilla to the tip of the fourth finger; hind limb length (HLL), from groin to the tip of the fourth toe. To remove potential effects of allometry, size was adjusted 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)—accessable in the R package GroupStruct (available at https://github.com/chankinonn/GroupStruct). The morphometrics of each species were adjusted separately and then concatenated so as not to conflate intra- with interspecific variation (Reists 1986). All data were then scaled to their standard deviation to insure they were analyzed on the basis of correlation and not covariance and were log-transformed to insure they were normally distributed.

An analysis of variance (ANOVA) was performed on a data set coded for species to search for the presence of statistically significant mean differences (p < 0.05) among characters across the selected subset of species in the luii and lichtenfelderi groups. Character means bearing statistical differences among species were subjected to a TukeyHSD test to ascertain which species pairs differed significantly from each other for those particular characters. A Student t-test was also performed on a second data set coded for only habitat preference (karst versus non-karst) to search for the presence of statistically significant mean differences (p < 0.05) among the same subsets of species coded for habitat. Violin plots with inserted boxplots were generated in order to visualize the range, frequency, mean, 50% quartile, and degree of differences between the dependent variables for both data sets bearing statistically different mean values.

The morphospatial clustering of the two separate data sets (species and habitat preference) were visualized using principal component analysis (PCA) along the ordination of the first two principal components (PC) using the Adegenet package in R (Jombart et al. 2010) and implemented by the prcomp () command. The data were log-transformed prior to analysis in order to normalize their distribution so as to ensure characters with very large or very low values could not over-leverage the results owing to intervariable nonlinearity. All statistical analyses were performed using R.3.1.2 (R Core Team 2018).

The ML, BI, and BEAST analyses recovered strong nodal support (UFB 98–100/BPP 1.00) for the monophyly of all four species groups with the kuroiwae group being the strongly supported (100/1.00) sister group to the remaining three groups (Fig. 3). The ML analysis weakly recovered (88) the lichtenfelderi and yingdeensis groups as sister lineages but the BEAST analysis recovered this relationship with strong support (1.00; Fig. 4). The BI analysis recovered the luii and yingdeensis groups as sister lineages, although the support is so low (0.51), the three groups effectively form a polytomy. The ML and BI analyses recovered the identical inter-specific relationships within the species groups but discordant relationships with the BEAST analysis regarding the luii group. The ML and BI analyses recovered a poorly supported (G. catbaensis (G. araneus, G. gezhi)) clade but the BEAST analysis recovered G. catbaensis as the strongly supported (0.99) sister species to the remainder of the luii group species (Figs 3, 4, respectively). All analyses recovered a polyphyletic Goniurosaurus luii with respect to G. kadoorieorum (not shown in the pruned tree of Fig. 4). The mito-nuclear data set of Liang et al. (2018) differed from all the above analyses in that their ML and BI analyses (79/0.99) placed the yingdeensis group as the sister group to a sister lineage comprised of the luii group and lichtenfelderi group (87/1.00).

The AICc scores for the three transition rate models of the SCM analysis were ARD = 34.547134 and SYM and ER = 32.099451. The SCM analysis using either the SYM or ER model suggests that a karst habitat preference is the most probable ancestral condition for Goniurosaurus (57.0% probability), the kuroiwae group (62.7%), the luii group (90.0%), and the yingdeensis group (95.7%; Fig. 4). The probable ancestral condition for the lichtenfelderi group is non-karst (55.4%). The karst habitat preference of G. kwanghua and G. zhoui of the lichtenfelderi group is considered to have evolved independently given that the ancestral condition of the lichtenfelderi group and that of the most recent common ancestor of the sister species G. lichtenfelderii and G. hainanensis was not karst-adapted (Fig. 4).

In both the species and habitat preference PCA analyses, PC1 accounted for 49.1% of the variation in the data set and loaded most heavily for limb length (FLL and HLL), snout length (SL), eye diameter (ED), interorbital distance (IO), head width (HW), and head length (HL). PC2 accounted for an additional 13.3% of the variation and loaded most heavily for body width (BW) and body height (BH) (Figs 5, 6; Table 4).

The PCA analysis of the karst-adapted Goniurosaurus catbaensis, G. huuliensis, and G. luii of the luii group demonstrates that their body shapes greatly overlap in morphospace despite there being several slight, but statistically significant mean differences among them (Fig. 5; Table 3). Additionally, none of the plots of the karst-adapted species overlap with that of the granite stream-adapted species G. lichtenfelderii along the ordination of PC1.

The PCA analysis using habitat preference as the dependent variable among the four species, showed that the karst-adapted and granite-stream-adapted species plot separately as before along taxonomic lines and that collectively, the former have significantly longer axilla-groin lengths (AG); longer, wider, and thicker heads (HL, HW, and CH); longer snouts (SL); longer limbs (FLL and HLL); wider interorbital distances (IO); larger eyes (ED) and larger ear openings (AD) (Fig. 6). Many of these characters—longer head and snout, larger eyes, longer trunk, longer limbs—occur in many other distantly related karst-adapted species of Cyrtodactylus (Grismer et al. 2016a, 2020b; Kaatz et al. 2021; Nielsen and Oliver 2017), indicating that these are convergent adaptions to a karstic life style within and between the gekkotan families.

The ML and BI analyses of Liang et al. (2018) and the BEAST analysis herein (Fig. 4) recovered Goniurosaurus catbaensis as the strongly supported sister species to the remainder of the luii group. Whereas the ML and BI analysis herein, recovered G. catbaensis as the very poorly supported (60/0.51) sister species of the G. araneus plus G. gezhi clade (Fig. 3). Given that three of the five analyses strongly supported the former relationship and two analyses poorly supported the latter, we prefer the placement of G. catbaensis as the sister species to the remainder of the luii group (Fig. 4). Given the very low nodal support of the latter, it essentially renders that portion of the tree a polytomy and as such, does not effectively contradict the strongly supported sister species position of G. catbaensis in the other trees.

Goniurosaurus kadoorieorum of the luii group (represented by only 16S) is nested within G. luii in both the ML and BI analyses, rendering G. luii polyphyletic (Fig. 3). The same relationship was recovered in the 16S phylogeny of Zhu et al. (2020a). This, and the lack of diagnostic characters separating G. kadoorieorum from G. luii (Yang and Chan 2015; Ngo et al. 2016), indicates the two species should be considered conspecific and as such, G. kadoorieorum Yang and Chan, 2015 is relegated here to a junior synonym of G. luii Grismer, Viets, and Boyle, 1999. In all analyses, G. huuliensis is consistently recovered as the sister species to G. luiisensu lato and its species status is not questioned (Figs 3, 4).

Wide-ranging more inclusive studies pertaining to ecosystems management are becoming commonplace in light of climate change and widespread habitat destruction. Such studies reconcile data from a broad range of disciplines in order to address issues that may bear on ecosystems management. Foundational to many of these studies is a basic understanding of species ecology and habitat preference—correlated here with ecomorphology (Cabral et al. 2009; Harfoot et al. 2014). Baseline information on habitat and microhabitat requirements of any species are paramount to understanding how they interact with, and navigate through, their environment (e.g. Grant and Grant 2008; Greene 2005; Losos 2010) and as such, the contextualization of ecosystem management may ultimately turn on these simple points (Meiri 2018; Sinervo et al. 2010).

Integrating the phylogenetic patterns of biodiversity and the morphological adaptations of habitat preference that, in part, underpin species radiations, can fundamentally contribute to conservation management programs (Grismer et al. 2020a, 2021; Erwin 1991; Vane-Wright et al. 1991; Williams et al. 1991; Vázquez and Gittleman 1998; Moritz et al. 2000; Forest et al. 2007; Sgro et al. 2010; Harvey et al. 2011; Rolland et al. 2012; Winter et al. 2012; Shaffer et al. 2015; Beaumont and Wang 2019; Fay et al. 2019; Holderegger et al. 2019)—especially in the karstic regions of northern Vietnam where anthropogenic impact is degrading the habitat and reducing the density of localized populations of Goniurosaurus (Ngo et al. 2019b). Northern Vietnam and many of its offshore islands in the Gulf of Tonkin, harbor large areas of fragmented karstic habitats scattered across their landscapes (Cerrano et al. 2006; Do 2001, 2014; Luo et al. 2016; Ngo et al. 2019a) that are inhabited by an exceptionally large number of endemic plants and animals (Do 2001; Sterling et al. 2006; Clements et al. 2008; Luo et al. 2016; von Oheimb et al. 2017). The obligate restriction of many species to fragmented karstic environments—such as all species of the luii and yingdeensis groups—functionally transforms these environments into habitat islands (Clements et al. 2006, 2008; von Oheimb et al. 2017), which in some cases, bear an unprecedented degree of range-restricted endemism (e.g. Sgro et al 2012; Harvey et al. 2011; Grismer et al. 2018a, 2021).

Unfortunately, Goniurosaurus species are particularly attractive (Fig. 2) and over-harvested for the illegal pet trade (Stuart et al. 2006; Yang and Chan 2015; Ngo et al. 2019b). This is an additional threat to these range-restricted endemics from imperiled karstic environments (Grismer et al. 1997; Orlov et al. 2008; Ziegler et al. 2008; Nakamura et al. 2014; Yang and Chan 2015; Honda and Ota 2017; Zhou et al. 2018; Ngo et al. 2019a; Qi et al. 2020a, b; Zhu et al. 2020a, b). In fact, in areas of China and Vietnam, many populations have suffered huge declines in numbers, or even extirpation at some localities, due to the illegal commercial pet trade (Stuart et al. 2006; Yang and Chan 2015; Ngo et al. 2019b, 2021). We hope that this study will bring more clarity to the plight of this genus and continue to serve ongoing conservation and management programs.