The Kuroiwa's eyelid gecko Goniurosaurus kuroiwae is an endangered species in a state of relict endemism in the Central Ryukyus, Japan, and is divided into five subspecies. We analyzed variations in sequence data for approximately 1900 base positions of mitochondrial 12S and 16S rRNA, and cytochrome b genes from samples representing all recognized subspecies of G. kuroiwae together with those from congeneric species in order to test the relevant previous phylogenetic hypotheses and discuss biogeographical implications in the degree and pattern of genetic divergence within G. kuroiwae. Our results, while confirming a previous molecular phylogenetic hypothesis proposed on the basis of much smaller data set, negate the relationships hypothesized on morphological grounds by explicitly supporting: 1) the primary dichotomy, with substantial genetic divergence, between G. k. splendens from the Amami Island Group and the remaining subspecies all from the Okinawa Island Group; and 2) the presence of at least six independent lineages within the latter, indicating non-monophyly for two of the subspecies, G. k. kuroiwae and G. k. orientalis, in the current taxonomic definitions. The marked genetic divergence between populations of the two island groups seems to have initiated in the middle Miocene, i.e., prior to formation of straits that have consistently been separating these two island groups since the early Pleistocene. All populations of G. kuroiwae are regarded as endangered from the viewpoint of conservation genetics.
INTRODUCTION
The Ryukyu Archipelago is composed of approximately 150 subtropical islands ranging between Kyushu of Japanese main islands and Taiwan (Fig. 1). The fauna of this region is characterized by a high ratio of endemic taxa and genetically diverged populations, most of which are thought to have been isolated from their relatives as a result of insularization (e.g., Ikehara, 1996; Ota, 1998, 2000a). The archipelago is biogeographically classified into three regions: Northern, Central, and Southern Ryukyus (Ota, 1998). Of these, the Central Ryukyus region, consisting of the Okinawa and Amami Island Groups, and a few southern islets of the Tokara Island Group, is known to accommodate an especially high proportion of endemic taxa in non-volant terrestrial vertebrates that are in a more or less relict state (Matsui, 1994; Ota, 1998, 2000a; Takahashi et al., 2007).
Generally, reptiles and amphibians are considered ideal materials for biogeographic studies (e.g., Vences et al., 2003), because of their relatively low vagility in the presence of physical (e.g., sea straits, large rivers, or high mountain ranges) or climatic (e.g., temperature or aridity) barriers (Souza et al., 2003; Amato et al., 2008). However, molecular-based phylogenetic and biogeographic studies on reptiles and amphibians in the Central Ryukyus are still meager (e.g., Matsui et al., 2005a, b; Honda et al., 2008). Moreover, most of these studies combined conspecific or consubspecific specimens from each island as a minimum unit for the analyses (Ota et al., 2002; Matsui et al., 2007; Honda et al., 2008; but see Ota et al., 1999), although such an approach may lead to underestimation of the actual genetic diversity of a given taxon on a single island. This seems to be particularly true when we assess genetic diversity of threatened species, because an island assemblage representing an endangered taxon often exhibits various effects of population fragmentation that may give rise to complex within-island genetic structures (e.g., Ota, 2000b; Honda et al., 2012).
The Kuroiwa's eyelid gecko, Goniurosaurus kuroiwae (Namiye, 1912) is an endangered eublepharid species endemic to the Central Ryukyus (Grismer et al., 1994) (Fig. 1). For this species, five subspecies are recognized on the basis of morphological characters: G. k. kuroiwae from Okinawajima Island and three adjacent islets, Sesokojima, Kourijima, and Yagachijima Islands; G. k. orientalis (Maki, 1931) from Tonakijima, Tokashikijima, Akajima, and lejima Islands; G. k. yamashinae (Okada, 1936) from Kumejima Island; G. k. toyamai Grismer et al., 1994 from Iheyajima Island; and G. k. splendens (Nakamura and Uéno, 1959) from Tokunoshima Island (Grismer et al., 1994; Maenosono and Toda, 2007; Kurita and Kawamura, 2011; but see Grismer et al. [1999] for a different taxonomic arrangement). The former four subspecies are distributed in the Okinawa Island Group, whereas the last subspecies exclusively occurs in the Amami Island Group. The other congeners are known from southeastern continental China, Hainan Island, and northern Vietnam, including several offshore islets (Grismer, 1987; Grismer et al., 1999, 2002; Ziegler et al., 2008; Orlov et al., 2008; Wang et al., 2010, 2013). There are considerable geographical and morphological gaps between G. kuroiwae and the other congeneric species, and the former is considered to be in an extremely relict state resulting from an initial continental-insular vicariance and subsequent extinction of ancestral populations that should have once occurred broadly in East Asia, most likely under the effects of paleogeographical and paleoenvironmental dynamics (Tanaka, 1996; Ota, 1998). It has also been suggested that inter-subspecific diversification of G. kuroiwae may be the result of more recent changes in land configuration around the current Central Ryukyus, i.e., formation and fragmentation of landbridges (Grismer et al., 1994; Tanaka, 1996; Ota, 1998, 2003).
Fig. 1.
Maps of East Asia and the Central Ryukyus, showing sampling localities of Goniurosaurus kuroiwae. Sample sizes are given in parentheses, following the numbers of each locality. 1, Tete; 2, Mikyo; 3, Asan; 4, Dana; 5, Koshidake; 6, Gakiya; 7, Yona; 8, Kudeken; 9, lejima A; 10, lejima B; 11, Tokashiki; 12, Aharen; 13, Tonakijima A; 14, Tonakijima B; 15, Uegusuku; 16, Kitahara; and 17, Uezu.
Based on morphological data, Grismer et al. (1994) hypothesized phylogeny among subspecies of G. kuroiwae on an a priori assumption of insular monophyly. Applying the parsimony method to their data set, they argued that G. k. yamashinae from Kumejima Island first diverged from the remainder, and that the latter further branched into two pairs of sister taxa, G. k. kuroiwae and G. k. orientalis, and G. k. toyamai and G. k. splendens. From analyses of mitochondrial DNA sequence data for several insular samples of G. kuroiwae, however, Ota et al. (1999) obtained relationships that were not congruent with the hypothesis of Grismer et al. (1994) at all. Ota et al. (1999) also indicated a possible non-monophyly of G. k. kuroiwae on Okinawajima Island due to the closer affinity of its southern Okinawajima population with the G. k. orientalis population from lejima Island, an islet immediately north of Okinawajima Island, than with the “consubspecific” northern Okinawajima population. This result strongly suggested that the original population of G. kuroiwae may have diversified within Okinawajima Island prior to divergence of its subpopulation with ancestral lejima population currently identified to G. k. orientalis. However, due to insufficient numbers and sizes of samples examined, the Ota et al. (1999) results remain inconclusive, requiring further verification on the basis of analyses using samples from the whole of the species' range.
In this study, we examined geographic pattern and degrees of genetic differentiation in G. kuroiwae by analyzing longer mitochondrial DNA sequence data for samples from most islands and within-island districts, from which the species had been recorded. Our purposes are: 1) to convincingly evaluate the two foregoing phylogenetic hypotheses (i.e., those proposed by Grismer et al. [1994] and by Ota et al. [1999]), and to propose an alternative hypothesis as needed; 2) to discuss historical biogeography of G. kuroiwae on the basis of molecularly estimated divergence times among its local populations; and 3) to examine taxonomic and conservation implications in the results of molecular analyses.
MATERIALS AND METHODS
Samples analyzed
Tissues were obtained from 42 specimens representing all five subspecies of Goniurosaurus kuroiwae (Fig. 1). Their sampling localities, 17 in number, were located on all islands from which the species had been recorded, except for Kourijima, Sesokojima, Yagachijima, Akajima, and Yoronjima Islands. Of the latter five islands, Akajima Island, located close to Tokashikijima Island of the Okinawa Island Group (Fig. 1), was once noted for the discovery of a single specimen identified as G. k. orientalis (see Toyama, 1983), but no additional records have been made to present, most likely due to complete predation of the population by exotic weasels. Likewise, Yoronjima Island, representing the south-westernmost extremity of the Amami Island Group ca. 30 km northeast of Okinawajima Island (Fig. 1), was also reported as an island on which G. kuroiwae occurred until recently, but its population has evidently disappeared during the last several decades without being observed in life (Nakamura et al., 2013). The remaining three islets, each having an apparently small population of G. k. kuroiwae, are located close to Okinawajima Island (Maenosono and Toda, 2007; Kurita and Kawamura, 2011). Because these three islets are separated from the latter only by narrow (< 3 km) and shallow (< 20 m) straits (Fig. 1), failure in inclusion of their samples to the present analyses did not seemingly much affect the results as long as samples from Okinawajima Island were incorporated.
Tail tip and/or one of the digital tips were taken from each individual gecko in the field. Animals were then released at the exact points where they had been captured. We also incorporated representative samples of two Chinese—Vietnamese congeners (G. luii Grismer et al., 1999 and G. huuliensis Orlov et al., 2008) and four other eublepharid genera (Aeluroscalabotes felinus (Günther, 1864), Coleonyx mitratus (Peters, 1863), Eublepharis macularius (Blyth, 1854), and Holodactylus africanus Boettger, 1893) into the phylogenetic analyses as outgroups. Of these, the Chinese-Vietnamese Goniurosaurus are expected to constitute the primary outgroup to G. kuroiwae (see Ota et al., 1999; Jonniaux and Kumazawa, 2008; Wang et al., 2013). Of the 48 individuals examined in the present study, data for 12S and 16S rRNAs of six G. kuroiwae and five outgroups were published in Ota et al. (1999).
DNA amplification, sequencing, and alignments
Extraction, amplification, and sequencing procedures of DNA are described elsewhere (Honda et al., 1999a, b). Parts [approximately 1,900 base pairs (bp) in total] of mitochondrial 12S and 16S rRNA, and cytochrome b genes were amplified using the polymerase chain reaction with the following primers: L1091 (5′-AAACTGGGATTAGATACCCCACTAT-3′) and H1478 (5′-GAGGGTGACGGGCGGTGTGT-3′) (Kocher et al., 1989) for 12S rRNA, L2206 (5′-GGCCTAAAAGCAGCCACCTGTAAAGACAGCGT-3′) and H3056 (5′-CTCCGGTCTGAACTCAGATCACGTAGG-3′) (Hedges et al., 1993) for 16S rRNA, and L14731 (5′-TGGTCTGAAAAACCATTGTTG-3′) (presented here), L14841 (5′-AAAAAGCTTCCATCCAACATCTCAGCATGATGAAA-3′), H15149 (5′-AAACTGCAGCCCCTCAGAATGATATTTGTCCTCA-3′) (Kocher et al., 1989) and H16064 (5′-CTTTGGTTTACAAGAACAATGCTT-3′) (Burbrink et al., 2000) for cytochrome b. Direct sequencing was performed on ABI 310 and 3130 DNA PRISM Sequencers (Applied Biosystems). The 29 new haplotypes of these genes are deposited in GenBank (Accession number: AB853424-AB853483).
Alignments for two rRNA genes were determined based on the maximum nucleotide similarity using Clustal W 2.1 (Larkin et al., 2007) with default gap penalties (gap opening penalty = 10).
Phylogenetic analyses
We performed phylogenetic analyses using maximum likelihood (ML), maximum parsimony (MP), and Bayesian inference (BI) methods. Model selection for ML and BI analyses was performed using the Akaike information Criterion (AIC: Akaike, 1973) in jModeltest 0.1.1 (Posada, 2008). For ML and BI analyses, we also performed a partitioning scheme following recent studies (e.g., Brandley et al., 2005; Wiens et al., 2010). The schemes for the data were single-partition strategy, three-partition strategy by genes, four-partition strategy by codon position of cytochrome b (1st + 2nd vs. 3rd) and rRNAs, and five-partition strategy of the coding genes separately by codon (1st, 2nd vs. 3rd) and rRNAs. These partitioning strategies were assessed using BI as implemented in BEAST 1.7.5 (Drummond et al., 2012) with an uncorrelated lognormal relaxed clock model and Yule tree prior run for 10 million generations and sampled every 1000 generations. The results were assessed using Bayes factors (Kass and Raftery, 1995) calculated in Tracer 1.5 (Rambaut and Drummond, 2007) using log-likelihood scores from the posterior distribution. The five-partition strategy was selected as an optimum for the data set. In this strategy, the general time reversible model (Lanave et al., 1984) with a certain proportion of invariant sites (Gu et al., 1995) and a gamma distribution of changes (Yang, 1994) (GTR + I + G) was selected as the best model for 16S rRNA gene and second position of cytochrome b, and the general time reversible model with a gamma distribution (GTR + G) for 12S rRNA gene and the other positions of cytochrome b.
ML analysis was performed using Treefinder October 2008 (Jobb, 2008) under the models selected in the above process. MP phylogeny was estimated using the heuristic search algorithm for each tree-building methodology using PAUP* 4.0b (Swofford, 2002). We used 100 random-taxon addition replicates for all analyses to minimize the effect of entry sequence on the topology of the resulting cladogram. We conducted the analyses with accelerated character transformation (ACCTRAN) optimization and tree-bisection-reconnection (TBR) with characters unordered and equally weighted.
The confidence of branches in ML and MP trees was assessed using non-parametric bootstrapping (Felsenstein, 1985) with 1000 pseudoreplicates in Treefinder and PAUP*, respectively. Tree topologies with bootstrap proportions (BPs) of 70% or greater were regarded as sufficiently resolved nodes (Huelsenbeck and Hillis, 1993; Shaffer et al., 1997).
BI using the Markov chain Monte Carlo (MCMC) technique was also performed using MrBayes 3.2.2 (Huelsenbeck and Ronquist, 2001). We used the same settings for the nucleotide substitution model selected in jModeltest. We initiated four independent analyses with a random starting tree that ran for 10 million generations. We used the program Tracer to determine when the log likelihood of sampled trees reached stationary distribution. Because apparent stationarity of MCMC runs was reached at no later than one million generations, we conservatively discarded the first two million generations from each run as burn-in and sampled one of every 100 generations from the remaining eight million generations to calculate posterior probability for each branch in the Bayesian tree. Bayesian posterior probabilities (BPPs) of 0.95 or greater were considered significant supports (Larget and Simon, 1999; Huelsenbeck et al., 2001).
To test the monophyly of the subspecies, we used Shimodaira-Hasegawa (SH: Shimodaira and Hasegawa, 1999) and approximately unbiased (AU: Shimodaira, 2002) tests as implemented in Treefinder using the models and partitioning strategy described above. SH test measures whether some trees are better than others under likelihood, whereas site likelihoods for the optimal and constrained phylogenetic trees were estimated in the AU test.
Genetic structure, divergence times and historical demography
We calculated both haplotype diversity (h) and nucleotide diversity values (π) in each of the samples. Population pairwise fixation indices (F ST) were calculated and their significances were tested with a nonparametric permutations approach with 1000 permutations of haplotypes among sampling localities. For comparisons of haplotype frequencies among samples, we conducted an exact test of population differentiation (Raymond and Rousset, 1995). These analyses were performed using Arlequin 3.5 (Excoffier et al., 2005).
The rate of approximately 2.0%/myr is broadly used for molecular clock of vertebrate mitochondrial cytochrome b as shown in bony fishes (e.g., Rocha et al., 2005), amphibians (e.g., Matsui et al., 2007), reptiles (e.g., Carranza and Arnold, 2006), birds (e.g., Fernandes et al., 2013), and mammals (e.g., Rapson et al., 2012). As there was no reliable temporal calibration point directly applicable to our data, we used the rate of 2.0%/myr assuming a molecular clock for cytochrome b. Rough estimation of divergence time based on molecular clocks may give some idea of the absolute date of colonization although their accuracy still remains disputable (e.g., Avise, 2000; Heads, 2005). 12S and 16S rRNA genes were observed to evolve 0.371 and 0.378 times slower than cytochrome b gene in the ingroup portion, respectively. Thus, we used a rate of 0.74%/myr for 12S and 0.76%/myr for 16S rRNA following Guiller and Madec (2010). To obtain Bayesian estimates of the timing of diversification events, the program BEAST was used with a relaxed clock model with a lognormal distribution and Yule process with default parameters. The models of sequence evolution and partition strategy were identical to the above process. The program ran for 100 million generations with sampling taking place every 1000 generations. A burn-in of 20% was applied to obtain the node age estimates.
To examine whether G. kuroiwae populations are at equilibrium, we calculated Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) implemented in Arlequin. We chose these tests because of the increased statistical power in detecting significant changes in population size when using small sample sizes (Ramos-Onsins and Rozas, 2002). Under the assumption of neutrality, negative values are expected in populations that have undergone recent expansion because rare alleles are more numerous than expected. Positive values occur if rare alleles are eliminated from populations following genetic bottlenecks (Tajima, 1989).
Although Tajima's D and Fu's Fs may provide insights regarding whether or not populations have undergone expansion, they do not provide clues to the shape of population growth over time. Both non-significant negative values are indications that the populations in problem have not undergone expansive growth relative to a null hypothesis of population stability. However, such values are agnostic as to whether populations are expanding slowly, contracting, or remaining relatively constant in size (Fontaneila et al., 2008). Therefore, to estimate the shape of population growth through time, we constructed Bayesian skyline plots (BSP) (Drummond et al., 2005) implemented in BEAST. For each BSP, the appropriate substitution model and rate were determined above process. Genealogies and model parameters for each lineage were sampled every 1000 iteration for 10 million generations under a strict molecular clock with uniformly distributed priors and 20% of burn-in. Demographic plots for each analysis were visualized in Tracer.
Morphological analysis
Morphological MP analysis was also conducted for Grismer et al.'s (1994) taxon-character matrix using PAUP*, because Grismer et al. (1994) conducted the cladistic analysis only manually and thus the resultant tree was not tested for statistical supports at all. We used 100 random-taxon addition replicates for all analyses with ACCTRAN and TBR. The confidence of branches in the trees was assessed, using BP with heuristic search of 1000 replicates. Based on the phylogenetic hypothesis available to that date (Grismer, 1987, 1988, 1989, 1991), Grismer et al. (1994) used 10 eublepharid species representing non-Ryukyu Goniurosaurus and three other genera, Eublepharis, Holodactylus, and Hemitheconyx. Of these, we selected the following four species as outgroups for our morphological analysis Goniurosaurus lichtenfelderi (Mocquard, 1897), Eublepharis macularius, Holodactylus africanus and Hemitheconyx caudicinctus (Dumèril, 1851). On the resultant tree, we evaluated evolution of the 10 characters by comparing their states in the ingroup and outgroup taxa manually. Terminology to refer to the morphological characters and character states follows Grismer (1988).
RESULTS
Sequence variation and divergence
For all 42 individuals of Goniurosaurus kuroiwae, we determined a total 1942 base positions of 12S and 16S rRNA, and cytochrome b genes. The fragment of these genes consisted of 402, 466, 1074 bp, of which 38, 42, and 269 sites were polymorphic in the species, respectively. Combined sequences of these three genes revealed a total of 23 unique haplotypes in the species.
There was a relatively large divergence between the Ryukyu and Chinese—Vietnamese species of Goniurosaurus in total (uncorrected p-distance = 20.1%). Inter-specific nucleotide divergence between the latter was 2.7%, whereas the divergences between island populations of G. kuroiwae ranged from 1.1% to 11.6%. The mean of the uncorrected p-distances was 10.8% (5.9%, 6.0%, and 14.6% in 12S and 16S rRNAs, and cytochrome b, respectively) between G. k. splendens from the Amami Island Group and other subspecific populations of G. kuroiwae from the Okinawa Island Group, 4.1–6.4% in pairwise comparisons of the latter, and 1.1–1.2% between G. k. orientalis from lejima Island and G. k. kuroiwae from the southern part of Okinawajima Island. Samples from Tonakijima, Iheyajima, Tokashikijima, and lejima Islands completely or almost completely lacked within-island haplotype and nucleotide diversity (Table 1).
Table 1.
Sample size (N), number of haplotype (Nh), number of polymorphic site (Np), haplotype diversity (h), nucleotide diversity (π), Tajima's D, Fu's Fs, and haplotypes (sampling locality (S) (haplotype [number of haplotype]) for samples of G. kuroiwae. Asterisk (*) indicates significant difference between samples (P < 0.05). See Figs. 1 and 2 for sampling localities and haplotypes, respectively.
The results of the tests for population differentiation are shown in Table 2. Both F ST and the result of the exact test showed significant population differentiation except for those among G. k. orientalis from lejima Island and G. k. kuroiwae.
Phylogenetic relationships
Of the 1942 sites in a total 29 haplotypes including those from the outgroup taxa, 924 were variable and 672 parsimony informative. ML analysis generated a topology with InL = -10,258. MP analysis yielded six equally most parsimonious trees (tree length = 2080, consistency index = 0.65, homoplasy index = 0.36, and retention index = 0.75). The mean InL score of BI was -10,062.
The ML tree (Fig. 2) was almost identical to the MP and BI trees (not shown) and there were no inconsistency among the three analyses in topology of ingroup nodes with significant supports (i.e., BP ≥ 70% or BPP ≥ 0.95). In all three analyses, monophyly was strongly supported for both an assemblage of the Chinese—Vietnamese congeners, and G. kuroiwae (100%, 100%, and 1.00 supports in MP and ML BPs, and in BPP, respectively), although monophyly of the genus Goniurosaurus did not get sufficient supports (< 70%, < 70%, and < 0.95).
Within the ingroup portion, G. kuroiwae was divided into two sufficiently supported groups in all analyses, one solely consisting of G. k. splendens from Tokunoshima Island of the Amami Island Group (100%, 100%, and 1.00), and the other of the remaining subspecies all from the Okinawa Island Group (99%, 100%, and 1.00). Of the latter, two of the four morphologically defined subspecies, G. k. kuroiwae and G. k. orientalis, failed to compose monophyletic groups (Fig. 2). For each of these subspecies monophyly was invariably rejected by the SH and AU tests (P < 0.001 and P < 0.001, respectively), thus strongly suggesting the presence of at least six lineages as follows: 1) G. k. toyamai from Iheyajima Island (99%, 100%, and 1.00), 2) G. k. kuroiwae from the northern part of Okinawajima Island (93%, 100%, and 1.00), 3) an assemblage of G. k. orientalis from lejima Island and G. k. kuroiwae from the southern part of Okinawajima Island (100%, 100%, and 1.00), 4) G. k. orientalis from Tokashikijima Island (100%, 100%, and 1.00), 5) G. k. orientalis from Tonakijima Island (100%, 100%, and 1.00), and 6) G. k. yamashinae from Kumejima Island (93%, 100%, and 1.00). Of these, G. k. toyamai was shown to have diverged first (89%, < 70%, and 0.95), whereas G. k. yamashinae was suggested to be closely related (76%, < 70%, and 0.96) to a cluster consisting of G. k. orientalis from Tokashikijima and Tonakijima Islands (77%, < 70%, and 0.96).
Table 2.
Pairwise fixation indices (F ST; above diagonal) and results of exact test (Raymond and Rousset [1995]; below diagonal) for samples of G. kuroiwae. Asterisk (*) indicates significant difference between samples (P < 0.05).
Fig. 2.
Maximum-likelihood (ML) phylogeny of Goniurosaurus kuroiwae derived from 12S and 16S rRNA, and cytochrome b sequence data. Numbers below each node are bootstrap proportions (BPs) in 1000 replicates of ML and maximum parsimony (MP) analyses, and Bayesian posterior probabilities (BPPs) of eight million generations in Bayesian inference (BI), respectively. BP and BPP values smaller than 50% and 0.70, respectively, are omitted.
Divergence times and historical demography
The divergence times were estimated to be 51.9 mya between G. kuroiwae and the Chinese—Vietnamese congeners, 14.5 mya between G. k. splendens and the common ancestor of six lineages from the Okinawa Island Group, 6.0 mya between G. k. toyamai and the other five lineages from the Okinawa Island Group, 3.4–5.0 mya among the latter five lineages, 1.4 mya between populations from lejima Island and the southern part of Okinawajima Island (Fig. 3).
Tajima's D for the Iheyajima population and Fu's Fs for the Kumejima population showed significant negative values, suggesting recent expansion of these populations (Table 1). However, there were no significant deviations in the values for other populations including Fs for Iheyajima and D for Kumejima populations. Individual BSP analyses, conducted on the two major lineages (the Amami (Tokunoshima Island) and Okinawa Island Group lineages), for which sufficient sizes of samples were available, revealed qualitatively different demographic histories (Fig. 3). The Amami lineage showed relatively stable demographic history, whereas the Okinawa lineage supposedly experienced gradual population decline over the past 0.2–1.0 myr.
Fig. 3.
Chronogram of Goniurosaurus species inferred from Bayesian analyses (upper) and Bayesian skyline plot (BSP) for two lineages of G. kuroiwae (bottom). Numbers below each branch are estimated divergence times. White boxes on the nodes indicate 95% of the node age posterior probability distribution. Relationships among outgroup taxa and estimated values within island populations are omitted. In the BSP, solid lines represent the median value for the log of the population size and the dashed lines represent the upper and lower 95% credible intervals. The x-axis measures time in substitutions per site per million years.
Morphological phylogeny
The MP analysis based on 10 morphological characters yielded 29 equally most parsimonious trees (tree length = 17, consistency index = 0.71, homoplasy index = 0.29, and retention index = 0.74). The consensus tree based on 50% majority rule was consistent with the phylogenetic tree proposed by Grismer et al. (1994). In the ingroup portion (Fig. 4A), G. k. yamashinae first diverged from the others, and then, two pairs of sister taxa, G. k. kuroiwae and G. k. orientalis, and G. splendens and G. toyamai formed the sister groups to each other. However, BP supports for clades were low (< 70%), except for that for the clade consisting of G. k. kuroiwae and G. k. orientalis (BP = 77%).
DISCUSSION
Diversification of Goniurosaurus kuroiwae
Based on their cladistic analysis using data for morphological characters (six of coloration, three of scalation, and one of body stature), Grismer et al. (1994) hypothesized the primary dichotomy of the ancestral Goniurosaurus kuroiwae into G. k. yamashinae from Kumejima of the Okinawa Island Group and the common ancestor of the remainder. They also treated G. k. splendens from Tokunoshima of the Amami Island Group as composing a small branch, sister to G. k. toyamai from Iheyajima of the Okinawa Island Group (Fig. 4A). Based on the analyses of partial sequence data for mitochondrial 12S and 16S rRNA genes, however, Ota et al. (1999) cast doubts as to the validity of such hypothetical relationships, e.g., by indicating non-monophyly of G. k. kuroiwae, although both the numbers of samples and of nucleotide sites therein used were too small to support any convincing arguments. The present results, while confirming Ota et al.'s (1999) views with substantially larger data set and methodological improvement, go further to negate relationships of G. kuroiwae hypothesized by Grismer et al. (1994). Our analyses, for example, indicate that the primary dichotomy in G. kuroiwae actually occurred between G. k. splendens from the Amami Island Group and the remainder, all from the Okinawa Island Group (Fig. 4B). We suspect that the analysis by Grismer et al. (1994) suffered from an insufficient number of characters for convincing phylogenetic resolution (10 vs. 349 variable characters in the present study) and also presumably inappropriate properties of characters employed (i.e., morphological characters that might have possibly been susceptive to environmental changes; molecular characters used here are supposedly much more neutral in response to environmental changes).
Fig. 4.
Area cladograms based on the morphological data listed in Grismer et al. (1994) (A) and on the mitochondrial DNA sequences in the present study (B). Numbers below each node are bootstrap proportions of the 1000 replications of MP (< 70% are omitted). In the DNA phylogeny, nodes without enough supports are depicted as multifurcating branches (see Fig. 2). Each letter denotes a character listed in Grismer (1988): a/A) dorsal tubercles, b/B) ventral scales, c/C) scales at base of digits, d/D) lineate middorsal pattern, e/E) banding pattern, f/F) hind limb banding, g/G) interspace mottling, h/H) juvenile coloration, i/I) eye color, and j/J) body stature. The letter case denotes the putative primitive (lower case) or derived (upper case) state of a given character. The common ancestor is presumed to have exhibited primitive conditions in these characters (a, b, c, d, e, f, g, h, i, and j).
Our results also indicate a relatively large genetic differentiation between the representatives of G. kuroiwae from the two island groups (6.0%, 6.3%, and 15.1% in 12S and 16S rRNAs, and cytochrome b, respectively). These values greatly exceed those in the corresponding sequences between other pairs of closely related terrestrial vertebrate taxa from the Amami and Okinawa Island Groups, e.g., 1.77% in 12S rRNA and 2.05% in 16S rRNA between Odorrana splendida from Amamioshima and O. ishikawae from Okinawajima (Kuramoto et al., 2011); 5.0% in 16S rRNA between Rana kobai from Amamioshima and Tokunoshima and R. ulma from Okinawajima and Kumejima (Matsui, 2011). The values are also larger than those between a pair of good species of the Chinese—Vietnamese Goniurosaurus (2.1% in 12S rRNA, 1.3% in 16S rRNA, and 3.4% in cytochrome b between G. huuliensis and G. luii). These strongly suggest that the genetic divergence between G. k. splendens from Tokunoshima of the Amami Island Group and the remaining subspecies, all from the Okinawa Island Group actually has proceeded to that at the species level (see below for further discussion).
With respect to the four subspecies from the Okinawa Island Group, the present results also clearly negate the hypothetical relationships proposed by Grismer et al. (1994) (Fig. 4), as our analyses reveal the presence of at least six distinct lineages within the Okinawan populations and because the statistical tests (SH and AU tests) rejected monophyly for both G. k. orientalis and G. k. kuroiwae. None of the foregoing reports including references to the taxonomic arrangement of Goniurosaurus from the Ryukyus, such as Nakamura and Uéno (1963), Ota (1989), and Grismer et al. (1994), have examined morphological variation in G. k. kuroiwae or G. k. orientalis in detail, at least partly due to insufficiency in sizes of available samples other than those from Okinawajima (for G. k. kuroiwae) or from Tokashikijima (for G. k. orientalis). Contrary to these, our analyses explicitly indicate that G. k. orientalis from each of the other three islands constitutes independent lineages. Moreover, population on lejima, currently assigned to G. k. orientalis, was shown to be genetically closest to the southern Okinawajima population of G. k. kuroiwae. The northern Okinawajima population of G. k. kuroiwae constitutes a lineage independent of that accommodating the lejima-southern Okinawajima assemblage. Although Ota et al. (1999) predicted such relationships, they failed to propose the whole phylogeny of the species because they analyzed no more than four populations (one G. k. orientalis, two G. k. kuroiwae and one G. k. splendens). On the other hand, our results, showing relationships among samples from almost entire range of the species, clearly show that the within-island divergence has proceeded to those among morphologically recognized subspecies. This view is circumstantially supported by the distant relationship observed between the southern and northern Okinawajima populations of other terrestrial vertebrates (e.g., Honda et al., 2012).
Grismer et al. (1994) rooted on the branch between G. k. yamashinae and the others because the former shared 10 primitive character states with the Chinese—Vietnamese congeners (Fig. 4A), of which three (absence of enlarged scales at base of digits [c], whitish dorsal pattern in juveniles [h] and yellow-brown iris [i]) were putative plesiomorphs possessed solely by G. k. yamashinae among subspecies of G. kuroiwae. However, when we consider character evolution on the phylogeny derived from sequence data including the Chinese—Vietnamese Goniurosaurus as outgroups (Fig. 4B), character states (enlarged scales at the base of digits [C], orange-pink dorsal pattern in juveniles [H] and red iris [I]) used to argue for monophyly of the non-yamashinae assemblage in Grismer et al.'s (1994) analysis were supposed to have changed in the common ancestor of the species initially and then evolved in G. k. yamashinae reversely. The lineate tendencies in dorsal pattern (D), which was assumed to be a synapomorphy of G. k. kuroiwae and G. k. orientalis, may also be attributable to the consequence of reversed character evolution from the common ancestor of G. k. kuroiwae, G. k. orientalis and G. k. yamashinae to the current G. k. yamashinae. We could depict the absence of interspace mottling (G) as parallelism in G. k. splendens and G. k. toyamai, but an alternative hypothesis, i.e., reversal evolution from the common ancestor of the whole species into the clade consisting of G. k. kuroiwae, G. k. orientalis, and G. k. yamashinae would also be equally likely as well. Likewise, the absence of dorsal banding pattern (E) might be a consequence of parallel derivation in two lineages of G. k. kuroiwae, although our interpretation should be carefully tested on the basis of much better-resolved tree in future.
Biogeography of G. kuroiwae
Several previous authors have argued that the Ryukyu Archipelago had landbridge connections to Taiwan and continental China twice since its separation from the Eurasian Continent by formations of East China Sea and Japan Sea at the beginning of the Cenozoic era (e.g., Kizaki and Oshiro, 1977, 1980; lijima and Tada, 1990). Although complete consensus has not yet been attained on the geohistorical process of the archipelago (e.g., Ujiié, 1990; Kimura, 1996, 2002), many researchers, including a number of biogeographers, have attempted to verify the validities of those two landbridges (e.g., Ota, 1998). According to Kizaki and Oshiro (1977, 1980), the first landbridge lasted from the middle to late (Langhian to Tortonian) Miocene (ca. 15-10 mya), and the second during the early (Calabrian) Pleistocene (ca. 1.5 mya). Many terrestrial animals supposedly dispersed into the Ryukyu region from the continent via these landbridges (e.g., Kizaki and Oshiro, 1977), and Tanaka (1996) considered G. kuroiwae had also migrated into the Central Ryukyus from the continent by passing through the first landbridge in the late (Tortonian) Miocene (ca. 10 mya). However, our results are not congruent with Tanaka's (1996) view, as they indicated a much longer divergence time (52 mya) between G. kuroiwae and the Chinese—Vietnamese congeners. We suspect that the common ancestor of the current Ryukyu and the Chinese—Vietnamese Goniurosaurus had been broadly distributed around the eastern Eurasian Continent including the area west to the current Ryukyus some time between the Cretaceous and the early Paleogene (i.e., the Paleocene period) (Kizaki and Oshiro, 1977, 1980; lijima and Tada, 1990) and that the divergence between the Ryukyu and the continental (= current Chinese—Vietnamese) lineages may have been involved by the primary separation of the Central Ryukyus from the continent around the early Eocene.
With regard to the historical biogeographical process of divergences in amphibian and reptile lineages between the Amami and Okinawa Island Groups, three types of events have been postulated on the basis of DNA data, i.e., 1) dispersal (as applied to Plestiodon marginatus [Honda et al., 2008]), 2) vicariance (as applied to Odorrana amamiensis and O. narina [Matsui, 1994]; Cynops ensicauda [Tominaga et al., 2010]; O. splendida and O. ishikawae [Matsui et al., 2005b; Kuramoto et al., 2011]; and Echinotriton andersoni [Honda et al., 2012]), and 3) combination of both dispersal and vicariance (as applied to Microhyla ornata sensu lato [Matsui et al., 2005a]). Based on their cladistic analysis, in which G. k. splendens was placed as sister to G. k. toyamai, Grismer et al. (1994) tentatively ascribed the endemic occurrence of the former on Tokunoshima to dispersals by the ancestor common with the latter from Iheyajima. However, the present results suggest that G. k. splendens is actually much more likely to have emerged through the primary vicariance of the common ancestor of the entire G. kuroiwae between the two island groups in the middle (Langhian) Miocene (15 mya), i.e., much earlier than the formation of straits in the early (Calabrian) Pleistocene (ca. 1.55 mya: Osozawa et al., 2011) which should have initiated division of the terrestrial fauna of the Central Ryukyus into the two island group elements.
Of the populations of the Okinawa Island Group, G. k. toyamai from Iheyajima appears to have diverged first, and this view seems concordant with the bathymetry-based idea that Iheyajima was isolated earlier than most other components of the Okinawa Island Group (Ota, 1998), although estimated divergence date of G. k. toyamai (6.0 mya) is obviously much earlier than the putative dates of separation of Iheyajima and most other islands of this group (Calabrian Pleistocene at most). The common ancestor of the other five lineages seems to have diverged into three clades that geographically correspond to populations of northern Okinawajima, southern—central Okinawajima + lejima and Kumejima + Tonakijima + Tokashikijima (5.0-4.7 mya). Within the last clade, G. k. yamashinae from Kumejima diverged from the other (3.9 mya), followed by subdivision between G. k. orientalis from Tonakijima and Tokashikijima (3.4 mya). This series of diversifications in the area west-northwest of Okinawajima is estimated to have occurred in the Pliocene, i.e., long before the middle and late (Ionian and Tarantian) Pleistocene sea level changes, although their statistical supports are not strong in MP analysis.
Interestingly, the southern part of Okinawajima is most likely to have been connected to Tokashikijima and adjacent islets of the Kerama Islands via a landbridge during the Last Glacial Maxima (LGM: ca. 0.07-0.01 mya), because the sea level is most likely to have lowered by ca. 120 m then (Fairbanks, 1989; Ohshima, 1990; Ota et al., 1993; Park et al., 1996), and the straits separating Okinawajima from those islands are much shallower than 120 m (Ota, 1996, 1998). Indeed, the terrestrial newt, E. andersoni has been shown by molecular analysis to have recently experienced gene flows between Tokashikijima and the southern part of Okinawajima (Honda et al., 2012). Nevertheless, landbridge connections of some islands during the LGM did not seemingly always offer good opportunities of gene flows to this forest-dwelling lizard populations, judging from differentiations between Tokashikijima and southern Okinawajima populations (Fig. 2).
Grismer et al. (1994) were seemingly aware that the assignment of lejima population to G. k. orientalis along with populations of Tokashikijima and adjacent islets west and northwest of Okinawajima is not concordant with the putative paleogeography of the Central Ryukyus at all, because they did not refer to even a single previous paper that has hypothesized with evidence an exclusive landbridge connection of lejima with other islands inhabited by G. k. orientalis. Instead, they interpreted the apparent similarity among lejima and other island populations identified to this subspecies as the consequences of parallel derivation (contra Shimojana [1979]). Considering the possible low-tolerance of G. kuroiwae against dehydration (Grismer et al., 1994), gene flows by direct oversea dispersals between lejima and other islands are also highly unlikely. Thus, the closest genetic affinity of the lejima population with the southern Okinawajima population demonstrated in the present study is biogeographically so puzzling, requiring further studies.
Implications for taxonomy and conservation
The subspecific classification of G. kuroiwae has been more or less confusing during the last several decades because it has almost exclusively depended on very few morphological characters, such as body colorations. Nakamura and Uéno (1963) recognized only three subspecies: G. k. kuroiwae from Okinawajima, G. k. orientalis from Tonakijima, Tokashikijima and Kumejima, and G. k. splendens from Tokunoshima. They regarded G. k. yamashinae as a junior synonym of G. k. orientalis. Grismer (1988) synonymized G. k. orientalis with G. k. kuroiwae on the one hand, but resurrected G. k. yamashinae. On the other hand, Ota (1989) referred to all insular populations west and northwest of Okinawajima including the Kumejima population as G. k. orientalis; he assigned populations from Okinawajima and two adjacent islets (Sesokojima and Kourijima) to G. k. kuroiwae, and retained the Tokunoshima population under the name of G. k. splendens. In the comprehensive revision of G. kuroiwae, Grismer et al. (1994) recognized five subspecies by resurrecting Kumejima population as G. k. yamashinae and by describing Iheyajima population as G. k. toyamai. Later, Grismer et al. (1999) proposed to elevate all these subspecies to distinct species without giving any reason to justify such treatment.
As discussed above, G. k. splendens from Tokunoshima of the Amami Island Group was obviously primarily separated from the remaining populations of the species, all from the Okinawa Island Group with substantial genetic divergence (10.8% in p-distance). The value is larger than that between a combination of other good species of the genus (2.7%: see above). Morphologically, Grismer et al. (1994) unequivocally distinguished individuals from Tokunoshima from those of the Okinawa Island Group by the canonical discriminant analysis of morphometric characters, although they took no account of the result into their taxonomic conclusion. We thus consider it appropriate to treat the Tokunoshima population as a good species, Goniurosaurus splendens (Nakamura and Uéno, 1959), rather than retaining it in the subspecific status of G. kuroiwae.
Present results also suggest non-monophyly of each of G. k. orientalis and G. k. kuroiwae. Considering the fact that the type locality of G. k. orientalis is Tonakijima and that the Tokashikijima population represents a unique phylogenetic lineage with a relatively large genetic distance from the other Okinawan populations examined so far, the Tokashikijima population should be treated as a distinct subspecies. Likewise, the southern Okinawajima population of G. k. kuroiwae is likely to deserve recognition as another undescribed subspecies, because the type locality of G. k. kuroiwae is located in the northern part of Okinawajima. However, the former population exhibited close genetic and cladistic affinities with the lejima population of G. k. orientalis. Further comprehensive studies are desired for elucidation of genetic structures the Okinawajima and lejima populations including descriptions of new taxa with additional data, those on morphological diagnostic characters in particular.
From a conservational viewpoint, the Central Ryukyus is an important area due to the numerical dominance of relict endemic and more or less threatened species, (e.g., Ota, 1998, 2000b; Environment Agency of Japan, 2012). Of these, all G. kuroiwae subspecies are assigned to one of those endangered categories by the Japan Ministry of Environment (CR, CR, EN, VU and EN for G. k. toyamai, G. k. yamashinae, G. k. orientalis, G. k. kuroiwae and G. k. splendens, respectively: Environment Agency of Japan, 2012), and by the Environment Sections of the Okinawa and Kagoshima Prefectural Governments (CR, EN, EN, VU and EN1 for G. k. toyamai, G. k. yamashinae, G. k. orientalis, G. k. kuroiwae and G. k. splendens, respectively: Ota, 2003; Tanaka, 2005). Results of the present study suggest that each island population almost invariably constitutes an independent phylogenetic lineage by itself and that individuals from Okinawajima are further divided into two unique lineages. These most likely reflect the current underestimation of taxonomic diversity in G. kuroiwae.
Also, in the present study, sample from each island or island district showed a relatively large genetic divergence with possession of unique haplotypes. It is noteworthy that samples from small islets, such as Tonakijima, Tokashikijima, lejima, and Iheyajima, had seemingly lost their genetic diversities completely or almost completely in DNA sequences examined. In general, an effective population size with low genetic diversity could be extremely smaller than that expected from the total number of individuals in the population (Hedrick, 1996; Kelly, 2001). In addition, our demographic analyses suggest that the population size of the Okinawa Island Group has decreased recently, although population expansion was weakly suggested for two relatively large islands, Iheyajima and Kumejima. Because the assignment of wildlife to threatened categories in the red list almost invariably depends on a morphologically defined taxonomie names, further effort should be made to describe those undescribed and apparently endangered lineages that were elucidated in the present analyses.
ACKNOWLEDGMENTS
We thank M. Toyama, S. Tanaka, F. Sato, T. Hamaguchi, and K. Nakata for helping and encouraging us throughout this project. Special thanks are due to H. Wada, K. Yahata, M. Nakamura, T. Denda, and members of the Graduate School of Life and Environmental Sciences, University of Tsukuba and Tropical Biosphere Research Center, University of the Ryukyus, for continuous support for our laboratory experiments. Two anonymous reviewers provided valuable comments on the manuscript.
Our research was partially supported by a Grant-in-Aid from the Japan Ministry of Education, Culture, Sports, Science and Technology (Scientific Research C No. 22510244 to M. Honda), and by a grant from the Zoshinkai Fund for Protection of Endangered Animals (to H. Ota). All subspecies of Goniurosaurus kuroiwae are protected by the laws of Kagoshima and Okinawa Prefectural Governments of Japan and Sampling was conducted under permissions from the Kagoshima Prefectural Board of Education (reference No. 17 to M. Honda, 2012) and from the Okinawa Prefectural Board of Education (reference Nos. 11 and 37 to M. Toda, 2009– 2013).
