利用祖先區域重建法探討東亞群島多稜攀蜥（Japalura polygonata）之歷史生物地理學

全文

(1)國立臺灣師範大學生命科學系碩士論文. 利用祖先區域重建法探討東亞群島多稜攀蜥（Japalura polygonata）之歷史生物地理學 Historical biogeography of Japalura polygonata in East Asian archipelago inferred from ancestral area reconstruction methods. 研究生：楊尚芳 Shang-Fang Yang 指導教授：林思民博士 Si-Min Lin 中華民國 102 年 1 月.

(2) 致謝經過了兩年半的實驗室生活，歷經兩次琉球跳島採集與無數次的台灣島內採集，加上日以繼夜的 Lab work 及討論寫作，終於將這些心血轉化成我的第一本科學著作。最感謝的當然是指導教授林思民老師這些時間的指導，對於論文的方向給了我很大的自由，讓我可以在兩棲爬蟲動物的研究領域中做自己喜歡的研究。而另一方面，也感謝呂光洋老師能夠讓我在大學時期進入生態與演化學的研究領域，讓我早些時候就開始接觸這些有趣的事物並且有所收獲。而對於為何會選擇生物地理學有關的題目則是要感謝徐堉峰老師的啟發，沒有了那些在課堂上得到的知識，我也無法接觸到這個研究題材，也感謝老師在口試時給的建議與提點。另外，也特別感謝林仲平老師在口試時給予的諸多詳細的意見，使我得以順利修改論文。在生活上，我很開心我所待的實驗室都充滿著歡笑，也要感謝其中的各個成員們，有了你們使得實驗室生活充滿了無數的樂趣：呂光洋老師實驗室的大砲，小蛇，BaGa，詩雯，亞融，文琪，以及林思民老師實驗室的彥博，致維，俊文，展蔚，書書，恰恰，龍哥，阿傑，阿平，李昱，閣桓，曾威，阿寶，阿如。最後要感謝我的家人，在這段期間內支持我，使我可以順利的完成碩士班學業。在這裡向在我研究當中被犧牲的攀蜥們獻上最高的敬意，而要向被剪尾巴的攀蜥們說聲抱歉，讓你們受苦了！ i.

(3) Table of contents Introduction………………………………………………………………1 Materials and methods……………..……………………………………….4 Results…………………………..…………………………………..……12 Discussion………………………………………………………….……18 Conclusion………………………………………………………….…...22 References………………………………………………………….……23 Table1……………………………..……………………………….…….33 Table2……………………………..………………………………….….35 Table3………………………………..……………………………….….36 Table4…………………………………..………………………….…….38 Table5……………………………………..………………………….….39 Table6……………………………………..………………………….….40 Table7……………………………………..………………………….….41 Figure1………………………………………….……………………….42 Figure2…………………………………………….…………………….44 Figure3…………………………………………….…………………….45 Figure4…………………………………………….…………………….46 Figure5…………………………………………….…………………….47 Figure6…………………………………………….…………………….48 Figure7…………………………………………….…………………….50 AppendixI………..……………………………………………………...51. ii.

(4) 中文摘要播遷作用（dispersal）與割據作用（vicariance）是生物地理學中解釋生物分布最重要的兩種假說，而大陸型群島（continental archipelagoes）是最適合用於探討此兩種假說之島嶼類型。過去認為生物的分布大多是由割據作用所導致，然而近年來有越來越多的證據顯示，生物現今的分布其實是由播遷作用所主導。琉球群島及台灣為一列位於東亞大陸沿岸之大陸型群島，探討其上生物經歷的歷史事件將有助於我們了解大陸型群島生物相之起源。本研究旨在探討多稜攀蜥（Japalura polygonata）於琉球群島及台灣的歷史生物地理（historical biogeography），以了解影響此種生物現今分布的機制。我們使用來自27個族群的246隻個體，涵蓋了本物種所有的分布區域，並定序粒線體cytochrome b、16S ribosomal RNA、染色體Bach-1、Rag-1等一共四個片段進行親緣地理學分析。利用 maximum parsimony、maximum likelihood以及Bayesian inference建構多稜攀蜥各基因型之親緣關係，並且利用 BEAST 、 dispersal-extinctioncladogenesis model（DEC model）及statistical dispersal-vicariance analysis （S-DIVA）方法進行分析，以重建各族群之祖先地理區域及分化時間。結果顯示多稜攀蜥在各島嶼之間呈現高度的遺傳分化，總共包含八個種內的主要系群，但是與現今認定的四個亞種並不完全吻合。多稜攀蜥的共同祖先最可能起源於4.56百萬年前（上新世早期; early Pliocene）之南琉球及台灣地區，而後逐漸往北擴散至中琉球地區。DEC model分析顯示 iii.

(5) 於上新世時期發生一次割據事件（3.67百萬年前）與兩次播遷事件（3.07 百萬年前），而於更新世時期發生一次割據事件及兩次播遷事件（1.32及 0.47百萬年前）。而S-DIVA分析結果則與DEC model結果類似，但存在較多的播遷事件（六次播遷事件）。總結來說，我們的結果顯示整個多稜攀蜥於過去近五百萬年的歷史主張着此物種起源於南琉球及台灣地區後經歷了早期的割據事件接著經由跳島播遷（island-hopping dispersal）拓殖至此物種分布的最東北。. 關鍵字歷史生物地理、祖先區域重建、琉球群島、台灣島、飛蜥科、多稜攀蜥。. iv.

(6) Abstract Continental archipelagoes are suitable objects for examination of alternative hypotheses in biogeography. Vicariance was once believed to be the predominant mechanism that influenced distribution of organisms. However, dispersal as the alternative one has been widely discussed in recent years. The Ryukyu archipelago and Taiwan are suitable continental archipelago to study these alternative hypotheses in East Asia. In this study, we aim to infer the historical biogeography of a dispersal-limited terrestrial vertebrate, the Okinawa tree lizard (Japalura polygonata). We obtained 246 individuals sampled from 10 localities of Ryukyu archipelago and 17 localities in Taiwan, covering the entire distribution range of this species. DNA sequences of the mitochondrial cytochrome b, 16S ribosomal RNA, nuclear Bach-1 and Rag-1 genes were obtained from these individuals. We employed maximum parsimony, maximum likelihood and Bayesian inference analyses to reconstruct phylogenetic relationships among haplotypes. We also performed BEAST molecular dating, dispersal-extinction-cladogenesis (DEC) model and statistical dispersal-vicariance analysis (S-DIVA) to reconstruct ancestral geographical ranges and the divergence time of the populations of J. polygonata. Phylogeny of J. polygonata haplotypes revealed the existence of eight major clades, which were not consistent with the four currently recognized subspecies. Biogeographic and divergence time reconstructions showed that J. polygonata originated from Taiwan and the southern Ryukyus 4.56 million years ago (early Pliocene). The DEC models analysis showed that there were one vicariance event and two dispersal events occurred in the middle Pliocene (3.67 MYA) and late Pliocene (3.07 MYA), respectively, and one vicariance event and two dispersal events occurred in the Pleistocene v.

(7) (1.32 and 0.47 MYA, respectively). The S-DIVA analysis revealed similar results but revealed more dispersal events than DEC model analysis (six dispersal events). These results support that both dispersal and vicariance mechanisms were major mechanisms shaping the current distribution of J. polygonata in East Asian archipelago and Taiwan. In conclusion, the biogeographic history of J. polygonata suggested that this species was originated from Taiwan and southern Ryukyu region, experienced historical vicariance events, and then colonized to their north most distribution by island-hopping dispersal during the past 4.56 million years.. Keywords Agamidae, ancestral area reconstruction, historical biogeography, Japalura polygonata, Ryukyu archipelago, Taiwan.. vi.

(8) Introduction According to Alfred Russel Wallace’s classification, continental islands are one of three types of ‘true’ islands, which are characterized for their close relationship either in geological history or geographic distance to nearest continent (Wallace, 1895). Continental archipelagoes provide lots of suitable objects for examining different biogeography scenario caused by landmass connection or disjunction (Hedges et al., 1992; Hisheh et al., 1998; Atkins et al., 2001; Ohdachi et al., 2001; Poulakakis et al., 2003; Bittkau & Comes, 2005). Thus, studies of continental archipelagoes may be a profitable progress for the central debates in biogeography since 19th century: vicariance versus dispersal for the origins of discontinuous populations (de Queiroz, 2005; Morrone, 2009). The vicariance hypothesis assumes that the formation of geographical barriers disconnected continuous distribution ranges of ancestral organisms, and disrupts gene flow between the separated subpopulations (de Queiroz, 2005; Crisp et al., 2011). In contrast, dispersal hypothesis supposes that ancestral organisms dispersed over geographical barriers to occupy their new distribution ranges, and sufficiently limited subsequent gene flow between the parent and daughter populations (de Queiroz, 2005; Crisp et al., 2011). The East Asian Islands, comprising the Ryukyu archipelago and Taiwan, have the greatest potential to study the mechanisms of vicariance and dispersal in the Asia-Pacific region. In this area, the archipelago composed approximately 140 subtropical islands and forming a chain of disjunctive landmasses in northeast-southwest direction between Kyushu of 1.

