The Ryukyu Islands have experienced a complex history of island connection and separation since the Late Miocene. Questions persist regarding how organisms have evolved through changes in island configurations over millions of years particularly in the Central Ryukyus, which is characterized by high species endemism. We conducted comparative phylogenetic and phylogeographic analyses among island populations targeting five mammalian species (Pentalagus furnessi, Diplothrix legata, and three Tokudaia species), all of which belong to genera endemic to the Central Ryukyus, employing genome-wide SNP and mitochondrial DNA variation analyses. The SNP and mitochondrial DNA analyses across these genera revealed distinct lineages on each island (Amami-Ohshima, Tokunoshima, and Okinawajima, except for P. furnessi absent from Okinawajima), with the Amami-Ohshima population showing the highest genetic diversity. Divergence times between island populations were estimated to be in the Middle Pleistocene for P. furnessi and D. legata and ranged from the Late Miocene to the Early Pleistocene for Tokudaia. These findings, along with those of previous studies on other terrestrial animals, emphasize the significant impact of pre-Pleistocene island isolation on genetic divergence in species with limited dispersal abilities, while some species exhibited inter-island migration by the Middle Pleistocene.
後期中新世以降の島嶼集団ダイナミクス:中琉球に固有の哺乳類 3 属（Pentalagus，Diplothrix，および Tokudaia）における比較系統地理学的研究．琉球列島は後期中新世以降，島間の連結と分断を繰り返してきた複雑な歴史を有している．特に多くの固有種が生息する中琉球では，過去数百万年にわたる島嶼構造の変化の下，生物がどのように進化してきたかについて，未だ多くの疑問が残されている．本研究では，中琉球に属レベルで固有の 5 種の哺乳類（アマミノクロウサギ Pentalagus furnessi，ケナガネズミ Diplothrix legata，およびトゲネズミ類 Tokudaia 属の 3 種）を対象に，ゲノムワイドな一塩基多型（SNP）およびミトコンドリア DNA 変異に基づいた比較系統学的・系統地理学的解析を行った．SNP およびミトコンドリア DNA の双方の解析結果から，上記の 3 属いずれでも，各島集団（奄美大島，徳之島，および沖縄島．ただしアマミノクロウサギは沖縄島に生息しない）はそれぞれ独立した系統であることが確認され，奄美大島の集団の遺伝的多様性が最も高いことが明らかとなった．島集団間の分岐時期は，アマミノクロウサギとケナガネズミでは中期更新世と推定されたのに対し，トゲネズミ類では後期中新世から前期更新世に遡ると推定された．以上の結果は，他の陸生動物の先行研究の結果も合わせて考えると，分散能力が限定されている種については更新世以前に生じた島嶼の隔離が遺伝的分岐に重要な役割を果たした一方で，一部の種においては中期更新世まで島間の移動が生じていたことを示している．
Fig. 1.
Location of the Central Ryukyus in Japan. The Central Ryukyus is isolated from the Northern Ryukyus by the Tokara Tectonic Strait and from the Southern Ryukyus by the Kerama Gap. The areas on the map where sea depth is >1000 m are represented in white. The straits between the three major islands in the Central Ryukyus—Amami-Ohshima, Tokunoshima, and Okinawajima—are deeper than 200 m.
The Ryukyu Islands consist of approximately 150 islands stretching along the eastern edge of the Asian continent, spanning ∼1200 km between Kyushu and Taiwan (Fig. 1) and harboring many endemic taxa. These islands form a curved chain known as the Ryukyu Arc, generated by the subduction of the Philippine Sea Plate beneath the Eurasian Plate along the Ryukyu Trench. The Ryukyu Islands are divided into three geographic regions by the Tokara Tectonic Strait and the Kerama Gap, which also serve as significant biogeographic boundaries: the Northern, Central, and Southern Ryukyus (Ota 1998; Motokawa 2000; Komaki 2021; Fig. 1). The terrestrial fauna of the Northern Ryukyus is similar to that of Kyushu Island, while that of the Southern Ryukyus resembles that of Taiwan Island. In contrast, many relict species inhabiting the Central Ryukyus are highly endemic (Ota 1998; Motokawa 2000; Okamoto 2017; Sato 2017). This biogeographic pattern is believed to reflect the geological history of the Islands.
The paleogeography of the Islands is controversial based on topographic surveys and biogeographical insight, but the temporal and spatial perspectives remain unresolved (Furukawa and Fujitani 2014). Many studies on the paleogeography and biogeography of the Ryukyu region have focused on the timing of separation and reconnection of the three Ryukyu regions to surrounding areas, such as Taiwan and the Eurasian continent. The history of the formation of the Islands has been considered in four main periods. The first period is the Late Miocene (5.333–11.63 million years ago [Mya]), when the Ryukyu region, which was part of the eastern edge of the continent, began to separate from the continent about 10–6 Mya (Iryu et al. 2006; Gungor et al. 2012; Wang et al. 2014). Some theories have suggested that the Central Ryukyus have been isolated by sea since that time (Kizaki and Oshiro 1977, 1980; Ujiie 1986; Wang et al. 2014; Su et al. 2016; Xu et al. 2016). The second period is the Early Pleistocene (2.58–0.77 Mya), during which land bridges may have formed in some areas, although the extent of these bridges is controversial (Furukawa and Fujitani 2014). The Okinawa Trough, a back-arc basin separating the Ryukyu Arc and the East China Sea shelf, is thought to have begun expanding during this period (Kizaki and Oshiro 1977; Letouzey and Kimura 1986; Ujiie 1994; Miki 1995; Sibuet et al. 1995, 1998). It has been proposed that land bridges extended from the continent and Taiwan to the Central Ryukyus (Kizaki and Oshiro 1977, 1980; Ujiie 1986; Kimura 1996). Osozawa et al. (2012) suggested that the Ryukyu Islands were part of the continent until 1.55 Mya. However, the high endemism of the biota in the Central Ryukyus indicates pre-Pleistocene isolation rather than Pleistocene colonization from the continent (Hikida and Ota 1997; Ota 1998; Matsui et al. 2005; Su et al. 2016; Xu et al. 2016; Okamoto 2017). The third period is the Middle Pleistocene (0.77–0.129 Mya). During this time, the formation of sedimentary layers indicating the development of coral reefs south of the Tokara Tectonic Strait, suggests island fragmentation and a decrease in the land area due to shallow seas (Kizaki and Oshiro 1977; Ujiie 1986). Fossils of terrestrial mammals suggest that the Central Ryukyus have been isolated from other areas since at least the Middle Pleistocene period to the present (Kawamura et al. 2016). However, the topographic development of the Islands due to Quaternary crustal deformation during this period is complex (Otsubo and Ogata 2024), and little is known about the exact timing of the connections and separations among islands (Furukawa and Fujitani 2014). The final period is the Late Pleistocene (0.129–0.0117 Mya), during which the island configuration was explained by fluctuations of sea level during the last glacial period (Furukawa and Fujitani 2014). The major islands in the Central and Southern Ryukyus separated by straits deeper than 200 m are believed not to have reconnected (see Fig. 1). However, it is also suggested that a continuous land area was formed in the Central Ryukyus (Oshiro and Nohara 2000; Oshiro 2002).
Although several hypotheses have been proposed, the high endemism in the Central Ryukyus suggests that this region may have been isolated from the surrounding regions since the Late Miocene. However, there has been little discussion on the history of how the islands within the Central Ryukyus have connected and separated to form their current configuration after initially separating from the continent. Thus, there is still considerable uncertainty about how the endemic species of the Central Ryukyus have persisted and evolved over millions of years through complex geomorphological and environmental changes. Comparative phylogeography of many species inhabiting such old island systems provides important insight into the geological dynamics of the region while elucidating the mechanisms of organismal evolution on a long time scale (Parent et al. 2008; Poulakakis et al. 2013; Hu et al. 2016; McCulloch and Waters 2019).
Phylogeographic studies of non-volant vertebrates in the Central Ryukyus have focused on amphibians and reptiles (Ota 1998; Okamoto 2017). These studies have revealed that genetic divergence among island populations within the Central Ryukyus is highly diverse from the Miocene to the Middle Pleistocene and has sparked discussions on the origin in geological history (Honda et al. 2014; Kurita and Hikida 2014; Tominaga et al. 2015; Kaito and Toda 2016; Yang et al. 2018). Many of the estimated divergence times predate the Pleistocene and do not align with the proposed formation periods for land bridges during the Early and Late Pleistocene. Low dispersal capabilities and habitat-specific preferences likely strongly affected early-stage divergence and speciation following colonization, potentially restricting gene flow between adjacent island populations during the recent land bridge (Parent et al. 2008; Papadopoulou and Knowles 2015; Cros et al. 2020). In contrast, genetic population structures of terrestrial mammals with high dispersal capabilities may have been affected by recent migration events during the land bridge period and gene flow between island populations (Hirata et al. 2013, 2017; Segawa et al. 2022; Kinoshita et al. 2024). Therefore, phylogeographical inferences of terrestrial mammals in the Central Ryukyus would offer new insight into the formation of highly endemic faunal assemblages and the geographical history of the island system.
The Ryukyu Islands harbor nine endemic non-volant terrestrial mammal species, seven of which are confined to the Central Ryukyus (Amami rabbit Pentalagus furnessi, Ryukyu long-furred rat Diplothrix legata, three species of Tokudaia spiny rats, Orii's shrew Crocidura orii, and Watase's shrew C. watasei). Particularly, the genera Pentalagus, Diplothrix, and Tokudaia harbor relict species endemic to two or three major islands (Amami-Ohshima, Tokunoshima, and Okinawajima) in the Central Ryukyus (Fig. 2a–c). Divergence times inferred from mainland relatives of these genera are estimated to predate the Pleistocene (Sato 2017). However, estimates of the phylogenetic relationships and divergence times among island populations in the Central Ryukyus have only investigated Tokudaia using a small number of mitochondrial and/or nuclear gene markers (Suzuki et al. 1999; Murata et al. 2010; Sato 2017). Notably, Tokudaia osimensis in Amami-Ohshima and T. tokunosimensis in Tokunoshima lack a Y chromosome, a rare trait among mammals worldwide, suggesting that a unique sex determination mechanism was established after divergence from T. muenninki in Okinawajima possessing the Y chromosome (Honda et al. 1977, 1978; Murata et al. 2010). Additionally, karyotypic differences have been reported between T. osimensis and T. tokunoshimensis (Nakamura et al. 2007), likely resulting from prolonged isolation and independent evolution of the Tokudaia island populations. On the other hand, lineage relationships and the timing of divergence between island populations for P. furnessi and D. legata remain unexplored based on genetic information.
In this study, we conducted comparative phylogeographic analyses focusing on the inter-island phylogeny of three endemic mammalian genera, such as Pentalagus, Diplothrix, and Tokudaia, to deepen our understanding of the mechanisms contributing to the high endemic biodiversity in the Central Ryukyus. We investigated their genetic diversity and divergence times to address the following questions: Does lineage divergence among island populations date back before the Pleistocene, as reported for various reptile and amphibian species? Are there genetic traces of inter-island migration or secondary contact coinciding with the formation of the Pleistocene land bridges? Did the timing of genetic divergence between island populations coincide among the three genera, or exhibit genus-specific patterns?
Materials and methods
Samples
All species of the three genera (Pentalagus, Diplothrix, and Tokudaia) were categorized as endangered or critically endangered on the International Union for Conservation of Nature (IUCN) Red List (Correia 2016; Ishii 2016a, 2016b, 2016c; Yamada and Smith 2016), and designated as natural monuments (D. legata and three species of Tokudaia) or a special natural monument (P. furnessi) of Japan. All samples used in this study ( Supplementary Table S1 (03kinoshita_2024-0045_suppl.pdf)) were collected as road-killed animals to the Amami Wildlife Center or captured for skin tissue or hair for other research projects with permission from the Ministry of the Environment and Agency for Cultural Affairs in Japan. No animals were sacrificed for this study. In total, 37 samples of P. furnessi (Amami-Ohshima, n = 20; Tokunoshima, n = 17), 31 samples of D. legata (Amami-Ohshima, n = 14; Tokunoshima, n = 3; Okinawajima, n = 14), and 44 samples of Tokudaia (Amami-Ohshima, n = 15; Tokunoshima, n = 9; Okinawajima, n = 20) were used for genomic DNA extraction in this study either by the standard phenol-chloroform extraction method (Sambrook et al. 1989) or the QIAamp DNA Mini Kit (Qiagen, Valencia, CA, USA).
MIG-seq experiment and data analysis
To identify genome-wide polymorphisms in the genomic DNA, we utilized the multiplexed inter-simple sequence repeat genotyping by sequencing (MIG-seq) method (Suyama and Matsuki 2015), which targets loci between two simple sequence repeats for amplification via polymerase chain reaction (PCR). The entire process of creating the MIG-seq library adhered to the protocol outlined by Suyama et al. (2022). During the first PCR round, we selected a primer set featuring repeated motifs and anchor sequences: (ACT) 4TG, (CTA) 4TG, (TTG) 4AC, (GTT) 4CC, (GTT) 4TC, (GTG) 4AC, (GT) 6TC, and (TG) 6AC. The library produced after the second round of PCR contained indices for sample identification and Illumina adapter sequences and was sequenced using the MiSeq platform (Illumina, San Diego, CA, USA) with a MiSeq Reagent Kit v.3 (150 cycles) (Illumina). The raw reads from each individual were processed in two steps using Fastp ver. 0.23.2 (Chen et al. 2018). The front 5 bp and tail 10 bp of the reads were trimmed during the first step. During the second step, the reads were filtered by setting the Phred quality to ≥ Q30, and the front 5 bp of the reads were trimmed further. All reads were restricted to be at least 50 bp in length.
After raw read processing, the data were divided into three sets corresponding to each genus, and the SNPs were detected. SNP calling of P. furnessi and D. legata was conducted de novo using Stacks 2.3 (Catchen et al. 2013; Rochette et al. 2019). The analytical parameters M (number of mismatches allowed between stacks within individuals) and n (number of mismatches allowed between stacks between individuals) were tested with values ranging from one to six to identify the settings that yielded the highest number of SNPs. Consequently, for P. furnessi, M = 2 and n = 2 were adopted, and for D. legata, M = 3 and n = 6 were adopted. The whole-genome sequence (GCA_027943515.1) of T. osimensis was obtained from the database for the three Tokudaia species, and only the autosomes were used as references for mapping with the Burrows-Wheeler Aligner (BWA)MEM algorithm in BWA program v.0.7.12 (Li and Durbin 2009), followed by detection of the SNPs with ref_map.pl module in Stacks v2.3. The SNPs in the sample set of each genus were filtered to ensure at least 50% of the individuals on each island and 80% of the individuals across all islands, with one SNP per locus, minor alleles detected in at least two individuals, and heterozygosity below 0.70.
We estimated nucleotide diversity (π) for each population using STACKS populations v.1.44 (Catchen et al. 2013) and the proportion of heterozygous genotypes for each individual (number of heterozygous genotypes/total number of genotypes for that individual) using Tassel v.5.2.79 (Bradbury et al. 2007). Then, the average proportion of heterozygous genotypes (PHG), which is the number of heterozygous genotypes/number of sites in an individual), was calculated for each population. A cluster analysis was performed using STRUCTURE ver. 2.3.4 (Pritchard et al. 2000) with K = 2–5 for each genus-specific dataset. The Tokudaia samples were divided into two clusters, regardless of the K value, corresponding to Okinawajima and Amami-Ohshima + Tokunoshima. Therefore, the Amami-Ohshima and Tokunoshima samples were analyzed again excluding the Okinawajima samples with K = 2–4. Independent runs of SNP calling and filtering were performed for the Tokudaia Amami-Ohshima + Tokunoshima dataset as described above. The number of iterations for each K was set to ten runs, and the most frequent clustering pattern among the 10 runs was adopted based on the summary by CLUMPAK (Kopelman et al. 2015). We calculated pairwise genetic distances between individuals in each genus dataset under the p-distance model using ngsDist (Vieira et al. 2016) with the option --avg_nuc_dist and constructed a NeighborNet tree using SplitsTree v.4.14.5 (Huson and Bryant 2006).
To assess possible genetic divergence scenarios among insular populations in each genus, we performed approximate Bayesian computation (ABC) with a random forest approach (Breiman 2001) using DIYABC-RF v.1.2.1 (Collin et al. 2021). As there were only two P. furnessi populations (Amami-Ohshima and Tokunoshima), we considered a simple single-divergence scenario and estimated changes in population size and the timing of divergence. We tested six scenarios for D. legata and Tokudaia by assuming simple divergence among the three islands without admixture (Scenarios 1, 2, and 3) or with a single admixture event for a certain population (Scenarios 4, 5, and 6) ( Supplementary Fig. S1 (03kinoshita_2024-0045_suppl.pdf)). We adopted previous settings of the parameters for effective population sizes and the numbers of generations listed in Supplementary Table S2 (03kinoshita_2024-0045_suppl.pdf). All previous values were drawn from uniform distributions and an additional condition was applied to specify t1 < t2. We assumed changes in the effective population size of the ancestral population before the split at t2 (i.e., Na) in all scenarios and the islands at t1 (i.e., Nb) in Scenarios 1, 2, and 3 to improve simulation efficiency. Using all summary statistics obtained from 100 000 simulations per scenario as well as the optional axes of linear discriminant analysis, a random forest estimate was conducted to determine the best-fit model and to estimate parameters under the chosen model. The number of trees in the random forest was set to 1000. The best-fit model was chosen based on the number of classification votes and then the posterior probability for that scenario as well as global and local error rates were estimated to assess the choice of scenario and the quality of prediction. Although there was no information about the P. furnessi generation time, we scaled the time using a generation time of two years in other Leporidae (Marboutin and Peroux 1995). On the other hand, we assumed a generation time of one year for D. legata and Tokudaia as their reproductive seasons were limited mainly to autumn and winter (Okano et al. 2015a; Kotaka et al. 2021).
MtDNA sequencing and analysis
We sequenced 1140 bp of the mitochondrial cytochrome b (Cytb) gene using the primer sets listed in Supplementary Table S3 (03kinoshita_2024-0045_suppl.pdf). PCR was carried out in a total volume of 10 µL containing 1 µL template genomic DNA, 0.05 µM of each set of primers, and 5 µL AmpliTaq Gold 360 Master Mix kit (Applied Biosystems, Foster City, CA, USA). The PCR conditions were 10 min of denaturation at 95°C, 35 cycles consisting of three stages, 30 s of denaturation (95°C), 30 s of annealing (50–55°C), 60 s of extension (72°C), and 7 min of extension at 72°C. The PCR products were sequenced using the BigDye Terminator Cycle Sequencing Kit v3.1 (ABI), followed by automated sequencing on an ABI3130 Genetic Analyzer (ABI). The Cytb sequences were aligned manually using ProSeq version 3.5 software (Filatov 2009). The P. furnessi and D. legata datasets were merged with database sequences (the three P. furnessi samples analyzed in Yamada et al. 2002, and the one D. legata sample analyzed in Suzuki et al. 2000). A median-joining network was reconstructed for each genus dataset using Network ver. 10.2 to infer the relationships between the observed haplotypes (Bandelt et al. 1999).
Phylogenetic trees for the mtDNA sequences were reconstructed using the Bayesian inference method and BEAST2 version 2.7.3 (Drummond and Rambaut 2007). According to previous studies that have inferred phylogenetically closely related taxa for each targeted genus, the database sequences of Oryctolagus cuniculus (Accession No. OCU07566) and Bunolagus monticularis (AY292718) were used for P. furnessi, Rattus norvegicus (AY172581) and R. rattus (NC012374) were used for D. legata, and Apodemus speciosus (AB164496), A. agrarius (AB096817), and Mus musculus (AB125774 and AB205312) were used as outgroup species for Tokudaia. Divergence times among the island populations were estimated by setting a calibration point obtained from TimeTree5 (Kumar et al. 2022) as the divergence time of two outgroup genera or species for each targeted genus ( Supplementary Table S4 (03kinoshita_2024-0045_suppl.pdf)). The best substitution model was selected based on the Bayesian Information Criterion (BIC) values given by the results of “Find best DNA/Protein Models” implemented in MEGA v. 5.2.2 (Tamura et al. 2011). The best-fit substitution models for the P. furnessi, D. legata, and Tokudaia datasets were TN93 + G, TN93 + G, HKY + I, respectively. Coalescent constant population was assumed as the tree prior. Bayesian searches were conducted in BEAST2 using the Markov chain Monte Carlo method for 50 × 106 generations with trees sampled every 5000 generations. Each log file was checked to confirm convergence to the stationary posterior distribution and sufficient effective sample size of each parameter using the Tracer ver. 1.7.2 program (Rambaut et al. 2018). The trees were summarized using TreeAnnotator in the BEAST2 package with the settings “Maximum clade credibility tree” and “Mean heights” and visualized with FigTree v. 1.4.4 software (Rambaut 2016).
Results
Phylogeography based on genome-wide SNPs
The number of MIG-seq raw reads per individual ranged from 47 143 to 135 956 reads in P. furnessi, from 109 306 to 295 462 reads in D. legata, and from 44 384 to 402 003 reads in Tokudaia. After processing the data, we retained 526 SNPs for 28 P. furnessi individuals, 399 SNPs for 24 D. legata individuals, and 1971 SNPs for 37 Tokudaia individuals (Table 1). Independent SNP calling and filtering runs for the Tokudaia Amami-Ohshima and Tokunoshima datasets yielded 3310 SNPs which were only used for clustering analysis. The genetic diversity estimated from the SNP data was highest in the Amami-Ohshima population in all genera (the Amami-Ohshima population of P. furnessi: π = 0.295 and PHG = 0.225, D. legata: π = 0.267 and PHG = 0.201, Tokudaia: π = 0.45 and PHG = 0.60) (Table 1). The genetic diversity between Tokunoshima and Okinawajima in the case of D. legata was not substantially different (π = 0.136 and PHG = 0.101 in Tokunoshima, π = 0.103 and PHG = 0.092 in Okinawajima), but genetic diversity of Tokunoshima in Tokudaia was notably lower (π = 0.018 and PHG = 0.018 in Tokunoshima, π = 0.041 and PHG = 0.046 in Okinawajima).
Table 1.
Summary of diversity indexes based on genome-wide SNPs data from MIG-seq analysis
The STRUCTURE clustering analysis showed different genetic composition patterns among the island populations in each genus (Fig. 2d–f). Individuals from the Amami-Ohshima and Tokunoshima populations of P. furnessi were assigned to different clusters across all K values. Subpopulation structures were inferred between north and south within each island with an increasing K value from 3 to 5, but no clusters were shared between the islands. The D. legata Tokunoshima population appeared to be admixed between the clusters detected in Amami-Ohshima and Tokunoshima when K = 2, but it was assigned to a different cluster when K = 3. Although unclear due to the small sample size, individuals from the northern and southern parts of each island were inferred to have different cluster proportions when increasing K = 4 or 5. When analyzing all Tokudaia species from the three islands with K = 2, they were separated into Okinawajima and Amami-Ohshima + Tokunoshima. The two species from Amami-Ohshima and Tokunoshima were divided into two distinct clusters regardless of the K value.
Similar divergence patterns were uncovered in the NeighorNet analysis (Fig. 2g–i). Individuals from the Amami-Ohshima and Tokunoshima populations of P. furnessi noticeably diverged, with individuals from the northern and southern parts of each island in distinct clusters. On the other hand, while inter-island divergence was clear for D. legata, geographical divergence within Amami-Ohshima or Tokunoshima was unclear, as in the clustering analysis. The Tokudaia species from the three islands diverged separately but the branch length between Okinawajima and Amami-Ohshima + Tokunoshima was notably longer than the branch length between Amami-Ohshima and Tokunoshima. Branches separating island clusters were considerably longer than branches within each island cluster, compared to the NeighborNet tree of P. furnessi or D. legata.
We conducted an ABC analysis for each genus to infer population phylogeny and estimate divergence times among the insular populations (Fig. 2j–l). Based on a simple divergence scenario, the expected divergence times between the Amami-Ohshima and Tokunoshima populations of P. furnessi were estimated to be 218 062 (95% confidence interval [CI]: 39 599–428 990) generations ago. If we assume a generation time of two years, the divergence time was c.a. 0.44 Mya (95% CI: 0.86–0.08 Mya) during the Middle Pleistocene. By comparing six demographic scenarios for D. legata, Scenario 4, which assumed that the Tokunoshima population was established as an admixture between the Amami-Ohshima and Okinawajima lineages, received the highest percentage of classification votes at 65.9% resulting in high posterior probability (PP = 0.630 with global and local errors of 0.243 and 0.130, respectively). The first divergence between the ancestral lineages of the Amami-Ohshima and Okinawajima populations was estimated to have occurred 0.70 Mya (95% CI: 1.06–0.39 Mya) if we assume a generation time of one year. Then, the admixture event in Tokunoshima occurred at 0.40 Mya (95% CI: 0.64–0.20 Mya). These events date back to the Middle Pleistocene. The best model for Tokudaia was Scenario 1, in which the Okinawajima population branched off and then the Amami-Ohshima and Tokunoshima populations were separated. The t1 divergence time was estimated to be 1.06 Mya (95% CI: 2.22–0.35 Mya) during the Early Pleistocene, and that of t2 was estimated to be 5.08 Mya (95% CI: 7.72–1.86 Mya), ranging between the Late Miocene to the Early Pleistocene, if we assume a generation time to be one year.
Phylogeography based on mitochondrial DNA sequences
In total, eight Cytb haplotypes were detected from P. furnessi and D. legata, while 14 haplotypes were detected from Tokudaia. In all three genera, the Amami-Ohshima population exhibited the largest number of haplotypes, ranging from 4 to 11, whereas the other insular populations showed only 1 to 3. Haplotype network analysis revealed a separation of the haplotypes into island-specific clusters across the three genera (Fig. 3a–c). Based on the Cytb sequences, we constructed Bayesian phylogenetic trees for each genus and estimated the divergence times among the insular populations. The divergence time between the Amami-Ohshima and Tokunoshima populations of P. furnessi was estimated to be about 0.312 Mya (95% height posterior density [HPD]: 0.856–0.049 Mya) by setting the calibration point for the divergence between O. cuniculus and B. monticularis to be a median time of 9.9 Mya (CI: 15.7–7.2 Mya) ( Supplementary Table S4 (03kinoshita_2024-0045_suppl.pdf)). Pentalagus furnessi was closer to B. monticularis than O. cuniculus in the tree, and they diverged around 7.00 Mya (95% HPD: 10.59–3.57 Mya). Divergence among the insular populations of D. legata first occurred between Amami-Ohshima and Okinawajima + Tokunoshima at about 0.50 Mya (95% HPD: 1.557–0.074 Mya), followed by divergence between Okinawajima and Tokunoshima of about 0.25 Mya (95% HPD: 0.67–0.012 Mya). The calibration point of the tree was set at the divergence between R. norvegicus and R. rattus to be a median time of 3.51 (CI: 5.9–2.3) Mya ( Supplementary Table S4 (03kinoshita_2024-0045_suppl.pdf)), and D. legata was closer to R. rattus with an estimated divergence time of about 2.73 Mya (95% HPD: 3.96–1.39 Mya). Compared to the divergence times of P. furnessi and D. legata during the Middle Pleistocene, those among the three Tokudaia species were considerably older. The Okinawajima lineage was first separated from the others at about 5.70 Mya (95% HPD: 7.04–4.45 Mya). Then, the Amami-Ohshima and Tokunoshima lineages diverged about 2.83 Mya (95% HPD: 3.61–2.07 Mya) around the transition from the Pliocene to the Early Pleistocene, which was slightly older than estimated based on the DYIABC analysis using SNP data. The calibration point of the tree was set at the divergence between Apodemus and Mus to be a median time of 10.7 (95% CI: 11.5–9.7) Mya ( Supplementary Table S4 (03kinoshita_2024-0045_suppl.pdf)), and Tokudaia was slightly closer to the two Apodemus species than M. musculus with an estimated divergence time of about 10.26 Mya (95% HPD: 11.54–8.52 Mya). As mentioned above, the estimated genetic divergences between island populations of P. furnessi or D. legata were significantly more recent than the calibration points that we set for the outgroups. This large gap in evolutionary timescale between calibration points and divergence events potentially negatively impacts the accuracy of time estimates (Ho et al. 2005). However, we emphasize that the age estimates based on mtDNA overlapped with those based on SNPs in the nuclear genome.
Fig. 2.
Genetic relationships based on the MIG-seq SNP datasets. The distribution ranges from the IUCN Red List version 2024-1 of P. furnessi (a), D. legata (b), and Tokudaia (c). N and S on the maps indicate the northern and southern distributions of each island population. Admixture panels produced by STRUCTURE for K values from 2 to 5 are shown for P. furnessi (d) and D. legata (e). Admixture panels for the three Toludaia island populations when K = 2 and for the Amami-Ohshima and Tokunoshima dataset when K = 2–4 (asterisks) are shown (f). Phylogenetic relationships were inferred using NeighborNet analysis for P. furnessi (g), D. legata (h), and Tokudaia (i). Schematic of the demographic scenario in the DIYABC-RF analysis is shown for P. furnessi (j), and the best models selected are shown D. legata (k) and Tokudaia (l). Generations are shown on the y-axis (t0, t1, and t2). The effective population size was set for some ancestral populations before the taxa split (Na and Nb). The admixture proportion from the parental populations is shown with ra and 1-ra in (k). Competing scenarios analyzed for D. legata and Tokudaia are shown in Supplementary Fig. S1 (03kinoshita_2024-0045_suppl.pdf). Before setting the time of divergence and admixture events (t), the admixture proportion during admixture events between two parental populations (ra), and the effective population size for each population are summarized in Supplementary Table S2 (03kinoshita_2024-0045_suppl.pdf), and the estimated values are given in Table 2.
Fig. 3.
Genetic relationships based on the mitochondrial Cytb sequences. Haplotype networks for P. furnessi (a), D. legata (b), and Tokudaia (c) according to the median-joining method. Node sizes are proportional to the haplotype frequencies in each network. Slashes on the branches between nodes indicate mutations. The Bayesian phylogenetic trees for P. furnessi (d), D. legata (e), and Tokudaia (f). Posterior probabilities in the Bayesian inference analysis are shown only for major nodes of divergence from outgroups and between island lineages. Shaded bars on major nodes indicate 95% confidence intervals of the divergence time estimates.
Table 2.
Estimates of demographic parameters by DIYABC-RF
Discussion
Isolation in the Central Ryukyus
We conducted a comparative phylogeographic analysis of three endemic mammalian genera (Pentalagus, Diplothrix, and Tokudaia) in the Central Ryukyus. Our results indicate that these three genera have been isolated in the Central Ryukyus since the initial formation of the Ryukyu Arc before the Pleistocene. However, the genetic relationships and divergence times among the island populations within the Central Ryukyus vary by genus, indicating that each genus has undergone a distinct evolutionary history under the complex geographic dynamics of the region.
The divergence times for P. furnessi and Tokudaia from continental relatives were estimated to be in the Late Miocene, and Pliocene (5.33–2.58 Mya) for D. legata. This agrees with the results summarized by Sato (2017) regarding the divergence times of endemic terrestrial mammals in the Japanese archipelago. Since the early formation of the Ryukyu Arc, which corresponds with the separation of land masses from the eastern margin of the continent (Iryu et al. 2006; Gungor et al. 2012; Wang et al. 2014), these genera have been isolated from the continent and maintained to the present. Considering that the divergence time of D. legata from its continental relatives was somewhat younger than the other two genera, the ancestral D. legata species may have dispersed to this region later than the other two genera. Among the amphibians, reptiles, and invertebrates that are endemic to the Central Ryukyus, divergence times from their phylogenetically closest species in East Asia vary widely, ranging from the pre-Miocene to the Pliocene epochs (Honda et al. 2014; Kurita and Hikida 2014; Tominaga et al. 2015; Kaito and Toda 2016; Yang et al. 2018). Nevertheless, these molecular phylogenetic estimates of divergence times strongly support the ancient origins of the terrestrial vertebrates in the Central Ryukyus, which predate the expansion of the Okinawa Trough since the Early Pleistocene (Kizaki and Oshiro 1977; Letouzey and Kimura 1986; Ujiie 1994; Sibuet et al. 1995, 1998) and the hypothesized land bridge with the continent through the Southern Ryukyus or Taiwan during the Early Pleistocene (Kizaki and Oshiro 1977; Ujiie 1986; Kimura 1996).
Phylogenetic relationships among the island populations
We obtained consistent estimates for the divergence times of the island populations in the genome-wide SNP and Cytb sequence analyses. It was estimated that the ancestral Amami-Ohshima + Tokunoshima population of Tokudaia diverged from the Okinawajima population about the end of the Miocene, with the Amami-Ohshima and Tokunoshima populations subsequently diverged between the end of the Pliocene and the Early Pleistocene. These results are consistent with previous estimates based on several DNA markers (Suzuki et al. 1999; Sato and Suzuki 2004; Murata et al. 2010; Sato 2017), and our study supports these findings with genome-wide SNP data, confirming the species relationships and divergence depths. The ancient divergence between the Okinawajima and Amami-Ohshima + Tokunoshima populations, dating back to the Pliocene to Late Miocene, has also been estimated for other amphibians and reptiles (Tominaga et al. 2010; Honda et al. 2012, 2014; Kurita and Hikida 2014; Tominaga et al. 2015; Kaito and Toda 2016; Shibata et al. 2016; Kaito et al. 2017; Makino et al. 2020). Furthermore, the divergence time between the Amami-Ohshima and Tokunoshima populations is estimated to be more recent for many of these taxa.
The divergence time between the P. furnessi Amami-Ohshima and Tokunoshima populations was estimated to be in the Middle Pleistocene, aligning with estimates for some amphibians and reptiles (Matsui et al. 2005; Honda et al. 2012; Kaito et al. 2017; Yang et al. 2018). Although P. furnessi is not currently found on Okinawajima, fossils from early Middle Pleistocene sediments have been discovered on Okinawajima, indicating that this species existed in the Central Ryukyus until at least that time (Kawamura et al. 2016). There may not have been any migration from the Amami-Ohshima + Tokunoshima to Okinawajima after the period. In this study, we estimated the divergence times among island populations of D. legata for the first time. Our results indicate that the populations of Okinawajima, Amami-Ohshima, and Tokunoshima diverged during the Middle Pleistocene. Additionally, mitochondrial DNA shows that the populations of Okinawajima and Tokunoshima are phylogenetically closer to each other than to Amami-Ohshima, and the nuclear genome data indicate that the Tokunoshima population formed through an admixture between the Okinawajima and Amami-Ohshima lineages. This pattern has not been observed in other endemic terrestrial vertebrates of the Central Ryukyus.
The Middle Pleistocene was a rapid subsidence period following the uplift of the Ryukyus Arc, which changed the island configuration (Otsubo and Ogata 2024). After the Early Pleistocene, the maximum elevation of Okinawajima decreased from more than 1500 m to the current 498 m (Kuroda et al. 2002), and pollen fossil analysis has revealed changes in vegetation from the Early to Middle Pleistocene (Fujiki and Ozawa 2008). There is no information on changes in vegetation for Amami-Ohshima and Tokunoshima but significant environmental changes that occurred in the Central Ryukyus during the Middle Pleistocene may have promoted inter-island migrations, genetic differentiation, and the extinction of some populations. Few species other than D. legata are believed to have migrated between Okinawajima and Amami-Ohshima + Tokunoshima during the Middle Pleistocene. One example is the Okinawa tree lizard (Diploderma polygonatum), which migrated among the three islands during the Middle Pleistocene, possibly by overseas dispersal on “floating logs and mats” carried by the Kuroshio Current, which may have changed its course during this period (Yang et al. 2018). Diplothrix legata, which also utilizes arboreal habitats, may have had similar opportunities. Whether long-distance over-water rafting could have occurred during the past range expansions of terrestrial mammals remains a topic of considerable debate (Ali and Vences 2019). Additionally, two Ryukyu-endemic cicada species (Nagata et al. 2021) and the Yakushima violet Viola iwagawae (Nakamura et al. 2015) also diverged between the Okinawajima and Amami-Ohshima + Tokunoshima populations in the Middle Pleistocene, suggesting that species with higher dispersal abilities were capable of migration.
The results of our study combined with previous phylogeographic findings of other organisms suggest that Amami-Ohshima and Tokunoshima were connected during the Middle Pleistocene, allowing for the migration of terrestrial animals. However, as hypothesized for some amphibians and reptiles, reproductive isolation or competitive exclusion may have prevented secondary admixture for Tokudaia (Kaito and Toda 2016). On the other hand, the Amami island group (Amami-Ohshima + Tokunoshima) and the Okinawa island group were probably separated before the Early Pleistocene. Nevertheless, overseas dispersal during the Middle Pleistocene may have been triggered in some highly dispersive species due to the narrowing of the straits or changes in ocean currents.
Although our study determined the molecular phylogenetic relationships among extant populations, these three genera are historically more widely distributed within the Central Ryukyus based on fossil records (Kawamura et al. 2016; Nakamura 2018). Moreover, it has been suggested that a currently submerged ancient island played a crucial role in the dispersal of endemic biota between the Central and Southern Ryukyus (Watanabe et al. 2023). Future advances in genomics analysis of Ryukyu biota, coupled with further integration of paleobiology and paleogeography, are expected to uncover the mechanisms underlying the biodiversity in the Ryukyu Islands.
Genetic diversity and population structure within each island
Across the three genera examined in this study, the genetic diversity of the Amami-Ohshima population was the highest, which may reflect the breadth of its current distribution range and population size. The genetic diversity of the D. legata Tokunoshima and Okinawajima populations was similar based on the SNP data but the genetic diversity of Tokudaia was lowest in Tokunoshima. These findings are counterintuitive considering the significant reduction in the distribution range and population size of Tokudaia species in Okinawajima, which raises concerns about extinction (Yamada et al. 2010). Only one haplotype was found from Okinawajima in the mitochondrial Cytb of the Tokudaia species, whereas three haplotypes were found from Tokunoshima. In addition, a previous study that estimated the genetic diversity of D. legata on Okinawajima based on mtDNA markers suggested a population bottleneck (Okano et al. 2015b). The impacts of recent population declines on genetic diversity should be investigated using larger genomic and sample datasets.
Despite the limited number of samples, signs of geographic population structure were found in P. furnessi and D. legata within the Amami-Ohshima and Tokunoshima populations. Ohnishi et al. (2017) and Ando et al. (2018) used microsatellite markers to investigate the genetic diversity of P. furnessi within Amami-Ohshima and Tokunoshima, respectively, and found that isolation by geographic distance forms the basis of genetic population structure within each island. Ohnishi et al. (2017) suggested that some genetic clusters within Amami-Ohshima may have occurred due to past genetic exchanges with the Tokunoshima population. However, in this study, genetic clusters within each island were not shared between the islands, forming genetic population structures unique to each island. Ando et al. (2018) estimated that the divergence time between subpopulations within Tokunoshima was ∼4320 years ago. This study revealed the possibility of similar genetic population structures in D. legata for the first time. While no signs of the formation of a genetic population structure within each island were observed in Tokudaia this time, the analytical method in this study may have masked the genetic diversity within each island due to the deep genetic divergence among the three Tokudaia species. Therefore, future studies should use more extensive D. legata and Tokudaia sample sets to investigate the detailed genetic diversity and geographic patterns within each island.
Comprehending the population structure within islands and the mechanisms maintaining populations among islands is essential for the sustainable management of endangered island species (Frankham 2003). In this study, it is estimated that the island populations of each genus have been isolated and persisted in their current state since at least the Middle Pleistocene. Just as Tokudaia has evolved unique sex-determination mechanisms on each island, P. furnessi and D. legata may have undergone distinct evolutionary processes on a per-island basis. Each island population of the three genera should be managed as an independent conservation unit. Small island populations that have been isolated for a long time may have low genetic diversity and vulnerable genomes (Hamabata et al. 2019). In addition, the genetic population structure within the islands may have been formed not only by long-term population history but also by recent anthropogenic activity. Thus, intra-island demographic history and population substructure should be verified through analysis based on more extensive sample sets, and higher-precision genome data are crucial from a conservation genetics perspective.
Conclusions
Our study revealed that the divergence among island populations of P. furnessi and D. legata dated back to the Middle Pleistocene, while the divergence of Tokudaia extended from the Late Miocene to the Early Pleistocene. The Central Ryukyus may have played a role in promoting genetic differentiation among island populations through repeated connections and isolation since its separation from continents during the Late Miocene, thus becoming an effective model for assessing long-term lineage diversification and speciation in island biogeography. This island system would provide valuable comparative data to other well-studied old island systems where highly endemic biological communities have developed through millions of years of geological dynamics (Ali 2018; Pellissier et al. 2018; Heads and Grehan 2021). Furthermore, while sharing a long geological history of the Central Ryukyus, our findings on the differences in the phylogeographic histories among P. furnessi, D. legata, and Tokudaia would offer important insights for developing effective conservation strategies for these endangered species.
Acknowledgments:
We express our appreciation to Shinichiro Fukumoto in Rakuno Gakuen University and Kimiyuki Tsuchiya for their support and providing valuable samples to start up this study. We are also grateful to the Amami Wildlife Conservation Center, Tokunoshima Nature Conservation Office, and Ministry of the Environment for granting permissions for fieldwork and sample collection, as well as for supplying valuable samples and information. This work was supported by a Grant-in-aid for Fellows of the Japan Society for the Promotion of Science (JSPS) to 23K05944 and 20K15862.
© The Mammal Society of Japan
Supplementary data
Supplementary data are available at Mammal Study online.
Supplementary Table S1 (03kinoshita_2024-0045_suppl.pdf). Sampling and sequencing information.
Supplementary Table S2 (03kinoshita_2024-0045_suppl.pdf). Prior distributions in DIYABC-RF analysis.
Supplementary Table S3 (03kinoshita_2024-0045_suppl.pdf). List of primers for PCR of the Cytb gene.
Supplementary Table S4 (03kinoshita_2024-0045_suppl.pdf). Divergence times obtained from TimeTree5, used as calibration points on the mitochondrial phylogenetic trees.
Supplementary Fig. S1 (03kinoshita_2024-0045_suppl.pdf). Evolutionary scenarios tested for D. legata and Tokudaia in the DIYABC-RF analysis. Populations are abbreviated AMA (Amami-Ohshima), TOK (Tokunoshima), and OKI (Okinawajima), respectively. We set the effective population size for some ancestral populations before the taxa split (Na and Nb). ra represents the admixture rate during an admixture event between two parental populations.
