Pleistocene events and present environmental factors have shaped the phylogeography of the intertidal limpet Cellana toreuma (Reeve, 1855) (Gastropoda: Nacellidae) in Southeast Asia and China

Abstract

High levels of biodiversity and phylogeographic structure in marine species in Southeast Asia are strongly linked to Quaternary sea-level fluctuations and complex oceanographic conditions. Cellana toreuma is a common limpet on intertidal rocky shores and is widely distributed in the Western Pacific. Analyses of partial mitochondrial COI gene sequences from Southeast Asia, combined with previously published sequences from East Asia and Indonesia, revealed the existence of five well-supported clades with high genetic divergences (between 1.4 and 7.6%), namely the East Asia clade, the eastern Southeast Asia clade, the western Southeast Asia clade, the Pelabuan Ratu (Java) clade and the Ogasawara clade. The geographical distribution of the five clades is likely related to the history of glaciations and rapid postglacial population expansions. Analyses of pairwise Φ_ST and hierarchical analysis of molecular variance shows significant population structure among collections in East and Southeast Asia. These results suggest that historical events have had strong effects on the phylogeographic structure of C. toreuma. In addition, present environmental factors, such as unsuitable habitats and ocean currents, have also affected the genetic footprints of past environments.

INTRODUCTION

It is well known that the Indo-Pacific region is a biodiversity hotspot with high marine biodiversity and deep intraspecific divergences (Briggs, 1999; Hughes, Bellwood & Connolly, 2002; Roberts et al., 2002; Carpenter & Springer, 2005; Tittensor et al., 2010). Demonstration of the roles of historical and contemporary restrictions on gene flow in creating such high biodiversity and intraspecific divergence is a challenge, given the intricate geographical history and highly dynamic hydrology of the region (Briggs, 1999). However, biogeographical and phylogeographical studies can promote the understanding of current patterns of both population divergence and speciation in marine ecosystems (Bouchet, 2006), especially as increased information becomes available on geomorphology and physical oceanography (DeBoer et al., 2014).

Areas and configurations of marginal seas in the Indo-Pacific region changed dramatically during the glacial-interglacial cycles of the Pleistocene (Wang, 1999; Voris, 2000). Low sea levels during the glacial intervals resulted in exposure of shallow continental shelves (e.g. Sunda Shelf) and connection of adjacent islands by land bridges (Voris, 2000). The lowest sea level (−130 m) occurred during the Last Glacial Maximum (LGM) (Voris, 2000; Clark et al., 2009; Fig. 1). Patterns of allopatric distribution in the region suggest that low sea level resulting in separation between Pacific and Indian Oceans, and formation of isolated marine basins (such as South China Sea (SCS), Sulu Sea, Celebes Sea and Flores Sea), was a driver of diversification (Briggs, 1999; Barber, Erdmann & Palumbi, 2006; Timm & Kochzius, 2008; Chan et al., 2014). Phylogeographic patterns of barnacles in the Indo-Pacific showed that some populations had been isolated in the LGM, but panmictic populations were also observed (Tsang et al., 2011; Cheang et al., 2012). Maggs et al. (2008) proposed four models of population structure (indicated by haplotype networks and allele/haplotype frequency distributions) following glaciation, by which to infer glacial refugia and understand postglacial colonization processes. In Model I (a null model), no difference in haplotype frequency is observed among populations; these are panmictic with high levels of contemporary gene exchange and originate from a single refugium. In Model II, populations show a latitudinal cline in allelic richness; these populations originate from a single refugium and have limited contemporary gene flow. In Model III, the distributions of monophyletic groups are closely correlated with multiple refugia and there is little or no gene flow between them. In Model IV, divergent populations from multiple refugia have become mixed, and there is an increase in diversity in the zone of secondary contact. In addition, significant effects of events before the LGM on population structure of marine species have been revealed in some regions, including the northeastern Pacific (Marko et al., 2010), northwestern Pacific (Ni et al., 2014) and in northern European seas (Laakkonen et al., 2013).

Figure 1.

Map of East Asia and Southeast Asia showing sampling sites and the coastal currents in summer (after Fang et al., 1998; Qiao & Lu, 2008). Black squares represent sampling localities of the present study; open squares show localities from previously published studies (for abbreviations see Table 1). Shaded sea areas are continental shelves that would have been exposed during periods of low sea level (–130 m). Abbeviations: 1, China Coastal Current; 2, Yellow Sea Warm Current; 3, Tsushima Warm Current; 4, Kuroshio Current; 5, Taiwan Warm Current; 6, SCS Warm Current; 7, SCS Branch of Kuroshio or Donsha Current; 8, NW Luzon Coastal Current; 9, SCS Southern Anticyclonic Gyre; 10, circulation in Gulf of Thailand; 11, Southeast Vietnam Offshore Current; SJ, Sea of Japan; YS, Yellow Sea; ECS, East China Sea; SCS, South China Sea.

Contemporary ocean currents are known to play important roles in shaping the present-day genetic structuring of coastal species (Tsang et al., 2012; Wee et al., 2014). Ocean circulation patterns in Indo-Pacific region are particularly complex and change seasonally (Wyrtki, 1961; Hu et al., 2000). For example, the surface circulation of the SCS, mainly driven by monsoon winds, is cyclonic (southward) in winter and anticyclonic (northward) in summer (Fang et al., 1998; Fig. 1). The surface circulation in the southern SCS does not always flow northward in summer along the Indo-China Peninsula as the northward Southeast Vietnam Offshore Current collides with the southward western branch of the cyclonic eddy to its north near the Indo-China Peninsula (Cai, Long & Wang, 2007; Fig. 1). Some studies, however, have suggested that phylogeographic patterns do not necessarily follow the major ocean currents and that ecological considerations can also be significant (Wyrtki, 1961; Barber et al., 2000; Reid et al., 2006). Along the Chinese coastline, there is a genetic break within certain widespread species at the Yangtze River estuary, due to the combined effects of the freshwater discharge of the Yangtze River, coastal currents and substrates (Dong et al., 2012; see also Yu et al., 2014). Clearly a variety of present-day and historical mechanisms can lead to the population subdivision and genetic divergence of marine species in this biodiversity hotspot.

