Morphological and molecular phylogenetic analyses of Zepedanulus ishikawai (Arachnida: Opiliones: Laniatores: Epedanidae) in the southern part of the Ryukyu Archipelago

Abstract

Harvestmen (Arachnida: Opiliones) are soil animals with extremely low dispersal abilities that experienced allopatric differentiation. To clarify the morphological and phylogenetic differentiation of the endemic harvestman Zepedanulus ishikawai (Suzuki, 1971) (Laniatores: Epedanidae) in the southern part of the Ryukyu Archipelago, we conducted molecular phylogenetic analyses and divergence time estimates based on CO1 and 16S rRNA sequences of mtDNA, the 28S rRNA sequence of nrDNA, and the external morphology. A phylogenetic tree based on mtDNA sequences indicated that individuals of Z. ishikawai were monophyletic and were divided into clade I and clade II. This was supported by the nrDNA phylogenetic tree. Although clades I and II were distributed sympatrically on all three islands examined (Ishigaki, Iriomote, and Yonaguni), heterogeneity could not be detected by polymerase chain reaction–restriction fragment length polymorphism of nrDNA, indicating that clades I and II do not have a history of hybridisation. Also, several morphological characters differed significantly between individuals of clade I and clade II. The longstanding isolation of the southern Ryukyus from the surrounding islands enabled estimation of the original morphological characters of both clades of Z. ishikawai.

Introduction

Continental archipelagos are suitable for examining causal connections between geographical barriers, disjunct distributions, and genetic differentiation (Hedges et al. 1992; Hisheh et al. 1998; Atkins et al. 2001; Ohdachi et al. 2001; Poulakakis et al. 2003; Bittkau and Comes 2005; Velo-Antón et al. 2012). 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 2002; Osozawa et al. 2011). 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 1996, 2000; Kizaki and Ohshiro 1977, 1980; Ujiie 1990; Ohshiro 2003; Keally 2005; Shinjo 2015). Plants and animals on the islands consist largely of relictual taxa that presumably differentiated from their relatives on the adjacent continent (Ota 1998; 2012). 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 1998; Maekawa et al. 1999; Seo et al. 2004; Chiang and Schaal 2006; Seki et al. 2007; Nakamura et al. 2015; Kaito and Toda 2016).

In general, the low dispersal ability of flightless organisms leads to a low rate of gene flow, suggesting accelerated differentiation among populations (Zera 1981; Hansen 1983; Jablonski 1986; Smith and Farrell 2006; Miura et al. 2012). For example, Ikeda et al. (2011) 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 2006a). In the Ryukyu Archipelago, for example, Kilungius insulanus (Hirst, 1911) (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. 2015). Zepedanulus ishikawai (Suzuki, 1971) (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 1973; Tsurusaki and Suzuki, 2015). 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 1963; Coyne and Orr 2004; Bolnick and Fitzpatrick 2007), 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 1937). 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.

Fig. 1. Distribution of Z. ishikawai in the Ryukyu Archipelago, Japan.

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 2000, 2009). Several phylogenetic and phylogeographic studies of harvestman have been based on mitochondrial DNA sequences (Giribet et al. 1999; Thomas and Hedin 2008; Sharma and Giribet 2009, 2011; Derkarabetian et al. 2010, 2011; Giribet et al. 2010; Hedin and Thomas 2010; Burns et al. 2012; Sharma 2010; Sharma et al. 2012; Kumekawa et al. 2014, 2015; Didodemico and Hedin 2016; Emata and Hedin 2016; Cruz-López et al. 2019). 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. 2003; Jinbo et al. 2011). 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 1991; Wendel and Doyle 1998; Maddison and Knowles 2006; Toews and Brelsford 2012). Of these, introgressive hybridisation and incomplete lineage sorting frequently lead to phylogenetic incongruence among different data sets at lower taxonomic levels (Avise et al. 1983; Cronn and Wendel 2004). Incongruence among genetic trees should be confirmed using data from nuclear DNA in addition to mitochondrial DNA. In fact, Schlick-Steiner et al. (2010) 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. 1996; Friedrich and Tautz 1997; Babcock et al. 2001; Schmidt et al. 2006; Taylor and Rogers 2015; Che et al. 2016), 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. 1999; Nijman et al. 1999; Mráz et al. 2005; Oleksa et al. 2011; Vaini et al. 2014; Kumekawa et al. 2019).

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).

Table 1. Summary of Zepedanulus ishikawai analysed in this study. Specimen number corresponds to our phylogenetic results (Figs. 3, 4, 5, and 6).

