Abstract

Littorina brevicula is one of the most common gastropods in the supralittoral zone around Japan. The northernmost population of this species is around Hokkaido and the determinant of this northern limit is likely seawater and air temperature. To reconstruct an evolutionary history of this species, we investigated genetic differentiation among 12 populations (three from Hokkaido, six from Honshu and three from Kyushu) using a mitochondrial DNA marker (partial sequence of the NADH dehydrogenase 6 gene). The haplotype network showed shallow genetic divergence within the species, suggesting a bottleneck followed by population expansion. One major haplotype that occurred in 70.5% of all individuals examined was the most frequent in every population sampled. A second major haplotype was abundant around Kyushu but not found in Hokkaido. This skewed haplotype distribution resulted in significant genetic differentiation along the north-south axis of Japan. The importance of the southern clade, which included the second major haplotype, was supported by population genetic analyses of datasets that excluded either the southern clade or the northern clades. The north-south differentiation remained when datasets that excluded the northern clades were used, but disappeared when datasets that excluded the southern clade were used. The combined evidence of shallow divergence and the north-south population structure suggests that the L. brevicula population around Japan once declined and then expanded and colonized northward. Although the time of population reduction and recolonization could not be precisely estimated, the observation that this species is absent further north in Japan suggests that it would have been unable to survive in northern Japan during the last glacial maximum (LGM) and therefore recolonization likely occurred after the LGM, probably from south to north.

INTRODUCTION

Many species in the northern hemisphere probably underwent changes to their distribution as a result of strong influence from climatic oscillation during glacial periods, particularly in the last glacial maximum (LGM) that lasted until around 20,000 years ago and in the subsequent rapid warming. Such distributional changes can be detected by phylogeographic investigations using suitable molecular markers (Hewitt, 2000). Coastal ecosystems were considerably influenced by the LGM and thus various species exhibit genetic evidence of rapid expansion and/or colonization, probably occurring in the postglacial period. Examples include the nearshore fish Syngnathus leptorhynchus (Wilson, 2006) and Xiphister mucosus (Marko et al., 2010), seaweed Pelvetia canaliculata (Neiva et al., 2014), gastropods Nassarius nitidus (Albaina et al., 2012), Littorina sitkana and L. scutulata (Marko et al., 2010), and starfish Pisaster ochraceus and Evasterias troschelii (Marko et al., 2010). Around Japan, glacial effects on intertidal molluscs have been suggested in Batillaria cumingii (Kojima et al., 2004) and Cellana nigrolineata (Nakano, Sasaki & Kase, 2010). However, the more stable temperature of the ocean compared with that of the land allowed some coastal species to maintain their genetic diversity during glacial maxima. Furthermore, the glacial effects sometimes varied within a single species, resulting in a complex population genetic structure, as in Littorina saxatilis, in which some populations seemed to have colonized following climate change, while others maintained diversity during a long continuous history (Panova et al., 2011). Such a spatially variable effect within a single species has also been suggested for Nucella ostrina and N. lamellosa (Marko et al., 2010). The variable genetic consequences of the LGM between or within species are probably due to specific ecological characteristics of each species, i.e. distribution range, habitat, life cycle and abundance, and also to area-specific physical effects of climatic oscillation such as temperature and sea-level change.

In addition to the effects of climate oscillations over palaeontological time scales, genetic population structure within a species is affected by contemporary gene flow. Populations of marine species are often weakly structured (i.e. show a low level of differentiation among local populations and/or weak correlation between geography and genetic differentiation) over the species' range (Palumbi, 1994), due to a high level of gene flow caused by active and passive dispersal. In invertebrates that generally show poor active mobility, gene flow has been thought to be largely affected by the developmental mode and the length of the planktonic larval period, and colonization and contemporary gene flow in marine species are generally mediated by larval dispersal along ocean currents (Palumbi, 1994; Tsang et al., 2008). Recently, larval dispersal has been revealed to be less crucial for determining the genetic structure of marine species than previously supposed (Marko, 2004; Weersing & Toonen, 2009), although it remains an important factor. A comparative study of two planktonic- and two direct-developing Littorina species suggested that planktonic larval dispersal reduced genetic differentiation between local populations (Kyle & Boulding, 2000). In this study the planktonic-developing species (Littorina scutulata and L. plena) showed very low levels of genetic differentiation among local populations, while the direct-developing L. subrotundata showed higher structuring with larger differentiation between localities. Subsequent research revealed that direct-developing Littorina species showed spatial rather than temporal differentiation and that planktonic-developing species showed the reverse, due to sweepstakes-like reproductive success (Lee & Boulding, 2009). These studies show that Littorina is a good model in which to investigate the relationships between genetic population structure, larval dispersal and historical factors in marine coastal species.

Littorina brevicula occurs in the littoral fringe of the temperate coast of the northwestern Pacific. Its range extends from Hong Kong to Peter the Great Bay along the continental coast of Asia and from Okinawa to Hokkaido along the Japanese Archipelago (Reid, 1996; Okutani, 2000). On the central to northern Japanese coast, L. brevicula is one of the most common snails, occurring in dense aggregations on rocky and boulder shores and on artificial constructions such as breakwaters, tetrapods and slipways. However, it is scarce along the southeastern coast of Hokkaido, which is under the influence of a cold current (Ohgaki, 1983). The spawning season is January to April around Japan and Korea (Kojima, 1957; Son & Hong, 1998). The egg capsule is pelagic and the veligers hatch about 7 d after spawning (Son & Hong, 1998). The total length of the planktonic egg and larval period is several weeks (Golikov, 1976). Such a long planktotrophic period is expected to lower the genetic differentiation between local populations of L. brevicula (Reid, 1996). Nevertheless, a spatial genetic structure could potentially be found if samples from a wide geographic range were investigated using highly variable molecular markers (Palumbi, 1994). The fine details of genetic structure across a wide area may provide knowledge about the effects of environmental change on coastal species on a palaeontological scale.

Several studies of population genetics in L. brevicula using allozyme variation appeared in the 1990s (e.g. Zaslavskaya, Sergievsky & Tatarenkov, 1992; Tatarenkov, 1995; Zaslavskaya & Takada, 1998; Park et al., 1999) and some of these studies reported genetic differentiation among samples. This was variously attributed to reduced gene flow between distant populations (the east and west coasts of the Sea of Japan; Zaslavskaya & Takada, 1998), the effect of heavy metal pollution (Park et al., 1999) and possible local selection on specific genes (Tatarenkov, 1995). Since these studies did not aim to clarify the total phylogeography of L. brevicula throughout the distribution range, they provided little information about hierarchical genetic structure among populations. Nevertheless, Zaslavskaya & Takada (1998) showed genetic differentiation between distant populations, which suggested the possibility of some structuring throughout the species' range.