Cellana toreuma (Reeve, 1855) is a common gastropod limpet on intertidal rocky shores with a recorded geographical distribution from Japan to Philippines (42—15°N; Powell, 1973; Nakano & Ozawa, 2007; Nakano & Espinosa, 2010). Though the pelagic larval duration of C. toreuma is still uncertain, it may last between 8 and 13 d as recorded in some congeners (Bird et al., 2007). This species has previously been included in several molecular phylogenetic studies (see Table 1) and distinct clades have been revealed on the coast of China and Indonesia (Nakano & Espinosa, 2010; Dong et al., 2012). Our aim was to extend this investigation of C. toreuma, to reveal its broad-scale phylogeographical pattern across East and Southeast Asia, and to investigate the role of historical and environmental factors in the formation of this pattern. To do so, we have generated a database of new and previously published mitochondrial COI sequences. This was not designed as a taxonomic study but, in view of the considerable taxonomic confusion surrounding Cellana species in the tropical Indo-Pacific (e.g. Powell, 1973; Nakano & Ozawa, 2007; Nakano et al., 2009; Nakano & Espinosa, 2010), we also tested the monophyly of the nominal species C. toreuma.

Table 1.

Sampling localities for Cellana toreuma complex, diversity indices and neutrality tests from present study and published studies.

Locality	Abbrev.	Coordinates	n	h (mean ± SD)	π (mean ± SD)	Tajima's D	Fu's FS	Reference
Qingdao, China	QD	36°09′N, 120°36′E	21	0.452 ± 0.105	0.00081 ± 0.00082	−0.33029	−0.27034	Dong et al. (2012)
Dagongdao, China	DGD	36°02′N, 120°13′E	30	0.515 ± 0.089	0.00097 ± 0.00090	−1.11781	−2.13140	Dong et al. (2012)
Shengshan, China	SS	30°07′N, 121°55′E	15	0.133 ± 0.112	0.00023 ± 0.00040	−1.15945	−0.64899	Dong et al. (2012)
Zhoushan, China	ZS	29°57′N, 122°12′E	52	0.184 ± 0.072	0.00033 ± 0.00046	−1.76473	−4.92636	Dong et al. (2012)
Nanjiliedao, China	NJ	27°33′N, 120°41′E	30	0.193 ± 0.095	0.00045 ± 0.00057	−1.88948	−2.67393	Dong et al. (2012)
Pingtan, China	PT	25°27′N, 119°48′E	29	0.136 ± 0.085	0.00024 ± 0.00039	−1.50906	−2.31203	Dong et al. (2012)
Nanridao, China	NRD	25°14′N, 119°27′E	32	0.182 ± 0.090	0.00032 ± 0.00047	−1.72954	−3.48995	Dong et al. (2012)
Meizhoudao, China	MZD	25°14′N, 119°09′E	24	0.163 ± 0.099	0.00028 ± 0.00044	−1.51469	−2.07839	Dong et al. (2012)
Chongwu, China	CW	24°53′N, 118°56′E	32	0.182 ± 0.090	0.00032 ± 0.00047	−1.72954	−3.48995	Dong et al. (2012)
Xiamen, China	XM	24°25′N, 118°08′E	32	0.238 ± 0.099	0.00043 ± 0.00055	−1.88876	−4.54522	Dong et al. (2012)
Dongshan, China	DS	23°43′N, 117°32′E	32	0.063 ± 0.058	0.00011 ± 0.00026	−1.14244	−1.26483	Dong et al. (2012)
Nanao, China	NA	23°26′N, 117°00′E	30	0.193 ± 0.095	0.00034 ± 0.00048	−1.73178	−3.38072	Dong et al. (2012)
Shenzhen, China	SZ	22°32′N, 114°13′E	14	0.264 ± 0.136	0.00045 ± 0.00059	−0.34144	0.18574	Dong et al. (2012)
Hong Kong, China	HK	22°17′N, 114°10′E	31	0.125 ± 0.077	0.00021 ± 0.00037	−0.77374	−0.46749	Dong et al. (2012)
Weizhoudao, China	WZD	21°26′N, 109°03′E	14	0.143 ± 0.119	0.00024 ± 0.00042	−1.15524	−0.59478	Dong et al. (2012)
Sanya, China	SY	18°09′N, 109°34′E	29	0.623 ± 0.086	0.02122 ± 0.01097	1.66698	12.79311	Present study
Da Nang, Vietnam	DN	16°08′N, 108°19′E	23	0.391 ± 0.125	0.01866 ± 0.00982	−0.55842	9.38326	Present study
Ko Sichang, Thailand	KS	12°40′N, 101°28′E	30	0.416 ± 0.112	0.00079 ± 0.00079	−1.98300	−5.92770	Present study
Muara, Brunei	MR	04°56′N, 115°02′E	24	0.159 ± 0.095	0.00027 ± 0.00043	−0.68111	−0.24889	Present study
Kuching, Sarawak, Malaysia	KC	01°45′N, 110°19′E	29	0.000 ± 0.000	0.00000 ± 0.00000	0.00000	0.00000	Present study
Lingao, China	LG	19°59′N, 109°38′E	2	–	–	–	–	Lin et al. (2015)
Oga, Akita Pref., Japan	OA	39°53′N, 139°50′E	1	–	–	–	–	Nakano & Ozawa (2007)
Pelabuan Ratu, Java, Indonesia	PR	06°58′S, 106°32′E	1	–	–	–	–	Nakano & Ozawa (2007)
Yamada-cho, Iwate Pref., Japan	YIP	39°28′N, 141°57′E	1	–	–	–	–	Nakano & Espinosa (2010)
Akashouzaki, Fukui Pref., Japan	AFP	35°32′N, 135°41′E	1	–	–	–	–	Nakano & Espinosa (2010)
Morozaki, Aichi Pref., Japan	MAP	34°41′N, 136°58′E	1	–	–	–	–	Nakano & Espinosa (2010)
Yukinoura, Mie Pref., Japan	YMP	34°42′N, 136°30′E	1	–	–	–	–	Nakano & Espinosa (2010)
Hayama, Kanagawa Pref., Japan	HKP	35°16′N, 139°16′E	1	–	–	–	–	Nakano & Espinosa (2010)
Ayukawa, Miyagi Pref., Japan	AMP	38°16′N, 141°31′E	1	–	–	–	–	Nakano & Espinosa (2010)
Tappizaki, Aomori Pref., Japan	TAP	41°15′N, 140°20′E	1	–	–	–	–	Nakano & Espinosa (2010)
Shoudoshima, Kagawa Pref., Japan	SKP	34°25′N, 134°15′E	1	–	–	–	–	Nakano & Espinosa (2010)
Mihama-cho, Wakayama Pref., Japan	MWP	33°49′N, 136°03′E	1	–	–	–	–	Nakano & Espinosa (2010)
Reihoku, Kumamoto Pref., Japan	RKP	32°30′N, 130°03′E	1	–	–	–	–	Nakano & Espinosa (2010)
Ping Lang Qiaq, Taiwan	PLQ	22°47′N, 121°11′E	1	–	–	–	–	Nakano & Espinosa (2010)
Shitiping, Taiwan	STP	23°28′N, 121°30′E	1	–	–	–	–	Nakano & Espinosa (2010)
Yeliu, Taiwan	YL	25°13′N, 121°42′E	1	–	–	–	–	Nakano & Espinosa (2010)
Ningo, Ghana	NG	05°44′N, 000°10′E	2	–	–	–	–	Nakano & Espinosa (2010)
Ogasawara Is, Japan	OG	27°04′N, 142°21′E	13	0.628 ± 0.143	0.00125 ± 0.00111	−1.43714	−2.53494	T. Nakano (unpubl.); Bird et al. (2011); Nakano et al. (2009)