(9) Japan and east coast of China (Fig. 1). There are four major island groups in this chain: (1) Northern Ryukyu (including Osumi and Tokara islands); (2) Central Ryukyu (including Amami and Okinawa islands); (3) Southern Ryukyu (Miyako and Yaeyama islands); and (4) Taiwan and adjacent islands (Fig. 1). These islands attracted lots of attention from a biogeographic perspective because of their history of land connections and fragmentations (Ota, 1998; Chiang & Schaal, 2006; Honda et al., 2008; Nakamura et al., 2009; Takada et al., 2010; Tominaga et al., 2010). The fauna of the Ryukyu archipelago is characterized by a high ratio of endemic taxa, most of which are highly differentiated from their relatives as a result of island formation and long-term isolation (Ota, 1998, 2000). As a group of less mobile terrestrial vertebrate, over three-fourths (76.6 %) of terrestrial amphibians and reptiles species are endemic to the Ryukyu archipelago (Ota, 2000). Most of these species are restricted to finite geographic area, only few among which are widely distributed across the entire archipelago. There are many researches about these species but generally focused on species level relationship (Toda et al., 2001; Matsui et al., 2005; Honda et al., 2008) and biogeography (Lin et al., 2002; Ota et al., 2002), or population genetics (Toda et al., 1997). However, none of these researches discussed about mechanisms shaping the distribution of organisms in this area. Because of the wide distribution in Ryukyu archipelago and Taiwan, the Okinawa Tree Lizard Japalura polygonata (Hallowell, 1861) is a suitable material to study inter-island dispersal or vicariance mechanisms.. 2.

(10) Japalura polygonata distributes in northern Taiwan and the most islands of Ryukyu archipelago south of the Tokara Gap (Ota, 1991, 2003; Biodiversity Center of Japan, 2010). It consists of four subspecies according to their morphological differences: J. p. polygonata, J. p. ishigakiensis, J. p. donan and J. p. xanthostoma (Hallowell, 1861; van Denburgh, 1912; Ota, 1991, 2003). Most Japalura species distribute over mainland Asia including southwestern China, western China, Indo-China peninsula and northern India (Uetz, 2006), yet J. polygonata is the only one species specialized in insular environment. All of the four subspecies were abundant in the margin of low-altitude evergreen broad-leaf forests on these islands (Shang, 2008; Uchiyama et al., 2009). Mechanism for amphibians and reptiles to disperse among landmasses is to stride across ocean through logs and mats of vegetation by ocean current (Censky et al., 1998; Calsbeek & Smith, 2003; Vences et al., 2003). There is a strong western boundary current, Kuroshio Current, flowing through the East Asia margin. The Kuroshio Current begins off the east coast of Taiwan and flows northeastward past Japan, where it merges with the North Pacific Current. Thus, to some degree J. polygonata could possibly stride across ocean by Kuroshio Current. On the other hand, the emerging and submerging of land bridge between Ryukyu archipelago and Taiwan may also influence the dispersal of organisms in this area. During the Pleistocene, land bridges between Ryukyu archipelago and Taiwan were emerged and submerged frequently due to the oscillations of sea level (Ota, 1998). Therefore, J. polygonata in 3.

(11) Ryukyu archipelago and Taiwan were anticipated to experience periodic cycles of contact and isolation between adjacent islands. According to all information that mention above, the special status of continental archipelago system and the unique pattern of distribution with ecological and geographical restriction make this species as an extraordinary material to examine the processes of historical biogeography in this area. In this study, we aim to infer the historical biogeography of the Okinawa tree lizard populations in Ryukyu archipelago and Taiwan using mitochondrial and nuclear DNA sequencing. Bayesian phylogenetic analysis with relaxed molecular clock models (Drummond et al., 2006) was applied to reconstruct the relationship and divergence time of the lizards among the islands. Furthermore, we employ event-based and model-based biogeographical inferences that so-called statistical dispersal-vicariance analysis (S-DIVA) (Ronquist, 1997; Nylander et al., 2008; Harris & Xiang, 2009) and dispersal-extinction-cladogenesis (DEC) models (Ree et al., 2005; Ree & Smith, 2008), respectively. The aims of this study are to answer the following questions: (1) Which area is the origin of the J. polygonata? (2) What is the major mechanism (vicariance, dispersal, or both) shaped the current distribution of J. polygonata?. Materials and methods Sample collection A total of 246 individuals of Japalura polygonata, including all four. 4.

(12) subspecies (J. p. polygonata, J. p. ishigakiensis, J. p. donan and J. p. xanthostoma), were collected from 10 localities in Ryukyu archipelago and 17 localities in Taiwan from 2010 to 2011 (Table 1 and Fig. 1). Specimens were euthanized, and tissues were preserved in 95% ethanol (skeletal muscles or tail samples). We used five Japalura (J. splendida, J. swinhonis, J. luei, J. makii, and J. brevipes) and two agamids (Calotes emma and Acanthosaura lepidogaster) as out groups to root the phylogenetic tree.. DNA extraction, PCR and sequencing Total genomic DNA was extracted from tissue samples using Qiagen DNeasy Blood & Tissue Kit (Qiagen Inc., 2009) according to the manufacturer’s instructions. Total genomic DNA was used for amplifying two mitochondrial genes (cytochrome b; cyt b and 16S ribosomal RNA; 16S rRNA) and two nuclear genes (BTB and CNC homology 1; Bach-1 and recombination activating gene 1; Rag-1) via polymerase chain reaction (PCR). Primers for amplifying all four genes were list in Table 2. Double-stranded polymerase chain reactions were performed in 20 μl reaction volume with the following thermal cycle: 1 cycle at 94°C for 3 minutes, followed by 35 cycles at 94°C for 30 seconds, each annealing temperature (Table 2) for 40 seconds, and 72°C for 60 seconds, and with a final cycle at 72°C for 10 minutes. PCR products were run on 1.5% agarose gels in 1X TBE buffer to ensure that the lengths of the certain fragments were correctly amplified. Automatic PCR products sequencing were carried out on ABI 3730 by Genomics BioSci & Tech Corp. (Taipei, Taiwan). Raw chromatograph data was checked and edited DNA sequences by using 5.

(13) Sequencher v.4.9 (GeneCode, Boston, MA, USA). Sequences were aligned by using the default setting in ClustalX (Jeanmougin et al., 1998). Mitochondrial sequencing of cytochrome b gene was applied to all the 246 samples in our collection. 16S rRNA and the two nuclear genes were amplified and sequenced for 57 individuals selected from 1-4 individuals each locality which comprising all the major clades in this species. According to our purposes, we aligned sequences into four datasets that differ from individual number and gene. All of the 246 cyt b sequences were aligned into Dataset 1. Dataset 2 consisted 57 cyt b and 16S rRNA sequences. Bach-1 and Rag-1 sequences were aligned into Dataset 3 and Dataset 4, respectively.. Phylogenetic analyses All of the sequence substitution models of gene regions were determined under Akaike information criterion (AIC) in jModelTest v.0.1.1 (Posada, 2008) before phylogenetic analyses. We used Dataset 1 to reconstruct single gene phylogenetic relationships of J. polygonata (primary tree), which were assessed using maximum parsimony (MP), maximum likelihood (ML) and Bayesian inference (BI) (primary tree). The MP analysis was performed in PAUP* 4.0b (Swofford, 2002). Analyses utilized heuristic searches by following steps: starting trees determined by 100 random taxon additions, tree bisection-reconnection (TBR) branch swapping. Each sequence was treated as an operational taxonomic unit and each nucleotide site was treated as a character. All characters were unordered and equally weighted. Non-parametric 6.

(14) bootstrapping with 1000 pseudo replicates was performed to obtain robustness of node support for the resulting trees (Felsenstein, 1985). The ML analysis was performed in PhyML v.3.0 (Guindon et al., 2010) with 1000 non-parametric bootstrap replicates. BI was implemented in MrBayes v.3.1.2 (Ronquist & Huelsenbeck, 2003) using two searches with 10 million generations each and sampling every 1000 generations. The 50% majority consensus tree with Bayesian posterior probabilities (BPPs) of clades was calculated to obtain the Bayesian estimate of phylogeny after removing the first 25% of sampled generations as burn-in. Based on this primary tree, we selected 57 individuals from the main clades of this tree for sequencing other three gene regions (mitochondrial 16S rRNA, nuclear bach-1 and rag-1). Two nuclear genes were used in haplotype data analysis below. The identical sequences of individuals of both mitochondrial gene regions were combined in a single matrix (Dataset 2), and the phylogenetic relationships were assessed using maximum parsimony (MP), maximum likelihood (ML) and Bayesian inference (BI) (secondary tree). Sequence data were partitioned into two partitions by gene region within Dataset 2 (only in ML and BI). The parameter settings of MP were same as above. The partitioned ML phylogenetic analysis was performed in RAxML v.7.2.6 (Stamatakis, 2006; Silvestro & Michalak, 2011) with independent GTR+I+G substitution models applied to 2 partitions and 1000 non-parametric bootstrap replicates. The partitioned BI was implemented in MrBayes v.3.1.2 (Ronquist & Huelsenbeck, 2003). The independent GTR+I+G and TVM+G substitution models were applied. 7.

(15) to cyt b and 16S rRNA regions within Dataset 2 respectively. All other parameter settings of partitioned BI were the same as non-partitioned BI above. Owing to the comparatively low genetic polymorphism, gene genealogies of the two nuclear genes (Bach-1 and Rag-1) were expressed as haplotype network using the concatenated unphased sequence dataset (combine Dataset 3 and Dataset 4). The sequence substitution models were determined under Akaike information criterion (AIC) in jModelTest v.0.1.1 (Posada, 2008). Phylogenetic reconstructions among haplotypes were estimated using a maximum likelihood (ML) approach, as implemented in the software PhyML v.3.0 (Guindon et al., 2010) with GTR+I+G substitution models and 1000 non-parametric bootstrap replicates. The generated trees were used to estimate network haplotype in Haploviewer program (Salzburger et al., 2011). With limited information, nuclear sequences were not used in the following analyses including molecular dating and ancestral range construction.. Molecular dating analyses For the purpose to understand the temporal scale to the biogeography of J. polygonata, we employed site-specific partitioned Bayesian phylogenetic analyses with relaxed molecular clock approach in BEAST v.1.7.2 (Drummond et al., 2006; Drummond & Rambaut, 2007), a software package that simultaneously estimates node ages and tree topology. Because different tree priors are appropriate for different phylogenetic scales (e.g., interspecific and intraspecific levels; Drummond et al., 2006), 8.