More recently, DNA markers were used to investigate population structure in L. brevicula around Korea (Kim et al., 2003b; Kim, Rodriguez-Lanetty & Song, 2003a). Kim et al. (2003b) failed to detect genetic differences among populations, even for populations 1,000 km apart along the coastline. Kim et al. (2003a) found significant population differentiation, which they attributed to a loss of genetic diversity caused by heavy metal pollution of some populations.

In this study, we investigated the genetic structure among populations of L. brevicula in its the main distributional area around Japan, from Hokkaido to Kyushu (Fig. 1). These populations span c. 2,500 km and experience considerably different climatic conditions; therefore, we expected to observe genetic structure. We chose the mitochondrial NADH dehydrogenase 6 gene (ND6) as a genetic marker, because this region has previously been shown to be more polymorphic than the commonly-used cytochrome-b in L. brevicula (Kim et al., 2003a). We found population genetic structure that revealed the evolutionary history of L. brevicula around Japan, and suggested patterns of population decline, expansion and migration during recent glacial and postglacial periods.

Figure 1.

Ocean currents around Japan, and sampling locations of Littorina brevicula in the present study.

Figure 1.

Ocean currents around Japan, and sampling locations of Littorina brevicula in the present study.

MATERIAL AND METHODS

Sampling, DNA sequencing and analysis of haplotype genealogy

We examined 540 individuals collected from 12 populations in three regions: three populations from Hokkaido (SSB, ABS, USJ), six from Honshu (FKU, HCH, OBM, CHO, URG, OSK) and three from Kyushu (SIK, AMK, MKR) (Table 1, Fig. 1). Genomic DNA was extracted from a piece of muscle (c. 1 mm3) of ethanol-fixed specimens using the Pure Gene Kit (Qiagen) according to the manufacturer's protocol and resuspended in 100 µl of Tris-EDTA buffer.

Table 1.

Geographical origin and genetic information estimated from partial ND6 sequence for Littorina brevicula populations analysed in the present study.

Region Sample code and name Coordinates Collection date Sample size No. of haplotypes HR h π 
Hokkaido 1. SSB 44°32′N; 141°45′E July 2010 48 10 5.5 0.343 0.0011 
 2. ABS 44°03′N; 144°15′E June 2010 39 5.6 0.371 0.0012 
 3. USJ 41°56′N; 140°56′E June 2013 48 3.5 0.199 0.0004 
Mean of Hokkaido      4.8 0.304 0.0009 
Honshu 4. FKU 40°46′N; 140°03′E Oct 2011 48 11 7.6 0.647 0.0018 
 5. HCH 40°32′N; 141°33′E Nov 2011 48 11 6.8 0.469 0.0015 
 6. OBM 35°32′N; 135°42′E Oct 2011 48 12 7.6 0.578 0.0018 
 7. CHO 35°42′N; 140°52′E Jan 2012 47 17 9.8 0.646 0.0021 
 8. URG 35°15′N; 139°44′E Aug 2013 48 14 8.4 0.629 0.0025 
 9. OSK 34°17′N; 132°54′E Nov 2011 47 4.4 0.311 0.0009 
Mean of Honshu      7.4 0.546 0.0018 
Kyushu 10. SIK 32°48′N; 131°53′E Sept 2011 46 3.6 0.423 0.0011 
 11. AMK 32°36′N; 130°29′E Sept 2011 24 0.409 0.0009 
 12. MKR 31°15′N; 130°18′E May 2013 49 14 7.9 0.677 0.0020 
Mean of Kyushu      4.8 0.503 0.0013 
Total    540 78 6.0* 0.497** 0.0015** 
Region Sample code and name Coordinates Collection date Sample size No. of haplotypes HR h π 
Hokkaido 1. SSB 44°32′N; 141°45′E July 2010 48 10 5.5 0.343 0.0011 
 2. ABS 44°03′N; 144°15′E June 2010 39 5.6 0.371 0.0012 
 3. USJ 41°56′N; 140°56′E June 2013 48 3.5 0.199 0.0004 
Mean of Hokkaido      4.8 0.304 0.0009 
Honshu 4. FKU 40°46′N; 140°03′E Oct 2011 48 11 7.6 0.647 0.0018 
 5. HCH 40°32′N; 141°33′E Nov 2011 48 11 6.8 0.469 0.0015 
 6. OBM 35°32′N; 135°42′E Oct 2011 48 12 7.6 0.578 0.0018 
 7. CHO 35°42′N; 140°52′E Jan 2012 47 17 9.8 0.646 0.0021 
 8. URG 35°15′N; 139°44′E Aug 2013 48 14 8.4 0.629 0.0025 
 9. OSK 34°17′N; 132°54′E Nov 2011 47 4.4 0.311 0.0009 
Mean of Honshu      7.4 0.546 0.0018 
Kyushu 10. SIK 32°48′N; 131°53′E Sept 2011 46 3.6 0.423 0.0011 
 11. AMK 32°36′N; 130°29′E Sept 2011 24 0.409 0.0009 
 12. MKR 31°15′N; 130°18′E May 2013 49 14 7.9 0.677 0.0020 
Mean of Kyushu      4.8 0.503 0.0013 
Total    540 78 6.0* 0.497** 0.0015** 

Abbreviations: HR, haplotype richness (number of haplotypes standardized by the smallest sample size); h, haplotype diversity; π, nucleotide diversity; *, average for all samples; **, calculated pooling all samples as a single population.

A partial sequence of ND6 was amplified by PCR in a 30-μl reaction mixture containing template DNA (c. 200 pg), dNTPs, a pair of primers [Lbnd-F (5′-AGG TAC ATA TTC CTG CGC TCT GAA A-3′) and Lbnd-R (5′-GTG TGC GCA TGA AAT GTA T-3′)] (Kim et al., 2003a,b) and ExTaq (Takara Bio, Inc.), according to the manufacturer's instructions. The thermal cycling profile included precycling denaturation at 95 °C for 1 min, followed by 35 cycles of denaturation at 95 °C for 30 s, annealing at 50 °C for 20 s, and extension at 72 °C for 45 s. The PCR products were then examined by electrophoresis on a 2% agarose gel, purified with magnetic beads (AMPure, Agencourt), cycle-sequenced using the above-mentioned forward primer and the BigDye® Terminator v. 3.1 Cycle Sequencing Kit (Applied Biosystems) and loaded onto an automated sequencer, ABI PRISMTM 3130 (Applied Biosystems). The obtained sequences were aligned and edited to 450 bp using DNASIS-Mac v. 3.5 (Hitachi) and ClustalX v. 1.81 software (Thompson et al., 1997) for defining haplotypes and deposited in the DDBJ/Genbank database with accession numbers AB968435-AB968512.