Abbreviations: n, number of individuals; h, haplotype diversity; π, nucleotide diversity; bold font indicates P < 0.05.

MATERIAL AND METHODS

Sampling and sequencing

From January 2013 to January 2014, 135 Cellana toreuma specimens were collected from five rocky shore localities (23–30 individuals per locality) across Southeast Asia (Fig. 1, Table 1). Samples were preserved in absolute ethanol before DNA extraction. Total genomic DNA was extracted from foot-muscle tissue (1 mg) using the TIANamp Marine Animals DNA Kit (Tiangen Biotech Co., Beijing).

A partial sequence of the COI gene was amplified from all samples with universal primers LCO1490 and HCO2198 (Folmer et al., 1994). PCR reactions were conducted in a 25-μl reaction volume containing 2.5 μl of 10× buffer (Mg²⁺ Plus), 2 μl of 2.5 mM dNTPs, 1 μl of each 10 mM primer, 0.25 μl (1.25 U) of Taq DNA polymerase and 200 ng DNA template. Amplification was initiated with denaturing at 95 °C for 3 min, followed by 35 cycles of 95 °C for 1 min, annealing at 40 °C for 1 min and 72 °C for 1 min, and then a final extension at 72 °C for 10 min. The target amplicon was visualized in 1.5% agarose gels and PCR products were sent to a commercial company for sequencing (Invitrogen Biotechnology Co., Shanghai).

Phylogenetic analyses

The COI haplotype dataset used for phylogenetic analysis was composed of 135 sequences (431 bp in length) of C. toreuma from the present study (GenBank accession nos KR132806–KR132940) and other Cellana sequences available in GenBank (). Published sequences from 25 other nominal species of Cellana were included in the analysis, in order to test the monophyly of C. toreuma. Trees were rooted with a species of Nacella as outgroup. Prior to maximum likelihood (ML) and Bayesian analyses, the best-fitting substitution model was evaluated using jModelTest v. 2.1.1 (Darriba et al., 2012) under the Bayesian information criterion (HKY + G, G = 0.122). Bayesian analyses were performed using MrBayes v. 3.2 (Ronquist, Huelsenbeck & Teslenko, 2011). Markov-chain Monte Carlo (MCMC) searches were run with four chains for 6,000,000 generations, sampling every 100 generations. Convergence was evaluated by the standard deviation of split frequencies (<0.01) and potential scale-reduction factor (close to 1.0 for all parameters). The total 60,000 sample trees were used to estimate the consensus tree and the posterior probabilities (PP) of node support after discarding the first 25% of sampled trees. The ML tree was constructed using MEGA6 (Tamura et al., 2013) with the corresponding substitution model selected by jModelTest. Statistical support for each node was obtained from 10,000 bootstrap (BS) replications.

Population structure analyses

For analysing the phylogeographic patterns of C. toreuma, a population dataset was constructed with 572 sequences of COI (length 595 bp) from C. toreuma, including the 135 new sequences and 437 previously published sequences from China (Dong et al., 2012; Lin, Kong & Li, 2015), Japan, Taiwan, Indonesia and Ghana (Nakano & Ozawa, 2007; Nakano & Espinosa, 2010; Fig. 1, Table 1). In addition, based on the results of the phylogenetic analysis, 13 published sequences identified as ‘C. radiata enneagona’ from the Ogasawara Islands, southeastern Japan, were included (Nakano et al., 2009; Bird et al., 2011; T. Nakano, unpubl.).