Fig. 2. (A) Spined type and (B) nonspined type, and characters of Zepedanulus ishikawai (collected from Yonaguni Island): A, BL – body length; B, CL – carapace length; C, CW – carapace width; D, Fe1L – femur length, first leg; E, Fe2L – femur length, second leg; F, Fe3L – femur length, third leg; G, Fe4L – femur length, fourth leg; H, SL – spine length; I, Chel.L – chelicerae length; and J, PTL – palpal–tarsus length, and the morphology of the second scutal area. Arrowhead indicates the spine. Scale bar = 1 mm.

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. 2014). The 16S rRNA and 28S rRNA gene fragments were amplified and sequenced using the primers of Xiong and Kocher (1991) and Mallatt and Sullivan (1998), respectively. Amplification of DNA followed the method of Kumekawa et al. (2013). 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. 1994) implemented in MEGA, version 6, software (Tamura et al. 2013). 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 2003). 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 2004) and PAUP * 4.0b10 (Swofford 2002) 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. 2014) 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. 2013). 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 2000). 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. 2019). 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 2014) to 3.34% per million years (Pons and Vogler 2005) for the CO1 gene, and 0.45% (Gomez-Zurita et al. 2000) to 1.9% per million years (Percy et al. 2004) 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. 2014), and FIGTREE, version 1.4.2 (Rambaut 2014), 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. 2014, 2015). 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.

Fig. 3. Phylogenetic tree of Zepedanulus ishikawai and related species with outgroups by the maximum-likelihood method based on the sequences of CO1 and 16S rRNA in mtDNA. Numbers above or below branches indicate bootstrap values. Symbols indicate the origin of samples: ○, Ishigaki Island; △, Iriomote Island; □, Yonaguni Island.

Fig. 4. Phylogenetic tree of Zepedanulus ishikawai and related species with outgroups by Bayesian inference based on the sequences of CO1 and 16S rRNA in mtDNA. Numbers above or below branches indicate posterior probabilities. Symbols indicate the origin of samples: ○, Ishigaki Island; △, Iriomote Island; □, Yonaguni Island.

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.

Fig. 5. Phylogenetic tree of Zepedanulus ishikawai and related species with outgroups by the maximum-likelihood method based on the sequence of 28S rRNA in nrDNA. Numbers above or below branches indicate the bootstrap values. Symbols indicate the origin of samples: ○, Ishigaki Island; △, Iriomote Island; □, Yonaguni Island.

Fig. 6. Phylogenetic tree of Zepedanulus ishikawai and related species with outgroups by Bayesian inference based on the sequence of 28S rRNA in nrDNA. Numbers above or below branches indicate the posterior probabilities. Symbols indicate the origin of samples: ○, Ishigaki Island; △, Iriomote Island; □, Yonaguni Island.

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. 2001). 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).

Fig. 7. A, Restriction sites in the 28S rRNA sequence determined by polymerase chain reaction–restriction fragment length polymorphism (M, HaeIII site) and B, polymerase chain reaction–restriction fragment length polymorphism profile of clades I and II of Zepedanulus ishikawai (M, size markers).

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.

Fig. 8. Estimated divergence time (over millions of years ago; MYA) for the two clades of Zepedanulus ishikawai. Horizontal bars represent the upper and lower 95% of highest posterior density values. Symbols indicate the origin of samples: ○, Ishigaki Island; △, Iriomote Island; □, Yonaguni Island.

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.

Table 2. The measurement (average ± standard deviation) of Zepedanulus ishikawai.

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.

Fig. 9. Histograms of the morphological characteristics of Zepedanulus ishikawai in Iriomote Island. White and grey colours indicate clade I and clade II, respectively. 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; Chel.L, chelicerae length; PTL, palpal–tarsus length.

Fig. 10. Histograms of the morphological characteristics of Zepedanulus ishikawai in Ishigaki Island. White and grey colours indicate clade I and clade II, respectively. 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; Chel.L, chelicerae length; PTL, palpal–tarsus length.

Fig. 11. Histograms of the morphological characteristics of Zepedanulus ishikawai in Yonaguni Island. White and grey colours indicate clade I and clade II, respectively. 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; Chel.L, chelicerae length; PTL, palpal–tarsus 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 (1971) 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.

Fig. 12. Distribution of morphological characteristics in a simplified phylogenetic tree of Zepedanulus ishikawai and related species.

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 2001; Fukuda et al. 2001; Griffin and Barrett 2004; Hewitt 2004; Vink et al. 2005; Matsumura et al. 2009; Minamiya et al. 2009; Fukuda et al. 2011; Maura et al. 2014; Mezzasalma et al. 2015). Kumekawa et al. (2014) 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 2002; Osozawa et al. 2012). 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 1996). 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. 2012). 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 1998; Kawana 2002; Osozawa et al. 2012). These gaps acted as strong barriers, even for flying insect species, and resulted in the generation of endemic species (e.g., Osozawa et al. 2013, 2015). Tsurusaki (2006b) 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 (2006b). 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 1992; Osozawa et al. 2012). 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 1975; Ota 1998). 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 1992; Osozawa et al. 2012). 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 1927; Suzuki 1981). 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. 2003; Kiyoshi 2008; Osozawa et al. 2013, 2015; Hosokawa et al. 2014; Osozawa et al. 2015). In the case of harvestmen, Kumekawa et al. (2015) 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.