The haplotype genealogy was resolved with a parsimony network using the TCS Network Program (Clement, Posada & Crandall, 2000) under a 95% connection limit. The age of lineage divergence was estimated using the divergence rate of 2.46% per million years (Myr). The rate was based on a calibration for mtDNA including the ND6 region in L. saxatilis (Quesada et al., 2007), which assumed a sequence divergence rate of 1.83 ± 0.21% per Myr. Since Quesada et al. applied this rate for the 1.8-kbp sequence including a not very variable region, the calibration may underestimate the evolutionary rate of our ND6 data. To correct the rate, we examined the GenBank data for L. saxatilis (accession nos AM500948-500955, deposited by Quesada et al., 2007) and compared sequence divergence within 1.8 kbp of the full length of these sequences and 0.45 kbp of the ND6 region, orthologous to our sequence data. In the 1.8-kbp sequence, the number of different nucleotides per site was 0.145, while it was 0.195 in the ND6 region in the L. saxatilis data. We therefore employed the corrected divergence rate of 2.46% (1.83 × 0.195/0.145) per Myr.

Population genetic analysis

Arlequin v. 3.1 (Excoffier, Laval & Schneider, 2005) was used to estimate haplotype (h) and nucleotide diversity (π) in each population and to calculate pairwise FST (Weir & Cockerham, 1984). Genetic distance between populations based on the pairwise FST was visualized on a two-dimensional surface by nonmetric multidimensional scaling (nMDS) plotting using statistical software R v. 2.9.0 (R Development Core Team, 2005). To assess geographic-genetic correlation, the isolation-by-distance (IBD) model (Wright, 1943) was tested using Arlequin v. 3.1. For the IBD test, the matrix of geographic distances along the coastline between sampling sites was compared with the FST matrix, and the significance of the correlations was evaluated using the Mantel test with 10,000 permutations.

To test the significance of the hierarchical population structure, analysis of molecular variance (AMOVA; Excoffier, Smouse & Quattro, 1992) was conducted with Arlequin v. 3.1, comparing the three categories of clustering suggested by geography and/or nMDS plotting based on FST analysis: (1) three regional group clustering—Hokkaido, Honshu and Kyushu, i.e. [SSB, ABS, USJ], [FKU, HCH, OBM, CHO, URG, OSK] and [SIK, AMK, MKR], respectively; (2) three groups suggested by nMDS plotting based on pairwise FST analysis—Hokkaido + HCH, Honshu (except HCH) and Kyushu, i.e. [SSB, ABS, USJ, HCH], [FKU, OBM, CHO, URG, OSK] and [SIK, AMK, MKR], respectively; and (3) two groups of clustering by coast—the Sea of Japan and the Pacific Ocean (including the Yellow Sea), i.e. [SSB, ABS, FKU, OBM] and [USJ, HCH, CHO, URG, OSK, SIK, AMK, MKR], respectively. Since Category (3) took into account ocean currents, ABS on the coast of the Sea of Okhotsk, where a branch of the Tsushima current flows in from the Sea of Japan (Fig. 1), was included in the Sea of Japan cluster.

When we detected significant genetic structure, we evaluated the effect of each clade (haplotype group derived from the most abundant haplotype shown in Fig. 3) on population structure by preparing virtual datasets for population analyses, performing an IBD test, estimating pairwise FST and conducting AMOVA using these datasets. Datasets were individually modified by excluding an area-restricted clade (haplotype groups connected to each other in a haplotype network; Fig. 3) and then used to test whether or not the structure appeared. If we found that the structure was lost when a certain clade was excluded from the modified dataset, we would know that the excluded clade was essential for the structure.

To test recent population expansions, neutrality tests and mismatch distribution (MMD) analysis with Arlequin v. 3.1 were added. When our MMD analysis detected population and range expansion, the expansion age was estimated based on the τ value, assuming a sequence-divergence rate of 2.46% per Myr. Neutrality tests of Fu's FS and Tajima's D were conducted to detect the possibility of population expansion.

Based on the genealogy and distribution of haplotypes, we detected the population structure and inferred the shift of the distribution area following a series of steps in the flow chart in Figure 2. The results of several analyses were employed to test the determinants in Figure 2: pairwise FST estimation for step 1; test of isolation by distance (IBD), nMDS plotting of populations based on pairwise FST and AMOVA, for steps 2 and 4; construction of the haplotype network and survey of the genealogy for step 3; and estimation of indices of genetic diversity in each population and each area for step 5. The first step, checking heterogeneity between populations, is also informative to evaluate the level of contemporary gene flow, which causes panmictic populations all through the distribution area in extreme cases. The second step, checking the genetic-geographic association, can reveal genetic drift and is useful for detecting population structures skewed by human impact as in L. brevicula (Kim et al., 2003b) and Neptunea arthritica around Japan (Azuma et al., 2015).

Figure 2.

Flow chart to evaluate glacial effects on population structure and postglacial colonization inferred from genetic data. The path of inference in the present study follows the circled determinations (present or absent).

Figure 2.

Flow chart to evaluate glacial effects on population structure and postglacial colonization inferred from genetic data. The path of inference in the present study follows the circled determinations (present or absent).

RESULTS

Genealogy and distribution of haplotypes

A total of 78 haplotypes were detected. LbN6AB03 was the most abundant haplotype, occurring in 381 individuals (70.5% of all individuals examined in the present study). This haplotype was also the most frequent in all populations, and was located at the centre of the haplotype network (Fig. 3). The second most frequent haplotype, LbND6OT02, occurred in 44 individuals (8.1% of all individuals examined) from Honshu and Kyushu, and the third, LbND6AO10, occurred in 10 individuals from Honshu and Hokkaido. Other haplotypes occurred in less than 9 individuals and 65 haplotypes were private, i.e. found exclusively in a single population (see Supplementary Material). Given that nucleotide divergence between LbN6AB03 and LbN6OT02 was 0.22% per site, the time of divergence was estimated as 89 thousand years ago (kya). The largest number of nucleotide substitutions observed was six, 1.33% per site.