All COI sequences were assembled and edited by comparing both DNA strands using DNAMAN v. 7.0 software (LynnonBioSoft, Quebec) and aligned with MUSCLE (Edgar, 2004) implemented in MEGA6 (Tamura et al., 2013) with default settings. Standard molecular diversity indices, including haplotype diversity (h), nucleotide diversity (π), Fu's F_s (Fu, 1997) and Tajima's D (Tajima, 1989) were calculated for each locality using Arlequin v. 3.5 (Excoffier & Lischer, 2010).

A network of COI haplotypes illustrating genealogical relationships between haplotypes was constructed with a 95% maximum parsimony threshold using TCS v. 1.21 (Clement, Posada & Crandall, 2000). The uncorrected pairwise genetic distances (p distances) between clades were calculated in the Distances module in MEGA6 (Tamura et al., 2013) to quantify genetic divergence between clades.

Analyses of genetic differentiation between populations were conducted using 21 localities (Table 1) for which sufficiently large samples (>10 individuals) were available. Pairwise Φ_ST estimates among localities were calculated to evaluate the level of genetic differentiation and the significance was estimated with 10,000 permutations using Arlequin v. 3.5 (Excoffier & Lischer, 2010) followed by Bonferroni correction for the 210 pairwise population comparisons (Rice, 1989).

All localities (Table 1), even those with only one individual, were included for the genetic structure analyses. A simulated annealing approach implemented in SAMOVA v. 2.0 (Dupanloup, Schneider & Excoffier, 2002) was used for a varying number of groups (K) to define groups of populations that were phylogeographically homogeneous and maximally differentiated from each other. The analysis was run for 2–12 groups (5 simulations repeated for each K), with 1,000 permutations. Then the grouping associated with the highest F_CT (i.e. the best number of groups and the best population configuration) was selected. Finally a hierarchical analysis of molecular variance (AMOVA) (Excoffier, Smouse & Quattro, 1992) was conducted for the grouping defined by the SAMOAV using Arlequin v. 3.5 (Excoffier & Lischer, 2010).

Demographic history

Owing to the absence of a clear fossil record or datable vicariant event, the a divergence rate for COI of 0.85–1.15% per myr (calculated from cowries and used for estimating divergence time in the limpet genus Eoacmaea by Kirkendale & Meyer, 2004, and for C. nigrolineata by Nakano, Sasaki & Kase, 2010) was applied to estimate divergence time of C. toreuma. A substitution rate of 0.5 × 10⁻⁸ per site per year was deduced from this divergence rate.

The change of effective population size of each clade over time was estimated using Bayesian skyline plot (BSP) in BEAST v. 1.7.4 and BEAUti (Drummond et al., 2012). Analyses were conducted under the corresponding model suggested by jModelTest with constant Bayesian skyline tree priors under a strict clock model. A substitution rate of 0.5% per myr was used to convert estimates into units of time in the present study with a generation time of 1 year for C. toreuma (Wang & Wu, 1999). Default priors were used for parameter settings. Three independent Metropolis-coupled MCMC searches were run for 50,000,000 generations, with the first 5,000,000 generations discarded as burn-in. For each clade, results from replicate runs were pooled with LogCombiner v. 1.7 (Drummond et al., 2012) for the final reconstruction of the BSP. Tajima's D (Tajima, 1989), Fu's F_s (Fu, 1997) and mismatch distribution were also calculated for each clade with Arlequin v. 3.5 (Excoffier & Lischer, 2010) to validate the result from the BSP. Population expansion parameters (τ) were transformed to estimates of real time since expansion (t) using the formula τ = 2μkt, where μ is the mutation rate and k is the sequence length.

RESULTS

Phylogenetic analyses

ML and Bayesian analyses revealed a monophyletic C. toreuma (although not all individuals had previously been identified as belonging to this species, see below), composed of five clades with high support: the East Asia (EA) clade (BS 98, PP 1.00), the eastern Southeast Asia (ESA) clade (86, 0.98), the western Southeast Asia (WSA) clade (91, 1.00), the Pelabuan Ratu (PR) clade (95, 1.00) and the Ogasawara (OG) clade (95, 0.96) (Fig. 2). The EA clade was dominated by individuals from China, Japan and Taiwan, including two individuals from the north side of Hainan Island, China (LG) and 22 from the south side of Hainan (SY), in the SCS. Other individuals collected in the SCS (others from SY and three from DN in Vietnam) had haplotypes that were placed in the ESA clade. All other individuals from the SCS (DN, KS and KC) were combined in the WSA clade. These three clades (EA, ESA, WSA) therefore showed partial geographical overlap. The PR clade comprised one individual from southwestern Java (PR) and two from Ghana, West Africa (NG; where the species was presumably introduced; Nakano & Espinosa, 2010). The last of the five clades, here named OG, contained individuals from the Ogasawara Islands, previously identified as an endemic subspecies, C. radiata enneagona (Nakano et al., 2009). Two further individuals from Madagascar, also previously identified as C. radiata enneagona, grouped with other Cellana species elsewhere in the phylogeny (Fig. 2).

Figure 2.

Phylogenetic analysis of Cellana species based on COI sequences of C. toreuma complex from present study and other species available from GenBank (Table 1; ). A. Maximum likelihood tree with bootstrap values (>50). Sample sizes (>1) shown in the brackets. B. Bayesian tree with posterior probabilities (>0.5). Filled symbols indicate samples from present study and open symbols those from previously published studies. Abbreviations: EA, East Asia; ESA, eastern Southeast Asia; OG, Ogasawara; PR, Pelabuan Ratu; WSA, western Southeast Asia.

