Phylogeny and evolution of habitat preference in Goniurosaurus (Squamata: Eublepharidae) and their correlation with karst and granite-stream-adapted ecomorphologies in species groups from Vietnam

Maximum likelihood (ML) and Bayesian inference (BI) analyses using two mitochondrial (16S and cyt b) and two nuclear (CMOS and RAG1) genes and 103 specimens recovered the first phylogenies of all 23 extant species of Goniurosaurus. The analyses strongly supported the recognition of four monophyletic species groups with identical inter-specific relationships within the kuroiwae, lichtenfelderi, and yingdeensis groups but discordant topologies at some nodes within the luii group. Both analyses recovered a polyphyletic G. luii with respect to G. kadoorieorum, and owing to the lack of diagnostic characters in the latter, it is considered a junior synonym of G. luii. A stochastic character mapping analysis of karst versus non-karst habitat preference suggested that karstic landscapes may have played a major role in the evolution and diversification of Goniurosaurus. A karst habitat preference is marginally supported as the most probable ancestral condition for Goniurosaurus as well as for the kuroiwae, luii, and yingdeensis groups. However, a nonkarst habitat preference is marginally supported as the most probable ancestral condition for the lichtenfelderi group. Multivariate and univariate ecomorphological analyses of the karst-adapted G. catbaensis, G. huuliensis, and G. luii of the luii group and the granite-stream-adapted G. lichtenfelderi of the lichtenfelderi group demonstrated that their markedly statistically different body shapes may be an adaptive response that contributes to habitat partitioning in areas of northern Vietnam where they are nearly sympatric.


Introduction
Eublepharid geckos of the genus Goniurosaurus Barbour, 1908 comprise 23 saxicolous specialists (Uetz et al. 2021) that extend from the Ryukyu Archipelago in Japan, southward through East Asia to northern Vietnam. Goniurosaurus is a well-defined monophyletic group (Grismer 1988) comprised of four monophyletic species groups: the kuroiwae group containing six species endemic to the Ryukyu Archipelago, Japan; the lichtenfelderi group with five species from insular and mainland China and northern Vietnam; the luii group with eight species from northern Vietnam, some of its offshore islands, and southern China; and the yingdeensis group consisting of four species from southern China (Kurita et al. 2008;Nguyen et al. 2009;Nguyen 2011;Wang et al. 2013;Honda and Ota 2017;Liang et al. 2018;Qi et al. 2020aQi et al. , 2020bZhu et al. 2020aZhu et al. , 2020b. Apart from these species, Goniurosaurus sinensis Zhou, Peng, Huo and Yuan, 2019 is likely a junior synonym of another species from Hainan Island, China and not included herein (Qi et al. in progress). Phylogenetic relationships within Goniurosaurus have never been strongly supported nor consistent among different studies (e.g. Wang et al. 2013;Liang et al. 2018;Qi et al. 2020aQi et al. , 2020bZhu et al. 2020aZhu et al. , 2020b. This protracted state of discordance results, in part, from researchers focusing on different species groups as opposed to the entire genus, as well as using different genes or different combinations of genes with varying combinations of ingroup and outgroup species-all variables that bear significantly on tree construction (Wiens 1998;Zwickl et al. 2002;Heath et al. 2008;Wiens and Morrill 2011;Wainwright and Price 2016). The most commonly used genetic markers have been the mitochondrial genes 12S and 16S rRNA and cytochrome b (cyt b). Liang et al. (2018) were the first to address the challenges of properly aligning rRNA (Pyron et al. 2013) and constructed a well-supported mito-nuclear data set using 16S, cyt b, and the nuclear genes oocyte maturation factor MOS (CMOS), and recombination activating 1 (RAG1). Zhu et al. (2020b) also used this mito-nuclear combination, but examined only relationships within the lichtenfelderi group.
In an effort to continue building a more global understanding of the phylogenetic relationships within Goniurosaurus, we expanded the mito-nuclear data set of Liang et al. (2018) to include 103 individuals as opposed to 31 and 23 as opposed to 17 species, which for the first time, includes all extant species of the genus (Table 1). We used this phylogeny in a stochastic character state mapping (SCM) analysis (Revell 2012) of habitat preference to explore the role karstic landscapes may have played in the evolution and diversification of Goniurosaurus and if    habitat preference coevolved with ecomorphology in near sympatric species of the luii and lichtenfelderi groups in Vietnam (Ngo et al. 2021;Figs. 1, 2).
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  to ensure that convergence and stationarity of all para meters 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 ) 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.

Ancestral state reconstruction
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.
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(Thorpe , 1983Turan 1999;Lleonart et al. 2000)-accessable in the R package GroupStruct (available at https://github.com/chankinonn/GroupStruct). The morpho metrics 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).

Phylogenetic relationships
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).

Ancestral state reconstruction
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. lichtenfelderi and G. hainanensis was not karst-adapted (Fig. 4).

Ecomorphology
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. lichtenfelderi 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. 2016aKaatz et al. 2021;Nielsen and Oliver 2017), indicating that these are convergent adaptions to a karstic life style within and between the gekkotan families.

Discussion
Geckos in general are particularly well-adapted to karstic landscapes (see Luu et al. 2016;Grismer et al. 2014 and references therein; Google Scholar search using key words "karst" and "Gekkonidae") and Goniurosaurus is no exception, being that 19 of its 23 species (83%) occupy karstic habitats (Grismer et al. 1994(Grismer et al. , 1999Orlov et al. 2008;Ziegler et al. 2008;Yang and Chan 2015;Honda and Ota 2017;Zhou et al. 2018;Ngo et al. 2019a;Qi et al. 2020aQi et al. , 2020bZhu et al. 2020aZhu et al. , 2020b. It is clear that karstic landscapes have played a significant role in the evolution and diversification of Goniurosaurus being that it is the probable ancestral habitat preference for the genus and three of the four species groups. Even the ancestor of the nonkarst adapted ancestor of the lichtenfelderi group was karst-adapted (Fig. 4). Furthermore, within the species groups, the limited data herein would suggest that the karst-adapted species are specialized, range-restricted endemics (Grismer et al. 1994(Grismer et al. , 1999Orlov et al. 2008;Ziegler et al. 2008;Yang and Chan 2015;Honda and Ota 2017;Zhou et al. 2018;Ngo et al. 2019a;Qi et al. 2020aQi et al. , 2020bZhu et al. 2020aZhu et al. , 2020b. With the exception of G. lichtenfelderi, all the non-karst-adapted species are restricted to islands in the Ryukyu Archipelago (kuroiwae group) or Hainan Island (lichtenfelderi group). It may be that the absence of competition and/or predators in these insular habitats widened the fundamental niches of their ancestors and allowed some species to become more generalized in their habitat preference, which should be tested using new techniques combining phylogenetic history, character evolution, and ecological reconstruction programs.

Systematics of the luii group
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 Table 4. Summary statistics and principal component analysis scores for the mensural characters for Goniurosaurus catbaensis, G. huuliensis, G. luii, and G. lichtenfelderi 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  Boyle, 1999. In all analyses, G. huuliensis is consistently recovered as the sister species to G. luii sensu lato and its species status is not questioned (Figs. 3, 4).

Conservation
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 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;Shaffer et al. 2015   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 2001Do , 2014Luo 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(Clements et al. , 2008von 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. 2018aGrismer et al. , 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. 2019bNgo et al. , 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.