Figure 3.

Parsimony network of the mtDNA ND6 haplotypes of Littorina brevicula. A solid line between circles indicates a single nucleotide substitution. Circle size reflects haplotype abundance (number of individuals that had the haplotype) and the colours in the circles indicate the proportions of individuals from Hokkaido, Honshu and Kyushu (see also Fig. 1). Dotted lines delimit each clade (a monophyletic group of multiple haplotypes). Colours of dotted lines indicate clade distribution (northern, southern or middle).

Figure 3.

Parsimony network of the mtDNA ND6 haplotypes of Littorina brevicula. A solid line between circles indicates a single nucleotide substitution. Circle size reflects haplotype abundance (number of individuals that had the haplotype) and the colours in the circles indicate the proportions of individuals from Hokkaido, Honshu and Kyushu (see also Fig. 1). Dotted lines delimit each clade (a monophyletic group of multiple haplotypes). Colours of dotted lines indicate clade distribution (northern, southern or middle).

Several clades appeared in the haplotype network. Six of the clades were found only in Hokkaido and Honshu (northern clades), and one only in Honshu and Kyushu (a southern clade). The rest of the clades appeared in Hokkaido, Honshu and Kyushu, or only in Honshu (middle clades). Some of the following analyses employed modified datasets that excluded the northern clades or the southern clade, to evaluate the importance of these clades in creating the significant genetic structure that we detected.

Genetic population structure

As shown in Table 1, h in each population ranged from 0.199 for USJ to 0.677 for MKR and π ranged from 0.0004 for USJ to 0.0025 for URG. The mean diversities of each population were 0.304 (h) and 0.0009 (π) in Hokkaido, 0.546 (h) and 0.0018 (π) in Honshu, and 0.503 (h) and 0.0013 (π) in Kyushu. Although these diversity indices seemed highest in Honshu followed by Kyushu and Hokkaido, no significant differentiation between these three areas appeared in the Kruskal-Wallis test for both h2 = 4.5256, df = 2, P = 0.104) and π2 = 3.2488, df = 2, P = 0.197).

The pairwise FST values between Hokkaido and Kyushu populations was always higher than 0.07 and always was significantly different from 0, whereas none of the pairwise values from population comparisons within Hokkaido and Kyushu were significant (α = 0.01; Table 2). These results indicated a north-south genetic differentiation in this species. The HCH population in Honshu, which did not include haplotypes from the southern clade, showed a significant difference from populations in Kyushu and the URG population. The nMDS plot based on FST indicated that genetic distance was correlated with a spatial distribution at the region (island) level (Fig. 4A). In the two-dimensional plot, Hokkaido populations were on the left, Kyushu on the right and Honshu in the middle. However, the relationship among populations within each region did not always follow this genetic-geographic trend.

Table 2.

Lower left: pairwise FST between Littorina brevicula populations based on partial sequence of ND6.

Population  10 11 12 
Hokkaido 1. SSB  − − − − − − − 
2. ABS 0.010  − − − − − − 
3. USJ 0.000 0.010  − − − − 
Honshu 4.FKU 0.033 0.016 0.034  − − − − − − − − 
5. HCH 0.009 0.000 0.010 0.009  − − − − 
6. OBM 0.014 0.017 0.021 0.005 0.014  − − − − − − 
7. CHO 0.007 0.002 0.005 0.000 0.001 0.000  − − − − 
8. URG 0.010 0.021 0.017 0.016 0.020 0.011 0.002  − − − − 
9. OSK 0.006 0.008 0.009 0.011 0.006 0.006 0.000 0.010  − − − 
Kyushu 10. SEK 0.089 0.093 0.128 0.017 0.078 0.027 0.030 0.027 0.041  − − 
11. AMK 0.073 0.077 0.138 0.003 0.060 0.013 0.013 0.011 0.030 0.000  − 
12. MKR 0.080 0.079 0.098 0.019 0.071 0.032 0.034 0.030 0.042 0.000 0.000  
Population  10 11 12 
Hokkaido 1. SSB  − − − − − − − 
2. ABS 0.010  − − − − − − 
3. USJ 0.000 0.010  − − − − 
Honshu 4.FKU 0.033 0.016 0.034  − − − − − − − − 
5. HCH 0.009 0.000 0.010 0.009  − − − − 
6. OBM 0.014 0.017 0.021 0.005 0.014  − − − − − − 
7. CHO 0.007 0.002 0.005 0.000 0.001 0.000  − − − − 
8. URG 0.010 0.021 0.017 0.016 0.020 0.011 0.002  − − − − 
9. OSK 0.006 0.008 0.009 0.011 0.006 0.006 0.000 0.010  − − − 
Kyushu 10. SEK 0.089 0.093 0.128 0.017 0.078 0.027 0.030 0.027 0.041  − − 
11. AMK 0.073 0.077 0.138 0.003 0.060 0.013 0.013 0.011 0.030 0.000  − 
12. MKR 0.080 0.079 0.098 0.019 0.071 0.032 0.034 0.030 0.042 0.000 0.000  

Bold font indicates P value <0.05. Upper right upper: significance level of FST. + indicates P value <0.01 after sequential Bonferroni correction.

Figure 4.

nMDS plot of Littorina brevicula with pairwise FST values based on a 450-bp sequence of partial mtDNA ND6 for unmodified dataset (A) and a dataset excluding the southern clade (B). While population structure along the Japanese Archipelago from north to south (from left to right in this plot) appeared in A, no such structure was observed in B, indicating that the north-south structure was attributable to the existence of the southern clade.

Figure 4.

nMDS plot of Littorina brevicula with pairwise FST values based on a 450-bp sequence of partial mtDNA ND6 for unmodified dataset (A) and a dataset excluding the southern clade (B). While population structure along the Japanese Archipelago from north to south (from left to right in this plot) appeared in A, no such structure was observed in B, indicating that the north-south structure was attributable to the existence of the southern clade.

The sequence of Honshu populations did not follow the spatial pattern. FKU, the northernmost among Honshu populations, appeared closest to the Kyushu populations, while HCH, the second northernmost population, was genetically closer to the Hokkaido populations than to the other Honshu populations. The OSK population, the southernmost in Honshu, was genetically close to populations in Hokkaido and not to those in Kyushu. The reason for this variance from an IBD model will be discussed later.