Shells of specimens of C. toreuma from the East China Sea and SCS are shown in and specimens from Japan, Indonesia and Ghana were illustrated by Nakano & Espinosa (2010). Shell characters, such as thickness, surface sculpture and shell colour are widely variable in this group, so that no consistent differences could be observed between the five clades defined by the phylogenetic analysis.

Population genetic analyses

A total of 59 haplotypes were defined among the 585 sequences in the population dataset (135 new sequences, 450 previously published; Table 1, ). Low levels of haplotype diversity were found in all five main locations of the present study (mean h ± SD: 0.000 ± 0.000 to 0.623 ± 0.086; Table 1), which was similar to the results of Dong et al. (2012) (0.063 ± 0.058 to 0.515 ± 0.089; Table 1). Regarding nucleotide diversity, two populations (SY and DN) were higher than the others by at least one order of magnitude (mean π ± SD: SY = 0.02122 ± 0.01097, DN = 0.01866 ± 0.00982, others: 0.00000 ± 0.00000 to 0.00125 ± 0.00111; Table 1).

The five clades of C. toreuma revealed in the phylogenetic analyses were also observed in the COI haplotype network, deeply separated from each other since at least 29–40 bp substitutions were required to connect them, except for the 8 bp substitution connecting the OG and ESA clades (Fig. 3). The p distances among the five clades ranged from 1.4% between OG and ESA to 7.6% between PR and EA + WSA (Table 2). These values were much larger than the genetic distances within clades (0.3–0.4%, Table 2).

Table 2.

Uncorrected pairwise genetic distances (p distances) for COI among clades of Cellana toreuma complex.

	EA clade	ESA clade	WSA clade	PR clade	OG clade
EA clade	0.004	0.047	0.075	0.076	0.046
ESA clade	0.051	0.004	0.067	0.068	0.014
WSA clade	0.078	0.071	0.003	0.076	0.066
PR clade	0.079	0.072	0.080	0.003	0.074
OG clade	0.049	0.018	0.069	0.077	0.003

Mean within-clade divergences on diagonal. The lower left portion shows average between-clade divergences and the upper right the net between-clade divergences.

Figure 3.

TCS network of COI haplotypes for Cellana toreuma complex. Each circle represents a single haplotype and sizes are proportional to haplotype frequencies. The open square represents the C. radiata enneagona haplotype. The small white circles and numbers indicate mutational steps between haplotypes. Colours within the circles denote localities (see Table 1 and Fig. 2 for abbreviations).

Genetic diversity analyses showed that the lowest levels of h and π occurred in the EA clade (h: 0.2247 ± 0.0272; π: 0.000431 ± 0.000533; Table 3). The other four clades had similar levels of h and π (h: 0.5775–0.6667; π: 0.001173–0.002241; Table 3). The most common haplotype H1 (65.30% of all individuals) was from the EA clade and was widely distributed in China, Japan and Taiwan (Fig. 3; ). Two private haplotypes in the EA clade with relatively high frequency were H18 (2.91% of individuals) in SY and H30 (2.05% of individuals) in QD and DGD (Fig. 3; ). H14 in MR was the highest-frequency haplotype in the ESA clade and accounted for 3.76% of the individuals (Fig. 3; ). H2 in DN and KS (7.01% of individuals) and H7 in KC (4.96% of individuals) were the two highest-frequency haplotypes in the WSA clade (Fig. 3; ). H55 in OG (1.37% of individuals) was the highest-frequency haplotype in the OG clade.

Table 3.

Estimates of genetic diversity, neutrality tests and mismatch distributions for each clade in Cellana toreuma complex.

	No. of individuals (n)	Haplotype diversity (h)	Nucleotide diversity (π)	Fu's FS (P)	Tajima's D (P)	Mismatch distribution estimate (τ)	Real expansion time (t, years)
EA clade	469	0.2247 ± 0.0272	0.000431 ± 0.000533	−3.403 × 1038 (0.000)	−2.553 (0.000)	3.000	504,421
ESA clade	34	0.5775 ± 0.0973	0.001956 ± 0.001438	−3.064 (0.023)	−1.196 (0.107)	0.000	–
WSA clade	79	0.6018 ± 0.0361	0.001173 ± 0.000994	−6.566 (0.000)	−1.596 (0.028)	0.877	147,395
PR clade	3	0.6667 ± 0.3143	0.002241 ± 0.002302	1.061 (0.640)	0.000 (0.093)	2.289	384,705
OG clade	13	0.6282 ± 0.1431	0.001250 ± 0.001108	−2.535 (0.005)	−1.437 (0.084)	2.008	337,479

Bold font indicates P < 0.05.

Analyses of population differentiation were conducted between all pairs of the 21 well-sampled populations. Before sequential Bonferroni correction (P < 0.05) there were 128 (60.9%) significant comparisons (Table 4). After sequential Bonferroni correction (P < 0.00024), 99 (47.1%) of the comparisons were significant (Table 4). Pairwise distances between SY and localities from the Chinese mainland (Φ_ST = 0.1921–0.3404) were much lower than p distances between SY and localities from Southeast Asia (Φ_ST = 0.7201–0.8702) (Table 4). The 4 localities (DN, KS, KC and MR) from Southeast Asia presented higher levels of genetic differentiation than 15 localities from the Chinese mainland (Φ_ST = 0.8566–0.9993; P < 0.00026; Table 4). Among localities in Southeast Asia, the distances between DN and KS + KC (Φ_ST = 0.1061–0.2135) were much lower than the others (Φ_ST = 0.7201–0.9983) (Table 4).

Table 4.

