Renewed classification within Goniurosaurus (Squamata: Eublepharidae) uncovers the dual roles of a continental island (Hainan) in species evolution

Accepted Manuscript Renewed classification within Goniurosaurus (Squamata: Eublepharidae) uncovers the dual roles of a continental island (Hainan) in species evolution Bin Liang, Run-Bang Zhou, Yan-Lin Liu, Bei Chen, L. Lee Grismer, Ning Wang PII: DOI: Reference:

S1055-7903(17)30825-4 https://doi.org/10.1016/j.ympev.2018.06.011 YMPEV 6194

To appear in:

Molecular Phylogenetics and Evolution

Received Date: Revised Date: Accepted Date:

18 November 2017 25 April 2018 6 June 2018

Please cite this article as: Liang, B., Zhou, R-B., Liu, Y-L., Chen, B., Lee Grismer, L., Wang, N., Renewed classification within Goniurosaurus (Squamata: Eublepharidae) uncovers the dual roles of a continental island (Hainan) in species evolution, Molecular Phylogenetics and Evolution (2018), doi: https://doi.org/10.1016/j.ympev. 2018.06.011

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Short running title: A Revised view on Evolution of Goniurosaurus

Renewed classification within Goniurosaurus (Squamata: Eublepharidae) uncovers the dual roles of a continental island (Hainan) in species evolution

Bin Lianga,b,*,1, Run-Bang Zhouc,1, Yan-Lin Liud, Bei Chene, L. Lee Grismer f, Ning Wangb,*

a

Forestry Research Institute of Hainan Province, Haikou 571100, Hainan, P. R. China

b

Department of Ecology and Evolutionary Biology, University of Michigan, Ann Arbor, MI

48109, USA c

Forestry Department of Hainan Province, Haikou 570203, China

d

Institute of Forestry Ecology, Environment and Protection, Chinese Academy of Forestry,

Beijing 100091, China e

Department of Transport of Hainan Province, Haikou 570204, China

f

Herpetology Laboratory, Department of Biology, La Sierra University, 4500 Riverwalk

Parkway, Riverside, CA, 92515, USA

1

These authors contributed equally to this work.

*

Corresponding author: Bin Liang, E-mail: [email protected] Ning Wang, E-mail: [email protected]

ABSTRACT Continental islands are often dynamic in regard to the origin and evolution of their biota. Although colonizations from mainland Southeast Asia to Hainan Island have been reported, the role of Hainan Island as a source for continental biota has not been considered. Goniurosaurus is a genus comprised of nocturnal ground geckos. We reexamined the evolutionary history of Goniurosaurus using both molecular phylogenetics and morphological comparisons. All phylogenetic trees recovered G. zhoui as sister to G. hainanensis + G. lichtenfelderi, which together are the sister lineage of G. bawanglingensis. The recovery of this “Hainan clade” contradicts previous classifications that placed G. bawanglingensis within the G. luii group. Moreover, ancestral trait reconstruction revealed that body band number might have decreased two or three times independently within Goniurosaurus from four to three. The divergence between the continental G. luii group and the Hainan clade was estimated at ~34.7 Mya (CI = 22.3–48.6), possibly correlating with the vicariance event between Hainan Island and the mainland. G. lichtenfelderi diverged from G. hainanensis very recently, which might be associated with a historical dispersal event from Hainan Island to Vietnam during glacial periods. Our study improves the understanding of Goniurosaurus systematics and reveals the important role of Hainan Island in bidirectional colonizations.

Keywords: two-way colonization; Goniurosaurus biogeography; body band evolution; Hainan Island

2

1. Introduction The importance of islands in revealing evolutionary processes has been appreciated since the beginning of the field of evolutionary biology (Darwin, 1859). Islands, as ideal systems to study the effects of geographical isolation (e.g., Wang et al., 2016) and long-distance dispersal (Gillespie et al., 2012), have provided significant insights into many aspects of evolution (Patiño et al., 2017), including adaptive radiation (Lamichhaney et al., 2015) and speciation (Emerson, 2008). Unlike remote isolated archipelagos (e.g., Hawaii and the Galapagos) where interchange is commonly one-way (from mainland to island systems) due to the reduced interaction between oceanic islands and continents (Cowie and Holland, 2008), continental islands (such as Madagascar, and the islands in the South China Sea) are expected to provide more dynamic means of evolution and biotic exchange with nearby continental systems (Bellemain and Ricklefs, 2008). For example, islands can potentially act as a source of continental flora and fauna instead of just serving as sinks (e.g., Bellemain and Ricklefs, 2008; Heaney, 2007) especially when they are large (Gillespie et al., 2012; Wikström et al., 2010). As the largest continental island in the northern region of the South China Sea, Hainan Island possesses a variety of habitats that facilitate species diversification. Previously, many new taxa have been documented on this island, including birds (Olsson et al., 1993), amphibians (Chou et al., 2007), and reptiles (Grismer et al., 2002). The continuous discovery of new species has broadened our insights into understanding the evolutionary processes on this island, although the interaction between Hainan Island and its nearby continental landmass has received limited attention in previous studies (but see Chen et al., 2015; Mao et al., 2010; Wang et al., 2016). Goniurosaurus (Eublepharidae), the focus of this study, is a genus comprising 18 primarily nocturnal lizard species associated with granitic and karst landscapes (Uetz et al., 2017;

3

Zhou et al., 2018). These species were separated into four major groups (Fig. 1) based on morphological and molecular analyses (Grismer et al., 2002; Pyron et al., 2013; Wang et al., 2014). In this classification, species from Hainan Island are clustered into different groups based on their different morphological traits, so that G. hainanensis Barbour, 1908 fell within the G. lichtenfelderi group that also contains G. lichtenfelderi (Mocquard, 1897) from north-eastern Vietnam, while G. bawanglingensis Grismer et al., 2002 fell into the G. luii group with the remaining seven species distributed in northern Vietnam and Guangxi Province, China (Grismer et al., 2002). However, because previous phylogenetic analyses based on mitochondrial or nuclear DNA did not include G. bawanglingensis, its relationships to G. hainanensis, G. lichtenfelderi, and G. luii are still unclear and need further examination. Additionally, the newly identified species G. zhoui (Zhou et al., 2018) – also from Hainan Island – has not been assigned to any group. In the current study, we aim to explore the evolutionary history of Goniurosaurus species, with a focus on the species endemic to Hainan Island. The distribution of Goniurosaurus species both on Hainan Island and the nearby continent also provides a model system to explore the function of continental islands as foci of species diversification. Moreover, the variations of body bands on different species of Goniurosaurus also provide an opportunity to study the evolutionary pattern of body bands and their potential contribution to species’ life histories. Using both morphological comparisons and molecular phylogenetic reconstructions, we test 1) whether Goniurosaurus species from Hainan Island form a monophyletic or polyphyletic groups, 2) whether Hainan Island is a sink as to the evolution of Goniurosaurus geckos, and 3) whether different body band numbers are a result of a single or multiple transitions during the evolutionary history of Goniurosaurus geckos.

4