The Mantel test showed a significant correlation between genetic (FST) and geographic distance (P < 0.01), indicating that the set of study populations of L. brevicula followed an IBD model.

As shown in Table 3, AMOVA supported hierarchical population structure clustering (1) comprising Hokkaido, Honshu and Kyushu and (2) comprising Hokkaido + HCH, Honshu (except HCH) and Kyushu (P = 0.001 and 0.000, respectively). However, the AMOVA did not support the structure of clustering (3) that comprised two groups—the Sea of Japan and the Pacific Ocean (P = 0.762)—indicating that the west and east sides of the Japanese Archipelago could not be separated by their genetic profiles. In each clustering, the variance component was the highest in the category of ‘within a population’ and much higher than ‘among groups' or ‘among populations within a group’, because of the relatively high genetic diversity in each population and the genetic similarity between populations due to the high frequency of LbN6AB03 in every population.

Table 3.

Results of five separate AMOVA of Littorina brevicula populations using K2P distance.

 Source of variation Variance component % of variation P value Fixation indices 
(1) Among groups 0.00994 2.86 0.001 FCT = 0.028 
Among population within groups 0.00124 0.36 0.154 FSC = 0.003 
Within populations 0.33538 96.78 0.000 FST = 0.032 
(2) Among groups 0.01035 2.99 0.000 FCT = 0.029 
Among population within groups 0.00059 0.17 0.260 FSC = 0.001 
Within populations 0.33538 96.84 0.000 FST = 0.031 
(3) Among groups −0.00128 −0.37 0.762 FCT = −0.003 
Among population within groups 0.00852 2.48 0.000 FSC = 0.021 
Within populations 0.33538 97.89 0.000 FST = 0.025 
(2a) Among groups 0.01111 5.02 0.000 FCT = 0.050 
Among population within groups −0.00048 −0.21 0.570 FSC = −0.002 
Within populations 0.21044 95.1 0.000 FST = 0.048 
(2b) Among groups −0.00015 −0.12 0.679 FCT = −0.000 
Among population within groups 0.00219 0.81 0.000 FSC = 0.007 
Within populations 0.28352 99.31 0.000 FST = 0.007 
 Source of variation Variance component % of variation P value Fixation indices 
(1) Among groups 0.00994 2.86 0.001 FCT = 0.028 
Among population within groups 0.00124 0.36 0.154 FSC = 0.003 
Within populations 0.33538 96.78 0.000 FST = 0.032 
(2) Among groups 0.01035 2.99 0.000 FCT = 0.029 
Among population within groups 0.00059 0.17 0.260 FSC = 0.001 
Within populations 0.33538 96.84 0.000 FST = 0.031 
(3) Among groups −0.00128 −0.37 0.762 FCT = −0.003 
Among population within groups 0.00852 2.48 0.000 FSC = 0.021 
Within populations 0.33538 97.89 0.000 FST = 0.025 
(2a) Among groups 0.01111 5.02 0.000 FCT = 0.050 
Among population within groups −0.00048 −0.21 0.570 FSC = −0.002 
Within populations 0.21044 95.1 0.000 FST = 0.048 
(2b) Among groups −0.00015 −0.12 0.679 FCT = −0.000 
Among population within groups 0.00219 0.81 0.000 FSC = 0.007 
Within populations 0.28352 99.31 0.000 FST = 0.007 

(1) 3 groups = [Hokkaido] [Honshu] [Kyushu]. (2) 3 groups = [Hokkaido + HCH] [Honshu (except HCH)] [Kyushu]. (3) 2 groups = [Pacific Ocean + Yellow Sea] [Sea of Japan + Sea of Okhotsk]. (2a) Population groups similar to (2), but all northern clades were excluded before analysis. (2b) Population groups similar to (2), but the southern clade was excluded before analysis.

When the significance of the northern clades and the southern clade was tested using modified datasets, it appeared that the southern clade was essential for population structure. While the unmodified original dataset, the six modified datasets that each excluded one of the six northern clades, and a dataset that excluded all northern clades all showed significant IBD correlations (P < 0.01), the genetic-geographic correlation was not significant (P = 0.840) when a modified dataset that excluded the southern clade was used. Further, the nMDS plot based on pairwise FST values showed no genetic structure for the modified dataset that excluded the southern clade (Fig. 4B). The significance of the hierarchical structure of [Hokkaido + HCH], [Honshu except HCH] and [Kyushu] was supported by the lowest P-value (0.000) in AMOVA with the unmodified dataset. Table 3 shows the results of analyses on modified datasets testing this structure where all six northern clades were excluded ((2a), significant at P = 0.000) and where only the southern clade was excluded ((2b), not significant at P = 0.679), indicating the importance of the southern clade as the determinant of the north-south structure of L. brevicula populations around Japan.

Demographic history

The MMD analysis with pooled populations showed that the observed frequencies of haplotype difference did not deviate significantly from the sudden-expansion model when either the sum of squared deviation (SSD; P = 0.995) or Harpending's raggedness index (P = 0.724) was used. Additionally, the deviation from spatial expansion assuming a constant deme (local population) size model was not significant using either the SSD (P = 0.990) or Harpending's raggedness index (P = 0.730). The τ values corresponding to the sudden expansion of a population and spatial expansion were estimated to be 0.72 and 0.70, respectively. Given τ = 2ut (2u = divergence rate, t = time from expansion) and 2u = 2.46% per Myr, the population and spatial expansion dates were estimated to be 65.0 and 63.2 kya, respectively (i.e. Upper Pleistocene).

Both Fu's FS (−29.709, P = 0.0000) and Tajima's D (−2.653, P = 0.0000) indicated that the observed haplotype frequency deviated from the one expected in neutrality, suggesting that the population expanded.

DISCUSSION

Characterization of population genetic structure of Littorina brevicula around Japan

