Phylogeography of the Rock Shell Thais clavigera (Mollusca): Evidence for Long-Distance Dispersal in the Northwestern Pacific

Abstract

The present-day genetic structure of a species reflects both historical demography and patterns of contemporary gene flow among populations. To precisely understand how these factors shape current population structure of the northwestern (NW) Pacific marine gastropod, Thais clavigera, we determined the partial nucleotide sequences of the mitochondrial COI gene for 602 individuals sampled from 29 localities spanning almost the whole distribution of T. clavigera in the NW Pacific Ocean (~3,700 km). Results from population genetic and demographic analyses (AMOVA, ΦST-statistics, haplotype networks, Tajima’s D, Fu’s FS, mismatch distribution, and Bayesian skyline plots) revealed a lack of genealogical branches or geographical clusters, and a high level of genetic (haplotype) diversity within each of studied population. Nevertheless, low but significant genetic structuring was detected among some geographical populations separated by the Changjiang River, suggesting the presence of geographical barriers to larval dispersal around this region. Several lines of evidence including significant negative Tajima’s D and Fu’s FS statistics values, the unimodally shaped mismatch distribution, and Bayesian skyline plots suggest a population expansion at marine isotope stage 11 (MIS 11; 400 ka), the longest and warmest interglacial interval during the Pleistocene epoch. The lack of genetic structure among the great majority of the NW Pacific T. clavigera populations may be attributable to high gene flow by current-driven long-distance dispersal of prolonged planktonic larval phase of this species.

About the Authors:

Xiang Guo

Contributed equally to this work with: Xiang Guo, Dan Zhao

Affiliation: Division of EcoScience, Ewha Womans University, 52 Ewhayeodae-gil, Seodaemun-gu, Seoul, 120–750, Republic of Korea

Dan Zhao

Contributed equally to this work with: Xiang Guo, Dan Zhao

Affiliation: The Key Laboratory of Mariculture, Ministry of Education, Ocean University of China, Qingdao, 266003, China

Daewui Jung

Affiliation: Program in Cell Biology and Genetics and Department of Parasitology, College of Medicine, Chungbuk National University, Cheongju, 361–763, Republic of Korea

Qi Li

* E-mail: [email protected] (J-KP); [email protected] (QL)

Affiliation: The Key Laboratory of Mariculture, Ministry of Education, Ocean University of China, Qingdao, 266003, China

Ling-Feng Kong

Affiliation: The Key Laboratory of Mariculture, Ministry of Education, Ocean University of China, Qingdao, 266003, China

Gang Ni

Affiliation: The Key Laboratory of Mariculture, Ministry of Education, Ocean University of China, Qingdao, 266003, China

Tomoyuki Nakano

Affiliation: Seto Marine Biological Laboratory, Kyoto University, 459 Nishimuro, Shirahama, Wakayama Prefecture, Japan

Akihiko Matsukuma

Affiliation: The Key Laboratory of Mariculture, Ministry of Education, Ocean University of China, Qingdao, 266003, China

Sanghee Kim

Affiliation: Korea Polar Research Institute, 26 Songdomirae-ro Yeonsu-gu, Incheon, 406–840, Republic of Korea

Chungoo Park

Affiliation: School of Biological Sciences and Technology, Chonnam National University, Gwangju, 500–757, Republic of Korea

Hyuk Je Lee

Affiliation: Department of Biological Science, College of Science and Engineering, Sangji University, Wonju, 220–702, Republic of Korea

Joong-Ki Park

* E-mail: [email protected] (J-KP); [email protected] (QL)

Affiliation: Division of EcoScience, Ewha Womans University, 52 Ewhayeodae-gil, Seodaemun-gu, Seoul, 120–750, Republic of Korea

Introduction

The present-day genetic structure of a species reflects both historical demography (i.e., population history) and patterns of contemporary gene flow among populations [1]. An understanding of precisely how these factors shape current population structure has been one of the main issues in ecology and evolutionary biology. Historical demographic events might have a particularly large influence on the genetic structure if a species has not yet reached mutation-drift equilibrium owing to insufficient evolutionary time since a population expansion, which is often the case for marine organisms [2].

A series of marginal seas separates Asia from the Pacific, straddling the world’s largest subduction zone of the Western Pacific [3]. During major glaciations, the large volume of accumulated ice reduced sea levels by approximately 120–140 m [4]. The East China Sea (ECS), the South China Sea (SCS), and the East Sea (ES) (Sea of Japan, SJ) are thought to have become partly continuous landmasses during the glacial period due to the sea level decline in the northwestern (NW) Pacific [5]. Marine organisms in the NW Pacific became extinct or survive only in glacial refuges [6]; however, rising sea levels in the interglacial period resulted in reconnection of these three marginal seas and range extension in marine organisms. Therefore, historical glaciation is one of the important factors that shaped present-day phylogeographical patterns of marine species. Significant genetic differentiation by historical isolation of sea basins in the NW Pacific has been well documented for a variety of marine animal species including crustaceans, molluscs, and fishes [7]. Each marginal sea has been proposed as an independent refuge, resulting in deep phylogeographic divergence among lineages in some marine taxa, such as the red lip mullet Chelon haematocheilus [8], the venus clam Cyclina sinensis [9], and the mitten crab Eriocheir sensu stricto [10]. Compared with records of historical population fragmentation, information on how oceanographic factors influence the genetic structure and distribution of natural populations remains scarce in this region.