2. Materials and methods 2.1. Taxon sampling Samples of Goniurosaurus species from Hainan Island were collected from November 2015 to April 2017 by Run-Bang Zhou. To conduct molecular analyses, we collected tail tissues from five samples of G. zhoui (voucher numbers: BL-RBZ [Herpetological Collection by Bin Liang and Run-bang Zhou]-001, 004, 006, 007, 008), four samples of G. bawanglingensis (BLRBZ-021–024 from Bawangling National Nature Reserve) and two samples of G. hainanensis (BL-RBZ-041 from Shishan, Haikou city, BL-RBZ-042 from Liulianling, Wanning County). All tissue samples were preserved in 95% alcohol and stored at -20℃ at the Forestry Research Institute of Hainan Province (Supplementary Table S1). 2.2. DNA extraction, sequencing, and alignment Whole genomic DNA was extracted using the Wolact® Blood Genomic DNA Purification Kit (Vicband Life Sciences company, Hong Kong, China). We amplified two commonly used mitochondrial genes (16S rRNA and partial cytochrome b [CYTB]), as well as two nuclear genes (recombination activating gene 1 [RAG1] and oocyte maturation factor MOS [CMOS]) using the primers listed in Supplementary Table S2. We conducted the polymerase chain reaction (PCR) by using 40 μl PCR reaction volumes that include ~200 ng DNA and 20 μl 2 × Taq PCR Premix (TIANGEN Biotech [Beijing] Co., Ltd.). The PCR procedure was performed with an initial denaturation at 94°C for 5 min, 35 cycles of 94°C for 30s, 55°C for 30s and 72°C for 1min, followed by a final extension at 72°C for 10 min. All PCR products were examined for size on 1% agarose gel and purified using PCR purification Kit (QIAGEN). The PCR products were sequenced in both directions using ABI BigDye® Terminator v.3.1 on an

5