The haplotype distribution of L. brevicula showed a population genetic structure with north-south differentiation. As expected from the planktotrophic developmental mode of L. brevicula (Reid, 1996), genetic differentiation between geographically close populations is not significant and the significant genetic differentiation observed between distant populations agreed with the results of Zaslavskaya & Takada (1998). The correlation of geographic and genetic distance was statistically significant in the IBD test and the north-south structure was supported by both AMOVA and FST analyses. The deviation from IBD within the Honshu group (Fig. 4A) was probably caused by a stochastic balance of private haplotypes resulting from genetic drift and contemporary gene flow, similar to the lack of genetic structure in the Korean populations (Kim et al., 2003a). In the case of OSK, however, it was possible to attribute the anomalous situation to recent human disturbance. OSK showed an unexpected deviation from IBD (Fig. 4A) and the lowest observed values of both haplotype and nucleotide diversity among all the Honshu populations, which could be the result of high genetic drift. Its small effective population size could be attributed to the effects of human impacts such as pollution, as has been argued for the Korean populations (Kim et al., 2003b). OSK is in the Seto Inland Sea, where pollution likely accumulates more than in the open sea. Overall, genetic-geographic correlation and a fair level of genetic diversity within each population indicated that the examined populations had not been as strongly affected by recent human activity as was the case of Neptunea arthritica (Azuma et al., 2015). The higher resistance of L. brevicula to anthropogenic effects compared with that of N. arthritica is attributable to much larger population sizes and higher gene flow in L. brevicula than in N. arthritica.

Population genetic structure can be greatly affected by ocean currents. Littorina brevicula distributed around Japan currently experiences warm currents (the Kuroshio and Tsushima currents). Comparison of Figures 1 and 5 shows that the proportion of the southern clade in each population declines as one travels along the direction of both these currents, until no haplotypes of this clade occur in Hokkaido and HCH. This correlation implies corresponding dispersal pathways for the L. brevicula veligers. The HCH population, which is in Honshu where it is less affected by these warm currents, did not include the haplotypes characteristic of the southern clade. We argue below that gene flow by passive dispersal along ocean currents was an important factor for the colonization of L. brevicula on a palaeontological time scale. Surprisingly, contemporary gene flow mediated by currents has not been strong enough to establish uniform or panmictic populations throughout the distribution of L. brevicula; the southern clade has not reached Hokkaido from Kyushu even though it has been 20 kyr since the LGM.

Figure 5.

Distribution of mtDNA ND6 haplotype groups in each sampling locality of Littorina brevicula. LbND6AB03 was the most abundant in all populations, while the southern clade was abundant in southern but not in northern populations.

Figure 5.

Distribution of mtDNA ND6 haplotype groups in each sampling locality of Littorina brevicula. LbND6AB03 was the most abundant in all populations, while the southern clade was abundant in southern but not in northern populations.

AMOVA failed to find east-west genetic differentiation, a result that may be attributable to the remarkable north-south differentiation; the high variation within the east and west groups, which corresponds to a north-south difference, hindered any statistical significance for east-west differentiation in AMOVA, based on the analysis of variance. Additionally, the frequency of the southern clade within each population was similar across Honshu Island, as shown in Figure 5, indicating that east-west differentiation was fundamentally small. In fact, pairwise FST was 0.000 between CHO and OBM, which are at a similar latitude on the east and west coasts of Honshu Island, respectively.

A glacial effect on L. brevicula populations inferred from population genetics

We inferred the glacial effect following the flow chart in Figure 2. First, the presence of a possible glacial effect was evaluated in steps 1–3 and then, if a glacial effect was suggested, the position of glacial refugia and post-glacial colonization was identified in steps 4 and 5. The determinants and resolutions in this chart are based on empirical and theoretical knowledge accumulated in phylogeography as reviewed by Avise (2000), Hewitt (2000) and Freeland (2005). The determinant keys of steps 1 and 2 correspond to gene flow and genetic drift. The IBD model is a powerful tool for detecting continuous genetic-geographic structure based on the balance of gene flow and genetic drift. The key of step 3 corresponds to coalescent theory and sequence mismatch distribution (MMD) analysis detailed by Avise (2000). A comparison of genetic diversity in each population can be a good way to detect the direction of colonization, as the sequential founder effect likely causes lower diversity in newly established populations. The existence of southern refugia is plausible for species in the temperate zone of the northern hemisphere (Hewitt, 2000) and such southern refugia often caused postglacial colonization to proceed northward and likely resulted in higher diversity in southern than in northern populations.

Testing the determinants in the flow chart in Figure 2, the presence of step 1 (genetic difference between local populations) and step 2 (populations with genetic-geographic correlation) were supported by pairwise FST, IBD and AMOVA, and we assumed that the genetic-geographic pattern of L. brevicula was basically shaped by natural evolution and distribution over the palaeontological time scale. Since step 3 (deep genealogy of haplotypes) was absent, we assumed that recent population expansion probably occurred after the LGM. Significant deviation from neutral evolution, indicated by Fu's F and Tajima's D as well as MMD analysis, confirmed recent population expansion. Step 4 (north-south differentiation) was present and although step 5 (higher diversity in the south) could not be proved by statistical analyses, the means of diversity indices are lowest in Hokkaido compared with Honshu and Kyushu.

Though the glacial effect was highly likely, northward postglacial colonization could not be directly supported by genetic evidence. The skewed frequency of the southern clade, however, which appeared as a key factor in the north-south structure, suggested northward colonization. The fact that the highest proportion of the southern clade was found at MKR (which faces the point of east-west separation) and that the proportion of the southern clade decreased in a similar pattern along both the east and west sides of the Japanese Archipelago (Fig. 5) indicated a gradual colonization of the southern clade towards the north in a stepwise manner. Since the northern population of L. brevicula likely declined with a decline in air and seawater temperatures, it can be assumed that a shift in the distribution of L. brevicula occurred during glacial and postglacial periods. Following the decline in temperatures, the northern range of this species shifted southwards and it seems that a bottleneck effect left only a single haplotype, LbN6AB03, around Honshu. An additional haplotype, LbN6OT02, survived in the southern area of the Yellow Sea, where temperatures were higher than around Hokkaido or Honshu. In the present study, genetic diversities in Hokkaido populations tend to be lower than those in Honshu populations, though the differences cannot be statistically significant due to the low number of examined Hokkaido populations. This tendency, plus the direction of the ocean currents and the ecological character of L. brevicula as described below, implies that this species was absent around Hokkaido during the glacial period and recovered during the postglacial colonization.

Time estimation corresponding to population expansion

Population and range expansion was estimated to have occurred at 65.0 and 63.2 kya, long before the end of the LGM at 20 kya. The estimated date might be earlier than the actual expansion time, however, as the divergence rate of 2.46% per Myr, a rate modified in this study from the molecular clock calibration by Quesada et al. (2007), might overestimate the divergence time within L. brevicula. The calibration by Quesada et al. (2007) may not be adequate for the time scale of the present study, as the evolutionary rate of 1.83 ± 0.21% per Myr was based on the split of the gastropod genera Hydrobia and Peringia 5.33 Mya, far earlier than the period in the scope of our study, and such calibration based on speciation is known to overestimate divergence time at the population level within a species (Ho et al., 2005). Thus, for precise dating, molecular clock calibration based on a more recent event should be established for L. brevicula or closely-related Littorina species.