Ocean currents and utilization of a wide range habitat are also correlated with long-distance dispersal, and some marine invertebrates with no geographical structure over large geographic scales have been reported [11–13]. The NW Pacific Sea is characterized by intricate hydrology and geography; two influential ocean current systems (Kuroshio Current [KSC] and China Sea Coastal Current [CSCC]) and freshwater outflows from a continental body (Changjiang River) may play an important role in limiting the dispersal of the planktonic larvae of marine species. Chinese coastlines are influenced by the CSCC, including two currents: the cold Subei Coastal Current (SBCC; flowing southward) and the warm China Coastal Current (CCC; flowing northeastward) (Fig 1). Coastal areas of Japan, western Korea, and eastern Taiwan are influenced by the warm KSC, which includes a strong main stream and three branches, namely, the Tsushima Warm Current (TSWC), Yellow Sea Warm Current (YSWC), and Taiwan Warm Current (TWC). The KSC carries salty water with a relatively high temperature; however, the SBCC carries low salinity surface water along the ECS coast line and the CCC flows into the ECS from the SCS in the summer. Evidence for the influence of different surface water circulation systems on the species distribution and genetic structure of marine organisms has been observed in the acorn barnacle Tetraclita japonica in the NW Pacific [14]. It was also proposed that freshwater outflow from the Changjiang River acted as an extrinsic barrier to the southward dispersal of larvae (Cellana toreuma) from the north in the NW Pacific [5]. More studies are needed to determine whether the effects of the Changjiang River on genetic structure are species-specific or common in marine species in this area. In general, the three marginal seas (ECS, SCS, and ES [SJ]) in the NW Pacific provide a natural setting to investigate the influences of population fluctuation during Pleistocene oscillations and contemporary gene flow on population structure of marine species.

[Figure omitted. See PDF.]

Fig 1. Map of East Asia showing the sampling sites of Thais clavigera and the summer ocean currents redrawn from [9].

Populations are labelled with numbers that correspond with those shown in Table 1. Shaded sea areas indicate regions 120 m in depth that would have been exposed during periods of low sea level. EKWC, East Korea Warm Current; TSWC, Tsushima Warm Current; KSC, Kuroshio Current; YSWC, Yellow Sea Warm Current; SBCC, Subei Coastal Current; CCC, China Coastal Current; TWC, Taiwan Warm Current. Dashed lines (——) represent a border of water bodies among South China Sea, East China Sea and Pacific Ocean.

https://doi.org/10.1371/journal.pone.0129715.g001

Most marine invertebrates remain in a pelagic larval stage for several days to months, and they may be able to disperse and metamorphose into sedentary adults [15]. It is generally thought that the propagule duration of marine species is significantly correlated with dispersal capacity [16]. Therefore, marine species with longer planktonic larval stages are expected to exhibit relatively low levels of population structure as a result of their increased opportunity for gene flow [11]. However, some surveys have shown that high levels of gene flow among populations are not always coupled with the duration of the pelagic larval stage. Small, but significant genetic structure has been observed in some marine invertebrates with high dispersal capability, including limpets [17], pen shells [18], and lobsters [19].

The rock shell Thais clavigera is the most common gastropod species in intertidal rocky shores of East Asia including China, Korea, Japan, and Taiwan. The life span of this species is a minimum of 7 years [20] and it is capable of adapting to eurythermic and euryhaline environments [21]. Although there is no precise estimate of the pelagic larval duration for T. clavigera, larvae of its congeneric species T. haemastoma [22] and T. chocolata [23] are known to last more than 60 days. Morphologically, T. clavigera exhibits a wide range of variation in shell sculpture, such as shape, the absence/presence of blotches on the nodules of shell surfaces, and shell apertures [24,25] In a previous DNA barcoding analysis of Korean Thais species, a high degree of genetic divergence was discovered in some T. clavigera populations [26]. In addition, some genetic surveys of T. clavigera populations performed in certain NW Pacific localities have resulted in inconsistent conclusions [27–29]: Huang [29] revealed significant population structure along the China coast based on the partial sequence of mt COI; however, based on sequence variation in mt16S and COI gene fragments, Wang [27] concluded that there is no genetic differentiation between the ESC and SCS. Moreover, taxon sampling was limited to a few local areas in these studies (Taiwan [28], China Sea [27,29]), insufficient to understand genetic relationships among NW Pacific populations.

In the present study, to better understand contemporary genetic structure and historical demography of the NW Pacific T. clavigera populations, we sequenced a partial fragment of the mt COI gene from a total of 602 individuals sampled at 29 localities across the NW Pacific Ocean, including both its upstream (Taiwan and China mainland) and downstream (Korea and Japan) ranges. Specifically, we tested whether potential glacial refuges, ocean current systems, or freshwater outflows of the Changjiang River significantly influenced the current population genetic structure of T. clavigera in the NW Pacific.

Materials and Methods

Sample Collection and Sequencing

A total of 602 T. clavigera specimens were sampled from 29 localities across the NW Pacific Ocean spanning a distance of approximately 3,700 km in East Asia from September 2008 to June 2014 (Fig 1; Table 1). At least 20 individuals were collected and genetically analysed from each locality except two Japanese populations (Iwate and Kanagawa), from which only 10 and 15 samples were obtained, respectively. T. clavigera is not an endangered or protected species, and therefore all collections were made from public access area without specific permits.

[Figure omitted. See PDF.]

Table 1. Sampling information, geographic coordinates, diversity indices, and neutrality tests for 29 Thais clavigera populations based on COI.

https://doi.org/10.1371/journal.pone.0129715.t001

Total genomic DNA was extracted using the E.Z.N.A. Mollusc DNA Kit (Omega Bio-Tek Inc., Norcross, GA, USA) following the instructions supplied by the manufacturer. A universal primer set (LCO1490: 5′-GGTCAACAAATCATAAAGATATTGG-3′, HCO2198: 5′-TAAACTTCAGGGTGACCAAAAAATCA-3′) [30] was used to amplify the partial fragment of the mitochondrial COI gene. Polymerase chain reaction (PCR) was performed in a 50 μL reaction volume containing ten units of Ex-Taq Polymerase (Takara, Shiga, Japan), 2.5 mM dNTP mixture, 2.5 mM MgCl2 and 20 pmole of each primer with the following amplification conditions: an initial 30 s denaturation at 94°C, 40 cycles of 10 s at 98°C, 30 s at 47°C, 30 s at 72°C and a final 10 min extension at 72°C. The sequencing reaction was performed using the BigDye Terminator Cycle Sequencing Kit (Applied Biosystems, Carlsbad, CA, USA) and all COI fragments were sequenced in two directions on an ABI 3730 XL (Applied Biosystems, USA) automatic sequencer.