ABI Prism® 3100-Avant genetic analyzer (Beijing Genomics Institute, Shenzhen). The sequences were visually compared with the original chromatograms and published data to check for sequencing errors. All sequences were examined using MEGA7 (Kumar et al., 2016) and no premature stop codons were found for coding genes. We followed the IUPAC (International Union of Pure and Applied Chemistry) nucleotide code for heterozygous site. These sequences were combined with published Genbank data from 2–12 Goniurosaurus species as well as three outgroups (two from the Eublepharidae [Holodactylus africanus Boettger, 1893 and Coleonyx mitratus (Peters, 1863)] and one from Gekkonidae [Gekko chinensis Gray, 1842]; Table S1). Each locus was aligned with MAFFT v7 (Katoh and Standley, 2013, maxiterate = 1000). Tcoffee (available at http://tcoffee.crg.cat, Notredame et al., 2000) was also used to align 16S for further comparison (see below). Newly generated sequences were deposited in GenBank (Table S1). 2.3. Phylogenetic analyses The pairwise p-distance between individual sequences of each gene was first examined in MEGA7 (Kumar et al., 2016). After evaluating alternative substitution models for each gene using AICc (Akaike Information Criterion with a correction for small sample sizes) in jModelTest 2 (Darriba et al., 2012, Table 1), we conducted maximum likelihood (ML) analyses for each gene using the GTR + Г + I model (for CYTB) or the GTR + Г model (for the others) in RAxML v8.0.19 (Stamatakis, 2014) with ten iterations. Clade support was assessed with 500 nonparametric bootstrap replicates (i.e., bootstrap support, BS) in RAxML with the same model setting. We also performed ML and bootstrap analyses in IQ-TREE v1.6.1 (Nguyen et al., 2015), which can conduct model choice and tree inference simultaneously. In addition, we generated three concatenated datasets, with the first including all four genes (hereafter All), the second

6

including the two nuclear coding genes (hereafter Nuc), and the third including the two mitochondrial genes (hereafter Mt). It is noteworthy that the concatenated sequences of the same species may not come from the same individual (Table S1). However, because little variation was observed within species, concatenation of different individuals is unlikely to influence inferred species-level relationships. Both ML and Bayesian inference (BI) methods were used to conduct phylogenetic analyses on the concatenated datasets. Specifically, ML trees (partitioned by locus) were built in RAxML with the GTR + Г model and ten search iterations. Clade support was assessed as above. IQ-TREE was also used for partitioned (by locus) ML analyses (Chernomor et al., 2016) and ultrafast bootstrap (1000 replicates, Hoang et al. 2017). BI analyses were conducted using MrBayes v3.2 (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003; Ronquist et al., 2012) with the GTR + Г model (partitioned by locus) and random starting trees. Following the default prior settings that perform well for most data (Ronquist et al., 2012), two independent Markov chain Monte Carlo (MCMC) analyses were performed, each containing one cold and three heated chains, for 100 million generations and sampling every 10,000 generations. The convergence between runs was examined using the average standard deviation of split frequencies and was also assessed in RWTY (Warren et al., 2017, Supplementary Fig. S1). We also examined the effective samples sizes (ESS > 7000) for each parameter estimate in Tracer 1.6.0 (Rambaut and Drummond, 2007). Finally, the first 25% of the trees from each run were discarded as burn-in, and the clade supports were estimated using Bayesian posterior probability (PP). 2.4. Dating analyses

7

Because Goniurosaurus lacks a fossil record, we used secondary calibration points to estimate divergence times among the major clades of Goniurosaurus using *BEAST 1.8.3 (Drummond and Rambaut, 2007). We used the likelihood ratio test (Table 1) implemented in MEGA 7 (Kumar et al., 2016) to evaluate each gene for clock-like behavior, based on the ML scores for each gene tree estimated in RAxML with and without the molecular clock constraints. We randomly chose one individual sequences from each species to conduct dating analyses. However, G. lichtenfelderi AB308458 was excluded from further analyses due to a potential misidentification (see discussion below). After assigning the best substitution and clock models for each gene (Table 1), we conducted BEAST runs with a Yule prior (Heled and Drummond, 2012) and a secondary calibration at the most recent common ancestor (MRCA) of Holodactylus africanus and Coleonyx mitratus following a normal distribution with mean = 79 Mya (Million years ago) and stdev = 9.7 (based on estimated time interval from Gamble et al., 2011). However, preliminary runs failed to converge even after 500 million generations, which could possibly be due to the complexity of the partitioned models, the large proportion of missing data in certain genes, and/or joint prior effects (Brown and Smith, 2017). Given these considerations, we chose to use the most taxon-rich locus, 16S, for the dating analyses. We conducted two replicate analyses to generate posterior results, each with 300 million generations, sampling every 30,000 generations. We assessed convergence, mixing and burn-in using Tracer 1.6.0 (Rambaut and Drummond, 2007) and RWTY (Warren et al., 2017, Fig. S1), which suggested sufficient Markov chain Monte Carlo sampling (Drummond et al., 2006). After discarding the first 25% of the runs as burn-in and combining the output of both runs, we built a consensus tree in Tree Annotator. In addition, we also assessed prior behavior by conducting these analyses without any data (i.e.,

8

sampling from the marginal prior, Supplementary Fig. S2). In total, two duplicates were run for one billion generations each, sampling every 100 thousand generations, allowing for a precise estimation of the breadth of the marginal prior distributions. In addition to BEAST analyses, we also used treePL (Smith and O'Meara, 2012) to estimate divergence time based on an autocorrelated penalized likelihood model, which allows for different rates on different branches and uses a smoothing parameter to penalize rate variation (with higher value representing more clock-like data, Smith and O'Meara, 2012). We set a fixed date (79 Mya, Gamble et al., 2011) to the most recent common ancestor (MRCA) of H. africanus and C. mitratus on a ML tree in treePL. The ML tree used for the dating analysis was built using the concatenated All dataset (partitioned by locus) in RAxML, but only included one individual per species as indicated above. Optimal parameter settings were obtained through a preliminary run with the “prime” option. We also conducted the Random Subsample and Replicate Cross Validation (RSRCV) to determine the best rate smoothing value (i.e., 10 in our case). 2.5. Reconstruction of ancestral body band number We mapped body band numbers (four-band = 1, three-band = 0) for each Goniurosaurus species and reconstructed the ancestral band number in Mesquite 3.03 (Maddison and Maddison, 2015) using maximum likelihood and proportional branch lengths on the ML tree from RAxML including one individual per species. The ancestral states were found to maximize the probability of arriving at the observed states in the terminal taxa under a stochastic model of evolution (Schluter et al., 1997; Pagel, 1999). We applied both the one-rate model (Mk1, in which changes between traits occur at a same rate) and the two-rate model (Asymm. Mk, in which gains and losses occur at different rates) of character evolution. A likelihood ratio test was used to determine the best fitting model. 9

2.6. Reconstruction of ancestral geographical ranges The ancestral geographic ranges were inferred for Goniurosaurus species on the timecalibrated tree from BEAST using the R package BioGeoBEARS (Matzke, 2013). This program implements the likelihood version of biogeographical models, including DECLIKE (dispersalextinction-cladogenesis, Ree and Sanmartín, 2009), DIVALIKE (dispersal-vicariance analysis, Ronquist, 1997), and BAYAREALIKE (Bayesian-based BayArea, Landis et al., 2013), each can also be modified by adding a founder-event speciation parameter (+J). The inclusion of J seems particularly important to insular lineages as it allows one daughter species inheriting the entire ancestral range, whereas the other dispersing to a new geographical area (Matzke, 2014). We fitted all six models on our tree, so that we can test different biogeographical assumptions using a statistical framework of AIC. Based on the potential distribution areas of each species (e.g., Grismer et al., 1994; Grismer et al., 2002; Wang et al., 2014; Zhou et al., 2018), we coded each Goniurosaurus species as present or absent from four general geographical ranges [Hainan (H); Vietnam and Guangxi, China (V); Guangdong, China (D); and the Ryukyu Archipelago of Japan (J); Fig. 1]. We limited the maximum number of areas inferred at a node to four given the close geographic distance among them. Moreover, there is also evidence of historical connection among these landmasses (Teng and Lin, 2004). 3. Results and discussion 3.1. Recovery of the “Hainan clade” updates the phylogeny of Goniurosaurus upon which a new classification can be inferred. The resulting 16S and CYTB matrices from MAFFT alignment contained 17 and 13 species, 530 and 396 base pairs, respectively, while RAG1 and CMOS include eight and nine

10

species, 1002 and 372 base pairs, respectively (Table 1). With more extensive taxon sampling, the 16S gene provided a more complete understanding of the relationships among Goniurosaurus species. However, aligning ribosomal RNA sequences is challenging and can variously influence phylogenetic analyses (Pyron et al., 2013). In order to avoid alignment uncertainty to phylogenetic results, we conducted alignments for 16S using both MAFFT and Tcoffee. Indeed, we noticed some variations in two gap-rich regions (not shown) between alignments, although these differences did not influence the tree topology in our study. Even so, alignment effects on tree topology are important to keep in mind (e.g., Wang et al., 2012), and we believe consistent results from multiple methods (i.e., identical tree topology based on both MAFFT and Tcoffee alignments for 16S, Supplementary Fig. S3) are evidence of a robust phylogenetic tree. In addition, when only considering the position of overlapping species among the four genes, ML trees based on individual loci recovered consistent species-level relationships within Goniurosaurus (Fig. S3). The concordant relationships of most Goniurosaurus species among gene trees lend credence that the concatenation analyses are powerful and robust. The partitioned ML (both from RAxML and IQ-TREE) and BI trees with either All genes or Mt genes recovered identical tree topologies for Goniurosaurus. Besides, trees built with concatenated Nuc genes, although including less taxa, also exhibited overall topological consistency with other trees (Fig. S3). In this case, we focused our phylogenetic discussion based on concatenated All genes. To be specific, the G. kuroiwae group (including five species that are distributed in the Ryukyu Archipelago of Japan, Fig. 1) clustered into a monophyletic group, which formed a sister relationship to all other Goniurosaurus species (Fig. 2A). The species relationships within the G. kuroiwae group, although marginally supported, are consistent with a previous study (Honda et al., 2014).

11

In addition, we also retrieved the monophyletic G. yingdeensis group, which only includes two species and diverged from its sister clade subsequent to the G. kuroiwae group (Wang et al., 2014; Fig. 2). Although our analyses also recovered the G. luii group with high support, it does not include G. bawanglingensis (Fig. 2). In fact, all species on Hainan Island formed a monophyletic group (Hainan clade) that received 100% bootstrap support (BS) for ML and 1.0 posterior probability (PP) for BI analyses (Fig. 2A and Fig. S3). In this clade, G. lichtenfelderi and G. hainanensis have a close relationship (same as morphological classification, Grismer, 2000), which together are sister to G. zhoui. The development of purple brown-like dorsal ground color in G. hainanensis and G. lichtenfelderi (dark purple) and G. zhoui (light purple) seems provide additional evidence for the close relationship among these three species. It is noteworthy that one individual, G. lichtenfelderi AB308458, grouped within G. hainanensis. This individual also exhibited small distances from G. hainanensis individuals (0.008-0.01, Supplementary Table S3). Based on correspondence with Yoshi Kumazawa, who previously submitted the sequence (Jonniaux and Kumazawa, 2008), this sample was identified as “Goniurosaurus lichtenfelderi hainanensis” in their original record, but mislabeled as “Goniurosaurus lichtenfelderi” in their paper as they treated it as a subspecies of G. lichtenfelderi; we thus treated the sample as G. hainanensis and excluded it from other analyses to avoid confusion. Finally, instead of clustering within the G. luii group, G. bawanglingensis fell into the Hainan clade by being sister to all remaining taxa, which conflicts with a previous morphological classification (Grismer et al., 2002). We note that all molecular analyses, either based on concatenated datasets or individual loci, strongly supported (Fig. 2A and Fig. S3) the monophyly of the G. luii group excluding G. bawanglingensis, as well as the monophyly of the Hainan clade with G. bawanglingensis sister

12

to the remaining taxa. In fact, some morphometric evidences are also congruent with our molecular data (e.g., Grismer et al., 2002; Zhou et al., 2018). For example, the maximum SVL (snout–vent length) of a species longer than 107mm was used to differentiate the G. luii group from other groups (Grismer et al., 2002). However, G. bawanglingensis tend to have a shorter SVL (holotype SVL = 104mm), supporting its phylogenetic position out of the G. luii group, which was tentatively implied by a previous study (Grismer et al., 2002). Moreover, Zhou et al. (2018) also mentioned that species in the Hainan clade tend to have a larger number of precloacal pores (PPs, 21–46) in males in comparison with members of the other groups (i.e., 16–29 PPs in the G. luii group; 9–13 PPs in the G. yingdeensis group; and no PP in the G. kuroiwae group; Wang et al., 2014). Goniurosaurus bawanglingensis has the largest number of PPs (Zhou et al., 2018) even within the Hainan clade, again implying its distinctiveness from the G. luii group. Besides, species in the Hainan clade also exhibit dark spots and dark tubercles within the light-colored body bands, which are largely absent in other groups (Wang et al., 2014). Taken as a whole, the recovery of the Hainan clade not only updates our understanding of the phylogeny and classification of Goniurosaurus, but also emphasizes the important role played by Hainan Island in Goniurosaurus diversification. Since the two species of the G. lichtenfelderi group (i.e., G. lichtenfelderi and G. hainanensis) clustered into this clade, we assigned the other lineages to this group as well. However, as G. lichtenfelderi is the only mainland species within the Hainan clade, it could hardly represent the major geographic distribution and origin of this group. Consequently, it seems reasonable to change the group name from “the G. lichtenfelderi group” to “the G. hainanensis group” to emphasize the unique geographic area occupied by this group. However, given that group names usually follow the earliest described species and that group name changes may cause historical confusions, group

13

renaming is open for future discussion and we will still follow the original group name in our paper. 3.2. New phylogeny provides insight on body band evolution in Goniurosaurus Body bands are noticeable traits in most eublepharid geckos, and the number and style of bands have been used as important characteristics for species delimitation (Grismer, 1988; Domingos et al., 2014). However, the evolution of body bands has never been explored in this group by previous studies. The relative consistency of bright-colored body bands in Goniurosaurus species offers a good opportunity to examine the evolutionary pattern of body band number between the nuchal loop and the caudal constriction of different species. In Goniurosaurus, only four species have three bands while all other species have four bands. Among the four three-banded species, two (G. hainanensis and G. lichtenfelderi) are in the the G. lichtenfelderi group, and two (G. splendens Nakamura and Uéno, 1959 and some G. toyamai Grismer et al., 1994) are in the G. kuroiwae group. In order to understand the pattern of body band evolution, we inferred the ancestral body band number on the ML tree of the concatenated All dataset. There was no significant difference between the Mk1 (LnL = 9.14) and the Assym. Mk model (LnL = 8.26, p > 0.1), both obtained identical results as to the ancestral number of body bands. Although we are aware of the limitations of this analysis in regard to the lack of outgroup information due to the complexity in quantifying various band patterns, the results indicated that the ancestor of Goniurosaurus likely possessed four body bands, with body bands decreasing perhaps twice – at the MRCA of the G. kuroiwae group (the probability of having three bands is 0.51) and the MRCA of G. hainanensis and G. lichtenfelderi (Fig. 2B). In this case, the expression of four body bands was also regained once in the G. kuroiwae group. However, because the probability of having three bands is not significantly higher than having four bands 14

either at the MRCA of the G. kuroiwae group (0.51) or at the MRCA of G. toyamai and G. orientalis (0.53), it is also possible that the ancestors of these two clades still possessed four bands, and the three-band trait evolved twice independently in G. splendens and G. toyamai. The transverse bands developed in many geckos seem to play an important role in a species’ life history, with major functions including mimicry and/or camouflage as a defense mechanism (Autumn and Han, 1989; Heatwole, 1968). It is likely that the disruptive coloration from bands can break up the body outline so that the gecko shape may not be recognized (Badger and Netherton, 2002). Additionally, the bright orange, yellow, or white body bands in adult and especially juvenile geckos (e.g., Grismer et al., 1999; Wang et al., 2014; Zhou et al., 2018) may be aposematic to avoid predators. For example, the juvenile Tibetan frog-eyed Gecko (Teratoscincus roborowskii Bedriaga, 1906) can exhibit Batesian mimicry of a scorpion in the genus Mesobuthus (Buthidae, Autumn and Han, 1989). In addition, the color of body bands tends to become dim or blended with the body background color with ages (e.g., G. zhelongi, Fig. 2a-b in Wang et al., 2014; G. hainanensis and G. bawanglingensis, Fig. 1D-E in Zhou et al., 2018), similar to some populations of Coleonyx variegatus (Baird, 1858) and Coleonyx reticulatus Davis and Dixon, 1958 (Grismer, 1988; Stebbins and McGinnis, 2012). These mottled spots and bands on the back of adult geckos may thus provide camouflage against their surroundings (Bauer and Russell, 1989; Heatwole, 1968). Given the multiple potential roles and the significant functional changes of body bands within the lifetime of a Goniurosaurus gecko (from advertising as a juvenile to disruptive/deceptive coloration as an adult), our exploration of body band evolution, although with some limitations, may provide novel insight into trait evolution in eublepharid geckos and offer an avenue for future study. 3.3. Biogeographic evolution implied two roles of Hainan Island

15

We examined the biogeographic evolution of Goniurosaurus geckos based on molecular dating and ancestral range reconstruction. The secondary calibration we used in the current study is located at the crown clade outside Goniurosaurus, which might cause under estimation of evolutionary rate (Ho et al., 2008). However, the divergence times among clades recovered by our BEAST and treePL analyses are relative consistent with previous studies. For example, the divergence between G. araneus and G. luii is approximately 12.7 Mya in BEAST (CI: 7.4–19.5, Fig. 3) and 15.1 Mya in treePL (Supplementary Fig. S4), falling into the range (~12.5 +/- 7 Mya) of a previous estimation based on nuclear data and less taxon sampling (Gamble et al., 2011). Moreover, the crown age of the G. kuroiwae group (i.e., 14.4 Mya [CI: 7.8–22.4] from BEAST and 15.5 Mya from treePL) are also consistent with those in Honda et al. (2014). Since the results from BEAST incorporate rate variation and dating uncertainty, we follow its results for discussion. With the confidence of dating consistency between the current and previous studies, the diversification of Goniurosaurus possibly started in Paleocene/Eocene after the CretaceousPaleogene (K-Pg) boundary or even earlier (~61.8 Mya [CI: 42.3–82.1], Fig. 3). According to AIC, the best fitting model for biogeographical range evolution is DIVALIKE+J, although DIVALIKE model also exhibits high likelihood (0.31, Supplementary Table S4). However, range reconstructions based on both models are almost identical (Fig. S4), thus, we follow the results from the best model in the discussion. Due to the controversial relationships among outgroup species (e.g., Honda et al., 2014; Pyron et al., 2013), we did not include distribution ranges from outgroups in order to avoid error inference. However, the lack of outgroup information might also cause uncertainty of range reconstruction at the base of Goniurosaurus groups. Being aware of this limitation, the ancestor of Goniurosaurus could be distributed among the landmasses including all assigned four areas (HVDJ, Fig. 3). This seems

16

probable given that all these areas were connected before Cenozoic and the East Asian margin comprised the region where the Ryukyu Archipelago is currently located (Teng and Lin, 2004). As a consequence of the subduction of the Philippine Sea Plate under the continental Eurasian plate ~50 Mya (Hall, 2002; Ota, 1998), the Ryukyu Archipelago has separated from the eastern Eurasian continent and undergone extensive changes in configuration, which could have possibly driven the divergence between the G. kuroiwae group and its continental sister clade. This vicariance-divergence hypothesis was also suggested by previous studies (e.g., Honda et al., 2014). Besides, speciation within the G. kuroiwae group could be driven by both dispersal and vicariance among small islands (Honda et al., 2014). The G. yingdeensis group and the remaining Goniurosaurs diverged ~47.4 Mya (CI: 31.9–65.3), while the G. lichtenfelderi group and the G. luii group diverged ~34.7 Mya (CI = 22.3–48.6) during early Eocene to late Oligocene (Fig. 3). The ancestor of the G. lichtenfelderi group and the G. luii group was distributed in both Hainan and Vietnam/Guangxi, China (HV, Fig. 3) according to the ancestral range reconstruction. Based on biogeographical and palaeomagnetic evidence, Hainan Island was thought to be adjacent to northern Vietnam and the Guangxi province, China in the Eocene (Zhu, 2016). Because of the indentation of the Indian Plate into Eurasia, Hainan Island drifted southeast to the present location during the large-scale left-lateral motion along the Red River Fault (Liu and Morinaga, 1999; Zhu, 2016). Thus, the separation of Hainan Island from the mainland possibly promoted the isolation of the G. lichtenfelderi group from its continental ancestral group, consistent with the results from the ancestral range reconstruction that also suggested a vicariance between these two groups (Fig. 3). Within the G. lichtenfelderi group, G. bawanglingensis diverged from its sister lineages ~ 15.8 Mya (CI: 9.2–23.2). Goniurosaurus bawanglingensis is generally found in the vicinity of

17

granitic in upland mountain forests (Grismer et al., 2002), while G. zhoui and the majority G. hainanensis are lowland dwellers based on our observation (i.e., G. zhoui is a typical karst dweller with the range of altitudes between 220 and 300 meters and the major population of G. hainanensis inhabits igneous habitat or mossy forests in lowland Hainan and up to some of its highest mountains [i.e. the type locality]). Thus, habitat type might have driven the divergence between them (Zhou et al., 2018). The divergence between G. hainanensis and G. zhoui is estimated to be ~10.9 Mya (CI: 6.0–16.9) in the middle to late Miocene. According to recent geochronology and geochemistry data, their divergence may be associated with the volcanism and formation of large areas of igneous rocks in the northern part of Hainan Island (i.e., volcanism in the Leiqiong area [including the Leizhou Peninsula and the northern part of Hainan Island] may have occurred in the late Oligocene and gradually increased in frequency toward the Miocene and Pliocene, Ho et al., 2000). Such habitat changes might also promote the body color evolution in G. hainanensis. For example, major populations of G. hainanensis inhabit the igneous area in the northern part of Hainan Island, where the igneous rocks cover more than 4000 km2. The unique dark purple brown dorsal ground color in G. hainanensis matches the background color of the igneous rocks, which might play a camouflage function against predation (e.g., Marshall et al., 2015). After the divergence between G. hainanensis and G. zhoui, the ancestor of G. hainanensis may have spread to Vietnam near the mouth area of the Red River through the exposed continental shelf between Hainan Island and the mainland during glacial periods (Figs. 1 and 3), a dispersal pathway that has been widely reported in plants (Zhu, 2016) and animals (e.g., Chen et al., 2015; Lin et al., 2010). This seems consistent with the fact that a population of G. lichtenfelderi is found in Kuinong Chao (Iles de Norway, the type locality of G. lichtenfelderi), a

18

small land bridge island located near the mouth of the Red River where its continental population is distributed (Grismer, 2000; Orlov et al. 2008). In this case, this dispersal and subsequent vicariance (during interglacial periods) possibly promoted the origin of G. lichtenfelderi through a reverse colonization from Hainan Island to the continental mainland around 3.5 Mya (Fig. 3), an evolutionary process indicated from our ancestral range reconstruction and being emphasized only recently (Bellemain and Ricklefs, 2008; Patiño et al., 2017). In fact, the idea of ‘unidirectional colonization’ from continents (or larger areas) to islands (or smaller areas) has long been a prevailing paradigm of ecology/biogeography (Wilson, 1961), and has been adopted by several theories including the equilibrium theory of island biogeography (MacArthur and Wilson, 1967) and the metapopulation theory (Hanski, 1999). In addition, one-way colonization has further been reinforced by the diversity–invasibility theory, which suggests that less diverse communities (i.e., islands) should be more susceptible to invasion by exotic species (Elton, 1958). However, the “reverse colonization” process (from islands to continents) is thought to be possible as island species are the direct descendants of successful colonists and they may preserve their dispersal ability (Fridley et al., 2007). In fact, many studies have revealed several reverse colonization events, highlighting its prevalence along biogeographical evolution for many groups (e.g., Bellemain and Ricklefs, 2008; Heaney, 2007; Raxworthy et al., 2002; Sturge et al., 2009). In the current study, we clearly identify bidirectional routes of species colonization between Hainan Island and the nearby continent (Fig. 3). Although previous studies emphasized the role of Hainan Island as a sink for species diversification (e.g., Wang et al., 2016), we identify its function as a source and provide an example of “reverse colonization” in the South China Sea for the first time. In fact, the direction and probability of colonization between islands

19

and continents could be influenced by many factors (Bellemain and Richlefs, 2008), among which, habitat conditions and the ability to adapt to new environments could be essential. It has been found that competitive interactions can be overruled by habitat suitability at larger scales (Fridley et al., 2007). Although closely related to G. hainanensis, G. lichtenfelderi typically inhabits the granite bed of valley streams (Orlov et al., 2008), a different habitat from the igneous rocks preferred by G. hainanensis. The utility of different habitats and the ability to switch habitat type might promote the diversification of Goniurosaurus (Zhou et al., 2018) and facilitate a two-way colonization between island and continent (Fig. 1). Because the large limestone karst formations with tropical forest in southern East Asia provide diverse habitats for unique flora and fauna, the ability to adapt to specific substrates is crucially important, and may accelerate the speciation process in geologic and topographic heterogeneous areas (Gamble et al., 2012). 4. Conclusion In this study, we uncover the evolutionary history of Goniurosaurus using both phylogenetic and morphological data. This improved understanding of Goniurosaurus phylogeny reveals the important role of Hainan Island as a source of continental biota. Goniurosaurus zhoui is sister to G. hainanensis + G. lichtenfelderi, which together are sister to G. bawanglingensis. Such a finding contradicts previous classifications that placed G. bawanglingensis within the G. luii group. After reconstructing the ancestral body band numbers, we found that the three-bodyband trait might have independently evolved two or three times within Goniurosaurus. Biogeographical analyses showed that the divergence between the continental G. luii group and the G. lichtenfelderi group occurred approximately 35 Mya, possibly correlating with the vicariance event between Hainan Island and the mainland. On the other hand, the island species 20

G. hainanensis diverged from G. lichtenfelderi (a mainland species) very recently and the latter might be associated with a historical dispersal event from Hainan Island to Vietnam during glacial periods. The recovery of this biogeographic evolutionary pattern reveals the dual roles of Hainan Island acting as both a sink and a source of novel species diversity with respect to the continental mainland in the South China Sea. Acknowledgments This study was supported by Hainan Key Research and Development Program (No. ZDYF2018143 to BL) and National Natural Science Foundation of China (No. 31301894 to BL and No. 31360510 to NW). We thank the constructive comments and edits from the two anonymous reviewers and editor James A. Schulte II. We thank Rebecca T. Kimball in University of Florida for constructive advice on biogeographical analysis. Dr. Qing Song assisted in making maps. Joseph Brown, Greg Stull, Stephen Smith, Simon Uribe-Convers, and Oscar Vargas have greatly improved the readability of this manuscript. Special thanks to Joseph Brown for the consulting with BEAST analyses. We are grateful to the collaboration of the Bawangling National Nature Reserve in sample collection. Appendix A. Supplementary material Supplementary tables and Figures associated with this article can be found, in the online version, at XXX. Table S1. GenBank accession number for each sequence used in this study. Table S2. Information of primers used in this study. Table S3. P-distance between individual 16S sequences in the Hainan clade.

21

Table S4. Model test for ancestral range reconstruction in BioGeoBEARS. Fig. S1. Assessing convergence between Bayesian runs from MrBayes and BEAST using RWTY. Fig. S2. Pairwise comparison between the marginal prior (blue) and posterior (black) distributions for the age of certain crown clades within Goniurosaurus. Fig. S3. Phylogenetic analyses based on individual genes and concatenated datasets. Fig. S4. Dating analyses from treePL with the ancestral range reconstruction under the DIVALIKE model. References Autumn, K., Han, B., 1989. Mimicry of scorpions by juvenile lizards, Teratoscincus roborowskii (Gekkonidae). Chin. Herpetol. Res., 2, 60–64. Badger, D., Netherton, J., 2002. Lizards: A natural history of some uncommon creatures—Extraordinary chameleons, iguanas, geckos and more. Minnesota: Voyageur Press Inc., pp. 78–79. Baird, S.F., 1858. Description of new genera and species of North American lizards in the museum of the Smithsonian Institution. Proc. Acad. Nat. Sci. Philadelphia 10, 253–256. Barbour, T., 1908. Some new reptiles and amphibians. Bull. Mus. Comp. Zool. 51, 315–325. Bauer, A.M., Russell, A.P., 1989. A systematic review of the genus Uroplatus (Reptilia: Gekkonidae), with comments on its biology. J. Nat. Hist. 23, 169–203. Bedriaga, J.V., 1906. Verzeichnis der von der Central-Asiatischen Expedition unter Stabs-Kapitän W. Roborowski in den Jahren 1893-1895 gesammelten Reptilien. Annu. Mus. Zool. Acad. Impér. Sci. St.-Petersbourg, 10, 159–200 [1905]. Bellemain, E., Ricklefs, R.E., 2008. Are islands the end of the colonization road? Trends Ecol. Evol. 23, 536–537. Boettger, O., 1893. Übersicht der von Prof. C. Keller anlässlich der Ruspoli'schen Expedition nach den Somaliländern gesammelten Reptilien und Batrachier. Zool. Anz. 16, 113–119.

22

Chen, D., Chang, J., Li, S.H., Liu, Y., Liang, W., Zhou, F., Yao, C.T., Zhang, Z.W., 2015. Was the exposed continental shelf a long-distance colonization route in the ice age? The southeast Asia origin of Hainan and Taiwan partridges. Mol. Phylogenet. Evol. 83, 167–173. Chernomor, O., von Haeseler, A., Minh, B.Q., 2016. Terrace aware data structure for phylogenomic inference from supermatrices. Syst. Biol. 65, 997–1008. Chou, W.H., Lau, M.W.N., Chan, B.P.L., 2007. New treefrog of the genus Rhacophorus (Anura: Rhacophoridae) from Hainan Island, China. Raffles B. Zool. 55, 157–165. Cowie, R.H., Holland, B.S., 2008. Molecular biogeography and diversification of the endemic terrestrial fauna of the Hawaiian Islands. Philos. T. R. Soc. B. 363, 3363–3376. Darriba D., Taboada G.L., Doallo R., Posada D. 2012. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods 9, 772. Darwin, C., 1859. On the origin of species. John Murray, London. Davis, W.B., Dixon, J.R., 1958. A new Coleonyx from Texas. Proc. Biol. Soc. Washington 71, 149–152. Domingos, F.M.C.B., Bosque, R.J., Cassimiro, J., Colli, G.R., Rodrigues, M.T., Santos, M.G., Beheregaray, L.B., 2014. Out of the deep: Cryptic speciation in a Neotropical gecko (Squamata, Phyllodactylidae) revealed by species delimitation methods. Mol. Phylogenet. Evol. 80, 113–124. Drummond, A.J., Ho, S.Y.W., Phillips, M.J., Rambaut, A., 2006. Relaxed phylogenetics and dating with confidence. PLoS Biol. 4, e88. Drummond, A.J., Rambaut, A., 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 7, 214. Elton, C.S., 1958. The ecology of invasions by animals and plants. Methuen. Emerson, B.C., 2008. Speciation on islands: what are we learning? Biol. J. Linn. Soc. 95, 47–52. Fridley, J.D., Stachowicz, J.J., Naeem, S., Sax, D.F., Seabloom, E.W., Smith, M.D., Stohlgren, T.J., Tilman, D., Von Holle, B., 2007. The invasion paradox: Reconciling pattern and process in species invasions. Ecology 88, 3–17. Gamble, T., Bauer, A.M., Colli, G.R., Greenbaum, E., Jackman, T.R., Vitt, L.J., Simons, A.M., 2011. Coming to America: multiple origins of New World geckos. J. Evolution. Biol. 24, 231–244.

23

Gamble, T., Greenbaum, E., Jackman, T.R., Russell, A.P., Bauer, A.M., 2012. Repeated origin and loss of adhesive toepads in geckos. PLoS One 7: e39429. Gillespie, R.G., Baldwin, B.G., Waters, J.M., Fraser, C.I., Nikula, R., Roderick, G.K., 2012. Long-distance dispersal: a framework for hypothesis testing. Trends Ecol. Evol. 27, 47–56. Gray, J.E., 1842. Description of some new species of reptiles, chiefly from the British Museum collection. Zool. Misc. 2, 57–59. Grismer, L.L., 1988. Phylogeny, taxonomy, classification, and biogeography of eublepharid geckos. In: Phylogenetic Relationships of the Lizard Families. Essays Commemorating Charles L. Camp, Estes, R., Pregill, G., Eds, Stanford University Press, Stanford, pp. 369–469. Grismer, L.L., 2000. Goniurosaurus murphyi Orlov and Darevsky: A junior synonym of Goniurosaurus lichtenfelderi Mocquard. J. Herpetol. 34, 486–488. Grismer, L.L., Shi, H.T., Orlov, N.L., Ananjeva, N.B., 2002. A new species of Goniurosaurus (Squamata: Eublepharidae) from Hainan Island, China. J. Herpetol. 36, 217–224. Grismer, L.L., Ota, H., Tanaka, S., 1994. Phylogeny, classification, and biogeography of Goniurosaurus Kuroiwae (Squamata, Eublepharidae) from the Ryukyu Archipelago, Japan, with description of a new subspecies. Zool. Sci. 11, 319–335. Grismer, L.L., Viets, B.E., Boyle, L.J., 1999. Two new continental species of Goniurosaurus (Squamata : Eublepharidae) with a phylogeny and evolutionary classification of the genus. J. Herpetol. 33, 382–393. Hall, R., 2002. Cenozoic geological and plate tectonic evolution of SE Asia and the SW Pacific: computer-based reconstructions, model and animations. J. Asian Earth Sci. 20, 353–431. Hanski, I., 1999. Metapopulation Ecology. Oxford University Press. Heaney, L.R., 2007. Is a new paradigm emerging for oceanic island biogeography? J. Biogeogr. 34, 753–757. Heatwole, H., 1968. Relationship of escape behavior and camouflage in anoline lizards. Copeia 1968, 109-113. Heled, J., Drummond, A.J., 2012. Calibrated tree priors for relaxed phylogenetics and divergence time estimation. Syst. Biol. 61, 138–149. Ho, K.S., Chen, J.C., Juang, W.S., 2000. Geochronology and geochemistry of late Cenozoic basalts from the Leiqiong area, southern China. J. Asian Earth Sci. 18, 307–324.

24

Ho, S.Y., Saarma, U., Barnett, R., Haile, J., Shapiro, B., 2008. The effect of inappropriate calibration: three case studies in molecular ecology. PLoS One 3, e1615. Hoang, D.T., Chernomor, O., von Haeseler, A., Minh, B.Q., Le, S.V., 2017. UFBoot2: Improving the ultrafast bootstrap approximation. Mol. Biol. Evol. 35, 518–522. Honda, M., Kurita, T., Toda, M., Ota, H., 2014. Phylogenetic relationships, genetic divergence, historical biogeography and conservation of an endangered gecko, Goniurosaurus kuroiwae (Squamata: Eublepharidae), from the central Ryukyus, Japan. Zool. Sci. 31, 309–320. Huelsenbeck, J.P., Ronquist, F., 2001. MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics 17, 754–755. Jonniaux, P., Kumazawa, Y., 2008. Molecular phylogenetic and dating analyses using mitochondrial DNA sequences of eyelid geckos (Squamata: Eublepharidae). Gene 407, 105–115. Katoh, K., Standley, D.M., 2013. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. Kumar, S., Stecher, G., Tamura, K., 2016. MEGA7: Molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874. Lamichhaney, S., Berglund, J., Almen, M.S., Maqbool, K., Grabherr, M., Martinez-Barrio, A., Promerova, M., Rubin, C.J., Wang, C., Zamani, N., Grant, B.R., Grant, P.R., Webster, M.T., Andersson, L., 2015. Evolution of Darwin's finches and their beaks revealed by genome sequencing. Nature 518, 371–375. Landis, M.J., Matzke, N.J., Moore, B.R., Huelsenbeck, J.P., 2013. Bayesian analysis of biogeography when the number of areas is large. Syst. Biol. 62, 789–804. Lin, L.H., Ji, X., Diong, C.H., Du, Y., Lin, C.X., 2010. Phylogeography and population structure of the Reevese's Butterfly Lizard (Leiolepis reevesii) inferred from mitochondrial DNA sequences. Mol. Phylogenet. Evol. 56, 601–607. Liu, Y.Y., Morinaga, H., 1999. Cretaceous palaeomagnetic results from Hainan Island in south China supporting the extrusion model of Southeast Asia. Tectonophysics 301, 133–144. MacArthur, R.H., Wilson, E.O., 1967. The Theory of Island Biogeography. Princeton University Press. Maddison, W.P., Maddison, D.R., 2015. Mesquite: a modular system for evolutionary analysis. Version 3.03, http://mesquiteproject.org.

25

Maki, M., 1931. A new Banded Gecko, Eublepharis orientalis sp. nov. from Riu Kyu. Annot. Zool. Japon. 13, 9–11. Mao, X.G., Zhu, G.J., Zhang, S.Y., Rossiter, S.J., 2010. Pleistocene climatic cycling drives intra-specific diversification in the intermediate horseshoe bat (Rhinolophus affinis) in Southern China. Mol. Ecol. 19, 2754–2769. Marshall, K.L.A., Philpot, K.E., Stevens, M., 2015. Conspicuous male coloration impairs survival against avian predators in Aegean wall lizards, Podarcis erhardii. Ecol. Evol. 5, 4115–4131. Matzke, N.J., 2013. BioGeoBEARS: BioGeography with Bayesian (and likelihood) evolutionary analysis in R Scripts. R package (updated May 2016). Matzke, N.J., 2014. Model selection in historical biogeography reveals that founder-event speciation is a crucial process in island clades. Syst. Biol. 61, 951–970. Mocquard, F., 1897. Notes herpétologiques. Bulletin du Muséum National d'Histoire Naturelle. Paris (Ser. 1) 3, 211–217. Nakamura, K., Uéno, S.I., 1959. The geckos found in the limestone caves of the Ryukyu Islands. Memoirs of the College of Science, University of Kyoto 26, 45–52. Namiye, M., 1912. The geckos from the Okinawa Islands. Zool. Mag. 24, 442–445. Nguyen, L.T., Schmidt, H.A., von Haeseler, A., Minh, B.Q., 2015. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274. Notredame, C., Higgins, D.G., Heringa, J., 2000. T-Coffee: A novel method for fast and accurate multiple sequence alignment. J. Mol. Biol. 302, 205–217. Orlov, N.L., Ryabov, S.A., Tao, N.T., Truong, N.Q., Cuc, H.T., 2008. A new species of Goniurosaurus (Sauria: Gekkota: Eublepharidae) from North Vietnam. Russ. J. Herpetol. 15, 229–244. Okada, Y., 1936. A new Cave-Gecko, Gymnodactylus yamashinae from Kumejima, Okinawa group. Proc. Imp. Acad. 12, 53–54. Olsson, U., Alstrom, P., Colston, P.R., 1993. A new species of Phylloscopus warbler from Hainan Island, China. Ibis 135, 3–7. Ota, H., 1998. Geographic patterns of endemism and speciation in amphibians and reptiles of the Ryukyu archipelago, Japan, with special reference to their paleogeographical implications. Res. Popul. Ecol. 40, 189–204.

26

Pagel, M., 1999. The maximum likelihood approach to reconstructing ancestral character states of discrete characters on phylogenies. Syst. Biol. 48, 612–622. Patino, J., Whittaker, R.J., Borges, P.A.V., Fernandez-Palacios, J.M., Ah-Peng, C., Araujo, M.B., Avila, S.P., Cardoso, P., Cornuault, J., de Boer, E.J., de Nascimento, L., Gil, A., Gonzalez-Castro, A., Gruner, D.S., Heleno, R., Hortal, J., Illera, J.C., Kaiser-Bunbury, C.N., Matthews, T.J., Papadopoulou, A., Pettorelli, N., Price, J.P., Santos, A.M.C., Steinbauer, M.J., Triantis, K.A., Valente, L., Vargas, P., Weigelt, P., Emerson, B.C., 2017. A roadmap for island biology: 50 fundamental questions after 50 years of The Theory of Island Biogeography. J. Biogeogr. 44, 963–983. Peters, W.C.H., 1863. Über einen neuen Gecko, Brachydactylus mitratus aus Costa Rica. Monatsber. königl. Akad. Wiss. Berlin. 19, 41–44. Pyron, R.A., Burbrink, F.T., Wiens, J.J., 2013. A phylogeny and revised classification of Squamata, including 4161 species of lizards and snakes. BMC Evol. Biol. 13, 93. Rambaut, A., Drummond, A.J., 2007. Tracer version 1.6. Available at: http://tree.bio.ed.ac.uk/software/tracer/. Raxworthy, C.J., Forstner, M.R.J., Nussbaum, R.A., 2002. Chameleon radiation by oceanic dispersal. Nature 415, 784–787. Ree, R.H., Sanmartın, I., 2009. Prospects and challenges for parametric models in historical biogeographical inference. J. Biogeogr. 36, 1211–1220. Ronquist, F., 1997. Dispersal–vicariance analysis: a new approach to the quantification of historical biogeography. Syst. Biol. 46, 195–203. Ronquist, F., Huelsenbeck, J.P., 2003. Mrbayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574. Ronquist, F., Teslenko, M., van der Mark, P., Ayres, D.L., Darling, A., Hohna, S., Larget, B., Liu, L., Suchard, M.A., Huelsenbeck, J.P., 2012. MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 61, 539–542. Schluter, D., Price, T., Mooers, A. O., Ludwig, D., 1997. Likelihood of ancestor states in adaptive radiation. Evolution 51, 1699–1711. Smith, S.A., O'Meara, B.C., 2012. treePL: divergence time estimation using penalized likelihood for large phylogenies. Bioinformatics 28, 2689–2690.

27

Stamatakis, A., 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313. Stebbins, R.C., McGinnis, S.M., 2012. Field guide to amphibians and reptiles of California: revised edition. Vol. 103. University of California Press. pp. 247. Sturge, R.J., Jacobsen, F., Rosensteel, B.B., Neale, R.J., Omland, K.E., 2009. Colonization of south America from Caribbean Islands confirmed by molecular phylogeny with increased taxon sampling. Condor 111, 575– 579. Teng, L.S., Lin, A.T., 2004. Cenozoic tectonics of the China continental margin: insights from Taiwan. Geol. Soc. London Spec. Publ. 226, 313–332. Uetz, P., Freed, P., Hošek, J., The Reptile Database, http://www.reptile-database.org/ (accessed February 2017). Wang, N., Braun, E.L., Kimball, R.T., 2012. Testing hypotheses about the sister group of the Passeriformes using an independent 30-Locus data set. Mol. Biol. Evol. 29, 737–750. Wang, N., Liang, B., Wang, J., Yeh, C.F., Liu, Y., Liu, Y., Liang, W., Yao, C.T., Li, S.H., 2016. Incipient speciation with gene flow on a continental island: Species delimitation of the Hainan Hwamei (Leucodioptron canorum owstoni, Passeriformes, Aves). Mol. Phylogenet. Evol. 102, 62–73. Wang, Y.Y., Jin, M.J., Li, Y.L., Grismer, L.L., 2014. Description of a new species of Goniurosaurus (Squamata: Eublepharidae) from the Guangdong Province, China, based on molecular and morphological data. Herpetologica 70, 309–322. Wang, Y.Y., Yang, J.H., Cui, R.F., 2010. A new species of Goniurosaurus (Squamata: Eublepharidae) from Yingde, Guangdong Province, China. Herpetologica 66, 229–240. Wang, Y.Y., Yang, J.H., Grismer, L.L., 2013. A new species of Goniurosaurus (Squamata: Eublepharidae) from Libo, Guizhou Province, China. Herpetologica 69, 214–226. Warren, D.L., Geneva, A.J., Lanfear, R., 2017. RWTY (R We There Yet): an R package for examining convergence of Bayesian phylogenetic analyses. Mol. Biol. Evol. 34, 1016–1020. Wilson, E.O., 1961. The nature of the taxon cycle in the Melanesian ant fauna. Am. Nat. 95, 169–193. Yang, J.H., Chan, B.P.L., 2015. Two new species of the genus Goniurosaurus (Squamata: Sauria: Eublepharidae) from southern China. Zootaxa 3980, 67–80.

28

Zhou, R.B., Wang, N., Chen, B., Liang, B., 2018. Morphological evidence uncovers a new species of Goniurosaurus (Squamata: Eublepharidae) from the Hainan Island, China. Zootaxa 4369, 281–291. Zhu, H., 2016. Biogeographical evidences help revealing the origin of Hainan Island. PLoS One 11: e0151941. Ziegler, T., Truong, N.Q., Schmitz, A., Stenke, R., Rosler, H., 2008. A new species of Goniurosaurus from Cat Ba Island, Hai Phong, northern Vietnam (Squamata: Eublepharidae). Zootaxa 1771, 16–30.

29

Table 1 Sequence information and results from model tests.

16S CYTB CMOS RAG1

No.

Lengths jModelTest

34 22 18 15

AICc GTR+G GTR+I+G HKY+G HKY+G

530bp 396bp 372bp 1002bp

Molecular clock test With -2669.38 -2610.303 -1087.203 -3066.389

Without -2135.563 -2599.447 -764.259 -2485.52

p 6.70E-182 0.99 1.90E-115 1.70E-228

Clock model UCLN Strict clock UCLN UCLN

No. – Number of sequences in each gene locus. AICc – Akaike information criterion corrected for small sample size. UCLN – uncorrelated lognormal clock model. Molecular clock tests were conducted in MEGA7 by comparing the ML scores for each gene tree estimated in RAxML with and without the molecular clock constraints.

30

Figure Legends: Fig. 1. Distribution map of Goniurosaurus species. The type locality of each species is labeled by a species specific number. Species from the same group are indicated by the same color panel. The G. kuroiwae group (blue) includes five species (1, G. kuroiwae Namiye, 1912 in Okinawajima Island; 2, G. splendens Nakamura and Uéno, 1959 in Tokunoshima Island; 3, G. toyamai Grismer et al., 1994 in Iheyajima Island; 4, G. orientalis Maki, 1931 in Iejima Island; 5, G. yamashinae Okada, 1936 in Kumejima Island); The G. yingdeensis group (orange) includes two species that distribute in Yingde City, Guangdong Province, China (6, G. yingdeensis Wang et al., 2010; 7, G. zhelongi Wang et al., 2014); The G. luii group (grey) includes seven species between northern Vietnam and Guangxi (8, G. liboensis Wang et al., 2013 in the Maolan National Nature Reserve, Libo County, Guizhou province; 9, G. luii Grismer et al., 1999 in Pingxiang, Guangxi Province; 10, G. araneus Grismer et al., 1999 in Cao Bang, Vietnam; 11, G. huuliensis Orlov et al., 2008 in Huu Lien Nature Reserve, Vietnam; 12, G. catbaensis Ziegler et al., 2008 in Cat Ba Island of Vietnam; 13–14, G. kwangsiensis and G. kadoorieorum Yang and Chan, 2015 from the southwestern Guangxi, China, exact locality withheld); The G. lichtenfelderi group (green) includes four species (15, G. bawanglingensis in Bawangling, Hainan Island, China, 16, G. hainanensis in Wuzhishan, Hainan Island, China; 17, G. lichtenfelderi in Kuinong Chao [Iles de Norway] and adjacent mainland, north-eastern Vietnam; 18, G. zhoui in the central area of Hainan Island, exact locality withheld).

Fig. 2. A) The partitioned ML tree from RAxML based on the concatenated “All” dataset. The bootstrap support of the RAxML analysis (BS, before slash) and the posterior probability of the BI analysis (PP, after slash) are shown at major nodes. B) The reconstruction of ancestral body

31

band number in Mesquite based on the Assym. Mk model. The probability of ancestral body band number within Goniurosaurus (between the nuchal loop and the caudal constriction) are shown in pie charts at nodes, with blue representing three bands and orange representing four bands.

Fig. 3. The inferred chronogram from BEAST showing the divergence times among clades of Goniurosaurus. A secondary calibration (79 Mya) following a normal distribution (mean = 79, stdev = 9.7) was placed at the most recent common ancestor of Coleonyx and Holodactylus africanus. The blue bar at each node indicates the CI range of time estimation. The most probable ancestral ranges under the DIVALIKE+J model are placed at nodes, with the pie charts showing the probability of each potential range. The ancestral ranges reconstructed by DIVALIKE model are almost identical and are shown in Supplementary Fig. S4.

32

Highlights Molecular phylogenetic reconstruction renewed traditional classification of Goniurosaurus G. bawanglingensis, G. zhoui, G. hainanensis, and G. lichtenfelderi form a monophyletic clade G. bawanglingensis belongs to the G. lichtenfelderi group instead of the G. luii group Body band number might have independently decreased multiple times within Goniurosaurus Reverse colonization from Hainan to Vietnam promoted the divergence of G. lichtenfelderi

33

Graphical abstract

34