Ecological and palaeontological perspectives supporting the hypothesis of postglacial colonization

The hypothesis of glacial effects and postglacial colonization inferred from phylogeography is also supported by ecological and palaeontological perspectives on L. brevicula around Japan. Fossil records suggest that the speciation event of L. brevicula occurred no later than the upper Pliocene, while molecular phylogeny suggests that the divergence between L. brevicula and L. mandshurica occurred by the lower Pleistocene at the latest (Reid, Dyal & Williams, 2012). The species should therefore have a long history, but the genetic genealogy of the extant populations was shallow, indicating that the population declined in the upper Pleistocene. Probably only two mitochondrial lineages (LbN6AB03 and LbN6OT02) remained after the latest bottleneck event, which was likely associated with glacial cooling. Generally, temperature is arguably the most critical factor determining the geographical distribution of species (e.g. Raffaelli & Hawkins, 1996). Considering the temperature range in which extant populations are currently found, the low temperatures during glacial periods can be expected to have had severe effects on survival and reproduction of L. brevicula around northern Japan. Though the lower temperature limit for the reproduction of L. brevicula has not been reported, it is probably 3–5 °C, given that this species is rare along areas of the Pacific coast under the influence of the cold Kuril Current but abundant on the Sea of Okhotsk coast of Hokkaido (Ohgaki, 1983; personal observations). Since low seawater temperature likely determines the present northern limit of L. brevicula in Hokkaido, it is safe to assume that there were no L. brevicula around Hokkaido during glacial periods and that postglacial colonization gave rise to the extant Hokkaido populations.

We did not find higher genetic diversity in Kyushu, the southern area, where the mean diversity indices in each population were estimated to be lower than in Honshu (Table 1). This might be attributable to population decline around Kyushu in recent years. Although both haplotype and nucleotide diversity at the southernmost sampling point (MKR) was the highest among all examined populations, L. brevicula is not common there at present, suggesting that conditions are no longer optimal and might possibly be too warm for this species. The upper limit for spawning was reported to be 13 °C (Son & Hong, 1998) and high seawater temperatures around Kyushu, ordinarily greater than 15 °C even in January, may also limit the reproduction of L. brevicula and may have caused population decline and a reduction in genetic diversity in this area after postglacial colonization. Since genetic diversity in Kyushu was probably higher in the past (long after the LGM) compared to the present, we assume a southern refugium existed around Kyushu or further south during the LGM and that subsequent northward postglacial colonization resulted in higher genetic diversity in the south than in the north.

A glacial effect was also suggested in Batillaria cumingii (Kojima et al., 2004) and Cellana nigrolineata (Nakano et al., 2010) around Japan, although no evidence of postglacial colonization was provided in those studies. The present study represents the first report of genetic evidence for such colonization in an intertidal mollusc species along the Japanese Archipelago.

SUPPLEMENTARY MATERIAL

Supplementary Material is available at Journal of Molluscan Studies online.

ACKNOWLEDGEMENTS

The authors gratefully acknowledge editors D. G. Reid and E. G. Boulding for the English language review and two anonymous reviewers for their constructive comments and suggestions. K. Imai, Z. Ishizuka, T. Kokita, J. Shoji, H. Takahara and T. Sonoda are thanked for their kind assistance in collecting samples. The present study was partly supported by Grant-in-Aid for Advanced Study Projects from the Tokyo University of Agriculture.

REFERENCES

ALBAINA
N.
,
OLSEN
J.L.
,
COUCEIRO
L.
,
RUIZ
J.M.
,
BARREIRO
R.
2012
.
Recent history of the European Nassarius nitidus (Gastropoda): phylogeographic evidence of glacial refugia and colonization pathways
.
Marine Biology
 ,
159
:
1871
1884
.
AVISE
J.C.
2000
.
Phylogeography: the history and formation of species
 .
Harvard University Press
,
Cambridge
.
AZUMA
N.
,
MIRANDA
R.M.
,
GOSHIMA
S.
,
ABE
S.
2015
.
Phylogeography of Neptune whelk (Neptunea arthritica) suggests sex-biased impact of tributyltin pollution and overfishing around northern Japan
.
Journal of Molluscan Studies
 ,
81
:
131
138
.
CLEMENT
M.
,
POSADA
D.
,
CRANDALL
K.A.
2000
.
TCS: a computer program to estimate gene genealogies
.
Molecular Ecology
 ,
9
:
1657
1659
.
EXCOFFIER
L.
,
LAVAL
G.
,
SCHNEIDER
S.
2005
.
Arlequin v. 3.0: an integrated software package for population genetics data analysis
.
Evolutionary Bioinformatics Online
 ,
1
:
47
50
.
EXCOFFIER
L.
,
SMOUSE
P.E.
,
QUATTRO
J.M.
1992
.
Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data
.
Genetics
 ,
131
:
479
491
.
FREELAND
J.R.
2005
.
Molecular ecology
 .
John Wiley
,
Chichester
.
GOLIKOV
A.N.
1976
.
Gastropoda
. In:
The animals and plants of Peter the Great Bay
 , pp.
79
92
.
Nauka
,
Leningrad
.
(in Russian)
HEWITT
G.M.
2000
.
The genetic legacy of the Quaternary ice ages
.
Nature
 ,
405
:
907
913
.
HO
S.Y.W.
,
PHILLIPS
M.J.
,
COOPER
A.
,
DRUMMOND
A.J.
2005.
Time dependency of molecular rate estimates and systematic overestimation of recent divergence times
.
Molecular Biology and Evolution
 ,
22
:
1561
1568
.
KIM
S.J.
,
RODRIGUEZ-LANETTY
M.
,
SONG
J.I.
2003a
.
Genetic population structure of Littorina brevicula around Korean waters
.
Hydrobiologia
 ,
505
:
41
48
.
KIM
S.J.
,
RODRIGUEZ-LANETTY
M.
,
SUH
J.H.
,
SONG
J.I.
2003b
.
Emergent effects of heavy metal pollution at a population level: Littorina brevicula a study case
.
Marine Pollution Bulletin
 ,
46
:
74
80
.
KOJIMA
Y.
1957
.
On the breeding of periwinkle, Littorina brevicula (Philippi)
.
Bulletin of the Biological Station of Asamushi, Tohoku University
 ,
