Life-history features and oceanography drive phylogeographic patterns of the chiton Acanthochitona cf. rubrolineata (Lischke, 1873) in the northwestern Pacific

Sampling, DNA extraction, PCR and sequencing

A total of 180 A. cf. rubrolineata individuals were sampled from 11 localities in Korea, Japan, and China, and about 20 specimens per population except for the Japanese populations were sequenced for mitochondrial COI and 16S (Fig. 1 and Table 1). Collected samples were fixed and kept in 95% ethanol until DNA extraction. Genomic DNA was extracted using the E.Z.N.A. Mollusc DNA Kit (OMEGA Bio-tek, Norcross, GA, USA) following the manufacturer’s protocol. Two fragments of mitochondrial COI and 16S were amplified using the universal primer sets of LCO1490 (5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′) and HCO2198 (5′-TAA ACT TCA GGG TGA CCA AAA AAT CA-3′) (Folmer et al., 1994), and 16Sar (5′-CGC CTG TTT ATC AAA AAC AT-3′) and 16Sbr (5′-CCG GTC TGA ACT CAG ATC ACG T-3′) (16S; Palumbi, 1996). The polymerase chain reaction (PCR) was performed in a 50 µl volume containing 1.5 µl of template DNA, 5 µl of 10x Ex Taq buffer (Takara, Shiga, Japan), 4 µl of 0.2 mM dNTP mixture, 1 µl of each primer (10 pmole), and 0.25 µl Ex Taq Polymerase. The amplification conditions were initial denaturation at 95 °C for 1 min, followed by 40 cycles of denaturation at 94 °C for 30 s, annealing at 46 °C for 30 s for COI and 52 °C for 30 s for 16S, extension at 72 °C for 30 s, and a final extension at 72 °C for 10 min. The PCR product was purified with a QIAquick Gel Extraction Kit (Qiagen, Valencia, CA, USA) and then sequenced using an ABI PRISM 3730xl DNA analyzer (Applied Biosystems, Foster City, CA, USA) in both directions.

Sequence analyses

Sequences were assembled and edited in Geneious (Kearse et al., 2012) and then deposited in GenBank (accession numbers: 16S: MN205571–MN205750; COI: MN205751–MN205930). We merged the two mitochondrial segments for each individual using the software FasParser (Sun, 2017), and performed a partition homogeneity test in PAUP* 4.0 b10 (Swofford, 2002) to check the congruence between COI and 16S sequences with the heuristic search option (number of replicates = 100). The result (P-value = 0.680) suggested there was no incongruence between the two segments, and therefore the concatenated sequences were used in subsequent analyses. Haplotypes were defined using DnaSP v5 (Librado & Rozas, 2009), and their relationships were inferred using the TCS network in Popart v1.7 (Leigh & Bryant, 2015). The software jModelTest 2 (Darriba et al., 2012) was used to determine the best-fit model for sequence evolution. HKY+I was selected as the best model based on Bayesian Information Criterion and used in the following sets. Genetic diversity indices, including haplotype diversity (Hd), nucleotide diversity (Nd), and mean number of pairwise differences (Pd) were calculated for each of the 11 populations in ARLEQUIN v3.5 (Excoffier & Lischer, 2010).

For population structure analyses, population clusters were first estimated for 10 populations (excluding the Nomomachi population in Japan, because only three individuals were analyzed; see Table 1) based on F_CT values using SAMOVA 2.0 (Dupanloup, Schneider & Excoffier, 2002). The test was performed with K-values of 2-7, and the value for which the F_CT was highest was chosen as the best grouping. Based on the cluster identified from the SAMOVA, analysis of molecular variance (AMOVA) was then conducted with 10,000 permutations in ARLEQUIN v3.5 to estimate the partitioning of genetic variation. Since the HKY model was not available in ARLEQUIN, the most similar model, that of Tamura-Nei (Tamura & Nei, 1993) was used to estimate population structure. Under the same model, pairwise population comparisons (F_ST values) were calculated with 10,000 permutations in ARLEQUIN v3.5.