(16) we infer the divergence times among Japalura species using external calibration age constraints, and then applied the resultant values to internal calibration to infer the divergence times among J. polygonata lineages. In the external calibration step, we included 12 and 9 sequences from Chamaeleonidae and Agamidae, respectively (Okajima & Kumazawa, 2010) and aligned with our sequences (including 13 clades of J. polygonata, 5 other Japalura species and 2 out group species) to infer the divergence time of among J. polygonata and other congeners. We chose two calibration points as follow, (1) the divergence of Chamaeleonidae and Agamidae (84 MYA) (Amer & Kumazawa, 2005; Hugall et al., 2007); (2) the divergence between African and Arabian chameleons within genus Chamaeleo (5<T<13 MYA) (Macey et al., 2008) (see [Amer & Kumazawa, 2005] for geological evidence for this time constraint). Divergence times were calibrated by using BEAST software (Drummond & Rambaut, 2007). The divergence time between Chamaeleonidae and Agamidae was set to a normal distribution with mean = 84, standard deviation = 4. Furthermore, we incorporated a calibration within Chamaeleo which was set a normal distribution with mean = 9 and standard deviation = 2. Moreover, the best-fitting substitution models of each site-specific partition were determined under AIC in jModelTest v.0.1.1 (Posada, 2008). The model testing results indicated that the GTR+I+G substitution model were the best-fitting model for all three codon partitions. The substitution rate variation was modeled using an uncorrelated lognormal distribution, and a Yule process was employed as tree prior. The MCMC analysis was run for 108 generations, with a sample frequency of 1000. 9.

(17) After external calibration step, we employed an internal calibration step including 57 J. polygonata individuals and 7 outgroup taxa using the date distribution of the most recent common ancestor of four Japalura species (including J. polygonata, J. luei, J. makii and J. brevipes) and the most recent common ancestor of J. polygonata that both inferred by external calibration step as internal calibration (normal distribution with mean = 10.60 and standard deviation = 1.35; mean = 5.26 and standard deviation = 0.75, respectively). The model testing results indicated that the GTR+G substitution model was the best-fitting model for all three codon partitions. This analysis used a coalescent-constant size tree prior, and run for 108 generations with a sample frequency of 1000. The effective sample size (ESS) value of all time intervals that we need in both calibration analyses were confirmed in Tracer v.1.5 (Rambaut & Drummond, 2007) to ensure this value above 150. We removed the first 10% of sampled generations as burn-in. Finally, we obtained maximum clade credibility (MCC) tree in TreeAnnotator v.1.4.8 and visualized in FigTree v.1.3.1.. Ancestral area reconstruction In order to infer biogeographic events of J. polygonata among East Asian islands, ancestral area reconstructions were coded in five areas based on the presence of barriers between island groups (Fig.1): Amami Group (A); Okinawa Group (O); Miyako Group (M); Yaeyama Group (Y) and Taiwan Island (T). We. used. the. dispersal-extinction-cladogenesis 10. (DEC). model.

(18) implemented in the latest snapshot archive of LAGRANGE (Ree & Smith, 2008) (version 20120508) to estimate the likelihood of range inheritance scenarios at node in phylogeny, and allow incorporation of information about alternative dispersal routes or probabilities among these five areas. This analysis used the MCC tree obtained in the BEAST analysis above. In our DEC model analyses, we construct two different models (M0 and M1) representing different dispersal probabilities between these areas through time. We set the root age according to the result of molecular dating analysis above. The unconstrained model (M0) allows geographic ranges to include any combination of adjacent areas that we defined above, and permits constant dispersal probability (dispersal probability = 0.1) between any pair of them. In M1, the limit of geographic ranges was unconstrained, but dispersal probability was set to vary according to historical connections among areas: it was maximum (dispersal probability = 1) during 2.58-0.02 MYA (Pleistocene) because of the frequently fluctuation of sea level during this particular period which provided a well environment for dispersal between islands; another, the probability was set to a lower value (dispersal probability = 0.1) during the time intervals before and after the Pleistocene period (2.58-0.02 MYA) that mentioned above. For. comparison,. we. also. performed. statistical. approach. to. dispersal-vicariance analysis (S-DIVA) (Ronquist, 1997; Nylander et al., 2008; Harris & Xiang, 2009) implemented in RASP v.2.1a (Yu et al., 2010, 2011), a software that statistically evaluates the alternative ancestral ranges at nodes based on a set of trees (Nylander et al., 2008; Harris & Xiang, 2009), to account for phylogenetic uncertainty and uncertainty in area 11.

(19) optimization. For this analysis, ancestral areas were reconstructed on 1000 randomly chosen post burn-in trees of the BEAST analysis. The maximum number of ancestral areas at each node was set to five (unconstrained). Relative frequencies of ancestral areas reconstructed for each node were plotted on the MCC tree from the BEAST analysis. As pointed out by Ronquist (1996), ancestral area reconstructions in DIVA become less reliable as we approach the root node. This was verified in DIVA as a tendency for the root distribution to be large and include all the areas we gave. In order to improve the reliability at the basal node of J. polygonata clade, we run the analyses based on the larger data set than DEC model analyses including other three species of Japalura (J. brevipes, J. makii and J. luei) that endemic in Taiwan, so that the basal node of J. polygonata is no longer the root node in our analyses.. Results Phylogenetic relationships The characteristics of the two mitochondrial and two nuclear DNA regions are summarized in Table 1 and Table 3. The fully aligned length of the Dataset 1, Dataset 2, Dataset 3 and Dataset 4 matrixes were 1132 bp, 2311 bp, 1235 bp and 1330 bp (Table 4), respectively. Among the four fragments, the two mitochondrial genes had polymorphic and parsimony informative sites roughly 10 times higher than the two nuclear genes, as also shown in Table 4. The pairwise distance comparison of Japalura. 12.

(20) polygonata populations using mitochondrial sequences is provided in Table 5. The maximum likelihood (ML) analysis tree of the Dataset 1 is shown in Figure 2 (primary tree). Results of maximum parsimony (MP) and Bayesian inference (BI) analyses yielded similar topologies (Fig. 2). All three analyses showed that Japalura polygonata formed a monophyletic clade with high support values (ML bootstrap = 100%; MP bootstrap = 100; BI posterior probability = 1). Totally eight clades, all with high statistic supports and congruent to certain island group(s), were defined as follows: (I) the Amami-Okinawa clade; (II) the Miyako clade; (III) the Itiomote clade; (IV) the Ishigaki clade; (V) the Yonaguni clade; and (VI), (VII), (VIII), representing three local lineages from eastern, western, and northern Taiwan, respectively. Among the four currently recognized subspecies, only J. p. polygonata (clade I) and J. p. donan (clade V) are monophyly. Japalura p. donan is clustered with the Eastern Taiwan clade (VI), although the statistic support was low. Phylogeny of Dataset 2 is represented in Figure 3 (secondary tree). All the three phylogenetic analyses yielded to a congruent result (Fig. 3) and were highly similar to the topology of the primary tree (Fig. 2). As in the primary tree, all the eight clades defined were strongly supported with high statistic values (Fig. 3). The only difference was in interrelationship between clades VII and VIII: in Fig. 3, these two clades form a monophyly, while in Fig. 2, clade VII was cluster with the other clades. Phylogeny of. J. polygonata haplotypes. indicated. a revised. interrelationship among the four currently recognized subspecies. Only J. p. 13.

(21) polygonata (clade I) and J. p. donan (clade V) were monophyletic. The other two subspecies, P. p. ishigakiensis and P. p. xanthostoma, are paraphyletic.. Genotype network of nuclear genes The genotype network of the two nuclear genes (Bach-1 and Rag-1, concatenated unphased sequences) is shown in Fig. 4. Four major haplotype groups, (a) to (d), were defined as shown in Fig. 4. The grouping of nuclear genotype network represent a pattern similar to mitochondrial gene phylogeny: group (a) represented individuals collected form Amami and Okinawa islands, and was closely related to group (b) members, which were collected from Miyako and Yaeyame islands. Some individuals from Yaeyama group were cluster with individuals from northern and eastern Taiwan. In group (a), eight individuals (8/18 or 44.4%) belonged to the most frequent haplotype (Fig. 4), while in other three groups, all individuals belong to unique haplotype of their own (Fig. 4). Genetic diversity of these populations (Table 3) indicated a trend of decreasing genetic diversity as latitude increased.. Divergence time estimation The result of external calibration analysis is provided in Appendix I. Estimated time of the internal nodes of J. polygonata is provided in Table 6. Ages are presented as the mean and 95% highest posterior densities (HPD) of the posterior distribution. The chronogram of internal calibration analysis is shown in Figure 5. The phylogeny calculated from the BEAST 14.

(22) analysis resulted in the same topology as that estimated by ML and BI analyses (Fig. 3; Fig. 5). The first lineage differentiation within J. polygonata clade occurred during the Pliocene [4.56 (3.61-5.57) MYA] (node 1 in Fig. 5), close to the other two major differentiations (node 2 and node 23) within J. polygonata clade [4.18 (3.20-5.19) and 4.14 (3.17-5.14) MYA, respectively] (Fig. 5; Table 6). Except for the Yonagunijima population of J. polygonata, the chronogram supported a crown-age for other Ryukyu archipelago populations (clade I to clade IV) around 3.07 (2.19-4.03) MYA (node 3 in Fig. 5). The most recent common ancestor for the three Taiwanese clades (the eastern, the northern, and the western clades) were 0.60 (0.35-0.94), 0.86 (0.53-1.25) and 3.23 (2.36-4.15) MYA, respectively (Fig. 5; Table 6).. Ancestral area reconstruction The two biogeographic models tested in LAGRANGE generated results with similar global likelihood values (Table 7). Model M0 received higher global log-likelihood units (-28.24), but model M1 fell within 2 log-likelihood units of higher one. Thus, there did not exist a significantly “best” model. For nodes where likelihood units were not significantly different between multiple reconstructions within the confidence window of a 2 log-likelihood unit difference, the relative probability of the global likelihood is shown in Table 7. The colonization histories of two models are identical (Fig. 6A and Fig. 6B). Both have two vicariance events and four dispersal events (Fig. 6). Model M0 resulted in a higher number of alternative scenario falling within 15.