Pairwise Φ_ST values for 20 populations of Cellana toreuma complex (below diagonal) and associated P-values (above diagonal).

	QD	DGD	SS	ZS	NJ	PT	NRD	MZD	CW	XM	DS	NA	SZ	HK	WZD	SY	DN	KS	MR	KC	OG
QD	*	0.99990	0.04802	0.00069	0.00772	0.00822	0.00653	0.01148	0.00594	0.00693	0.00356	0.00624	0.02228	0.00178	0.05712	0.00119	0.00000	0.00000	0.00000	0.00000	0.00000
DGD	−0.0372	*	0.07237	0.00010	0.00832	0.00832	0.00277	0.01287	0.00228	0.00386	0.00208	0.00792	0.03495	0.00158	0.07098	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
SS	0.1167	0.0870	*	0.67795	0.81437	0.72518	0.69330	0.77755	0.79745	0.87011	0.54985	0.82081	0.22077	0.24651	0.73913	0.00356	0.00000	0.00000	0.00000	0.00000	0.00000
ZS	0.1622	0.1361	−0.0035	*	0.22255	0.39333	0.31413	0.41065	0.30858	0.23819	0.48183	0.29680	0.08356	0.16771	0.64895	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
NJ	0.1226	0.0992	−0.0103	0.0059	*	0.84348	0.45352	0.79923	0.70349	0.46391	0.21790	0.75240	0.09494	0.26938	0.79537	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
PT	0.1518	0.1183	−0.0005	0.0020	−0.0003	*	0.65716	0.61895	0.79349	0.87209	0.46045	0.94060	0.09583	0.80071	0.69884	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
NRD	0.1427	0.1112	−0.0057	0.0039	0.0004	−0.0005	*	0.58489	0.32630	0.57470	0.99990	0.62667	0.08712	0.23562	0.65588	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
MZD	0.1353	0.1062	−0.0027	0.0027	−0.0019	0.0008	−0.0002	*	0.72221	0.81576	0.38848	0.74745	0.13306	0.24354	0.75319	0.00287	0.00000	0.00000	0.00000	0.00000	0.00000
CW	0.1427	0.1156	−0.0057	0.0039	−0.0089	−0.0005	0.0000	−0.0006	*	0.56608	0.88754	0.82823	0.08583	0.23730	0.77527	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
XM	0.1276	0.1073	−0.0099	0.0052	0.0001	−0.0009	0.0000	−0.0020	0.0000	*	0.94664	0.75339	0.08732	0.23691	0.85180	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
DS	0.1851	0.1316	0.0150	−0.0002	0.0013	0.0013	−0.0159	0.0052	0	−0.0001	*	0.27453	0.08712	0.23621	0.51896	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
NA	0.1371	0.1114	−0.0064	0.0043	−0.0096	−0.0139	0.0001	−0.0007	−0.0108	−0.0002	0.0011	*	0.10207	0.71567	0.80507	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
SZ	0.1250	0.1009	0.0547	0.0625	0.0396	0.0701	0.0565	0.0584	0.0566	0.0426	0.1104	0.0529	*	0.11078	0.26215	0.02930	0.00000	0.00000	0.00000	0.00000	0.00000
HK	0.1653	0.1284	0.0179	0.0147	0.0114	−0.0177	0.0134	0.0163	0.0134	0.0110	0.0230	−0.0132	0.0873	*	0.25116	0.00000	0.00000	0.00000	0.00000	0.00000	0.00000
WZD	0.1121	0.0835	0.0003	−0.0028	−0.0106	0.0008	−0.0052	−0.0022	−0.0052	−0.0101	0.0189	−0.0061	0.0514	0.0191	*	0.03039	0.00000	0.00000	0.00000	0.00000	0.00000
SY	0.2308	0.2650	0.1981	0.3404	0.2667	0.2630	0.2727	0.2413	0.2742	0.2744	0.2762	0.2662	0.1925	0.2699	0.1921	*	0.00000	0.00000	0.00000	0.00000	0.00000
DN	0.8754	0.8925	0.8604	0.9242	0.8954	0.8951	0.8998	0.8847	0.8997	0.8992	0.9012	0.8961	0.8566	0.8987	0.8573	0.7353	*	0.03346	0.00000	0.00000	0.00000
KS	0.9907	0.9897	0.9929	0.9942	0.9927	0.9939	0.9936	0.9934	0.9936	0.9929	0.9949	0.9934	0.9920	0.9942	0.9928	0.8661	0.1061	*	0.00000	0.00000	0.00000
KC	0.9959	0.9941	0.9991	0.9975	0.9972	0.9986	0.9979	0.9985	0.9979	0.9973	0.9993	0.9979	0.9983	0.9987	0.9991	0.8702	0.2135	0.8100	*	0.00000	0.00000
MR	0.9901	0.9874	0.9951	0.9941	0.9929	0.9952	0.9942	0.9947	0.9943	0.9931	0.9966	0.9940	0.9936	0.9954	0.9949	0.7201	0.8588	0.9926	0.9983	*	0.00000
OG	0.9808	0.9793	0.9861	0.9900	0.9848	0.9892	0.9884	0.9877	0.9885	0.9869	0.9915	0.9879	0.9835	0.9898	0.9831	0.6735	0.8182	0.9877	0.9945	0.9661

Bold values are significant at P < 0.00024 after sequential Bonferroni correction. See Table 1 for abbreviations.

The SAMOVA analyses showed that the F_CT values increased with K and reached a plateau at K = 4 (F_CT = 0.96650; ). This grouping strategy (K = 4) was determined, because for K ≥ 5 at least one group had only one population, indicating that the group structure was disappearing. The AMOVA conducted with this grouping of populations revealed that 95.76% of the total variance was found among groups and was significant (Φ_CT = 0.9576; P < 0.001; Table 5). The percentage of variation was 1.18% among populations within groups (P < 0.001) and 3.06% within populations (P < 0.001) (Table 5). The four SAMOVA groups correspond to the EA clade, the WSA clade, clades OG + ESA and clade PR.

