Introduction
Continental archipelagos are suitable for examining causal connections between geographical barriers, disjunct distributions, and genetic differentiation (Hedges et al. Reference Hedges, Hass and Maxon1992; Hisheh et al. Reference Hisheh, Westerma and Schmitt1998; Atkins et al. Reference Atkins, Presto and Cronk2001; Ohdachi et al. Reference Ohdachi, Dokuchaev, Hasegawa and Masuda2001; Poulakakis et al. Reference Poulakakis, Lymberakis, Antonioub, Chalkia, Zourosb, Mylonasa and Valakos2003; Bittkau and Comes Reference Bittkau and Comes2005; Velo-Antón et al. Reference Velo-Antón, Zamudio and Cordero-Rivera2012). The Ryukyu Archipelago was the eastern margin of continental East Asia in the middle to late Miocene (15.97–5.333 million years ago). In the late Miocene to the early Pleistocene (5.333–0.774 million years ago), it fragmented into large islands and now includes approximately 140 islands scattered in an arc between Kyushu Island of Japan and Taiwan (Kimura Reference Kimura and Kimura2002; Osozawa et al. Reference Osozawa, Pavlis, Flower, Wakabayashi and Dilek2011). The islands are the product of the repeated formation and division of a land bridge to continental Asia that extended through Taiwan and Kyushu. Expansion and regression of the land bridge have been attributed to climatic oscillations in the Pliocene and Quaternary periods (Kimura Reference Kimura1996, Reference Kimura2000; Kizaki and Ohshiro Reference Kizaki and Oshiro1977, Reference Kizaki, Oshiro and Kizaki1980; Ujiie Reference Ujiie and Ujiie1990; Ohshiro Reference Ohshiro, Nishijima, Nishida, Shikatani and Shokita2003; Keally Reference Keally2005; Shinjo Reference Shinjo, Fujita, Arakaki, Denda, Hidaka, Hirose and Reimer2015). Plants and animals on the islands consist largely of relictual taxa that presumably differentiated from their relatives on the adjacent continent (Ota Reference Ota1998; Reference Ota2012). A number of biogeographic and evolutionary studies on the organisms of the Ryukyu Archipelago have been conducted and have enhanced our understanding of the biogeographic and evolutionary histories of organisms in continental archipelagos (e.g., Ota Reference Ota1998; Maekawa et al. Reference Maekawa, Lo, Kitade, Miura and Matsumoto1999; Seo et al. Reference Seo, Watabanabe, Hotta and Murakami2004; Chiang and Schaal Reference Chiang and Schaal2006; Seki et al. Reference Seki, Sakanashi, Kawaji and Kotaka2007; Nakamura et al. Reference Nakamura, Denda, Kokubugata, Huang, Peng and Yokota2015; Kaito and Toda Reference Kaito and Toda2016).
In general, the low dispersal ability of flightless organisms leads to a low rate of gene flow, suggesting accelerated differentiation among populations (Zera Reference Zera1981; Hansen Reference Hansen1983; Jablonski Reference Jablonski1986; Smith and Farrell Reference Smith and Farrell2006; Miura et al. Reference Miura, Torchin, Bermingham, Jacobs and Hechinger2012). For example, Ikeda et al. (Reference Ikeda, Nishikawa and Sota2011) reported that the loss of flight accelerates allopatric speciation among carrion beetles (Coleoptera: Silphidae). Flightless soil animals have also been found to exhibit extensive differentiation among allopatric populations (Tsurusaki Reference Tsurusaki2006a). In the Ryukyu Archipelago, for example, Kilungius insulanus (Hirst, Reference Hirst1911) (Opiliones: Epedanidae) is a species of harvestman that is genetically differentiated between Amami-Oshima and Okinawa Islands in the middle of the archipelago (Kumekawa et al. Reference Kumekawa, Ito, Miura, Yokoyama, Tebayashi, Arakawa and Fukuda2015). Zepedanulus ishikawai (Suzuki, Reference Suzuki1971) (Opiliones: Epedanidae) is another species of harvestman, belonging to the same family as K. insulanus. It is approximately 5 mm in body length, lives in humid places under forest litter, and is distributed in the southern Ryukyu Archipelago, namely Miyako, Ishigaki, Iriomote, Yonaguni, and Senkaku Islands (Fig. 1; Suzuki Reference Suzuki1973; Tsurusaki and Suzuki, Reference Tsurusaki, Suzuki and Aoki2015). This scattered distribution of Z. ishikawai might promote the differentiation of populations on different islands. In fact, allopatric speciation, where geographic isolation results in speciation, is widely accepted (Mayr Reference Mayr1963; Coyne and Orr Reference Coyne and Orr2004; Bolnick and Fitzpatrick Reference Bolnick and Fitzpatrick2007), and the accidental differentiation such as a genetic drift found between allopatric populations has led to substantial postzygotic isolation it is also known to bring about (Dobzhansky Reference Dobzhansky1937). The presence of physical barriers, such as oceans, rivers, or high mountains, will cause populations isolated by these boundaries to undergo different accumulations of variation over time from isolation. If reproductive isolation is acquired through this, it is likely that even if the isolation is subsequently resolved, the populations will not be able to interbreed with each other, and each population will have its own history.
Recently developed molecular techniques allow phylogenies to be reconstructed and the systematic relationships among species or populations within a given species to be evaluated. Application of mitochondrial DNA as a marker for phylogenetic and phylogeographic surveys enables evaluation of intraspecific population structures and differentiation (e.g., Avise Reference Avise2000, Reference Avise2009). Several phylogenetic and phylogeographic studies of harvestman have been based on mitochondrial DNA sequences (Giribet et al. Reference Giribet, Rambla, Carranza, Baguna, Riutort and Ribera1999; Thomas and Hedin Reference Thomas and Hedin2008; Sharma and Giribet Reference Sharma and Giribet2009, Reference Sharma and Giribet2011; Derkarabetian et al. Reference Derkarabetian, Steinmann and Hedin2010, Reference Derkarabetian, Ledford and Hedin2011; Giribet et al. Reference Giribet, Vogt, Perez, Sharma and Kury2010; Hedin and Thomas Reference Hedin and Thomas2010; Burns et al. Reference Burns, Hedin and Shultz2012; Sharma Reference Sharma2010; Sharma et al. Reference Sharma, Buenavente, Clouse, Diesmos and Giribet2012; Kumekawa et al. Reference Kumekawa, Ito, Tsurusaki, Hayakawa, Ohga and Yokoyama2014, Reference Kumekawa, Ito, Miura, Yokoyama, Tebayashi, Arakawa and Fukuda2015; Didodemico and Hedin Reference Didodemico and Hedin2016; Emata and Hedin Reference Emata and Hedin2016; Cruz-López et al. Reference Cruz-López, Monjaraz-Ruedas and Francke2019). It is also widely accepted that selecting suitable genetic markers is important in phylogeny and phylogeography studies, and a number of studies using mitochondrial DNA data have demonstrated the utility of the mitochondrial genes cytochrome c oxidase subunit I (CO1) and 16S ribosomal RNA (rRNA). Particularly for shallow-level animal molecular phylogeny and phylogeography, these genes enable accurate identification at the species (or lower) level through comparative analyses of sequence variations in short, standardised fragments of the genome (Hebert et al. Reference Hebert, Cywinska, Ball and deWaard2003; Jinbo et al. Reference Jinbo, Kato and Ito2011). We thus assessed the differentiation among Z. ishikawai populations on different islands using the sequences of CO1 and 16S rRNA in mitochondrial DNA.
The congruence of multiple genetic trees enables assessment of the reliability of a hypothesis regarding the evolutionary relationships of organisms. Indeed, even phylogenetic incongruence may provide a window into evolutionary processes that cannot be accessed using a single gene data set. Although a mitochondrial DNA phylogenetic tree allows clarification of the evolutionary relationships of a given species, it cannot clarify certain evolutionary processes because of its maternal inheritance. Theoretically, incongruence between a genetic tree and species phylogeny, or among different genetic trees, can result from a variety of evolutionary processes, such as conflation of orthology and paralogy, incomplete lineage sorting of ancestral polymorphism, and introgressive hybridisation (Rieseberg and Soltis Reference Rieseberg and Soltis1991; Wendel and Doyle Reference Wendel, Doyle, Soltis, Soltis and Doyle1998; Maddison and Knowles Reference Maddison and Knowles2006; Toews and Brelsford Reference Toews and Brelsford2012). Of these, introgressive hybridisation and incomplete lineage sorting frequently lead to phylogenetic incongruence among different data sets at lower taxonomic levels (Avise et al. Reference Avise, Shapiro, Daniel, Aquadro and Lansman1983; Cronn and Wendel Reference Cronn and Wendel2004). Incongruence among genetic trees should be confirmed using data from nuclear DNA in addition to mitochondrial DNA. In fact, Schlick-Steiner et al. (Reference Schlick-Steiner, Steiner, Seifert, Stauffer, Christian and Crozier2010) claimed that it is necessary to analyse various types of information to clarify the relationships among and differentiation within, populations within of a species. Among nuclear DNA sequences, 28S rRNA is particularly useful for differentiating taxa among groups of organisms (Nunn et al. Reference Nunn, Theisen, Christensen and Arctander1996; Friedrich and Tautz Reference Friedrich and Tautz1997; Babcock et al. Reference Babcock, Heraty, De Barro, Driver and Schmidt2001; Schmidt et al. Reference Schmidt, Driver and De Barro2006; Taylor and Rogers Reference Taylor and Rogers2015; Che et al. Reference Che, Wang, Shi, Du, Zhao, Lo and Wang2016), suggesting that this gene could be used to resolve the phylogenetic relationships of Z. ishikawai. In addition, the analysis of polymerase chain reaction–restriction fragment length polymorphisms has been used to show whether hybridisation and introgression occur in a variety of species (Hollingsworth et al. Reference Hollingsworth, Bailey, Hollingsworth and Ferris1999; Nijman et al. Reference Nijman, Bradley, Hanotte, Otsen and Lenstra1999; Mráz et al. Reference Mráz, Chrtek, Fehrer and Plačková2005; Oleksa et al. Reference Oleksa, Chybicki, Tofilski and Burczyk2011; Vaini et al. Reference Vaini, Grisolia, Prado and Porto-Foresti2014; Kumekawa et al. Reference Kumekawa, Miura, Fujimoto, Ito, Arakawa, Yokoyama and Fukuda2019).
Therefore, we evaluated the differentiation of Z. ishikawai populations among islands using 28S rRNA in nuclear DNA and compared the results with gene fragments of mitochondrial DNA. We also analysed the morphological characters of Z. ishikawai to determine whether the populations could be distinguished morphologically.
Materials and methods
Sample collection
We collected Z. ishikawai by hand from three major islands of the Yaeyama Islands, in the southernmost part of the Ryukyu Archipelago (Ishigaki, Iriomote, and Yonaguni Islands; Table 1; Fig. 2).
DNA extraction, amplification, sequencing, and polymerase chain reaction–restriction fragment length polymorphism analysis
We extracted DNA from a leg using DNeasy kits (Qiagen, Valencia, California, United States of America), according to the manufacturer’s protocol for animal tissue samples. The isolated DNA was resuspended in an appropriate volume of Tris-ethylenediaminetetraacetic acid buffer (pH 8.0) and stored at −20 °C until use. For all specimens, we amplified the CO1 and 16S rRNA genes in mitochondrial DNA and the 28S rRNA gene in nuclear DNA. The CO1 gene fragments were amplified and sequenced using primers published previously (Kumekawa et al. Reference Kumekawa, Ito, Tsurusaki, Hayakawa, Ohga and Yokoyama2014). The 16S rRNA and 28S rRNA gene fragments were amplified and sequenced using the primers of Xiong and Kocher (Reference Xiong and Kocher1991) and Mallatt and Sullivan (Reference Mallatt and Sullivan1998), respectively. Amplification of DNA followed the method of Kumekawa et al. (Reference Kumekawa, Murjoko, Hayakawa, Ohga, Mori and Miyazaki2013). The reaction solution (50 µL) contained 50 ng of total DNA, 10 mM Tris (hydroxymethyl) aminomethane hydrochloride buffer (pH 8.3) with 50 mM potassium chloride and 1.5 mM magnesium dichloride, 0.2 mM of each deoxyribonucleotide triphosphate, 1.25 U Taq DNA polymerase (TaKaRa, Tokyo, Japan), and 0.5 µM of each primer. The following thermal cycle was applied: incubation at 94 °C for 10 seconds, followed by 45 cycles of incubation at 94 °C for 1.5 minutes, 48 °C for 2 minutes, and 72 °C for 3 minutes, with a final extension at 72 °C for 15 minutes. After amplification, the reaction mixtures were subjected to electrophoresis in 1% low melting–temperature agarose gels, and the products were purified using QuickSpin kits (Qiagen) according to the manufacturer’s specifications. We sequenced the purified polymerase chain reaction products using a BigDye Terminator Cycle Sequencing Kit (Applied Biosystems, Tokyo, Japan) and ABI PRISM 3100-Avant Genetic Analyzer (Applied Biosystems) according to the manufacturers’ instructions.
We also carried out polymerase chain reaction–restriction fragment length polymorphism analysis, because autapomorphic characters (clade II, see Results) of 28S rRNA are involved in the restriction site HaeIII (GGCC). We designed the following primers to assess the length of the restriction fragments and conducted polymerase chain reaction as above, using the primers Zepe28S-RFLP-F: 5'-GCC CTC GAA GTT CGG AGG CG-3' and Zepe28S-RFLP-R: 5'-CAC CCG GCT ACC TAC CGG TT-3'. The amplified products were digested by HaeIII at 37 °C for more than 1 hour. The digested DNAs were separated in a 1.5% agarose gel and the size of each band was determined.
Data analysis
To construct phylogenetic trees for Z. ishikawai and related species (Table 1), 58 nucleotide sequences were aligned using the ClustalW program (Thompson et al. Reference Thompson, Higgins and Gibson1994) implemented in MEGA, version 6, software (Tamura et al. Reference Tamura, Stecher, Peterson, Filipski and Kumar2013). The positions containing gaps and missing data were eliminated manually. A Bayesian analysis was conducted using MrBayes, version 3.1.2, software (Ronquist and Huelsenbeck Reference Ronquist and Huelsenbeck2003). We considered the model with the lowest log-likelihood score to be the most accurate description of the substitution pattern. The Akaike information criterion value of each model was calculated. The best-fit model (mitochondrial DNA, GTR + I + G model; nuclear DNA, F81 model) selected by hLRT in MrModeltest 2.3 (Nylander Reference Nylander2004) and PAUP * 4.0b10 (Swofford Reference Swofford2002) software was used. We ran four Markov chains for 50 million generations and sampled Markov chains every 100 generations; 20% of the runs were discarded as burn-in. The effective sample size (ESS) of each parameter indicated by Tracer version 1.6, software (Rambaut et al. Reference Rambaut, Suchard, Xie and Drummond2014) confirmed that the number of Markov chain generations was sufficient (ESS > 200). The maximum-likelihood method was conducted using MEGA, version 6.06, software (Tamura et al. Reference Tamura, Stecher, Peterson, Filipski and Kumar2013). The model test function in MEGA was used to select the appropriate model for maximum-likelihood analysis by considering the model with the lowest Bayesian information criterion scores as the most accurate description of the substitution pattern. For each model, the Akaike information criterion–corrected value, maximum-likelihood value (lnL), and the number of parameters (including branch lengths) are presented (Nei and Kumar Reference Nei and Kumar2000). The nonuniformity of evolutionary rates among sites may be modelled using a discrete gamma distribution (+ G) with five rate categories and a parameter for invariable sites (+ I). The relative values of instantaneous r should be considered when evaluating them. For simplicity, the sum of r values is made equal to 1 for each model. To estimate maximum-likelihood values, a tree topology was automatically computed. Nodal support in maximum-likelihood trees was measured by bootstrapping (1000 samples).
The divergence times of the clades were estimated for all individuals used for phylogenetic analyses. Divergence times were inferred using an uncorrelated relaxed lognormal clock, implemented in BEAST, version 2.6.2, software (Bouckaert et al. Reference Bouckaert, Vaughan, Barido-Sottani, Duchêne, Fourment and Gavryushkina2019). Because there are no mitochondrial DNA molecular clocks available for Zepedanulus or related harvestman species, we conservatively adopted a wide range of estimated divergence rates for arthropods from 0.21% (Clouse and Wheeler Reference Clouse and Wheeler2014) to 3.34% per million years (Pons and Vogler Reference Pons and Vogler2005) for the CO1 gene, and 0.45% (Gomez-Zurita et al. Reference Gómez-Zurita, Juan and Petitpierre2000) to 1.9% per million years (Percy et al. Reference Percy, Page and Cronk2004) for the 16S gene. Prior distributions for mean substitution rates were sampled from a log normal distribution (clock.rate = 0.006, sigma = 0.8 for the CO1 gene; clock.rate = 0.005, sigma = 0.5 for the 16S gene). We used the same substitution model as the Bayesian analysis and used the Yule process as a tree prior. We constrained the monophyly of two clades. The analysis was run for 500 million generations, sampled every 50 000 steps, and the first 1000 samples were discarded as burn-in. We used Tracer, version 1.6 (Rambaut et al. Reference Rambaut, Suchard, Xie and Drummond2014), and FIGTREE, version 1.4.2 (Rambaut Reference Rambaut2014), software to check for convergence and to visualise the results.
Morphological analysis
For morphological analysis, the body length, carapace length, carapace width, femur lengths of the first to fourth legs, spine length, chelicerae length, and palpal–tarsus length were measured for each Z. ishikawai specimen (Fig. 2). We measured the morphological characters of 139 Z. ishikawai individuals (39 from Ishigaki Island, 34 from Iriomote Island, and 66 from Yonaguni Island). Measurements were performed using cellSens standard imaging software (Olympus Co., Tokyo). Statistical analysis was conducted by Tukey–Kramer test (P < 0.05).
Results
Molecular analysis of Zepedanulus ishikawai and related species
We reconstructed phylogenetic relationships using 43 samples (27 from Ishigaki Island, eight from Yonaguni Island, and six from Iriomote Island) of Z. ishikawai and related species (Table 1), including previously published sequences of the CO1 and 16S rRNA genes in mitochondrial DNA. Of these, even though only one sample of Epedanellus tuberculatus was used in the analysis of the 16S rRNA and CO1 sequences, four samples of Penaeus japonicus and two of K. insulanus were analysed because previous studies indicated that these species are divided into four and two groups, respectively (Kumekawa et al. Reference Kumekawa, Ito, Tsurusaki, Hayakawa, Ohga and Yokoyama2014, Reference Kumekawa, Ito, Miura, Yokoyama, Tebayashi, Arakawa and Fukuda2015). The combined length of the CO1 and 16S rRNA sequences in Z. ishikawai was 765 bp (16S rRNA, 287 bp; CO1, 478 bp). We reconstructed phylogenetic trees based on the maximum-likelihood (Fig. 3) and Bayesian inference (Fig. 4) methods. The branching patterns of the Z. ishikawai maximum-likelihood phylogenetic tree were identical to those of the tree that was based on Bayesian inference (Figs. 3 and 4). The phylogenetic trees based on the CO1 and 16S rRNA genes indicated that individuals of Z. ishikawai composed a monophyletic group with a high bootstrap value (96%) and high posterior probability (1.00). The species Z. ishikawai was divided into two clades, henceforth referred to as clade I and clade II (Fig. 3). Clade I was not strongly supported in the maximum-likelihood tree (68%), but Bayesian analysis indicated a high posterior probability for the clade (1.00). Clade I involved three monophyletic subgroups clades IA, IB, and IC. Clade IA consisted of individuals from Ishigaki Island and was the sister group of clade IB, which comprised from Iriomote Island. Clade IC was the sister group to the combined clade of IA + IB, and comprised samples from Yonaguni Island. These three groups had high bootstrap values (70–100%) and posterior probabilities (0.99–1.00). Clade II, which had high support values (100% maximum-likelihood bootstrap and 1.00 Bayesian posterior probability values), also comprised three monophyletic subgroups clades IIA, IIB, and IIC. The branching patterns of clade II were identical to those of clade I, and individuals of Yonaguni Island, Ishigaki Island, and Iriomote Island constituted monophyletic groups of clades IIA, IIB, and IIC with high support values (100% maximum-likelihood bootstrap and 1.00 Bayesian posterior probability values), respectively.
The phylogenetic tree based on the 28S rRNA gene in nuclear DNA was constructed using Z. ishikawai and Tokunosia tenuipes, which was revealed to be a sister species of Z. ishikawai based on the mitochondrial DNA phylogenetic tree (Figs. 5 and 6). The length of the Z. ishikawai 28S rRNA sequence was 978 base pairs. As for mitochondrial DNA data sets, we reconstructed phylogenetic trees based on the maximum-likelihood (Fig. 5) and Bayesian inference (Fig. 6) methods. The major branching patterns of the phylogenetic trees based on 28S rRNA were similar to those of the mitochondrial DNA trees. The phylogenetic trees indicated that individuals of Z. ishikawai comprised a monophyletic group and were divided into two groups. These groups corresponded to clades I and II in the mitochondrial DNA phylogenetic tree.
Our phylogenetic analyses of mitochondrial DNA and 28S rRNA sequences both indicated the presence of two identical clades. Clades I and II were sympatrically distributed on all three islands examined, and no incongruence of mitochondrial DNA and nuclear DNA was observed. Therefore, whether hybrids of individuals of clades I and II on each of the three islands could form is unclear because assessing heterogeneity at polymorphic sites in nuclear DNA using the results of Sanger sequencing can be problematic (Simsek et al. Reference Simsek, Tanira, Al-Baloushi, Al-Barwani, Lawatia and Bayoumi2001). Based on the above analyses, we found that the 28S rRNA sequences of clade II had HaeIII restriction sites, whereas those of clade I did not (Fig. 7). We conducted polymerase chain reaction–restriction fragment length polymorphism to evaluate hybridisation and introgression between clade I and clade II of Z. ishikawai. No hybridisation or introgression between clade I and clade II of Z. ishikawai was detected on the three islands (Fig. 7).
The dated tree of Z. ishikawai is shown in Figure 8. The tree exhibited that the divergence of clade I and clade II occurred about 28 million years ago (47–13 million years ago: 95% 4-hydroxyphenylpyruvate dioxygenase). In addition, the divergence time of the three groups of clade I occurred about 25 million years ago, and the divergence time of the three groups of clade II occurred about 10 million years ago.
Morphological analysis of Zepedanulus ishikawai
All examined morphological characters of clade I were significantly larger than those of clade II (Table 2; Figs. 9–11). No significant differences within the same clade were found among the islands. The spine length could not be measured for individuals of clade I because they lacked an eye-spine and spine on the second suctal area. Therefore, the spine length was measured only for individuals of clade II.
Abbreviations: Fe1L, femur length, first leg; Fe2L, femur length, second leg; Fe3L, femur length, third leg; Fe4L, femur length, fourth leg; BL, body length; CL, carapace length; CW, carapace width; SL, spine length; PTL, palpal–tarsus length; Chel.L, chelicerae length.
Discussion
A phylogenetic tree of the mitochondrial DNA sequences of Z. ishikawai indicated that the species consisted of two phylogenetic groups clades I and II which was supported by the phylogenetic tree based on the 28S rRNA gene. Our morphological analyses indicated that individuals of clade I could be distinguished from those of clade II not only quantitatively (significant differences in body length, carapace length, carapace width, chelicerae length, palpal–tarsus length, spine length, and the femur lengths of the first to fourth legs) but also qualitatively (presence of an eye-spine and a spine in the second suctal area). Both groups contain individuals on three islands (Ishigaki, Iriomote, and Yonaguni), suggesting that both groups are distributed sympatrically throughout the southern Ryukyu Archipelago. This has consequences for the diversification of Z. ishikawai, because two morphologically and genetically differentiated groups were apparent at the species level. We did not detect incongruence of the branching patterns of the phylogenetic trees based on mitochondrial DNA and nuclear DNA sequences (Figs. 3, 4, 5, and 6) or on nuclear DNA heterogeneity by polymerase chain reaction–restriction fragment length polymorphism, despite collecting samples from sympatric populations on each island (Fig. 7). In addition, the body colour differed between clades on all the islands (clade I, yellow; clade II, dark brown on Ishigaki and Iriomote Islands, and clade I, brown; clade II, black on Yonaguni Island; Fig. 12). We therefore infer that these two groups represent reproductively isolated species. Suzuki (Reference Suzuki1971) reported that the total body length of Z. ishikawai was 3.65 mm in males and 2.81–3.24 mm in females and that the legs of the type specimen of Z. ishikawai were comparatively long. Although a new taxonomic treatment is needed, based on Suzuki’s description of legs and total body length, clade I is Z. ishikawai and a new name is needed for clade II.
The phylogeographic patterns of organisms on the islands have been influenced by the formation of the continental island systems, which offer relatively simple arenas for the evolutionary dynamics of community assembly, because harvestmen are generally small and spatially isolated beyond their vagility. Prior phylogeographic studies reported that the current distribution of biological diversity cannot be understood without information about how organisms responded to historical changes in geological and climatic conditions (e.g., Arbogast and Kenagy Reference Arbogast and Kenagy2001; Fukuda et al. Reference Fukuda, Yokoyama and Ohashi2001; Griffin and Barrett Reference Griffin and Barrett2004; Hewitt Reference Hewitt2004; Vink et al. Reference Vink, Thomas, Paquin, Hayashi and Hedin2005; Matsumura et al. Reference Matsumura, Yokoyama, Fukuda and Maki2009; Minamiya et al. Reference Minamiya, Yokoyama and Fukuda2009; Fukuda et al. Reference Fukuda, Song, Ito, Hayakawa, Minamiya, Kanno, Hayakawa, Minamiya and Yokoyama2011; Maura et al. Reference Maura, Salvi, Bologna, Nascetti and Canestrelli2014; Mezzasalma et al. Reference Mezzasalma, Dall’Asta, Loy, Cheylan, Lymberakis and Zuffi2015). Kumekawa et al. (Reference Kumekawa, Ito, Tsurusaki, Hayakawa, Ohga and Yokoyama2014) suggested that Pseudobiantes japonicus comprised four allopatric groups with morphological and genetic differences, and based on the current distribution and historical changes of suitable habitats, these groups had undergone independent northward spread by various routes from allopatric refugia distributed along the coasts of the Pacific Ocean in the western part of the main island of Japan. It is intriguing to consider how Z. ishikawai colonised the southern part of the Ryukyu Archipelago despite the presence of a barrier to expansion. The southern part of the Ryukyu Archipelago formed the eastern margin of continental East Asia in the mid to late Miocene (15.97–5.333 million years ago) and then fragmented into large islands, breaking the direct land connection to China, in the late Miocene to the early Pleistocene (5.333–0.774 million years ago; Kimura Reference Kimura and Kimura2002; Osozawa et al. Reference Osozawa, Shinjo, Armid, Watanabe, Horiguchi and Wakabayashi2012). This suggests that Z. ishikawai colonised the archipelago region from continental East Asia during this period. Subsequently, sea levels ranged from 200 m higher to 200 m lower than present-day sea level, likely due to climatic oscillations during the Quaternary period. Lower islands were then either largely or completely submerged during the interglacial age (Kimura Reference Kimura1996). During this period, subsidence created two deep-water passages through the island arc, namely the Tokara Gap and the Kerama Gap, which divide the Ryukyu area into three groups: the northern, central, and southern Ryukyus (Osozawa et al. Reference Osozawa, Shinjo, Armid, Watanabe, Horiguchi and Wakabayashi2012). Although the islands were repeatedly fragmented into smaller islands and reconnected in the middle to late Pleistocene (0.774–0.0017 million years ago), the existence of land connections across the two gaps after their opening is unlikely, even during the Quaternary glacial sea level minima, because the ocean depth is > 1000 m between the constituent island groups (Ota Reference Ota1998; Kawana Reference Kawana and Kimura2002; Osozawa et al. Reference Osozawa, Shinjo, Armid, Watanabe, Horiguchi and Wakabayashi2012). These gaps acted as strong barriers, even for flying insect species, and resulted in the generation of endemic species (e.g., Osozawa et al. Reference Osozawa, Su, Oba, Yagi, Watanabe and Wakabayashi2013, Reference Osozawa, Takáhashi and Wakabayashi2015). Tsurusaki (Reference Tsurusaki2006b) reported that Opiliones fauna is divided into three categories along the above island groups, based on cluster analysis using the similarity index. Therefore, Z. ishikawai could not cross the Kerama Gap and is included in the third category of Tsurusaki (Reference Tsurusaki2006b). Alternatively, the land connection between Yonaguni Island, which is about 130 km from Taiwan Island and is the westernmost area in the distribution of Z. ishikawai, and Taiwan Island ceased in the early Pleistocene (2.58 million years ago), roughly contemporaneous with disconnection of the southern Ryukyus from continental East Asia (Koba Reference Koba1992; Osozawa et al. Reference Osozawa, Shinjo, Armid, Watanabe, Horiguchi and Wakabayashi2012). The distance between the two islands could not prevent the movement of flying organisms, such as winged insects and birds, nor the spread of plants with long-distance dispersal abilities (Hatusima Reference Hatusima1975; Ota Reference Ota1998). However, the flightlessness and low vagility of Z. ishikawai could restrict its movement between Yonaguni and Taiwan Islands, and this species likely experienced the long-standing isolation of the southern Ryukyus from the surrounding islands. The deep genetic divergences we observed within the Z. ishikawai species complex suggest that the geological history of the Ryukyus played key roles in driving geographic isolation of populations since their divergence from their common ancestor in the archipelago.
The endemism of Z. ishikawai in the southern Ryukyus is supported by the historical isolation of the area from surrounding regions. However, our phylogenetic results indicated that Z. ishikawai comprises two groups and three subgroups; the former are considered different species, and the latter are considered differentiation among islands. Based on the geographic history of the area and the divergence time, we considered the following hypotheses for the formation of the current distribution of Z. ishikawai: (1) single colonisation by a common ancestor and speciation into clades I and II and (2) independent colonisation of two differentiated ancestors. The former hypothesis seems unlikely, based on the estimated divergence time of clade I and clade II of > 20 million years ago (Fig. 8), when the southern part of the Ryukyu Archipelago was connected to Taiwan and the eastern margin of continental East Asia (Koba Reference Koba1992; Osozawa et al. Reference Osozawa, Shinjo, Armid, Watanabe, Horiguchi and Wakabayashi2012). In this geographical setting, speciation within the limited prospective area of the southern part of Ryukyu archipelago was unlikely even with the limited dispersal ability of harvestmen, thereby supporting the latter hypothesis. The ancient divergence between clade I and clade II suggests that each clade of Z. ishikawai had accumulated variations independently, leading to considerable morphological differentiation. We estimate that the colonising ancestors of both Z. ishikawai clades originated from continental Southeast Asia, because the genus Zepedanulus is also recorded in, for instance, northern Thailand and Malacca (Roewer Reference Roewer1927; Suzuki Reference Suzuki1981). However, it is unclear whether the ancestors of Z. ishikawai colonised the Ryukyu Archipelago region via Taiwan or invaded the archipelago directly from continental Asia. Systematic studies and phylogenetic analyses using a range of Epedanidae species are required to assess the history of diversification of Z. ishikawai and related species.
Our phylogenetic results suggest that Z. ishikawai did not frequently colonise islands in the southern part of the Ryukyu Archipelago after divergence of the two clades, because all monophyletic subgroups within each clade of this species comprised samples from an island, and individuals from the same island did not show paraphyletic or polyphyletic assemblies. Therefore, each subgroup of Z. ishikawai has an independent history of differentiation on each island, without additional introductions from other islands. Indeed, geographical variations among islands are also seen in other organisms in the Ryukyu Archipelago (Maki et al. Reference Maki, Yamashiro and Matsumura2003; Kiyoshi Reference Kiyoshi2008; Osozawa et al. Reference Osozawa, Su, Oba, Yagi, Watanabe and Wakabayashi2013, Reference Osozawa, Takáhashi and Wakabayashi2015; Hosokawa et al. Reference Hosokawa, Nikoh and Fukatsu2014; Osozawa et al. Reference Osozawa, Takáhashi and Wakabayashi2015). In the case of harvestmen, Kumekawa et al. (Reference Kumekawa, Ito, Miura, Yokoyama, Tebayashi, Arakawa and Fukuda2015) reported that K. insulanus in the central part of the Ryukyu Archipelago exhibited genetic differentiation between Amami-Oshima and Okinawa Islands based on the sequence of the CO1 in mitochondrial DNA. Although we detected morphological differences between clades of Z. ishikawai, whether differentiation occurred among subgroups of each clade is unclear because few individuals were collected from each island. Further analyses using larger numbers of samples from the islands are needed to clarify the morphological differences within each clade of Z. ishikawai.
Acknowledgements
The authors thank Drs K. Ito, Y. Minamiya, and H. Hayakawa for discussion in this study, and C. Uemoto, Y. Ozaki, K. Yoshioka, Y. Kamakura, Y. Mori, Y. Nibuno, T. Miyashita, and F. Maekawa for their field and laboratory assistance. The authors also thank Drs. N. Tsurusaki and T. Suguro for providing samples. This study was partly supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Science and Culture of Japan.