(23) 2 log-likelihood units of the optimal reconstruction (Table 7). Due to there is no significantly best model, we selected the highest likelihood value model (model M0) as the optimal reconstruction (Fig. 6A). In model M0, there are eight uncertainty nodes, including root node (node 1) (Table 7) that alternative reconstruction might influence the interpretation of biogeographic scenarios. The dispersal and vicariance events (indicated by arrows and black bar in Figure 6, respectively) of model M0 were inferred as follows under optimal reconstruction (Fig. 6A): (1) the J. polygonata originated in Yaeyama Group (Y) and Taiwan Island (T); (2) a vicariance event occurred between Yaeyama Group and Taiwan Island between node 2 and 20; (3) two dispersal events occurred from Yaeyama Group to Miyako Group (M) and Okinawa Group (O) between node 3 and 4; (4) a vicariance event occurred between Miyako Group and Okinawa Group between node 4 and 5; (5) a dispersal event occurred from Okinawa Group to Amami Group (A) between node 8 and 9; (6) same as (5) but between node 8 and 12. Result of S-DIVA analysis was mostly congruent with DEC model inferences. As in DEC model reconstructions, the range of the most basal node is undetermined (Table 7). There are three possible ancestral area types on root node (node 1): T: 35.04% marginal probability, YT: 35.04% marginal probability and OMYT: 29.89% marginal probability. Because of the vicariance event have zero cost in DIVA method (Ronquist, 1997) the optimal area reconstruction of node will be the widespread ancestor. Thus, the optimal area reconstruction of root node (node 1) is the widespread one (YT) (Fig. 7). Based this result, the dispersal and vicariance events 16.

(24) (indicated by arrows and black bar in Figure 7, respectively) were inferred as follows (Fig. 7): (1) the J. polygonata originated in Yaeyama Group (Y) and Taiwan Island (T); (2) a dispersal event occurred from Yaeyama Group to Taiwan Island between node 1 and 2; (3) a dispersal event occurred from Taiwan Island to Yaeyama Group and followed by a vicariance event between node 2 and 20; (4) two dispersal events occurred from Yaeyama Group to Miyako Group (M) and Okinawa Group (O) between node 3 and 4; (5) a vicariance event occurred between Miyako Group and Okinawa Group between node 4 and 5; (6) a dispersal event occurred from Okinawa Group to Amami Group (A) between node 8 and 9; (7) same as (6) but between node 8 and 12. In our reconstructions, the DEC model and S-DIVA methods generated two biogeographic scenarios with only slight difference (Fig. 6A and Fig. 7). The DEC model suggested that there were four dispersal events and two vicariance events, and the S-DIVA model suggested six dispersal and two vicariance events. Except the first two dispersal events in S-DIVA result, all other events are identical between two methods. The first dispersal event (between node 1 and node 2; Fig. 7) is interpreted from a single area ancestor Y (Yaeyama Group) to a widespread descendant YT (Yaeyama Group and Taiwan Island). The second dispersal event (between node 2 and node 20; Fig. 7B) is interpreted from a single area ancestor T to a widespread descendant YT.. 17.

(25) Discussion Early divergence of J. polygonata lineages Given that Japalura polygonata is strongly supported as a monophyletic group, biogeographic history of this insular species could be reliably inferred. Our biogeographic reconstructions are conclusive on the ancestral area of J. polygonata from the two ancestral area reconstruction methods, which clearly support a south-originated (Yaeyama Group and Taiwan Island) ancestor (Fig 6A; Fig 7). The crown-age of J. polygonata based on our dating analysis supports an origin of this species around 4.56 (3.61, 5.57) MYA (late Miocene to middle Pliocene) (Fig. 5). There is only one vicariance event in the early history of J. polygonata based on DEC model result (Fig. 6A). The event on node 20 suggested a clear vicariance scenario between Y (Yaeyama Group) and T (Taiwan Island) occurred 3.67 (2.71, 4.69) MYA (middle Pliocene) (Fig. 6A). This result could be considered as the separation between populations in Yayeyama and Taiwan caused by the formation of Yonaguni Strait (between Taiwan Island and Yonagunijima) due to collision of the Luzon arc in Pliocene (about 5 MYA) (Lallemond et al., 2001). In contrast, the result of S-DIVA revealed a different interpretation of biogeographic scenarios in the early history of J. polygonata. The results of S-DIVA on node 1, node 2 and node 20 all suggested a widespread ancestor YT separated into Y and T descendants (Fig. 7). Under this reason, the interpretation of biogeographic scenario between node 1 and node 2 is dispersal from Y to T, and between node 2 and node 20 is also dispersal but. 18.

(26) from T to Y (Fig. 7). These dispersal events occurred in late Miocene to middle Pliocene [4.56 (3.61, 5.57) and 4.18 (3.20, 5.19) MYA; Fig. 5]. There is a strong northward current, Kuroshio Current, flowing through the East Asia margin. The Kuroshio Current is probably a major force that helped these lizards stride across ocean through logs and mats of vegetation during late Miocene to middle Pliocene. Because there are a series of geological events (opening of the Okinawa Trough and formation of Yonaguni Strait) (Kimura, 2000; Lallemond et al., 2001) in this time period allowing the Kuroshio Current pass through the northeast of Taiwan Island, this ocean current have the opportunity become a media of species dispersal.. Colonization of Miyako, Okinawa and Amami Group Both methods agree on inferring that J. polygonata colonized northeastward from Yaeyama Group to Miyako Group and Okinawa Group in the late Pliocene [3.07 (2.19, 4.03) MYA; Fig. 5]. The Okinawa Group population finally colonized to Amami Group, which deduced to have occurred twice in the middle Pleistocene [0.47 (0.29, 0.67) MYA; Fig. 5] (Fig. 6A; Fig.7). During the late Pliocene period (3.3 to 3 MYA), the global temperature was 2 to 3°C higher than present (Robinson et al., 2008), and the global sea level also 25 meters higher than present (Dwyer & Chandler, 2009) suggested that there were no connections between islands in East Asian archipelago during this time period. These evidences response to that the origin of Miyako Group population was not cause by range expansion via 19.

(27) land bridges but cross-ocean dispersal from Yaeyama Group in the late Pliocene [3.07 (2.19, 4.03) MYA; Fig. 5]. Owing to the extreme depth of Kerama Gap (sill depth 1050 meters; Choi et al., 2002), the land mass between Okiwanajima and Miyakojima islands have never connected with each other even during the coldest glacial era, which has lowered the sea level for approximately 130 to 140 meters (Millien-Parra & Jaeger, 1999; Kimura, 2000). As a result, the vicariance event between Miyako Group and Okinawa Group might not cause by the sea level fluctuation in the early Pleistocene [1.32 (0.91, 1.80) MYA; Fig. 5]. In early Pleistocene, present-day Amamioshima, Tokunoshima, and Okinawajima islands formed a large super-island (Ota, 1998) until the super-island was beginning separate to several small super-islands during the late Pleistocene (Ota, 1998). Thus, the two dispersal events from Okinawa Group to Amami Group in the middle Pleistocene [0.47 (0.29, 0.67) MYA; Fig. 5] might not due to the cross-ocean dispersal but range expansion inside the super-island.. Island-hopping across ocean The whole biogeographic history of J. polygonata suggested that this species originate from Taiwan and southern Ryukyu region, and arrived at their current distribution in the northern Ryukyu by ‘island-hopping dispersal’ during almost past five million years. This model has been applied to explain the current distribution of many taxa on Canary Islands, e.g. the wild olive tree Olea (Hess et al., 2000), Gallotia lizards (Thorpe et 20.

(28) al., 1994), Chalcides skinks (Brown & Pestano, 1998), and Hegeter beetles (Juan et al., 1996, 1998). All of these studies showed a consistent stepwise dispersal from the east to the west. Flightless taxa could stride across ocean barrier by rafting on vehicle-like floating vegetation. Considering the formation of these vehicles, tropical storms are known to generate the large, floating tree ‘islands’, as well as associated precipitation, that might make a successful ocean voyage of this type possible (Thiel & Haye, 2006). On the other hand, release and recapture experiments with buoys have demonstrated that flotsam can drift across large distances in the Caribbean Sea and may take a meandering course affected by dominant currents and storm conditions (Molinari et al., 1979). Thus, short to medium distance dispersal events accumulate to island-hopping dispersal in these taxa may be more than previously thought, and J. polygonata provides a continental archipelago example.. Taxonomy versus phylogeny incongruence This study provides the first full-scale phylogeny of J. polygonata with a complete taxon sampling (from Taiwan to Amamioshima). Molecular phylogeny indicated that not all the four currently recognized subspecies are monophyletic (Fig. 3). The only one precisely congruent to current taxonomy is the clade from Amamioshima and Okinawajima populations (J. p. polygonata). The Yonagunijima population, currently recognized as J. p. donan, form another clade but clustered with the East Taiwan population, which is recognized as J. p. xanthostoma. The remainders (J. p. ishigakiensis and J. p. xanthostoma) are paraphyletic respective to the 21.