Molecular Diversity

Forward and reverse sequences of the COI target fragment were edited, assembled, and merged into consensus sequences using the Geneious software program [31]. Consensus sequences were aligned with ClustalX 1.81 using default settings [32]. Standard molecular diversity indices including haplotype diversity (h), nucleotide diversity (π), and the mean number of pairwise differences (k) were calculated using the Arlequin 3.5 software package [33]. The number of haplotypes was determined by a Bayesian coalescent-based program implemented in DnaSP V5 [34]. All haplotype sequences were deposited in GenBank with the accession numbers KP116312–KP116913.

Phylogenetic Analyses and Genealogical Network Construction

To infer the phylogenetic relationships among haplotypes, Bayesian inference (BI) and neighbour-joining (NJ) analyses were performed with two congeneric species (T. bronni and T. luteostoma) as outgroups using MrBayes 3.1.2 [35] and Mega 6.0 [36], respectively. The Markov-chain Monte Carlo (MCMC) search was run with four chains for 10 million generations with a sampling frequency of 1/1,000 trees. The best-fit model of nucleotide substitutions for Bayesian analyses was estimated using jModeltest v.0.1.1 [37], and selected based on the Akaike Information Criterion ([38]. HKY+G+I was identified as the most appropriate model and used for subsequent analyses. A network showing the genetic relationships among haplotypes was constructed using a median-joining algorithm in Network 4.612 [39]. The maximum parsimony (MP) option was run on the output file to delete all superfluous median vectors and links [40].

Population Genetic Structure

A hierarchical analysis of molecular variance (AMOVA) was conducted using Arlequin to assess the population genetic structure of T. clavigera in the NW Pacific. Three independent AMOVA analyses were carried out based on our a priori expectation for each of the focal factors (potential glacial refuges, Changjiang River, and ocean current systems). Our groupings (see Table 2 for details) were as follows (Table 1 & Fig 1 also show sampling sites and location codes, respectively): (1) for glacial refuges, PC (Pacific Ocean) included populations 1 and 3–11, ES (SJ) included 2, 14, and 15, ECS included 12–26, and SCS included 27–29; (2) for Changjiang River, CJN (Changjiang River North) included 18–22 and CJS (Changjiang River South) included 23–24 and 26–29; (3) for ocean systems, KSC included 1–17 and 25, and CSCC included 18–24 and 26–29. These AMOVA analyses partitioned the total molecular variance among groups (ΦCT), among populations within groups (ΦSC), and among populations whatever the groups (ΦST) and tested if those Φ-statistic values were statistically significant using 10,000 random permutations. The HKY model, which fit the data best according to jModeltest, cannot be implemented in Arlequin, so the Tamura-Nei model was selected to correct for multiple substitutions.

[Figure omitted. See PDF.]

Table 2. Analysis of molecular variance (AMOVA) results of population structure.

Significant P-values are indicated in bold.

https://doi.org/10.1371/journal.pone.0129715.t002

Pairwise genetic differentiation between all 29 populations was further assessed with Φ-statistics in Arlequin 3.5 [33]. The significance of each pairwise comparison was estimated using 10,000 permutations, and corrections for multiple tests were made following a sequential Bonferroni procedure [41]. To estimate isolation by distance (IBD), a nonparametric Mantel test was performed online to evaluate the association between matrices of pairwise comparisons among sampling locations of shortest geographical distances (estimated in km using Google Earth version 4.3) and genetic distances (log-transformed) with 10,000 randomizations of the data using Isolation By Distance Web Service 3.23 [42].

Demographic Analyses

Historical demographic analysis was conducted using two different methods. Tajima’s D [43] and Fu’s Fs [44] statistics were calculated to test for neutrality using DnaSP. The significance levels of Tajima’s D and Fu’s Fs were evaluated under 10,000 permutations. The parameter τ obtained from the mismatch distribution was used to estimate the time elapsed since the sudden population expansion using the equation τ = 2μt, where μ is the mutation rate of the marker (per locus per generation) and t is the number of generations (Rogers & Harpending, 1992). In order to convert parameters into quantitative estimates of time, we used a mutation rate of 7.9 × 10−9 substitutions/site/year, the genus Nucella-specific average substitution rate for COI [45], because Nucella is ever known the most closely relatives of T. clavigera (belongs to the same family Muricidae) in which fossil-based molecular clock calibration is thus far available. The generation time was assumed to be 1 year.

For comparison with the mismatch distribution analysis, a Bayesian Skyline Plot (BSP) analysis was executed to examine changes in population size across time in Beast v1.7 software [46]. Three independent MCMC samplings were performed to assure the consistency of the results. Chains were run for 100 million generations and sampled every 1,000 generations, with the first 10% of generations discarded as “burn-in” under the HKY+G+I model determined by jModeltest v.0.1.1 with a constant skyline model and Bayesian skyline tree priors. All operators were optimized automatically. In all runs, the effective sample size yielded by MCMC chains for the parameters of interest was greater than 300. Finally, the results were visualized with the Tracer 1.6 software program [47].

Results

Genetic Diversity and Phylogenetic Relationships

A 658 bp sequence of the mt COI gene fragment was determined for 602 individuals of 29 T. clavigera populations sampled across the NW Pacific Ocean (Table 1; Fig 1). A total of 165 polymorphic sites were identified and 298 haplotypes were encountered. No deletion or insertion mutations were detected. Of these 165 variable sites, 44 were single-nucleotide polymorphic, and the remains were parsimony informative. The COI data in NW Pacific T. clavigera populations showed a high level of genetic diversity; haplotype diversity (h) for all studied populations was 0.982, ranging from 0.944 to 1.000. Nucleotide diversity (π) for all examined samples was 0.0081 and ranged from 0.0064 to 0.0102 for individual populations. The average number of nucleotide differences (k) between haplotypes was 5.34 (Table 1).

The phylogenetic trees for T. clavigera mtDNA COI sequence data inferred using NJ and BI methods were both unresolved, i.e., all nodes received low bootstrap support (all nodes ≤ 50%), suggesting that the 298 haplotypes sampled in all study populations lacked phylogeographic structure (S1 Fig). For three different groupings based on either potential glacial refuges, the Changjiang River, or the present-day ocean current systems, the centrality of the haplotype networks was occupied by the two most dominant haplotypes; however, these accounted for less than 20% of all sampled individuals (Fig 2). Furthermore, regardless of how we divided these populations into groups, there were no significant genealogical branches or geographic clusters detected in the three haplotype networks. For potential refuges, the two most dominant haplotypes accounted for 16.8% (101/602) of the total individuals (Fig 2A). The most dominant haplotype (9.5%, 57/602) occurred in all populations except Chungcheongnam-do (CN), Pohang-si (PH), Zhoushan (ZS), and Haikou (HK). We detected the second most dominant haplotype (7.3%, 44/602) in 22 out of 29 populations. The great majority of detected haplotypes (76%, 227/298 haplotypes) occurred in only a single individual (i.e., singletons).

[Figure omitted. See PDF.]

Fig 2. Network of Thais clavigera using COI data.

The size of circles is proportional to haplotype frequency and median vectors are indicated with a red pie chart. The different shadings of the haplotypes refer to the respective locations in terms of the three population groupings: (a) potential refuges (PC, Pacific Ocean; ECS, East China Sea; SCS, South China Sea; ES, East Sea [SJ, Sea of Japan]); (b) ocean currents (KSC, Kruoshio Current with the branch currents; CCC, China Coastal Current and Subei Coastal Current); (c) Changjiang River (CJS, the South of Changjiang River; CJN, the North of Changjiang River).

https://doi.org/10.1371/journal.pone.0129715.g002

Population Genetic Structure

The three-level hierarchical AMOVA analyses revealed small but statistically significant genetic structure among populations of T. clavigera, irrespective of whether populations are grouped with respect to potential refuges, Changjiang River, or ocean currents (Table 2). We found that 1.47–5.14% of molecular variation occurred among groups (ΦCT; P < 0.01), whereas most variation was observed within populations (ΦST; 93.84–96.92% of the total variation, P < 0.001). Among the three factors analysed, however, the Changjiang River had the largest influence on genetic structure, indicating it acts as a physical barrier to larval dispersal; the proportion of variation among groups was approximately five times greater than the proportion of variation among populations within groups (ΦCT = 0.051, ΦSC = 0.011; Table 2).

Pairwise comparisons across sampling locations showed a near absence of population structure across a wide geographical range (3,700 km). Pairwise ΦST values between populations varied from -0.027 to 0.155, most of which were non-significant after sequential Bonferroni correction (P > 0.05). Interestingly, however, pairwise ΦST values from eight populations exclusively north or south of the Changjiang River (Suzaki [SU; 4 (location code)], Pohang-si [PH; 14], Qingdao [QD; 19], Rizhao [RZ; 20], Zhoushan [ZS; 23], Wenzhou [WZ; 24], Shenzhen [SZ; 27] and Haikou [HK; 29]) were small, but statistically highly significant (P < 0.05) (Table 3). These results substantiate the existence of a barrier to gene flow across the Changjiang River; however, there was no evidence of IBD estimated as a correlation between the shortest geographic distance by sea and genetic distance by the Mantle test (r = 0.46, P = 1.00) (Fig 3).

[Figure omitted. See PDF.]

Fig 3. Isolation by distance plot for all Thais clavigera populations from the northwestern Pacific Ocean.

https://doi.org/10.1371/journal.pone.0129715.g003

[Figure omitted. See PDF.]

Table 3. Pairwise ΦST among Thais clavigera populations based on COI (below diagonal) and associated P-values (above diagonal) after the sequential Bonferroni correction.

https://doi.org/10.1371/journal.pone.0129715.t003

Demographic History

We detected significant population expansion of T. clavigera based on multiple lines of evidence from neutrality tests and the mismatch distribution (Table 1). Tajima’s D rejected neutrality (P < 0.05) for 23 out of 29 populations. Similarly, Fu’s Fs statistic was significantly different from zero (P < 0.05) for all populations except for the Taiwan samples. Both tests of neutrality were significant (i.e., indicated significant deviations from neutrality) when samples were grouped as a single data set. The mismatch distribution was unimodal for all examined samples, supporting a model of sudden expansion (Fig 4). The expansion time parameter τ was estimated from mismatch distribution analysis to be 4.005. Assuming a substitution rate of 7.9 × 10−9 for COI data, the time since expansion was estimated to be approximately 385 kyr. The population expansion was further validated by the results of BSP analysis, which revealed that population sizes began to expand approximately 400 ka (Fig 5). These estimates indicate a population expansion of T. clavigera during marine isotope stage (MIS) 11, the longest and warmest interglacial interval of the past 500 kyr.

[Figure omitted. See PDF.]

Fig 4. Mismatch distribution of Thais clavigera based on COI.

The dotted line with circles represents the observed distribution, whereas the solid line shows the expected value under the sudden demographic expansion model.

https://doi.org/10.1371/journal.pone.0129715.g004

[Figure omitted. See PDF.]

Fig 5. Demographic history of Thais clavigera estimated using Bayesian skyline plots from COI data.

The black line represents the median population estimates, while the blue lines are the upper and lower bounds of the 95% confidence interval.

https://doi.org/10.1371/journal.pone.0129715.g005

Discussion

The genetic structure of a species reflects the effect of historical demography as well as contemporary gene flow among populations. The factors that contribute to the present-day genetic structure can be inferred by genetic analysis of molecular data. In this study, we utilized sequence variation of the partial mtDNA COI gene for 602 individuals sampled from 29 localities spanning an extremely long coastline of approximately 3,700 km. This sampling represents almost the whole distribution of T. clavigera in the NW Pacific Ocean. We found that genetic variation among the NW Pacific populations was generally low, perhaps owing to a combination of high contemporary gene flow and recent common ancestry of haplotypes. Nevertheless, both AMOVA and pairwise ΦST analyses indicated weak, but significant genetic structure across the Changjiang River (Tables 2 & 3), suggesting the presence of geographical barriers to continuous larval dispersal at this locality; however, it should be noted that no signal of IBD was detected, nor did the haplotype networks show distinct genealogical branches or geographic clusters (Figs 2 & 3). These results support the hypothesis that the populations of T. clavigera that were examined have a high level of gene flow throughout the NW Pacific Ocean.

Although there are some exceptions [e.g., the marine clam genus Lasaea [48]], marine invertebrate species with a long-lived planktotrophic larval stage are generally capable of long-distance dispersal, and their offspring are spread several hundred to thousand kilometres away from their origin by the prevailing surface flow of the ocean current [16,49,50]; hence, the long-distance dispersal associated with a prolonged pelagic larval stage and the present-day oceanic current have been regarded as the most influential factors contributing to continuous gene flow over a wide geographic scale in many marine invertebrates, including molluscan species [11–13,18]. A reproductive ontogenetic study has reported that T. clavigera undergoes indirect development with a planktonic veliger larval stage [20], and its pelagic larval stage lasts up to approximately 2 months. Additionally, T. clavigera is perennial with an average lifespan of 7 years. The species can reach sexual maturity during their second year, and accordingly have a reproductive lifespan of approximately 6 years [20]. It should also be noted that in the NW Pacific, there are two influential current systems in surface water circulation (Fig 1), namely, the KSC (Kuroshio Current) and the CSCC (China Sea Coastal Current), with the KSC flowing northward year-round [51,52], and the CCC (China Coastal Current) entering the ECS (East China Sea) from the SCS (South China Sea) in the summer [7]. These prevailing currents transport a great number of warm-water marine species from their tropical centre to the north and expand their ranges [53]. The long planktonic larval stage in T. clavigera may facilitates gene flow by current-driven dispersal of pelagic larvae and consequently decreases genetic structure among distant populations spanning over 3,700 km in the NW Pacific coastline. Aside from long-distance dispersal ability and ocean currents, we hypothesize that its utilization of a wide range of habitats is a key factor for successful colonization in a new environments. T. clavigera is abundant in the intertidal zone over a wide range of environmental conditions including different temperatures and salinities [21]; thus it ability to inhabit a wide range of eurythermic and euryhaline environments may also indicate its potential to colonize new environments.

In this study, the AMOVA analysis and pairwise ΦST values (Tables 2 & 3) suggest that the Changjiang River poses a weak but significant barrier to gene flow among some T. clavigera populations, indicating that the larval pool is not well mixed geographically across this area despite the long planktonic larval stage. Nonetheless, it is evident that there is generally a low level of genetic structure among populations, and we did not detect geographic clusters in the haplotype networks, consistent with previous results in marine invertebrates with long-lived pelagic larvae [18,19]. In addition to potential glacial refuges in the NW Pacific region, the Changjiang River and ocean circulation systems are two potentially important geographical barriers shaping current population structure [5,14]. We observed genetic divergence (ΦCT = 0.051, P = 0.003) between the northern and southern populations of the Changjiang River in our AMOVA analysis (Table 2). As documented earlier, freshwater outflows from the Changjiang River may act as physical barriers that limit northward dispersal of planktotrophic larvae from southern populations [5]. Genetic subdivision may also be attributed to the habitat of T. clavigera near the mouth of the Changjiang River. T. clavigera is most commonly found in shadowy crevices in intertidal rocky shores. There are relatively well developed mudflat areas formed by the deposition of sediments near the mouth of the Changjiang River. These conditions provide a relatively insufficient rocky shore substratum, and are consequentially unsuitable for the settlement of T. clavigera larvae. When specimens were sampled near the northern mouth of the Changjiang River (e.g., Nantong), only a few T. clavigera individuals were found; by contrast, near the southern mouth, T. clavigera were very abundant in the rocky seashore of the Zhoushan archipelago. This habitat discontinuity may have reduced effective gene flow between the northern and southern populations of the Changjiang River.

It has been reported that oceanographic patterns play an important role in maintaining genetic and phenotypic differentiation in the acorn barnacle Tetraclita japonica in the NW Pacific [14]. Moreover, in southern Australia, the major ocean currents influence the phylogeography and population structure of the intertidal barnacle Catomerus polymerus [54]; however, in the present study of T. clavigera, we found a lack of genetic structure across major ocean current systems, we found very low, but statistically significant genetic structure between the two major ocean circulation systems (CSCC and KSC) by AMOVA (ΦCT = 0.015, P = 0.002), but ΦSC (i.e., structure among populations within groups) was equivalent and statistically significant (ΦSC = 0.018, P < 0.001). The lack of genetic structure between the two current systems may be attributable to high gene flow owing to the long spawning time and prolonged planktonic larval phase of T. clavigera. Our findings that the most common two haplotypes occur at almost every site support this high gene flow hypothesis. Some additional observations suggest that although no water mass from a sub-branch of the KSC reaches the SCS in the summer, in some years this does occur [55]. Moreover, in other seasons, a south-westward current from Kuroshio flowing into the SCS has been observed [55]. In southern China and Taiwan, T. clavigera spawning could occur from spring to summer (February to August) [20,28]; furthermore, its pelagic larval duration lasts up to approximately 2 months. Therefore, T. clavigera larvae are likely to enter and mix into the CSCC system from the Taiwan coastline, which may increase gene flow between populations in the two circulation systems to some extent, therefore resulting in a lack of genetic differentiation.

From haplotype network analysis for COI data, we found a complicated network pattern that suggested that T. clavigera populations underwent a demographic expansion (Fig 2). Also, the observed pattern of mtDNA variation in T. clavigera further supports the hypothesis of non-equilibrium historical processes such as population range expansion. We observed very high COI haplotype diversity due to an excess of singleton variants (76% of the 298 detected haplotypes) coupled with relatively low nucleotide diversity. The retention of a surplus of rare COI variants may indicate a recent population expansion of T. clavigera; otherwise, these rare variants are predicted to be eliminated by genetic drift [56]. This suggests that mutation-drift equilibrium has not yet been attained in T. clavigera in the NW Pacific [57], an interpretation consistent with the significantly negative neutrality test statistics, a clear unimodal mismatch distribution, and BSP analysis. Furthermore, both mismatch distribution and BSP revealed an MIS 11 population expansion approximately 400 ka (Figs 4 & 5). This stage is the longest and the warmest during the Pleistocene epoch [12] and has been described as a super-interglacial period because of its long duration of 25–30 kyr [58]. The sequence of land mollusc species fossils in the Chinese loess-soil shows that the summer monsoon was particularly strengthened during MIS 11, which is typical of warmer climates [59]. Paleontological and palaeoecological estimates of MIS 11 deposits from Japan, Hawaii, Bermuda, and the Bahamas suggest a global sea level had risen during this stage [60–62]. Such climatic conditions are necessary to allow warm-water species to reach the northern Pacific and expand their range.

The Pleistocene glacial age, and particularly the last glacial maximum (LGM) approximately 20,000 years ago, had an important influence on the evolution and genetic structure of marine organisms. Many species in various marine realms appeared to arise at the beginning of LGM [2,63]; however, population expansion of T. clavigera is assumed to have occurred pre-LGM. These results differ from the traditional view of demographic expansion that it occurs during the period of LGM. Nevertheless, population expansion that occurred pre-LGM has been reported in species such as the cold-water barnacle Chthamalus challengeri [13] and the marine snail Concholepas concholepas [12]. In the NW Pacific, Ni et al. [7] also estimated a period of 120–140 ka that corresponded with dramatic population expansion of various species, including molluscs [13,64], fishes [65,66], and crustaceans [67,68]. These earlier reports indicate the importance of pre-LGM events in determining the demography of marine populations and should be considered in the future. Nonetheless, it is still unclear why the glacial events did not significantly impact the current population structure and demographic history of these species.

Conclusions

To better understand contemporary genetic structure and historical demography of the NW Pacific T. clavigera populations, we determined the partial sequence of the mt COI gene from 602 individuals sampled from 29 localities across the NW Pacific Ocean. We observed a high level of genetic diversity within each of sampled populations, and no significant genealogical branches or geographic clusters, suggesting high levels of gene flow among populations throughout the NW Pacific. Nevertheless, we detected low, but significant genetic differentiation that corresponds to habitat conditions and freshwater discharge from the Changjiang River. Since we used only a single mtDNA marker, further studies of T. clavigera using multiple nuclear markers are required to validate the observed genetic structure. Also, population genetic studies of other marine species with a planktonic larval phase in the region will provide additional insight into the phylogeographic patterns of NW Pacific organisms.

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. Neighbour-joining (NJ) tree of COI haplotypes for populations of Thais clavigera from the northwestern Pacific.

Both NJ and Bayesian inference (BI) analyses yield the same topology. Bootstrap values for NJ (the former number) and the posterior probabilities for BI (the latter number) analyses are indicated at the nodes.

https://doi.org/10.1371/journal.pone.0129715.s001

(TIF)

Acknowledgments

The authors thank Sang-Hwa Lee (Chungbuk National University), Hideki Takami, Morihiko Tomatsuri, Koetsu Kon, Shinji Iseda, Keiju Okano, Hikaru Hamaguchi, Takeo Kurihara, Takenori Sasaki, Shinichi Ikebe, Ryoji Takada, Yoshitake Takada, Tamotsu Hamada (Japan), and Kwen-Shen Lee (National Museum of Natural Science, Taiwan) for providing Thais specimens.

Author Contributions

Conceived and designed the experiments: JKP QL. Performed the experiments: XG DZ. Analyzed the data: XG DZ. Contributed reagents/materials/analysis tools: JKP QL. Wrote the paper: XG JKP. Interpreted the data and edited the manuscript: HJL. Collected the samples: JKP DZ DJ SK CP LFK TN GN AM.

Citation: Guo X, Zhao D, Jung D, Li Q, Kong L-F, Ni G, et al. (2015) Phylogeography of the Rock Shell Thais clavigera (Mollusca): Evidence for Long-Distance Dispersal in the Northwestern Pacific. PLoS ONE 10(7): e0129715. https://doi.org/10.1371/journal.pone.0129715

References

1. Avise JC (2004) Molecular markers, natural history and evolution. 2nd edition. Geographic population strucutre and gene flow. Massachusetts: Sinauer Associates, Inc. Publishers. pp. 248–282.

2. Marko PB, Hoffman JM, Emme SA, McGovern TM, Keever CC, Cox LN, et al. (2010) The ‘Expansion–Contraction’model of Pleistocene biogeography: rocky shores suffer a sea change? Molecular Ecology 19: 146–169. pmid:20092033

3. Wang PX (1999) Response of Western Pacific marginal seas to glacial cycles: paleoceanographic and sedimentological features. Marine Geology 156: 5–39.

4. Lambeck K, Esat TM, Potter E-K (2002) Links between climate and sea levels for the past three million years. Nature 419: 199–206. pmid:12226674

5. Dong Y-W, Wang H-S, Han G-D, Ke C-H, Zhan X, Nakano T, et al. (2012) The impact of Yangtze River discharge, ocean currents and historical events on the biogeographic pattern of Cellana toreuma along the China coast. PLoS One 7: e36178. pmid:22563446

6. Hewitt GM (1996) Some genetic consequences of ice ages, and their role in divergence and speciation. Biological Journal of the Linnean Society 58: 247–276.

7. Ni G, Li Q, Kong L, Yu H (2014) Comparative phylogeography in marginal seas of the northwestern Pacific. Molecular Ecology 23: 534–548. pmid:24600706

8. Liu JX, Gao TX, Wu SF, Zhang YP (2007) Pleistocene isolation in the Northwestern Pacific marginal seas and limited dispersal in a marine fish, Chelon haematocheilus (Temminck & Schlegel, 1845). Molecular Ecology 16: 275–288. pmid:17217344

9. Ni G, Li Q, Kong L, Zheng X (2012) Phylogeography of bivalve Cyclina sinensis: testing the historical glaciations and Changjiang River outflow hypotheses in northwestern Pacific. PLoS One 7: e49487. pmid:23145178

10. Xu J, Chan T-Y, Tsang LM, Chu KH (2009) Phylogeography of the mitten crab Eriocheir sensu stricto in East Asia: Pleistocene isolation, population expansion and secondary contact. Molecular phylogenetics and Evolution 52: 45–56. pmid:19236929

11. Addison JA, Ort BS, Mesa KA, Pogson GH (2008) Range-wide genetic homogeneity in the California sea mussel (Mytilus californianus): a comparison of allozymes, nuclear DNA markers, and mitochondrial DNA sequences. Molecular Ecology 17: 4222–4232. pmid:19378402

12. Cárdenas L, Castilla JC, Viard F (2009) A phylogeographical analysis across three biogeographical provinces of the south-eastern Pacific: the case of the marine gastropod Concholepas concholepas. Journal of Biogeography 36: 969–981.

13. Cheang CC, Tsang LM, Ng WC, Williams GA, Chu KH, Chan BKK (2012) Phylogeography of the cold-water barnacle Chthamalus challengeri in the north-western Pacific: effect of past population expansion and contemporary gene flow. Journal of Biogeography 39: 1819–1835.

14. Tsang LM, Chan BK, Ma KY, Chu KH (2008) Genetic differentiation, hybridization and adaptive divergence in two subspecies of the acorn barnacle Tetraclita japonica in the northwestern Pacific. Molecular Ecology 17: 4151–4163. pmid:19238711

15. Purcell JF, Cowen RK, Hughes CR, Williams DA (2006) Weak genetic structure indicates strong dispersal limits: a tale of two coral reef fish. Proceedings of the royal society B: biological sciences 273: 1483–1490. pmid:16777742

16. Shanks AL (2009) Pelagic larval duration and dispersal distance revisited. The biological bulletin 216: 373–385. pmid:19556601

17. Hoffman JI, Clarke A, Clark MS, Fretwell P, Peck LS (2012) Unexpected fine-scale population structure in a broadcast-spawning Antarctic marine mollusc. PLoS One 7: e32415. pmid:22403655

18. Xue D-X, Wang H-Y, Zhang T, Liu J-X (2014) Population Genetic Structure and Demographic History of Atrina pectinata Based on Mitochondrial DNA and Microsatellite Markers. PLoS One 9: e95436. pmid:24789175

19. Stamatis C, Triantafyllidis A, Moutou K, Mamuris Z (2004) Mitochondrial DNA variation in Northeast Atlantic and Mediterranean populations of Norway lobster, Nephrops norvegicus. Molecular Ecology 13: 1377–1390. pmid:15140084

20. Tong LK (1988) The reproductive biology of Thais clavigera and Morula musiva (Gastropoda: Muricidae) in Hong Kong. Asian Marine Biology 5: 65–75.

21. You Z-J, Chen Z-Y (2010) Systematic taxonomy of Thais (Gastropoda: Muricidae) along Zhejiang Coast. Journal of Zhejiang Ocean University (Natural Science) 29: 306–317.

22. Scheltema RS (1971) Larval dispersal as a means of genetic exchange between geographically separated populations of shallow-water benthic marine gastropods. The biological bulletin 140: 284–322.

23. Soledad Romero M, Gallardo CS, Bellolio G (2004) Egg laying and embryonic-larval development in the snail Thais (Stramonita) chocolata (Duclos, 1832) with observations on its evolutionary relationships within the Muricidae. Marine Biology 145: 681–692.

24. Hayashi T (1999) Genetic differentiation between the two forms of Thais clavigera (Küster, 1858)(Mollusca, Gastropoda) in Tanabe Bay, central Japan. Zoological science 16: 81–86.

25. Tan K, Liu L (2001) Description of a new species of Thais (Mollusca: Neogastropoda: Muricidae) from Taiwan, based on morphological and allozyme analyses. Zoological science 18: 1275–1289. pmid:11911084

26. Lee S-H, Kim T, Lee J, Lee J-R, Park J-K (2011) Mitochondrial DNA sequence variations and genetic relationships among Korean Thais species (Muricidae: Gastropoda). Animal Systematics, Evolution and Diversity 27: 1–17.

27. Wang X (2013) Research on genetic diversity of Thais luteostoma and Thais clavigera: Master's thesis: Zhejiang Ocean University.

28. Hsieh Y-C (2001) Genetic structure of populations of oyster drill (Thais clavigera) along the west coast of Taiwan. Master's thesis: National Sun Yat-sen University.

29. Huang J-R (2013) Study on the taxonomy of Thais and population genetic structure of Thais clavigera in China coast. Master's thesis: Xiamen University.

30. Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R (1994) DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Molecular Marine Biology and Biotechnology 3: 294–299. pmid:7881515

31. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. (2012) Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28: 1647–1649. pmid:22543367

32. Thompson JD, Gibson TJ, Plewniak F, Jeanmougin F, Higgins DG (1997) The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic acids research 25: 4876–4882. pmid:9396791

33. Excoffier L, Lischer HE (2010) Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Molecular ecology resources 10: 564–567. pmid:21565059

34. Librado P, Rozas J (2009) DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452. pmid:19346325

35. Huelsenbeck JP, Ronquist F (2001) MRBAYES: Bayesian inference of phylogenetic trees. Bioinformatics 17: 754–755. pmid:11524383

36. Tamura K, Stecher G, Peterson D, Filipski A, Kumar S (2013) MEGA6: molecular evolutionary genetics analysis version 6.0. Molecular biology and evolution 30: 2725–2729. pmid:24132122

37. Posada D (2008) jModelTest: phylogenetic model averaging. Molecular biology and evolution 25: 1253–1256. pmid:18397919

38. Posada D, Buckley TR (2004) Model selection and model averaging in phylogenetics: advantages of Akaike information criterion and Bayesian approaches over likelihood ratio tests. Systematic biology 53: 793–808. pmid:15545256

39. Bandelt H-J, Forster P, Röhl A (1999) Median-joining networks for inferring intraspecific phylogenies. Molecular biology and evolution 16: 37–48. pmid:10331250

40. Polzin T, Daneshmand SV (2003) On Steiner trees and minimum spanning trees in hypergraphs. Operations Research Letters 31: 12–20.

41. Rice WR (1989) Analyzing tables of statistical tests. Evolution: 223–225.

42. Jensen JL, Bohonak AJ, Kelley ST (2005) Isolation by distance, web service. BMC genetics 6: 13. pmid:15760479

43. Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123: 585–595. pmid:2513255

44. Fu Y-X (1997) Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147: 915–925. pmid:9335623

45. Nicole Cox L, Zaslavskaya NI, Marko PB, Rocha L (2014) Phylogeography and trans-Pacific divergence of the rocky shore gastropod Nucella lima. Journal of Biogeography 41: 615–627.

46. Drummond AJ, Rambaut A (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC evolutionary biology 7: 214. pmid:17996036

47. Rambaut A, Suchard M, Drummond A (2013) Tracer v1.6. Available: http://tree.bio.ed.ac.uk/software/tracer. Accessed 11 December 2013.

48. O'Foighil D (1989) Planktotrophic larval development is associated with a restricted geographic range in Lasaea, a genus of brooding, hermaphroditic bivalves. Marine Biology 103: 349–358.

49. Kinlan BP, Gaines SD (2003) Propagule dispersal in marine and terrestrial environments: a community perspective. Ecology 84: 2007–2020.

50. Nunes FL, Norris RD, Knowlton N (2011) Long distance dispersal and connectivity in amphi-Atlantic corals at regional and basin scales. PLoS One 6: e22298. pmid:21799816

51. Ito T, Kaneko A, Furukawa H, Gohda N, Koterayama W (1995) A structure of the Kuroshio and its related upwelling on the East China Sea shelf slope. Journal of Oceanography 51: 267–278.

52. Lee H-J, Chao S-Y (2003) A climatological description of circulation in and around the East China Sea. Deep Sea Research Part II: Topical Studies in Oceanography 50: 1065–1084.

53. Liu JY (2013) Status of marine biodiversity of the China seas. PLoS One 8: e50719. pmid:23320065

54. York KL, Blacket MJ, Appleton BR (2008) The Bassian Isthmus and the major ocean currents of southeast Australia influence the phylogeography and population structure of a southern Australian intertidal barnacle Catomerus polymerus (Darwin). Molecular Ecology 17: 1948–1961. pmid:18363669

55. Hu J, Kawamura H, Hong H, Qi Y (2000) A review on the currents in the South China Sea: seasonal circulation, South China Sea warm current and Kuroshio intrusion. Journal of Oceanography 56: 607–624.

56. Aris-Brosou S, Excoffier L (1996) The impact of population expansion and mutation rate heterogeneity on DNA sequence polymorphism. Molecular Biology and Evolution 13: 494–504. pmid:8742638

57. Zane L, Marcato S, Bargelloni L, Bortolotto E, Papetti C, Simonato M, et al. (2006) Demographic history and population structure of the Antarctic silverfish Pleuragramma antarcticum. Molecular Ecology 15: 4499–4511. pmid:17107479

58. Candy I, Schreve DC, Sherriff J, Tye GJ (2014) Marine Isotope Stage 11: Palaeoclimates, palaeoenvironments and its role as an analogue for the current interglacial. Earth-Science Reviews 128: 18–51.

59. Guo Z, Biscaye P, Wei L, Chen X, Peng S, Liu T (2000) Summer monsoon variations over the last 1.2 Ma from the weathering of loess-soil sequences in China. Geophysical Research Letters 27: 1751–1754.

60. Kariya C, Hyodo M, Tanigawa K, Sato H (2010) Sea-level variation during MIS 11 constrained by stepwise Osaka Bay extensions and its relation with climatic evolution. Quaternary Science Reviews 29: 1863–1879.

61. Hearty P, Kindler P, Cheng H, Edwards R (1999) A+ 20 m middle Pleistocene sea-level highstand (Bermuda and the Bahamas) due to partial collapse of Antarctic ice. Geology 27: 375–378.

62. Hearty PJ (2002) The Ka'ena highstand of O'ahu, Hawai'i: further evidence of Antarctic ice collapse during the middle Pleistocene. Pacific Science 56: 65–81.

63. Maggs CA, Castilho R, Foltz D, Henzler C, Jolly MT, Kelly J, et al. (2008) Evaluating signatures of glacial refugia for North Atlantic benthic marine taxa. Ecology 89: S108–S122. pmid:19097488

64. Kang J, Park J, Choi T (2012) Genetic differentiation of octopuses from different habitats near the Korean Peninsula and eastern China based on analysis of the mDNA cytochrome C oxidase 1 gene. Genetics and Molecular Research 11: 3988–3997. pmid:23212336

65. Xiao Y, Zhang Y, Gao T, Yanagimoto T, Yabe M, Sakurai Y (2009) Genetic diversity in the mtDNA control region and population structure in the small yellow croaker Larimichthys polyactis. Environmental biology of fishes 85: 303–314.

66. Song N, Zhang XM, Sun XF, Yanagimoto T, Gao TX (2010) Population genetic structure and larval dispersal potential of spottedtail goby Synechogobius ommaturus in the north-west Pacific. J Fish Biol 77: 388–402. pmid:20646163

67. Liu Y, Liu R, Ye L, Liang J, Xuan F, Xu Q (2009) Genetic differentiation between populations of swimming crab Portunus trituberculatus along the coastal waters of the East China Sea. Hydrobiologia 618: 125–137.

68. Guo E, Liu Y, Cui Z, Li X, Cheng Y, Wu X (2012) Genetic variation and population structure of swimming crab (Portunus trituberculatus) inferred from mitochondrial control region. Molecular biology reports 39: 1453–1463. pmid:21625862