Table 5.

Results from analysis of molecular variance (AMOVA) of Cellana toreuma group.

Source of variation	d.f.	Sum of squares	Variation components	Percentage of variation	Fixation index	P-value
Among groups	3	3988.921	20.66010	95.76	0.9576	0.0000
Among locations within groups	34	148.919	0.25438	1.18	0.2780	0.0000
Within locations	539	356.065	0.66060	3.06	0.9694	0.0000
Total	576	4493.904	21.57508

The grouping refers to the best SAMOVA population grouping: group 1 [QD, DGD, SS, ZS, NJ, PT, NRD, MZD, CW, XM, DS, NA, SZ, HK, WZD, SY, LG, OA, YIP, APP, MAP, YMP, HKP, AMP, TAP, SKP, MWP, RKP, PLQ, STP, YL]; group 2 [DN, KS, KC]; group 3 [MR, OG]; group 4 [PR, NG]. See Table 1 for locality abbreviations.

Demographic history

Neutrality testes were negative and significant in both the EA clade (Fu's F_S = −3.403 × 10³⁸, P < 0.001; Tajima's D = −2.553, P < 0.001) and WSA clade (Fu's F_S = −6.566, P < 0.05; Tajima's D = −1.596, P < 0.05) (Table 3), suggesting recent demographic expansion or the action of selection. The significance of Fu's F_S and Tajima's D was inconsistent in the ESA clade (Fu's F_S = −3.064, P = 0.023; Tajima's D = −1.196, P = 0.107; Table 5) and OG clade (Fu's F_S = −2.535, P = 0.005; Tajima's D = −1.437, P = 0.084; Table 5). Skewed unimodal distributions were observed for all five clades in the mismatch distribution analysis (Fig. 4). The timing of population expansion could be estimated with the tau values (τ). Using a substitution rate of 0.5% per myr, the approximate expansion times for the EA, WSA, PR and OG clades were calculated as 504 kyr, 147 kyr, 384 kyr and 337 kyr BP, respectively. BSP analysis revealed a gentle growth in population size from 270 kyr BP for the EA clade and 155 kyr BP for the WSA clade, before a prolonged period of stability (). This also revealed that both the ESA and OG clades were at a stable population size (). Because there were only two haplotypes in the PR clade, no BSP could be constructed. The overall flat nature of the skylines of each clade could result from low nucleotide divergence among the haplotypes within each clade.

Figure 4.

Mismatch distributions for the five COI clades of the Cellana toreuma complex. The histograms are the observed frequencies of pairwise divergences and the grey line shows the expectation under the sudden population expansion model.

DISCUSSION

Phylogeny and systematics of Cellana toreuma

Historically, the taxonomy of Cellana species has depended to a large extent on the morphology of the shell, but this shows wide intraspecific variation, leading to taxonomic uncertainty (Powell, 1973). While recent molecular studies have begun to clarify the limits and relationships of phylogenetic species of Cellana (e.g. Nakano & Ozawa, 2007; Nakano et al., 2009; Nakano & Espinosa, 2010; Dong et al., 2012), problems remain in terms of morphological identification and nomenclature. While the present study revealed the individuals identified nominally as C. toreuma to be a monophyletic group (hereafter termed the C. toreuma complex), this clade also included some identified in a previous study as C. radiata enneagona (Nakano et al., 2009). These latter individuals were widely separated in the phylogeny from other taxa defined by Powell (1973) as ‘subspecies’ of C. radiata. Clearly, the C. radiata complex requires taxonomic revision (Nakano & Sasaki, 2011).

Within the C. toreuma complex, we found no consistent shell morphological differences among the individuals we sampled (see also similar shells figured by Nakano & Espinosa, 2010), although the endemic taxon in the Ogasawara Islands (‘C. radiata enneagona’) is apparently distinctive (see Powell, 1973). Our phylogenetic analysis revealed five highly divergent and well-supported clades within the C. toreuma complex, here referred to as the EA, WSA, PR, ESA and OG clades (Fig. 2). These clades differed by 1.4–7.6% in their COI sequences. Hebert et al. (2004) suggested a ‘ten times threshold rule’ to identify likely cases of species-level divergence. In the present study, the genetic distance between the OG clade and the ESA clade (1.4%) was only about three times the distance within clades. However, for the other three clades the minimum distance between clades (4.6%) was more than ten times the maximum divergence within the clades (0.4%). This suggests that at least four clades could be recognized as distinct phylogenetic species within the C. toreuma complex (EA, WSA, PB, ESA + OG). Previous phylogenetic studies of Cellana species have also discovered a sister relationship between ‘C. toreuma’ from Japan and ‘C. radiata enneagona’ from the Ogasawara Islands (Nakano & Ozawa, 2007; Nakano et al., 2009) and have suggested that C. enneagona could be recognized at specific level (Nakano & Ozawa, 2007).

In the Indo-West Pacific many rocky shore species which were thought to be widely distributed have been shown by molecular analysis to be composed of a number of morphologically cryptic species. Examples include Echinolittorina malaccana (as E. trochoides; Reid et al., 2006), Tetraclita squamosa (Chan, Tsang & Chu, 2007a, b) and Chthamalus malayensis (Tsang et al., 2012). Our results further confirm that Southeast Asia is a hotspot of marine biodiversity and suggest that there are potentially many more cryptic species to be discovered.

Historical events driving the phylogeographical pattern in Southeast Asia and east China