(29) former. Therefore, this result provides an example of subspecies boundary contradicting molecular phylogenetic data. Nowadays, an increasing number of empirical examples have shown incongruence for the definition of subspecies between traditional taxonomy and molecular phylogeny, no matter in birds (Phillimore & Owens, 2006), mammals (Tursi et al., 2012) or reptiles (Daniels et al., 2005). The widespread use of integrative taxonomic approaches may increase the objectivity of taxon descriptions in all groups of organisms. Thus, the taxonomic status between four subspecies of J. polygonata will confirm after collecting more data or using more approaches in future studies.. Conclusion In this study, the whole biogeographic history of J. polygonata suggested that this species was originate from Taiwan and southern Ryukyus, experienced historical vicariance events between the two regions driven by geological events, and then colonized to the northern islands by island-hopping. dispersal.. Furthermore,. dispersal. is. the. dominant. mechanism that shaped the current distribution of J. polygonata, and this result probably cause by the ocean current through rafting on vehicle-like floating vegetation.. 22.

(30) References Amer, S.A.M. & Kumazawa, Y. (2005) Mitochondrial DNA sequences of the Afro-Arabian spiny-tailed lizards (genus Uromastyx; family Agamidae):. phylogenetic. analyses. and. evolution. of. gene. arrangements. Biological Journal of the Linnean Society, 85, 247-260. Atkins, H., Preston, J. & Cronk, Q.C.B. (2001) A molecular test of Huxley's line: Cyrtandra (Gesneriaceae) in Borneo and the Philippines. Biological Journal of the Linnean Society, 72, 143-159. Biodiversity Center of Japan. (2010) The National Survey on the Natural Environment Distributional Survey of Japanese Animals: Animal Distribution Atlas of Japan. Biodiversity Center of Japan, Tokyo, Japan. Bittkau, C. & Comes, H.P. (2005) Evolutionary processes in a continental island system: molecular phylogeography of the Aegean Nigella arvensis alliance (Ranunculaceae) inferred from chloroplast DNA. Molecular Ecology, 14, 4065-4083. Brown, R.P. & Pestano, J. (1998) Phylogeography of skinks (Chalcides) in the Canary Islands inferred from mitochondrial DNA sequences. Molecular Ecology, 7, 1183-1191. Calsbeek, R. & Smith, T.B. (2003) Ocean currents mediate evolution in island lizards. Nature, 426, 552-555. Censky, E.J., Hodge, K. & Dudley, J. (1998) Over-water dispersal of lizards due to hurricanes. Nature, 395, 556.. 23.