We used the coalescent-based program MIGRATE-N v. 3.6 (Beerli & Felsenstein, 1999; Beerli & Palczewski, 2010) to estimate mutation-scaled migration rate (M) among the groups defined in SAMOVA and effective population size (Θ) of each group based on maximum-likelihood estimates. Several short runs were performed to check for convergence of chains, and appropriate upper bounds were set for each parameter after the test runs. Three long chains were run with 2,000,000 steps and sampled every 100 steps, and the first 20,000 steps were discarded as burn-in. A static heating scheme was set at four temperatures of 1, 1.5, 3, and 100,000. Five replicates were run to check the consistency of the estimates, and stationarity of the Markov chain Monte Carlo (MCMC) was assessed by examining the effective sample size (ESS). We transformed the migration rate to the number of migrants per generation (Nm) using the formula Nm = Θ*M/x (where x is the inheritance parameter, and here x = 1 for mitochondrial gene) (Beerli & Palczewski, 2010).

Neutrality tests were applied to each group defined in SAMOVA and the entire population (including all the sampled individuals) to infer historical demography. Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) tests were conducted with 10,000 replications in ARLEQUIN v3.5. In addition, a mismatch distribution analysis was performed using DnaSP v5 to test the signal of historical expansion. Bayesian skyline plot (BSP) analysis was performed for the COI data set using BEAST v2.4.8 (Bouckaert et al., 2014) to infer historical demography of effective population size for the entire population and individual groups as defined in SAMOVA. We selected the HKY model with four rate categories of heterogeneity as the site evolution model. For the clock model, we used a log-normal relaxed clock. There is no reliable molecular clock yet for mitochondrial COI sequences of this chiton species. However, a molecular clock for the mtDNA COI of other two Acanthochitona species was estimated, ranging from 1.0 to 2.0% per million years (myr) based on the minimum splitting time between Caribbean and Eastern Pacific presumed geminate species pairs of the Isthmus of Panama (Ríos et al., 2014). Here we applied the lower bound of 1.0%/myr for the COI sequences and a generation time of two years (Avila-Poveda & Abadia-Chanona, 2013; Cherns, 1999) to convert the parameters to real demographic times. As the priors, the Coalescent Bayesian Skyline model was chosen while setting the free standard deviation (ucldStdev) parameter to Exponential. The analysis was performed with an MCMC chain length of 20 million iterations for the dataset including all individuals, and 10 million for each group, sampling every 1,000 generation to estimate the effective population size change. Tracer 1.6 was used to check and visualize the results (available from http://beast.bio.ed.ac.uk/Tracer). The 16S data set was not used here because it is uninformative and there is no reported molecular clock yet for this species or any other Acanthochitona species.

Significant population structure among three regional groups

A most striking pattern in A. cf. rubrolineata is the existence of significant population structure among the populations suggested by the various analyses. The best grouping method in SAMOVA divided all study populations into three groups, which geographically correspond to the coastal regions of Korea, China, and Japan (Table 2). Pairwise F_ST estimates also showed that all significant values came from the comparison of populations between different groups (Table 3). Estimations in MIGRATE-N revealed asymmetric gene flow among the three groups: among the estimations of six possible directions of gene flow, only the directions from China to the other two regional groups were significant (Fig. 3). These results indicate that genetic connectivity (i.e., gene flow) is limited or at least partially impeded among populations of the three regions, contributing to significant population structure among them. A comprehensive review of the correlation between pelagic larval duration and dispersal distance of marine species has also shown that species with a larval stage of less than a week usually disperse a short span of 1–10 km (Shanks, 2009). Larval behavior (especially larval stage duration time) has been shown to play a critical role in determining dispersal distance for many marine species (Cowen & Sponaugle, 2009; Hellberg et al., 2002; Pineda, Hare & Sponaugle, 2007). Thus, it is not surprising to observe here significant population structure in A. cf. rubrolineata populations, considering its limited dispersal capacity (inferred from the life-history traits of many chiton species). Fertilized chiton eggs are known to hatch out after some period that mostly depends on temperature into trochophore larvae (ciliated, free-swimming larvae) and then become increasingly benthic and soon metamorphose into a juvenile form within a few days to about 12 days after fertilization depending on chiton species (Pearse, 1979). From the developmental study of a congeneric species (A. crinitus) that undergoes metamorphosis within 3-4 days after fertilization (Fritsch et al., 2015), A. cf. rubrolineata is assumed to have a short planktonic larval duration (albeit variable depending on temperature), likely preventing long-distance dispersal (i.e., gene flow) of A. cf. rubrolineata populations across our three studied regions. A short larval stage would likely affect population structure among the three regions, especially considering their geographic separation and habitat discontinuity, and the general lack of stepping stone habitats between them.

Contribution of the warm current system to population structure among the regional groups

Although the magnitude of gene flow during the planktonic larval stage seems limited for chitons (Ríos et al., 2014; Todt et al., 2008), they can still be transported via rafting on floating macroalgae or terrestrial vegetation like trees (Eernisse, Draeger & Pilgrim, 2018). A long-distance dispersal (up to 600 km) by kelp rafting has been reported for two species of subantarctic chitons that dispersed their propagules from subantarctic islands to New Zealand (Fraser, Nikula & Waters, 2010). However, dispersal by rafting in A. cf. rubrolineata may be limited by the warm current system in the YS. The trifurcated branches of the Kuroshio Current (including the Yellow Sea Warm Current, the Tsushima Warm Current, and Kuroshio Extension) flow northeastward into their downstream pathways that separate the three regional groups. This oceanographic surface flow may act as a physical barrier that, albeit not completely, limits the opportunity for gene flow via larval dispersal and/or occasional rafting of the individuals across the three regions. Note that, as inferred from the MIGRATE-N results, this barrier effect may be weak for the directions from China to both Korea and Japan. Similar population structuring resulting from hydrographic patterns has also been reported in other marine species in the NWP (Ni et al., 2017) and other ocean systems, including the snakeskin chiton Sypharochiton pelliserpentis in New Zealand (Veale & Lavery, 2011) and coral reef fishes in the Caribbean (Taylor & Hellberg, 2006).

Genetic homogeneity within the regional groups

In contrast to inter-regional genetic structuring proposed to be driven by the three branches of the Kuroshio Current, substantial genetic homogeneity was revealed within each regional group, with no significant population structure detected (Table 3). Genetic connectivity in each group is assumed to be facilitated by the present-day surface currents along the coastline of each region (Fig. 1). These coastal currents may play a significant role in maintaining genetic connectivity that contrasts with the warm currents: they are expected to enhance gene flow by transporting planktonic larvae and/or rafting adults among local populations, reducing genetic differentiation among the within-group populations. Previous studies have revealed a close correlation between the coastal current system and population homogeneity in many organisms in the YS (e.g., Li et al., 2010; Ni et al., 2015b; Xue et al., 2014). Relatively short distances between populations and habitat continuity within each region are also likely potential factors underlying their genetic homogeneity.

Demographic history and the impact of Last Glacial Period

In estimating the historical population demography of A. cf. rubrolineata, results both from the neutrality test and BSP analysis supported a sudden historical population expansion. By reconstructing effective population size change in BSP, the estimated timing of population expansion for the entire population and Korea group were calculated to be 60–70 kya (the entire population) and 100 kya (Korea group) (Figs. 4A, 4B), after and/or approximately corresponding to the onset of Last Glacial Period (LGP, 110 to 15 kya), the most recent glacial period within the Quaternary Ice Age (Severinghaus & Brook, 1999). Abrupt climate change and subsequent sea level fluctuation during the LGP is believed to have greatly influenced the demographic history of contemporary marine species worldwide (Hewitt, 2004). Considering the topography of the YS, which is characterized by a flat, shallow, and extensive continental shelf (with an average depth of 55 m), population size changes resulting from the LGP are assumed to have been very significant for marine species in this region (reviewed in Ni et al., 2014 and references therein). The China group showed a significant population expansion more recently (after ∼25 kya), likely reflecting the impact of the Last Glacial Maximum about 19–21 kya (Yokoyama et al., 2000). Different population demography, however, was observed for the Japan group, which showed a gradual increase since 550 kya, which is inconsistent with all other groups. We consider this result may be an artifact from the small number of individuals (3–7 individuals/population) included in the analysis, or else this population has its own historical demography. Note that due to the absence of a reliable molecular clock for mitochondrial genes of A. cf. rubrolineata, the estimated population expansion times of different groups in this study remain speculative. In addition, the expansion times could be more recent if a molecular clock several times faster was applied under the hypothesis of time-dependent molecular rate (Ho et al., 2011). Population genetic analysis of this study is based on a single locus of mtDNA sequence data, and therefore further analysis using multi-locus data from more extensive sampling of the species (particularly additional sampling from Japanese populations) would be required to resolve this uncertainty.