Limpets of the C. toreuma complex have a phylogeographic structure similar to that of Model IV proposed by Maggs et al. (2008). There is a likely zone of secondary contact at Hainan Island, with representatives of both the ESA and EA clades, where recolonists from two or more refugia came into contact. Similarly, in northern Vietnam we found individuals from both the ESA and WSA clades. During the Pleistocene (2.6 Ma–11.7 ka), barriers that formed during times of low sea level influenced population connectivity in tropical coastal species by restricting larval dispersal (and thus gene flow), resulting in high population-genetic structuring (Ludt & Rocha, 2015). Historical isolation of marginal seas and refugia as a result of sea-level changes could have imposed barriers for C. toreuma and resulted in divergent allopatric clades, which have since come into secondary contact. However, this conclusion is tentative, as it is impossible to distinguish the secondary contact hypothesis from the alternative that some individuals simply retain ancestral mtDNA polymorphisms, or that the haplotype sharing represents secondary contact and introgression. Without additional loci or distinguishing morphological traits for these clades, we cannot determine if the mtDNA clades represent species with some secondary contact (with or without introgression) or geographically isolated species that retain some ancestral polymorphisms. Thus, these hypotheses need further study in the future.

The impacts of lowered sea level on the changing topography of this region are well understood (Woodruff, 2010; Keyse et al., 2014; Ni et al., 2014). During intervals of low sea level in the Pleistocene, the East China Sea (ECS) was reduced to a narrow trough and the SCS became an inland sea, the ECS and SCS being separated by the Taiwan land bridge (Fig. 1). This land bridge had significant vicariant effects on biogeographic and phylogeographical patterns of coastal species (including fishes, molluscs and crustaceans) in the northwestern Pacific (Ni et al., 2014) and could also have brought about the divergence of the EA and ESA clades of C. toreuma. During interglacial periods, the EA clade experienced a rapid demographic expansion by large-scale larval dispersal via coastal currents, leading to its wide distribution in EA. The ESA clade is distributed in Hainan, Vietnam and Brunei (SY, DN and MR), localities which were located along the historic coastline of the SCS during glacial intervals, when the sea level was up to 130 m lower than at present (Fig. 1). We hypothesize that the ESA clade simply colonized the nearby available habitats as sea level subsequently rose and did not undergo significant demographic expansion. The results of neutrality tests, mismatch distribution and BSP analyses support this hypothesis.

The formation of the WSA and PR clades could be due to isolation in refugia during glacial intervals. During the Pleistocene glaciations, enclosed seas such as the SCS, Sulu Sea, Celebes Sea and Flores Sea (Fig. 1) are believed to have played important roles in the speciation and population-level divergences of marine fauna and flora (Briggs, 1999; Barber et al., 2006; Timm & Kochzius, 2008; Chan et al., 2013; Ludt & Rocha, 2015). These isolated refugia could account for the formation of the WSA and PR clades of C. toreuma. The wide distribution of the WSA clade was perhaps achieved by northward expansion from a southern refugium since 147 kyr BP when sea level rose. However, the fine-scale patterns of the WSA and PR clades cannot be determined because there is a large sampling gap in Indian Ocean. As Nakano et al. (2009) suggested, the formation of the OG clade may be associated with the southward extension of an ancestral Cellana species using the new islands of the Ogasawara chain as stepping-stones and their subsequent isolation and divergence caused by the rising sea level and/or tectonic activity. Further study is necessary to resolve the phylogeographic history of limpets in this region.

Local environmental factors maintaining genetic subdivision

Current oceanographic conditions may promote, maintain or homogenize the footprint of genetic divergence created by historical events. The genetic similarity between populations in Thailand and Vietnam (KS and DN) could be explained by the influence of the South China Sea Warm Current (SCSWC), which is believed to promote gene flow in the barnacle C. malayensis along the coastline of southern China, Vietnam and eastern Thailand (Tsang et al., 2012). Localities SY and DN have been recognized as two possible locations of secondary contact under the scenario of Model IV (of Maggs et al., 2008). Results of pairwise Φ_ST analysis suggest that there is limited contemporary gene flow between SY and DN, even though both of them are under the influence of SCSWC, indicating that with little or no contemporary gene exchange they have been able to maintain their genetic record of historical recolonization.

Habitat availability also appears to be important for the formation of genetic structure in the C. toreuma complex. Although two localities in northwestern Borneo (MR and KC) are close in geographic distance and both of them are under the influence of the NW Luzon Coastal Current (Fig. 1), MR shows haplotypes from the ESA clade while KC has a unique haplotype from the WSA clade. Between these two localities there is a stretch of about 600 km of sandy beach, which could interrupt gene flow. In southern Australia, it has been shown that a sandy beach without rocky substrate causes a phylogeographical break in the barnacle Catomerus polymerus (York, Blacket & Appleton, 2008). Similarly, in the Coral Triangle habitat availability has a role in shaping the high biodiversity of this area (Lourie & Vincent, 2004; Williams & Reid, 2004; Lourie, Green & Vincent, 2005; Reid et al., 2006).

In conclusion, while historical events have likely played the dominant role in shaping the present phylogeographical pattern of the C. toreuma complex in Southeast Asia and East Asia, some contemporary factors, including ocean currents and habitat availability, are also important in modifying or maintaining the genetic footprint of past environments.

ACKNOWLEDGEMENTS

This work was supported by grants from National Basic Research Program of China (2013CB956504), National Natural Science Foundation of China (41276126, 41476115), Program for New Century Excellent Talents in University of Fujian Province and New Century Excellent Talents of Ministry of Education. Authors appreciate the help of Prof. Gray Williams (The University of Hong Kong), Dr David Marshall (University of Brunei Darussalam) and Prof. Shichun Sun (Ocean University of China) in collecting samples. The authors also thank two anonymous reviewers and Associate Editor Peter Marko for helpful suggestions that have improved this paper.