We conducted phylogenetic analyses using mitochondrial COI gene sequences of Tylos granuliferus, a semiterrestrial coastal isopod in East Asia, to clarify possible phylogeographic patterns and to assess relationships between present and past marine environments and genetic population structures. Our results strongly suggest the presence of four clades of T. granuliferus, one of which consists of three subclades. The distribution pattern of clades and subclades seems to have been affected by ocean current activities. Our results also suggest that historical changes in oceanic environments and the subsequence bifurcation of current streamlines affected the first and second divergences of T. granuliferus during the late Miocene and near the beginning of the Pliocene, respectively.
INTRODUCTION
The processes underlying marine biodiversity constitute a topic of considerable scientific focus, as both phyletic and genetic diversity are higher among marine animals than among terrestrial animals (Gray, 1997). In coastal areas, species diversify along coastlines, showing unique distribution patterns (e.g., Burton and Lee, 1994). Because the distributions of their populations can be treated as one-dimensional (Pielou, 1979) and each population is expected to disperse through restricted routes, we can analyze the relationships between these distributions and phylogenetic, ecological, and geographic factors with comparative ease. Thus, phylogeographic research into coastal animals provides fruitful information not only about the individual biogeographic studies themselves, but also about the development of the marine biota.
The distributions of populations are affected by both environmental and historical factors (Avise, 2000). The effects of environmental factors, such as ocean currents, on coastal animals must be strong, especially because of their strong transportation capacities, and some coastal species show current-related distribution patterns. Furthermore, organisms that have relatively low migration capacities or frequencies, such as those with short planktonic periods or that inhabit areas far from offshore environments, show complex genetic structures within a single species, which correspond to past oceanic conditions (Dawson, 2001).
In contrast, some coastal species that inhabit semiterrestrial regions do not show current-related distribution patterns. Populations of these coastal species seem to show geohistorical structures rather than structures affected by currents and other oceanic conditions (e.g., Thoropa spp.: Fitzpatrick et al., 2009). Phylogenetic analyses of mitochondrial DNA data have revealed that some semiterrestrial coastal isopods display noncurrent-related and complex distribution patterns in several areas (Ligia spp.: Itani, 2000; Hurtado et al., 2010).
The genus Tylos, a semiterrestrial coastal isopod, is an ideal model organism with which to assess the relationships between present/past marine environments and genetic population structures. It inhabits sandy beaches in tropical and temperate regions (Kensley, 1974), not far from the high-water mark. It can only withstand seawater for several hours (Brown and Odendaal, 1994). Lacking a planktonic stage, Tylos also has low migratory ability, and individuals hide under the sand for most of the day (Brown and Odendaal, 1994). Thus, members of the genus are expected to show both current-related and noncurrent-related distributional patterns.
The Japanese Archipelago is an important area in which to study the ecological and evolutionary processes of coastal species because it has a long, branched shape and two warm ocean currents, the Tsushima and Kuroshio Currents, flow around it into the surrounding Sea of Japan and Pacific Ocean, respectively. Two cold currents also exist in the area: the Oyashio Current, located along the northeastern side of the archipelago; and the Liman Current, flowing along the northeastern coast of continental Asia (Fig. 1). The presence and collision of these currents influence the oceanic environments and contribute to the formation of the marine climatic zones of Japan (Nishimura, 1981; Fujikura et al., 2010).
Fig. 1.
Sampling localities for Tylos granuliferus, with population numbers. The numbers of the localities are as in Table 1. Closed circles indicate populations of specimens of C-I, open pentagons indicate C-II, open squares indicate C-III, and open triangles indicate C-IV. Open lozenges indicate localities at which individuals of C-I and C-III coexist. The Japanese islands (a: Hokkaido Island; b: Honshu Island; c: Shikoku Island; d: Kyushu Island; RA, Ryukyu Archipelago) are shown with the water currents and seas around the islands (TC: Tsushima Current; KC: Kuroshio Current; OC: Oyashio Current; LC: Liman Current; YS: Yellow Sea; ECS: East China Sea; SIS: Seto Inland Sea; TsS: Tsushima Strait; ToS: Tokara Strait; and KrS: Kerama Strait).
In this study, we examined the East Asian species T. granuliferus, which occurs throughout Japan and adjacent areas, including Korea and Russia (Nunomura, 1990; Schmalfuss, 2003). Our aims were to analyze the genetic diversity and genetic structure of T. granuliferus around the Japanese Archipelago, and to discuss the paleogeographic and paleoenvironmental implications of the phylogeny established.
MATERIALS AND METHODS
Sampling
In total, 222 specimens of T. granuliferus, including both adults and juveniles, were collected by hand or sift-sorting from either sandy or gravel beach at 56 sites throughout Japan and at two sites on the Korean Peninsula (Table 1 and Fig. 1). Although we failed to select sampling sites in a random manner due to the rarity of the species, designation of 58 sites as representative of all populations from the Japanese Archipelago and adjacent areas should minimize any substantial error in a gross estimation of phylogenetic relationships.
All the specimens were stored in a freezer at -30°C or in 100% ethanol at room temperature before analysis. As outgroups, we included a Tylos sp. which obviously belongs to an undescribed species from Hachijojima Island, Japan, in the analyses, together with previously published data for T. neozelanicus from New Zealand (EU364624), Ligia oceanica (DQ442914), and Armadillidium vulgare (GU130251).
DNA extraction, PCR, and DNA sequencing
Total DNA was extracted from leg samples from each specimen with the DNeasy Blood and Tissue Kit (Qiagen, CA). A fragment of the mitochondrial gene encoding cytochrome oxidase subunit I (COI) was amplified by PCR using the primers LCO1490 and HCO2198 (Folmer et al., 1994). The cycling program used consisted of an initial denaturation step at 94°C for 5 min, followed by 35 cycles of 94°C for 1 min, 48°C for 1 min, and 72°C for 1 min, before a final extension at 72°C for 5 min.
The amplified products were purified with exonuclease and alkaline phosphatase (ExoProStar; Amersham Biosciences, NJ) and directly sequenced with the primers mentioned above using an ABI BigDye Terminator v3.1 Cycle Sequencing Kit and Applied Biosystems 3730 DNA Analyzer (Applied Biosystems, CA). The 115 new COI haplotypes of T. granuliferus and one haplotype of Tylos sp. were deposited in GenBank (accession numbers: AB763432– AB763451, AB763453–AB763460, AB763462–AB763464, AB763466– AB763530, AB763532–AB763534, AB763536–AB763537, AB763539– AB763553).
Phylogenetic analyses
The partial COI sequences were aligned using MAFFT v.6 (Katoh and Toh, 2008). Before their analysis, the p-distances of each position of codons were plotted against those of the total sequences in order to confirm no evidence of mutation saturations by eye.
The schemes used to partition the data were a single-partition strategy and a three-partition strategy (1st, 2nd vs. 3rd). Each partitioning scheme was assessed using a Bayesian analysis implemented in BEAST 1.7.5 (Drummond et al., 2011) with an uncorrelated lognormal relaxed clock and a Yule tree model prior with 10 million generations, sampling every 1000 iterations. The three-partition strategy of the coding gene by codon position was selected as the optimal strategy for the present data using Bayes factors (Kass and Raftery, 1995) calculated in TRACER 1.5 (Rambaut and Drummond, 2007).
The phylogenetic relationships of T. granuliferus were reconstructed using the maximum likelihood method (ML) with RAxML (Stamatakis, 2006), with the three-partition strategy and the GTR + gamma model for each codon, which was determined to be the best model with the Akaike information criterion using Kakusan4 (Tanabe, 2011). The confidence of branches on the ML tree was assessed with 1,000 bootstrap pseudoreplicates. Tree topologies with 70% or greater bootstrap proportions (BPs) were regarded as having sufficiently resolved nodes, whereas those with BPs between 50% and 70% were considered to be weakly supported (Huelsenbeck and Hillis, 1993; Shaffer et al., 1997).
Bayesian inference (BI) using the Markov chain Monte Carlo (MCMC) technique was also performed with MrBayes 3.2 (Ronquist et al., 2012), using the partitioning schemes and substitution models described above. We initiated four independent analyses with a random starting tree that ran for 10 million generations. We used the program Tracer 1.5 (Rambaut and Drummond, 2007) to determine when the log likelihood of the sampled trees reached a stationary distribution. Because the apparent stationary of the MCMC runs was reached at no later than one million generations, we conservatively discarded the first 25% generations from each run as “burn-in” and sampled one of every 100 generations to calculate the posterior probability for each branch on the Bayesian tree. Bayesian posterior probabilities (BPPs) of 0.95 or greater were considered significant support for a given branch (Larget and Simon, 1999; Huelsenbeck et al., 2001).
Parsimony-based haplotype networks of the COI region for each assemblage were constructed using TCS 1.2 (Clement et al., 2000), with the connection limit set at 95%.
Table 1.
Sampling sites, number of specimens (N), haplotype information, number of polymorphic site (Np), haplotype diversity (h), nucleotide diversity (π), Tajima's D, and Fu's Fs of Tylos granuliferus. Populations with N < 2 were omitted. Asterisk (*) indicates significant difference between samples (P < 0.05). See Figs. 1 and 2 for sampling localities and haplotypes, respectively.
Divergence times, genetic structure, and historical demography
Estimates of the divergence times of T. granuliferus were calculated with BEAST 1.7.5 (Drummond et al., 2011). It has been reported that the evolutionary rates of the COI genes in the Crustacea vary from 1.66%/Ma to 3.0%/Ma (Ketmaier et al., 2003), although a rate of approximately 2.0%/Ma is broadly used for the molecular clock of vertebrate mitochondrial cytochrome b (e.g., Rocha et al., 2005; Fernandes et al., 2012; Rapson et al., 2012). Because there is no reliable temporal calibration point directly applicable to our data, we used a rate of 2.5%/Ma as a representative value, which was reported to be the COI divergence rate in a stenasellid isopod (Ketmaier et al., 2003). Rough estimates of divergence times based on molecular clocks may give some idea of the absolute date of colonization, although their accuracy remains disputable (e.g., Avise, 2000; Heads, 2005).
The analysis was run for 60 million generations, with the first six million generations discarded as “burn-in” and the parameter values were sampled every 1000 generations, imposing a strict clock on the substitution rates (Drummond and Rambaut, 2007) and coalescent constant population size because divergences among closely related lineages follow coalescent process (e.g., Ho et al., 2005; Sota et al., 2005). The parameter estimates and convergences were checked using Tracer 1.5 (Rambaut and Drummond, 2007).
We calculated both haplotype diversity (h) and nucleotide diversity values (π) in each of the samples. Population pairwise fixation indices (FST) were calculated and their significances were tested with a non-parametric permutations approach with 1,000 permutations of haplotypes among sampling localities. To examine whether populations are at equilibrium, we calculated Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997). We chose these tests because of the increased statistical power in detecting significant changes in population size when using small sample sizes (Ramos-Onsins and Rozas, 2002). In addition, demographic parameters were estimated based on mismatch distribution analysis. These analyses were performed using Arlequin 3.5 (Excoffier et al., 2005). In these analyses, we tentatively treated clades and subclades (see below) as populations in the gross estimation of diversity of each lineage.
RESULTS
Sequence divergence
We determined 627 bp of the partial COI sequences without any insertions or deletions for 222 individuals, and detected a total of 115 haplotypes within the species T. granuliferus, of which 232 positions of sequences displayed variability. Interspecific nucleotide replacements of Tylos ranged from 119 bp (T. granuliferus vs. Tylos sp.) to 135 bp (Tylos sp. vs. T. neozelanicus).
The genetic distances (p-distances) among the clades were relatively high (7.7%–16.6%), whereas the intra-clade genetic distances were low (Clade I [C-I]: 2.6%; C-II: 2.6%; C-III: 0.9%; and C-IV: 0.4%). The observed haplotype (h) and nucleotide diversities (π) of C-Ic were higher than those of C-Ia and C-Ib (Table 3).
Fig. 2.
Maximum likelihood (ML) tree derived from mitochondrial COI gene sequences of Tylos granuliferus (InL = -5584). Node support is given as bootstrap proportions (BPs) for the ML analysis, followed by the Bayesian posterior probabilities (BPPs) for BI method. BP and BPP values smaller than 50% and 0.70, respectively, as well as those within clades and subclades, are omitted.
Phylogenetic analysis
The ML tree (Fig. 2) was nearly identical to the BI tree (not shown), and there was no inconsistency between the two analyses in the topology of the ingroup nodes, with significant support (i.e., BP ≥ 70% or BPP ≥ 0.95). The results showed the monophyly of T. granuliferus, with complete support (100% BP and 1.00 BPP; Fig. 2). The specimens of T. granuliferus were divided into four major clades, which corresponded to their geographic distributions: the main island of Japan and the Korean peninsula (C-I; 82%, 1.00), the Ryukyu Archipelago (C-II; 88%, 1.00), Kyushu, Shikoku, and the northeastern of Honshu (C-III; 90%, 1.00), and a narrow area in central eastern Honshu (C-IV; 100%, 1.00). Both analyses indicated the monophyly of an assemblage of C-I and C-II (83%, 0.98) and that of C-III and C-IV (98%, 1.00). Clade-I was further divided into three subclades: samples from the northern part of Kyushu, northwestern Honshu, Hokkaido, and west part of the Korean Peninsula (C-Ia; 79%, 1.00); along the Seto Inland Sea (C-Ib; 97%, 1.00); and from southern Kyushu and east Shikoku (C-Ic; 78%, 1.00). Clade-Ia and C-Ib seemed to constitute a monophyletic group (71%, 0.83), and leaving C-Ic as the most basal lineage within C-I. All of the clades and subclades were distributed allopatrically, except for the samples of C-I and C-III from Shikoku and Kyushu Islands. The distributions of C-I and C-III seemed to overlap in Shikoku Island, and the coexistence of these two clades was observed along the southeastern coast of Kyushu Island (sites 24 and 28 in Fig. 1 and Table 1).
Haplotype network analysis
The haplotype networks and their geographic distributions are shown in Figs. 3 and 4. A star-like network was found for C-Ia, with the most common haplotype H001, which was shared by 32 individuals collected from 10 different sites, and several other haplotypes that were shared between the Pacific coast and the Sea of Japan coast. A haplotype group (H017, H028–029, H031–033, H035–036, and H039–042) in C-Ia was confined to Kyushu and the western part of the main islands of Japan.
All members of C-II occurred in the Ryukyu Archipelago clade, which showed relatively high genetic divergence among specimens. The analysis resolved three separate networks within C-II with no shared haplotypes among the islands, and two independent haplotypes.
Although some individuals within C-III shared common haplotypes, such as H085, H091 and H097, many others did not share. Although the network of C-IV had four haplotypes, the number and size of samples examined are too small to detect any definite patterns.
Estimation of divergence times and demographic population histories
Figure 5 gives divergence times and their 95% credibility intervals (CIs). The initial divergence of T. granuliferus into two major lineages (C-I–II and C-III–IV) was estimated to have occurred during the Tortonian (late) Miocene (9.0 Ma, CI: 7.3–10.9 Ma). The divergence times between C-I and C-II and between C-III and C-IV were estimated to have occurred near from boundary of the Pliocene and Miocene (ca. 5.3 Ma), approximately at 5.0 Ma (CI: 4.0–6.3 Ma) and 4.8 Ma (CI: 3.6–6.1 Ma), respectively, whereas those among the subclades of C-I occurred during Calabrian (early) Pleistocene—Piacenzian (late) Pliocene, C-Ia–Ib and C-Ic approximately at 2.6 Ma (CI: 2.0–3.3 Ma), and C-Ia and C-Ib at 1.7 Ma (CI: 1.2–2.3 Ma). The other intra-clade divergences were assumed to have begun during the Calabrian Pleistocene—Piacenzian Pliocene (C-II: 2.3 Ma, CI: 1.7– 3.0 Ma), Ionian—Calabrian Pleistocene (C-III: 0.9 Ma, CI: 0.6– 1.3 Ma), and Tarantian—Ionian Pleistocene (C-IV: 0.2 Ma, CI: 0.1–0.4 Ma), whereas intra-subclade divergences of C-I occurred in the Ionian—Gelasian Pleistocene (C-Ia, 0.9 Ma, CI: 0.6–1.2 Ma; C-Ib: 0.6 Ma, CI: 0.3–1.0 Ma; and C-Ic, 1.8 Ma, CI: 1.3–2.3 Ma).
The results of Fst showed significant population differentiation (Table 2). Tajima's D and Fu's Fs showed significant negative values in each clades except for C-Ic, C-II and C-IV (Table 3), whereas only a few populations demonstrated significant negative values (Table 1). Mismatch distribution analyses showed all of clades and subclades did not differ from sudden expansion and spatial expansion models significantly.
Fig. 3.
Statistical parsimony networks of COI haplotypes of Tylos granuliferus. Haplotype groups were separated from one another at steps exceeding the 95% confidence connection limit. Attached letters adjacent to the symbols indicate haplotype names, consistent with those used in Table 1. Symbol size is proportional to the sample size.
Fig. 4.
Geographic distributions of COI haplotypes for each population. Symbol size is proportional to the sample size. Patterns for each symbol indicate those used in Fig. 3.
Fig. 5.
Chronogram of Tylos inferred from BI of the COI gene. Numbers above each branch are estimated divergence times. The white boxes on the nodes indicate the 95% of the node age posterior probability distribution. Relationships within clades or subclades and their estimates are omitted. Se: Serravallian, To: Tortonian, Me: Messinian, Za: Zanclean, Pi: Piacenzian, Ge: Gelasian, Ca: Calabrian, Io: Ionian, and Ta: Tarantian.
Table 2.
Pairwise fixation indices (FST) for clades and subclades of Tylos granuliferus. Asterisk (*) indicates significant values (P < 0.05).
Table 3.
Results of haplotype diversity (h), nucleotide diversity (π), Tajima's D and Fu's Fs, and Harpending's Raggedness index (RI) and probablities of RI for sudden expansion and spatial expansion models. Asterisk (*) indicates significant values (P < 0.05).
DISCUSSION
Distribution patterns of populations and their relationships to current activities
Two hypotheses exist regarding the distribution patterns of littoral animals. Some authors have suggested a close relationship between the distribution patterns of populations and marine environmental factors, such as current activity, seawater temperature, and salinity (e.g., Hansen et al., 2007). However, several semiterrestrial coastal species do not show clear relationships between their distributions and local ocean currents (e.g., coastal tiger beetles: Satoh et al., 2004, Ligia spp.: Itani, 2000; Hurtado et al., 2010).
Our results from a substantially larger data set analyzed with methodological improvements show that the distributions of two lineages of T. granuliferus are strongly related to ocean currents, as reported for other coastal species (e.g., anadromous fish Leucopsarion petersii: Kokita and Nohara, 2011). Two well-supported clades, C-I and the assemblage of C-III and C-IV, are distributed along the Tsushima and Kuroshio Currents, respectively. This strongly suggests that the present and past warm ocean currents have played a major role in the transportation and subsequent divergence of this semiterrestrial isopod. These distribution patterns also correspond to marine climatic zones, because the ocean currents determine the marine and coastal environments around the Japanese Archipelago (Nishimura, 1981; Fujikura et al., 2010).
Although our results show that T. granuliferus displays current-related distribution patterns, some other semiterrestrial coastal species display noncurrent-related patterns, as mentioned above. Members of Ligia are adapted to inland environments (Itani, 2000; Taiti et al., 2003; Lee, 1994), which suggests that Ligia has a tolerance for dry conditions and has diversified with inland dispersal. Inland dispersal has also been reported in some aquatic freshwater insects (Kovats et al., 1996). Indeed, species of the genus Tylos could have invaded inland regions and migrated to other beaches via sandy substrates (Brown and Odendaal, 1994). However, there are no large sand-hills in the Japanese Archipelago. Thus, the activities of warm currents along the continents and islands are considered to be a substantial factor in the formation of the distribution pattern of populations of this species. Although T. granuliferus lacks a planktonic stage and inhabits intertidal zone of sandy beaches seemingly avoiding water immersion (Imafuku, 1976), some members seem to have gradually extended their distributional ranges along coasts by ocean currents if powerful waves or rafting substrata were available (e.g., Kensley, 1974; Thiel et al., 2003), and such a transport mechanism is suggested in order to explain the wide geographic distribution of some terrestrial isopods (Taiti et al., 1992) or littoral amphipods (Myers, 1993).
Our results also suggest that the opening of the strait and the direction and bifurcation of the ancient ocean currents affected the divergence of these populations. The distribution of the Ryukyu Archipelago lineage (C-II) is restricted to the area south of the Tokara Strait where the Kuroshio Current branches off from the Tsushima Current and flows towards the Pacific Ocean at present, whereas the sister taxon (C-I) is distributed on the northwestern coasts, north of the strait, including on the main islands of Japan. Although the time of formation of this strait is still disputed (Nakagawa et al., 2001; Yamamoto et al., 2003; Odawara et al., 2005; see Iryu et al. [2006] for further discussion), Kizaki and Oshiro (1977, 1980) calculated that it dates back to no later than the Pliocene (2.6–5.3 Ma). Our gross estimate of the divergence time of C-I and C-II (5.0 Ma, CI: 4.0–6.3 Ma) is congruent with this period. This suggests that the divergence between C-I and C-II was initiated by the opening of the strait and the subsequent bifurcation of the Kuroshio Current.
The Tsushima Strait is located between Kyushu and the Korean Peninsula, and is presumed to have been established in the early (Calabrian) Pleistocene (ca. 1.7 Ma: Kitamura et al., 2001). Besides this establishment, inflows of the Tsushima Current into the Sea of Japan via another strait are suggested based on the presence of warm-water mollusks and planktonic foraminifers in the Pliocene and Gelasian Pleistocene (ca. 3.2, 2.9, 2.4, and 1.9 Ma; Kitamura and Kimoto, 2006). The divergence times among the subclades of C-I (2.6 and 1.7 Ma) are presumed to be consistent with these inflows of the Tsushima Current into the Sea of Japan.
The estimated divergence time between the common ancestor of C-I and C-II and that of C-III and C-IV (9.0 Ma, CI: 7.3–10.9 Ma) also seems corresponds to late (Tortonian) Miocene when a land bridge formed between the continent and the Japanese Archipelago through the Ryukyu Archipelago (Kizaki and Oshiro, 1977). This may have caused dramatic changes in paleogeographic and paleoenvironmental factors. In contrast, the occurrence of C-Ib along the Seto Inland Sea is attributed to complex seawater movements, including the ebb and flow of tides.
Historical biogeography
Our results suggest that the initial divergence within T. granuliferus corresponded to paleogeographic and paleoenvironmental changes during the late (Tortonian) Miocene (9.0 Ma, CI: 7.3–10.9 Ma). Based on our results, together with evidence that the genus Tylos is common within tropical/temperature regions (Kensley, 1974), avoids low-salinity conditions, and is inactive at low temperatures (Imafuku, 1976), the ancestor of T. granuliferus is assumed to have originated in a relatively low-latitude marine environment, somewhere on the subtropical coasts of East Asia, and to have then spread into both the Ryukyu and Japanese Archipelagos during the Miocene when these two archipelagos were connected to continental Asia (Kizaki and Oshiro, 1977). The paleo-East China Sea at this time was not a suitable habitat for marine and coastal species because it was dry land, so the ancestral populations must have been distributed only along the Pacific Ocean coast. The initial divergence of T. granuliferus might have been affected by the geological changes induced by the first phase of the formation of the Okinawa Trough (6–10 Ma; Shinjo, 1999).
At the beginning of the Pliocene (5.3 Ma), changes in the oceanic environments and the fragmentation of the Ryukyu Archipelago (Kizaki and Oshiro, 1977, 1980) could have affected the divergence between C-I and C-II (5.0 Ma, CI: 4.0–6.3 Ma), and also that between C-III and C-IV (4.8 Ma, CI: 3.6–6.1 Ma). The latter is presumed to have occurred allopatrically, with C-III in Kyushu, Shikoku and Honshu, and C-IV in the northeastern part of Honshu (Fig. 1). Interestingly, we observed a unique distribution pattern that has not been reported for other coastal species in Japan (C-IV in Fig. 1). Clade IV is distributed in a narrow restricted area located along the Pacific coast of northern Honshu, where the Kuroshio and Oyashio Currents meet. Pelc et al. (2009) have suggested that the areas in which ocean currents mix become barriers that inhibit gene flow among populations of coastal species. Therefore, it is probable that the integrity of C-IV has been maintained by the collision of two ocean currents. The opening of the strait between the Pacific Ocean and the Sea of Japan and the subsequent bifurcation of the ancient Kuroshio Current (ca. 3.2, 2.9, 2.4, and 1.9 Ma; Kitamura and Kimoto, 2006) and that of the present Tsushima Strait in the Pleistocene (ca. 1.7 Ma) should have enhanced the divergences among subclades of C-I (2.6 and 1.7 Ma). Considering relatively high genetic divergence of C-Ic (Table 3) and the direction of the Tsushima Current (Fig. 1), C-Ic distributed along the southern coast of Kyushu is supposed to be the ancestral lineage of C-I. This view is circumstantially supported by the findings that recent population expansions are expected in C-Ia and C-Ib, because these two subclades showed significant negative values in Tajima's D and Fu's Fs and because they did not differ from sudden expansion model in mismatch distribution analyses significantly. In order to confirm of validity of our scenario, relationships between ocean activities and divergences within each subclade should be seriously tested in future in the light of independent evidence with accumulation of detailed paleogeographic and climatic data.
ACKNOWLEDGMENTS
We thank T. Kakui, S. Shimizu, Y. Uchiyama, T. Nishimura, Y. Yoshida, M. Shimomura, K. Ishii, S. Karasawa, M. Masumoto, K. Miyazaki, T. Suguro, and H. Takasu for providing samples and locality information. We extend special thanks to H. Wada and members of the Graduate School of Life and Environmental Sciences, University of Tsukuba, for their continuous support for our laboratory experiments. Two anonymous reviewers provided valuable comments on the manuscript. This study was partly supported by the Research Institute of Marine Invertebrates Foundation and a Grant-in-Aid for JSPS Fellows.