(31) Chiang, T.Y. & Schaal, B.A. (2006) Phylogeography of plants in Taiwan and the Ryukyu Archipelago. Taxon, 55, 31-41. Choi, B.H., Kim, K.O. & Eum, H.M. (2002) Digital bathymetric and topographic data for neighboring seas of Korea. Journal of Korean Society of Coastal and Ocean Engineers, 14, 41-50 (in Korean). Crisp, M.D., Trewick, S.A. & Cook, L.G. (2011) Hypothesis testing in biogeography. Trends in Ecology and Evolution, 26, 66-72. Daniels, S.R., Heideman, N.J.L., Hendricks, M.G.J., Mokone, M.E. & Crandall, K.A. (2005) Unraveling evolutionary lineages in the limbless fossorial skink genus Acontias (Sauria: Scincidae): are subspecies equivalent systematic units? Molecular Phylogenetics and Evolution, 34, 645-654. de Queiroz, A. (2005) The resurrection of oceanic dispersal in historical biogeography. Trends in Ecology and Evolution, 20, 68-73. Drummond, A.J. & Rambaut, A. (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology, 7, 214. Drummond, A.J., Ho, S.Y.W., Phillips, M.J. & Rambaut, A. (2006) Relaxed phylogenetics and dating with confidence. PLoS Biology, 4, e88. Dwyer, G.S. & Chandler, M.A. (2009) Mid-Pliocene sea level and continental. ice volume based. on. coupled benthic Mg/Ca. palaeotemperatures and oxygen isotopes. Philosophical Transactions of the Royal Society A: Mathematical, Physical & Engineering Sciences, 367, 157-168. Felsenstein, J. (1985) Confidence limits on phylogenies: an approach using the bootstrap. Evolution, 39, 783-791. 24.

(32) Guindon, S., Dufayard, J.F., Lefort, V., Anisimova, M., Hordijk, W. & Gascuel, O. (2010) New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Systematic Biology, 59, 307-321. Hallowell, E. (1861) Report upon the Reptilia of the North Pacific Exploring Expedition, under command of Capt. John Rogers, U. S. N. Proceedings of the Academy of Natural Science of Philadelphia, 12, 480-510. Harris, A.J. & Xiang, Q.Y. (2009) Estimating ancestral distributions of lineages with uncertain sister groups: a statistical approach to Dispersal-Vicariance Analysis and a case using Aesculus L. (Sapindaceae) including fossils. Journal of Systematics and Evolution, 47, 349-368. Hedges, S.B., Hass, C.A. & Maxson, L.R. (1992) Caribbean biogeography: molecular evidence for dispersal in West Indian terrestrial vertebrates. Proceedings of the National Academy of Sciences of the United States of America, 89, 1909-1913. Hess, J., Kadereit, J.W. & Vargas, P. (2000) The colonization history of Olea europaea L. in Macaronesia based on internal transcribed spacer 1 (ITS-1) sequences, randomly amplified polymorphic DNAs (RAPD), and intersimple sequence repeats (ISSR). Molecular Ecology, 9, 857-868. Hisheh, S., Westerman, M. & Schmitt, L.H. (1998) Biogeography of the Indonesian archipelago: mitochondrial DNA variation in the fruit bat, Eonycteris spelaea. Biological Journal of the Linnean Society, 65, 25.

(33) 329-345. Honda, M., Okamoto, T., Hikida, T. & Ota, H. (2008) Molecular phylogeography of the endemic five-lined skink (Plestiodon marginatus) (Reptilia: Scincidae) of the Ryukyu Archipelago, Japan, with special reference to the relationship of a Northern Tokara population. Pacific Science, 62, 351-362. Hugall, A.F., Foster, R. & Lee, M.S. (2007) Calibration choice, rate smoothing, and the pattern of Tetrapod diversification according to the long nuclear gene RAG-1. Systematic Biology, 56, 543-563. Jeanmougin, F., Thompson, J.D., Gouy, M., Higgins, D.G. & Gibson, T.J. (1998) Multiple sequence alignment with Clustal X. Trends in Biochemical Sciences, 23, 403-405. Juan, C., Oromi, P. & Hewitt, G.M. (1996) Phylogeny of the genus Hegeter (Tenebrionidae, Coleoptera) and its colonization of the Canary Islands deduced from Cytochrome Oxidase I mitochondrial DNA sequences. Heredity, 76, 392-403. Juan, C., Ibrahim, K.M., Oromi, P. & Hewitt, G.M. (1998) The phylogeography of the darkling beetle, Hegeter politus, in the eastern Canary Islands. Proceedings of the Royal Society B: Biological Sciences, 265, 135-140. Kimura, M. (2000) Paleogeography of the Ryukyu Islands. Tropics, 10, 5-24. Lallemand, S., Font, Y., Bijwaard, H. & Kao, H. (2001) New insights on 3-D plates interaction near Taiwan from tomography and tectonic implications. Tectonophysics, 335, 229-253. 26.

(34) Lin, S.M., Chen, C.A. & Lue, K.Y. (2002) Molecular phylogeny and biogeography of the grass lizards genus Takydromus (Reptilia: Lacertidae) of East Asia. Molecular Phylogenetics and Evolution, 22, 276-288. Macey, J.R., Kuehl, J.V., Larson, A., Robinson, M.D., Ugurtas, I.H., Ananjeva, N.B., Rahman, H., Javed, H.I., Osman, R.M., Doumma, A. & Papenfuss, T.J. (2008) Socotra Island the forgotten fragment of Gondwana: unmasking chameleon lizard history with complete mitochondrial genomic data. Molecular Phylogenetics and Evolution, 49, 1015-1018. Matsui, M., Ito, H., Shimada, T., Ota, H., Saidapur, S.K., Khonsue, W., Tanaka-Ueno, T. & Wu, G.F. (2005) Taxonomic relationships within the Pan-Oriental narrow-mouth toad Microhyla ornata as revealed by mtDNA analysis (Amphbia, Anura, Microhylidae). Zoological Science, 22, 489-495. Millien-Parra, V. & Jaeger, J.J. (1999) Island biogeography of the Japanese terrestrial mammal assemblages: an example of a relict fauna. Journal of Biogeography, 26, 959-972. Molinari, R.L., Atwood, D.K., Duckett, D., Spillane, M. & Brooks, I. (1979) Surface currents in the Caribbean Sea as deduced from satellite tracked drifting buoys. Proceedings of the Gulf Caribbean Fisheries Institute, 1979, 106-113. Morrone, J.J. (2009) Evolutionary biogeography: an integrative approach with case studies. Columbia University Press, New York. Nakamura, K., Suwa, R., Denda, T. & Yokota, M. (2009) Geohistorical and 27.

(35) current environmental influences on floristic differentiation in the Ryukyu Archipelago, Japan. Journal of Biogeography, 36, 919-928. Nylander, J.A., Olsson, U., Alstrom, P. & Sanmartín, I. (2008) Accounting for phylogenetic uncertainty in biogeography: a Bayesian approach to dispersal-vicariance analysis of the thrushes (Aves: Turdus). Systematic Biology, 57, 257-268. Ohdachi, S., Dokuchaev, N.E., Hasegawa, M. & Masuda, R. (2001) Intraspecific phylogeny and geographical variation of six species of northeastern Asiatic Sorex shrews based on the mitochondrial cytochrome b sequences. Molecular Ecology, 10, 2199-2213. Okajima, Y. & Kumazawa, Y. (2010) Mitochondrial genomes of acrodont lizards: timing of gene rearrangements and phylogenetic and biogeographic implications. BMC Evolutionary Biology, 10, 141. Ota, H. (1991) Taxonomic redefinition of Japalura swinhonis. nther. (Agamidae: Squamata), with a description of a new subspecies of J. polygonata from Taiwan. Herpetologica, 47, 280-294. 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. Researches on Population Ecology, 40, 189-204. Ota, H. (2000) Current status of the threatened amphibians and reptiles of Japan. Population Ecology, 42, 5-9. Ota, H. (2003) A new subspecies of the Agamid lizard, Japalura polygonata. (Hallowell,. 1861). (Reptilia:. Squamata),. from. Yonagunijima island of the Yaeyama group, Ryukyu archipelago. 28.

(36) Current Herpetology, 22, 61-71. Ota, H., Honda, M., Chen, S.L., Hikida, T., Panha, S., Oh, H.S. & Matsui, M. (2002) Phylogenetic relationships, taxonomy, character evolution and biogeography of the lacertid lizards of the genus Takydromus (Reptilia: Squamata): a molecular perspective. Biological Journal of the Linnean Society, 76, 493-509. Phillimore, A.B. & Owens, I.P.F. (2006) Are subspecies useful in evolutionary and conservation biology? Proceedings of the Royal Society B: Biological Sciences, 273, 1049-1053. Posada, D. (2008) jModelTest: phylogenetic model averaging. Molecular Biology and Evolution, 25, 1253-1256. Poulakakis, N., Lymberakis, P., Antoniou, A., Chalkia, D., Zouros, E., Mylonas, M. & Valakos, E. (2003) Molecular phylogeny and biogeography of the wall-lizard Podarcis erhardii (Squamata: Lacertidae). Molecular Phylogenetics and Evolution, 28, 38-46. Ree, R.H. & Smith, S.A. (2008) Maximum likelihood inference of geographic range evolution by dispersal, local extinction, and cladogenesis. Systematic Biology, 57, 4-14. Ree, R.H., Moore, B.R., Webb, C.O. & Donoghue, M.J. (2005) A likelihood framework for inferring the evolution of geographic range on phylogenetic trees. Evolution, 59, 2299-2311. Robinson, M.M., Dowsett, H.J. & Chandler, M.A. (2008) Pliocene role in assessing future climate impacts. Eos, Transactions American Geophysical Union, 89, 501-502. Ronquist, F. (1996) DIVA, version 1.1. Computer program and manual 29.

(37) available by anonymous FTP from Uppsala University, Uppsala, Sweden. Available at: ftp.uu.se or ftp.systbot.uu.se. Ronquist, F. (1997) Dispersal-vicariance analysis: a new approach to the quantification of historical biogeography. Systematic Biology, 46, 195-203. Ronquist, F. & Huelsenbeck, J.P. (2003) MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics, 19, 1572-1574. Salzburger, W., Ewing, G.B. & von Haeseler, A. (2011) The performance of phylogenetic algorithms in estimating haplotype genealogies with migration. Molecular Ecology, 20, 1952-1963. Shang, G. (2008) A Field Guide to Lizards in Taiwan. Bookzone Ltd., Taiwan. Silvestro, D. & Michalak, I. (2012) raxmlGUI: a graphical front-end for RAxML. Organisms Diversity & Evolution, 12, 335-337. Stamatakis, A. (2006) RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics, 22, 2688-2690. Swofford, D.L. (2002) Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4. Sinauer Associates, Sunderland, Massachusetts. Takada, M., Tachihara, K., Kon, T., Yamamoto, G., Iguchi, K., Miya, M. & Nishida, M. (2010) Biogeography and evolution of the Carassius auratus-complex in East Asia. BMC Evolutionary Biology, 10, 7. Thiel, M. & Haye, P.A. (2006) The ecology of rafting in the marine environment. III. Biogeographical and evolutionary consequences. 30.

(38) Oceanography and Marine Biology, 44, 323-429. Thorpe, R.S., McGregor, D.P., Cumming, A.M. & Jordan, W.C. (1994) DNA evolution and colonization sequence of island lizards in relation to geological history: mtDNA RFLP, cytochrome b, cytochrome oxidase, 12s rRNA sequence, and nuclear RAPD analysis. Evolution, 48, 230-240. Toda, M., Hikida, T. & Ota, H. (2001) Discovery of sympatric cryptic species within Gekko hokouensis (Gekkonidae: Squamata) from the Okinawa Islands, Japan, by use of allozyme data. Zoologica Scripta, 30, 1-11. Toda, M., Nishida, M., Matsui, M., Wu, G.F. & Ota, H. (1997) AIIozyme variation among East Asian populations of the Indian rice frog, Rana limnocharis (Amphibia: Anura). Biochemichal Systematics and Ecology, 25, 143-159. Tominaga, A., Ota, H. & Matsui, M. (2010) Phylogeny and phylogeography of the sword-tailed newt, Cynops ensicauda (Amphibia: Caudata), as revealed by nucleotide sequences of mitochondrial DNA. Molecular Phylogenetics and Evolution, 54, 910-921. Townsend, T.M., Alegre, R.E., Kelley, S.T., Wiens, J.J. & Reeder, T.W. (2008) Rapid development of multiple nuclear loci for phylogenetic analysis using genomic resources: an example from squamate reptiles. Molecular Phylogenetics and Evolution, 47, 129-142. Tursi, R.M., Hughes, P.T. & Hoffman, E.A. (2012) Taxonomy versus phylogeny: evolutionary history of marsh rabbits without hopping to 31.

(39) conclusions.. Diversity. and. Distributions,. DOI:. 10.1111/j.1472-4642.2012.00915.x. Uchiyama, R., Maeda, N., Numata, K. & Seki, S. (2009) A Photographic Guide; Amphibians and Reptiles in Japan. Heibonsha Ltd., Japan. Uetz,. P.. (2006). The. Reptile. Database.. Available. at:. http://www.reptile-database.org (accessed May 23 2006). van Denburgh, J. (1912) Concerning certain species of reptiles and amphibians from China, Japan, the Loo Choo Islands, and Formosa. Proceedings of the California Academy of Sciences, 3, 187-258. Vences, M., Vieites, D.R., Glaw, F., Brinkmann, H., Kosuch, J., Veith, M. & Meyer, A. (2003) Multiple overseas dispersal in amphibians. Proceedings of the Royal Society B: Biological Sciences, 270, 2435-2442. Wallace, A.R. (1895) Island Life, 2nd edn. Macmillan, London. Yu,. Y.,. Harris,. A.J.. &. He,. X.. (2010). S-DIVA. (Statistical. Dispersal-Vicariance Analysis): A tool for inferring biogeographic histories. Molecular Phylogenetics and Evolution, 56, 848-850. Yu, Y., Harris, A.J. & He, X. (2011) RASP: reconstruct ancestral state in phylogenies. v.2.1. beta.. http://mnh.scu.edu.cn/soft/blog/RASP.. 32. Available. at:.

(40) Table 1. Sample localities, localities abbreviation (Abb), GPS coordination, number of sequences (N), number of haplotypes (H), haplotype diversity (h), nucleotide diversity (π), Tajima’s D, Fu’s Fs, and clade number (Cno) for the Japalura polygonata populations from mitochondrial cytochrome b sequences. Subspecies. Locality. Abb. GPS coordination. N. H. h. π. Tajima’s D. Fu’s Fs. J. p. polygonata. Amami Group. A. -. 12. 7. 0.909. 0.00467. 1.38936. 0.12472. Amamioshima. AM. 28.278232, 129.425574. 6. 4. 0.867. 0.00106. -0.44736. -1.45444*. I. Tokunoshima. TKN. 27.825605, 128.931667. 6. 3. 0.733. 0.00077. -0.05002. -0.42679. I. Okinawa Group. O. -. 41. 27. 0.973. 0.01349. 0.66963. -4.03105. North Okinawajima. ONN. 26.733281, 128.190334. 11. 10. 0.982. 0.00437. -0.66252. -4.56515*. I. South Okinawajima. ONS. 26.171679, 127.822033. 8. 3. 0.607. 0.00069. 0.06935. -0.22360. I. *. I I. J. p. ishigakiensis. Cno. Iheyajima. IH. 27.052009, 127.974468. 7. 5. 0.905. 0.00143. -1.02379. -2.01916. Kumejima. KM. 26.379277, 126.763192. 15. 9. 0.924. 0.00305. -0.25643. -2.31464. Miyako Group. M. -. 12. 3. 0.318. 0.00056. -1.17901. -0.18049. Miyakojima. MK. 24.799842, 125.318561. 12. 3. 0.318. 0.00056. -1.17901. -0.18049. Yaeyama Group. Y. -. 52. 30. 0.971. 0.04033. 2.13477. 3.04818. Ishigakijima. IG. 24.453392, 124.200583. 28. 14. 0.931. 0.00176. -1.17817. -8.94940***. III. ***. IV. Iriomotejima. IM. 24.395731, 123.824815. 13. 9. 0.923. 0.00140. -1.48471. J. p. donan. Yonagunijima. YG. 24.446451, 122.979414. 11. 7. 0.873. 0.00254. -0.67508. -1.79007. J. p. xanthostoma. Taiwan Island. T. -. 129. 55. 0.968. 0.05207. 1.19267. 5.67282. Taiwan (North). NT. -. 74. 30. 0.949. 0.01369. -1.38650. -0.83074. Gongzihliao. KL. 25.145869, 121.783804. 13. 4. 0.679. 0.00181. 0.22498. 1.12800. VIII. Yangmingshan. YM. 25.136224, 121.540165. 18. 11. 0.935. 0.00500. -1.21860. -1.75920. VIII. Hushan. XY. 25.029305, 121.583618. 11. 3. 0.618. 0.00594. 2.25711. 6.61003. VIII. 33. -6.39246. II. V.

(41) Xindian. XD. 24.933822, 121.536696. 9. 1. 0.000. 0.00000. NA *. NA. VIII. 2.52818. VIII. Wulai. WL. 24.838942, 121.528620. 10. 7. 0.911. 0.01722. -1.84698. Yingge. JH. 24.972859, 121.348305. 11. 3. 0.564. 0.00199. 1.22875. 2.47881. VIII. Pingxi. PS. 25.023586, 121.721210. 1. 1. NA. NA. NA. NA. VIII. Ruifang. HT. 25.092143, 121.833778. 1. 1. NA. NA. NA. NA. VIII. Taiwan (East). ET. -. 33. 15. 0.813. 0.00810. 0.85309. -0.08656. Jiaoxi. JS. 24.807719, 121.711112. 2. 2. 1.000. 0.00088. 0.00000 *. 0.00000. VI. -2.38086. VI. Dong’ao. TA. 24.482977, 121.839255. 13. 9. 0.910. 0.00365. -2.01914. Gulu. DT. 24.605150, 121.684661. 2. 2. 1.000. 0.00353. 0.00000. 1.38629. VI. Nanshan. NS. 24.406296, 121.360187. 1. 1. NA. NA. NA. NA. VI. Fuxing. FX. 24.650647, 121.427288. 15. 3. 0.257. 0.00024. -1.49051. -1.54636*. VI. Taiwan (West). WT. -. 22. 11. 0.896. 0.02509. -0.37051. 6.39697. Qingquan. XJ. 24.575833, 121.102696. 6. 4. 0.800. 0.00200. -0.78648. -0.27174. VII. Simaxian. ML. 24.410883, 120.946013. 1. 1. NA. NA. NA. NA. VII. Daxueshan. HP. 24.235605, 120.903942. 13. 5. 0.744. 0.00125. -0.99929. -0.90677. VII. Lienhuachih. NT. 23.920337, 120.884527. 2. 1. 0.000. 0.00000. NA. NA. VII. Note: NA, not available; *, p < 0.05; **, p < 0.01; ***, p < 0.001. 34.

(42) Table 2. PCR primer of cytochrome b, 16S rRNA, Bach-1 and Rag-1 sequences developed for this study. Gene. Primers. Annealing Temp.. Sequences. Aligned length (bp). Citations. 54℃. 5’ ACCCCT ACTT AAAAAYC 3’ 5’ AYTY AAAAAYCA C TT 3’ 5’ CCC TCTTTASTTTACAA TC 3’ 5’ AATA TA A ATAA TTT 3’. 1132. This study This study This study This study. Ja16SF1 Ja16SR1. 52℃. 5’ AC CAATA A AAA TACT C 3’ 5’ A ATATA ACC ACCT ATT 3’. 1179. This study This study. Bach-1. BACH1F1 BACH1R2. 59℃. 5’ ATTT AHCCYTTRCTTCA TTT C 3’ 5’ ACCTCACATTCYT TTCYCTR C 3’. 1235. Townsend et al. (2008) Townsend et al. (2008). Rag-1. JaRAG1F1 JaRAG1R1. 57℃. 5’ CCA TTCAT ACCA TA AT A 3’ 5’ ACTAACATCTCCCATTCCATCAC 3’. 1330. This study This study. Cyt b. 16S. JaCBF1 JaCBF2 JaCBR1n JaCBRinter. 35.

(43) Table 3. Sample localities, localities abbreviation (Abb), number of sequences (N), number of haplotypes (H), haplotype diversity (h), nucleotide diversity (π), Tajima’s D, Fu’s Fs, and clade number (Cno) for the Japalura polygonata populations from nuclear Bach-1 and Rag-1 phased-sequences.. Subspecies. Locality. Abb. J. p. polygonata. Amami Group. H. h. π. Tajima’s D. Fu’s Fs. N. H. h. π. Tajima’s D. Fu’s Fs. A. 12. 1. 0. 0.00000. 0.00000. NA. 12. 1. 0.000. 0.00000. 0.00000. NA. Amamioshima. AM. 6. 1. 0.000. 0.00000. 0.00000. NA. 6. 1. 0.000. 0.00000. 0.00000. NA. Tokunoshima. TKN. 6. 1. 0.000. 0.00000. 0.00000. NA. 6. 1. 0.000. 0.00000. 0.00000. J. p. xanthostoma. O. 24. 8. 0.656. 0.00082. -1.44001. -4.45834. I. NA *. Cno. I. ***. 24. 5. 0.493. 0.00048. -1.51776. *. 6. 1. 0.000. 0.00000. 0.00000. NA. I. -2.07552. *. North Okinawajima. ONN. 6. 3. 0.600. 0.00054. -1.13197. -0.85842. South Okinawajima. ONS. 6. 3. 0.733. 0.00076. 0.31063. -0.30414. 6. 3. 0.600. 0.00075. -1.23311. -0.18945. I. Iheyajima. IH. 6. 3. 0.733. 0.00070. -0.05002. -0.42679. 6. 2. 0.600. 0.00045. 1.44510. 0.79518. I. Kumejima. KM. 6. 2. 0.333. 0.00054. -1.13197. 0.95213. I. 6. 2. 0.533. 0.00040. 0.85057. 0.62543. *. 6. 3. 0.733. 0.00085. 1.39259. 0.02028. 6. 3. 0.733. 0.00085. 1.39259. 0.02028. Miyako Group. M. 6. 5. 0.933. 0.00189. 0.36689. -1.78570. Miyakojima. MK. 6. 5. 0.933. 0.00189. 0.36689. -1.78570. Yaeyama Group. J. p. donan. Rag-1. N. Okinawa Group. J. p. ishigakiensis. Bach-1. Y. 18. 13. 0.948. 0.00553. -1.16648. -3.10507. 18. 11. 0.889. -4.81440. II. **. 0.00211. -1.51663 *. -0.49899. III. Ishigakijima. IG. 6. 6. 1.000. 0.00480. -0.19835. -1.83372. 6. 4. 0.800. 0.00150. -1.36732. Iriomotejima. IM. 6. 5. 0.933. 0.00691. -0.45176. 0.55466. 6. 4. 0.800. 0.00150. -1.36732. -0.49899. IV. Yonagunijima. YG. 6. 3. 0.600. 0.00081. -1.23311. -0.18945. 6. 4. 0.867. 0.00105. 0.33839. -1.15958. V. *. -11.73308. Taiwan Island. T. 54. 29. 0.951. 0.00525. -1.43702. Taiwan (North). NT. 24. 15. 0.953. 0.00482. -0.66294. Gongzihliao. KL. 4. 3. 0.833. 0.00324. Yangmingshan. YM. 4. 4. 1.000. Hushan. XY. 4. 3. 0.833. **. 54. 33. 0.972. 0.00531. -1.07479. -3.81893. 24. 18. 0.978. 0.00514. -0.66995. -6.80435**. 0.46705. 1.16315. 4. 3. 0.833. 0.00414. 0.08338. 1.65464. VIII. 0.00513. -0.32685. -0.21932. 4. 4. 1.000. 0.00639. -0.26532. 0.14168. VIII. 0.00297. 1.16799. 1.03139. 4. 3. 0.833. 0.00752. 0.18932. 2.62417. VIII. 36. -16.36062. ***.

(44) Xindian. XD. 4. 2. 0.500. 0.00081. -0.70990. 1.09861. 4. 3. 0.833. 0.00376. 1.44213. 1.50590. VIII *. Wulai. WL. 6. 3. 0.600. 0.00194. -0.49605. 1.36614. 6. 6. 1.000. 0.00311. -0.33767. -2.55169. Yingge. JH. 2. 2. 1.000. 0.00972. 0.00000. 2.48491. 2. 2. 1.000. 0.00150. 0.00000. 0.69315. VIII. Pingxi. PS. 0. NA. NA. NA. NA. NA. 0. NA. NA. NA. NA. NA. VIII. Ruifang. HT. 0. NA. NA. NA. NA. NA. 0. NA. NA. NA. NA. NA. VIII. Taiwan (East). ET. 16. 8. 0.842. 0.00256. -0.48384. -1.32336. 16. 9. 0.892. 0.00254. -0.99265. -2.12164. Jiaoxi. JS. 4. 2. 0.500. 0.00202. -0.79684. 2.59805. 4. 4. 1.000. 0.00388. -0.52807. -0.48049. VI. Dong’ao. TA. 4. 3. 0.833. 0.00162. -0.78012. 0.13353. 4. 3. 0.833. 0.00188. -0.79684. 0.46110. VI. Gulu. DT. 4. 2. 0.500. 0.00081. -0.70990. 1.09861. 4. 3. 0.833. 0.00150. -0.78012. 0.13353. VI. Nanshan. NS. 0. NA. NA. NA. NA. NA. 0. NA. NA. NA. NA. NA. VI. Fuxing. FX. 4. 3. 0.833. 0.00243. 0.95621. 0.73089. 4. 3. 0.833. 0.00138. 1.08976. 0.00617. VI. -2.26869. 14. 8. 0.901. 0.00260. -1.27999. -1.49169. *. VIII. Taiwan (West). WT. 14. 9. 0.912. 0.00310. -1.76581. Qingquan. XJ. 4. 3. 0.833. 0.00148. 1.08976. 0.00617. 4. 3. 0.833. 0.00075. -0.70990. -0.88730. VII. Simaxian. ML. 2. 2. 1.000. 0.01053. 0.00000. 2.56495. 2. 1. 0.000. 0.00000. 0.00000. NA. VII. Daxueshan. HP. 4. 4. 1.000. 0.00189. 0.65010. -1.62218*. 4. 4. 1.000. 0.00539. -0.62393. -0.06549. VII. Lienhuachih. NT. 4. 4. 1.000. 0.00216. -0.21249. -1.41422. 4. 1. 0.000. 0.00000. 0.00000. NA. VII. Note: NA, not available; *, p < 0.05; **, p < 0.01; ***, p < 0.001. 37.

(45) Table 4. Length after alignments, variable and parsimony-informative sites, and GC content of the 4 gene sequences from Japalura polygonata.. Aligned length (bp) Variable sites Parsimony-informative sites Mean % G+C content. Cyt b. 16S. Bach-1. Rag-1. 1132 317 284 43.5. 1179 136 130 44.4. 1235 31 14 43.3. 1330 21 14 45.7. 38.

(46) Table 5. Pairwise distance comparison of Japalura polygonata populations. Upper: 16S rRNA. Lower: cytochrome b. AM AM. ONN. TKN. ONS. IH. KM. MK. IM. IG. ET. NT. YG. WT. 0.0019. 0.0049. 0.0040. 0.0079. 0.0063. 0.0118. 0.0094. 0.0220. 0.0401. 0.0340. 0.0423. 0.0396. 0.0060. 0.0051. 0.0069. 0.0074. 0.0127. 0.0106. 0.0231. 0.0395. 0.0352. 0.0414. 0.0390. 0.0009. 0.0068. 0.0049. 0.0106. 0.0097. 0.0188. 0.0369. 0.0326. 0.0411. 0.0372. 0.0060. 0.0040. 0.0097. 0.0088. 0.0180. 0.0361. 0.0317. 0.0403. 0.0363. 0.0083. 0.0153. 0.0148. 0.0240. 0.0403. 0.0385. 0.0439. 0.0394. 0.0120. 0.0111. 0.0220. 0.0364. 0.0348. 0.0420. 0.0381. 0.0066. 0.0199. 0.0380. 0.0303. 0.0405. 0.0399. 0.0191. 0.0366. 0.0294. 0.0391. 0.0388. 0.0360. 0.0242. 0.0350. 0.0358. 0.0401. 0.0404. 0.0400. 0.0412. 0.0346. ONN. 0.0056. TKN. 0.0078. 0.0082. ONS. 0.0096. 0.0096. 0.0044. IH. 0.0129. 0.0136. 0.0113. 0.0131. KM. 0.0198. 0.0194. 0.0176. 0.0195. 0.0217. MK. 0.0233. 0.0226. 0.0216. 0.0225. 0.0268. 0.0283. IM. 0.0243. 0.0252. 0.0233. 0.0250. 0.0270. 0.0281. 0.0295. IG. 0.0533. 0.0538. 0.0501. 0.0524. 0.0551. 0.0575. 0.0582. 0.0552. ET. 0.0579. 0.0594. 0.0596. 0.0605. 0.0607. 0.0659. 0.0644. 0.0624. 0.0600. NT. 0.0706. 0.0704. 0.0704. 0.0711. 0.0738. 0.0741. 0.0697. 0.0737. 0.0689. 0.0748. YG. 0.0729. 0.0718. 0.0711. 0.0709. 0.0753. 0.0774. 0.0699. 0.0752. 0.0705. 0.0742. 0.0790. WT. 0.0803. 0.0805. 0.0798. 0.0803. 0.0805. 0.0857. 0.0832. 0.0846. 0.0814. 0.0851. 0.0864. 0.0473 0.0942. Note: AM: Amamioshima, ONN: Okinawajima (North), TKN: Tokunoshima, ONS: Okinawajima (South), IH: Iheyajima, KM: Kumejima, MK: Miyakojima, IM: Iriomotejima, IG: Ishigakijima, ET: eastern Taiwan, NT: northern Taiwan, YG: Yonagunijima, WT: western Taiwan.. 39.

(47) Table 6. Mean values and 95% highest posterior density (HPD) intervals of divergence time estimated by relaxed molecular clock analysis using concatenated sequences of cytochrome b and 16S in BEAST. Serial numbers of nodes are defined in Figure 5. Node. Crown age (Ma). 95% HPD interval. 1 2 3 4 5. 4.56 4.18 3.07 1.32 1.14. 3.61-5.57 3.20-5.19 2.19-4.03 0.91-1.80 0.79-1.56. 6 7 8 9 10 11 12 13 14. 0.76 0.57 0.47 0.23 0.08 0.10 0.22 0.07 0.05. 0.52-1.05 0.38-0.80 0.29-0.67 0.12-0.36 0.02-0.17 0.03-0.19 0.09-0.37 0.01-0.15 0.01-0.11. 15 16 17 18 19 20 21 22 23. 0.09 0.17 0.05 0.08 0.11 3.67 0.22 0.60 4.14. 0.02-0.18 0.07-0.31 0.01-0.12 0.02-0.19 0.03-0.24 2.71-4.69 0.10-0.39 0.35-0.94 3.17-5.14. 24 25. 0.86 3.23. 0.53-1.25 2.36-4.15. 40.

(48) Table 7. Results for the biogeographic method tested in DEC and S-DIVA models. Statistics of DEC models include global log-likelihood units (lnL), estimates of dispersal (D), extinction (E) rates (events per million year) and maximum-likelihood scenarios of range inheritance for the eight nodes. If a node has multiple scenarios within 2 log-likelihood units of the optimal reconstruction, the two most likely scenarios are shown, and the relative probability of each is indicated in brackets. Statistics of the two most likely scenarios of eight nodes from S-DIVA analysis are shown, with marginal probability listed in brackets.. Method Constraints DEC. M0. M1. S-DIVA. Max areas = 5. lnL -28.24. -28.34. -. D 0.070. 0.091. -. E 0.000. 0.000. -. The reconstructions for ancestral taxon of the following nodes: Node 1 (root) a. Node 2. Node 3. Node 4. Node 5. Node 6. Node 7. Node 8. YT∣T (0.60),. Y∣YT (0.68),. Y∣Y (0.68),. OM∣Y (0.80),. O∣M (0.91),. O∣O (0.89),. O∣O (0.83),. O∣O (0.66),. Y∣T (0.14). Y∣Y (0.12). MY∣Y (0.18). M∣Y (0.12). AO∣M (0.09). AO∣O (0.11). AO∣O (0.16). AO∣A (0.09). YT∣T (0.91),. Y∣YT (0.93),. Y∣Y (0.94),. OM∣Y (0.77),. O∣M (0.96). O∣O (0.93),. O∣O (0.87),. O∣O (0.69),. Y∣T (0.03). Y∣Y (0.03). MY∣Y (0.04). M∣Y (0.20). AO∣O (0.07). AO∣O (0.12). AO∣A (0.08). YT (35.04),. YT (27.58),. Y (55.13),. OMY (49.98),. OM (99.95),. O (99.95),. O (99.95),. O (99.97),. T (35.04). Y (27.58). OMY (44.83). MY (49.98). AOM (0.05). AO (0.05). AO (0.05). AO (0.03). a. In the eight nodes results, areas shown indicate the range inherited by the daughter lineages. When a bar separates two ranges, the first range is inherited by the upper branch in Figure 6 and the second. range is inherited by the lower branch. M0: area constraint at nodes = 5, dispersal probability = 0.1 constant through all time. M1: area constraint at nodes = 2, dispersal probability = 0.1 constant through all time. M2: area constraint at nodes = 5, dispersal probability = 1 during the Pleistocene, 0.1 before and after Pleistocene. M3: area constraint at nodes = 2, dispersal probability = 1 during the Pleistocene, 0.1 before and after Pleistocene. Note: A: Amami Group, O: Okinawa Group, M: Miyako Group, Y: Yaeyama Group, T: Taiwan Island.. 41.

(49) Figure 1A. Sample localities (black circles) of Japalura polygonata populations used in this study. Pie charts represented the lineage composition of each J. polygonata population. Precise collecting localities are provided in Table 1. The colors used in the pie charts correspond to the clade definition in Figure 2, 3 and 5: Purple: Amami Group and Okinawa Group; Light steelblue: Miyako Group; Light green: Yaeyama Group (Ishigakijima); Light salmon: Yaeyama Group (Iriomotejima); Yellow: Yaeyama Group (Yonagunijima); Blue: Taiwan (North); Red: Taiwan (East); Dark green: Taiwan (West). 42.

(50) Figure 1B. Sample localities of Japalura polygonata populations used in this study. Pie charts represented the lineage composition of each J. polygonata clade on Taiwan Island. Precise collecting localities are provided in Table 1. The colors used in the pie charts correspond to the clade definition in Figure 2, 3 and 5: Blue: Taiwan (North); Red: Taiwan (East); Dark green: Taiwan (West).. 43.

(51) Figure 2. Maximum likelihood (ML) tree based on cytochrome b sequences of Japalura polygonata. The support values of each node are bootstrap confidence (BS) results for maximum parsimony and maximum likelihood analyses, and Bayesian posterior probabilities (PP) for Bayesian inference. BS values lower than 50% and PP values lower than 95% are indicated by “-”. Asterisks (*) indicate topological incongruence between MP, ML and BI topologies. Eight defined clades with associated distributional areas are listed beside the tree. The four currently recognized subspecies are list on the right. 44.

(52) Figure 3. Maximum likelihood tree (ML) based on cytochrome b and 16S concatenated sequences of Japalura polygonata. The support values of each node are bootstrap confidence (BS) results for maximum parsimony and maximum likelihood analyses, and Bayesian posterior probabilities (PP) for Bayesian inference analysis. BS values lower than 50% and PP values lower than 95% are indicated by “-”. Asterisks (*) indicate topological incongruence between MP, ML and BI topologies. Eight defined clades with associated distributional areas are listed beside the tree. The four currently recognized subspecies are list on the right. 45.

(53) Figure 4. Genotype network of Japalura polygonata using two nuclear genes (Bach-1 and Rag-1), concatenated and unphased sequences. Lines represent single nucleotide substitutions; dots indicate absent haplotypes. Dash line groups all haplotypes into four major groups (a) to (d). A: Amami Group, O: Okinawa Group, M: Miyako Group, Y: Yaeyama Group, T: Taiwan Island. 46.