8
:
59
62
.
KOJIMA
S.
,
HAYASHI
I.
,
KIM
D.
,
IIJIMA
A.
,
FUROTA
T.
2004
.
Phylogeography of an intertidal direct-developing gastropod Batillaria cumingi around the Japanese Islands
.
Marine Ecology Progress Series
 ,
276
:
161
172
.
KYLE
C.J.
,
BOULDING
E.G.
2000
.
Comparative population genetic structure of marine gastropods (Littorina spp.) with and without pelagic larval dispersal
.
Marine Biology
 ,
137
:
835
845
.
LEE
H.J.
,
BOULDING
E.G.
2009
.
Spatial and temporal population genetic structure of four northeastern Pacific littorinid gastropods: the effect of mode of larval development on variation at one mitochondrial and two nuclear DNA markers
.
Molecular Ecology
 ,
18
:
2165
2184
.
MARKO
P.B.
2004
.
‘What's larvae got to do with it?’ Disparate patterns of post-glacial population structure in two benthic marine gastropods with identical dispersal potential
.
Molecular Ecology
 ,
13
:
597
611
.
MARKO
P.B.
,
HOFFMAN
J.M.
,
EMME
S.A.
,
MCGOVERN
T.M.
,
KEEVER
C.
,
COX
L.M.
2010
.
The expansion-contraction model of Pleistocene demography: rocky shores suffer a sea change?
Molecular Ecology
 ,
19
:
146
169
.
NAKANO
T.
,
SASAKI
T.
,
KASE
T.
2010
.
Color polymorphism and historical biogeography in the Japanese patellogastropod limpet Cellana nigrolineata (Reeve) (Patellogastropoda: Nacellidae)
.
Zoological Science
 ,
27
:
811
820
.
NEIVA
J.
,
ASSIS
J.
,
FERNANDES
F.
,
PEARSON
G.A.
,
SERRÃO
E.A.
2014
.
Species models and mitochondrial DNA phylogeography suggest an extensive biogeographical shift in the high-intertidal seaweed Pelvetia canaliculata
.
Journal of Biogeography
 ,
41
:
1137
1148
.
OHGAKI
S.
1983
.
Distribution of the family Littorinidae (Gastropoda) in Hokkaido, with special emphasis on the distribution in Akkeshi Bay
.
Nanki Seibutu
 ,
25
:
173
180
(in Japanese).
OKUTANI
T.
2000
.
Marine Mollusks in Japan
 .
Tokai University Press
,
Kanagawa
.
PALUMBI
S.R.
1994
.
Genetic divergence, reproductive isolation, and marine speciation
.
Annual Review of Ecology and Systematics
 ,
25
:
547
572
.
PANOVA
M.
,
BLAKESLEE
A.M.H.
,
MILLER
A.W.
,
MÄKINEN
T.
,
RUIZ
G.M.
,
JOHANNESSON
K.
,
ANDRE
C.
2011
.
Glacial history of the North Atlantic marine snail, Littorina saxatilis, inferred from distribution of mitochondrial DNA lineages
.
PLoS One
 ,
6
:
e17511
.
PARK
K.S.
,
SONG
J.I.
,
CHOE
B.L.
,
KIM
S.J.
1999
.
Amylase polymorphism of Littorina brevicula from polluted and unpolluted sites, Korea
.
Bulletin of Environmental Contamination and Toxicology
 ,
63
:
633
638
.
QUESADA
H.
,
POSADA
D.
,
CABALLERO
A.
,
MORÁN
P.
,
ROLÁN-ALVAREZ
E.
2007
.
Phylogenetic evidence for multiple sympatric ecological diversification in a marine snail
.
Evolution
 ,
61
:
1600
1612
.
R DEVELOPMENT CORE TEAM
.
2005
.
R: a language and environment for statistical computing
.
R Foundation for Statistical Computing
,
Vienna, Austria
.
ISBN 3-900051-07-0
, .
RAFFAELLI
D.
,
HAWKINS
S.
1996
.
Intertidal ecology
 .
Chapman & Hall
,
London
.
REID
D.G.
1996
.
Systematics and evolution of Littorina
 .
Ray Society
,
London
.
REID
D.G.
,
DYAL
P.
,
WILLIAMS
S.T.
2012
.
A global molecular phylogeny of 147 periwinkle species (Gastropoda, Littorininae)
.
Zoologica Scripta
 ,
41
:
125
136
.
SON
M.H.
,
HONG
S.U.
1998
.
Reproduction of Littorina brevicula in Korean waters
.
Marine Ecology Progress Series
 ,
172
:
215
223
.
TATARENKOV
A.N.
1995
.
Genetic heterogeneity in populations of Littorina brevicula (Philippi) (Mollusca: Gastropoda) in the northern part of Peter the Great Bay (Sea of Japan)
.
Veliger
 ,
38
:
85
91
.
THOMPSON
J.D.
,
GIBSON
T.J.
,
PLEWNIAK
F.
,
JEANMOUGIN
F.
,
HIGGINS
D.G.
1997
.
The ClustalX windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools
.
Nucleic Acids Research
 ,
24
:
4876
4882
.
TSANG
L.M.
,
CHAN
B.K.K.
,
MA
K.Y.
,
CHU
K.H.
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
.
WEERSING
K.
,
TOONEN
R.J.
2009
.
Population genetics, larval dispersal, and connectivity in marine systems
.
Marine Ecology Progress Series
 ,
393
:
1
12
.
WEIR
B.S.
,
COCKERHAM
C.C.
1984
.
Estimating F-statistics for the analysis of population structure
.
Evolution
 ,
38
:
1358
1370
.
WILSON
A.B.
2006
.
Genetic signature of recent glaciation on populations of a near-shore marine fish species (Syngnathus leptorhynchus)
.
Molecular Ecology
 ,
15
:
1857
1871
.
WRIGHT
S.
1943
.
Isolation by distance
.
Genetics
 ,
28
:
114
138
.
ZASLAVSKAYA
N.I.
,
SERGIEVSKY
S.O.
,
TATARENKOV
A.N.
1992
.
Genetic biochemical comparison of Atlantic and Pacific species of Littorina (Gastropoda: Littorinidae)
.
Genetika
 ,
28
:
89
98
(in Russian)
ZASLAVSKAYA
N.I.
,
TAKADA
Y.
1998
.
Allozyme variation and behavioural dimorphism among populations of Littorina brevicula (Philippi) from Japan
.
Hydrobiologia
 ,
378
:
53
57
.