Phylogeography of bivalve Meretrix petechialis in the Northwestern Pacific indicated by mitochondrial and nuclear DNA data

The marine clam Meretrix petechialis is an important economic shellfish species in Northwestern Pacific, but little is known about its phylogeographical pattern. Here, we analyzed 311 samples from 22 locations along the northwestern Pacific using combined profiling of one mitochondrial gene (the first subunit of cytochrome coxidase, COI) and one nuclear DNA marker (the internal transcribed spacer region 1, ITS-1) to investigate contemporary genetic structure and reconstruct phylogenetic history of this species. The results revealed that two distinct phylogeographic lineages dominated marginal seas—the East China Sea (ECS) and the South China Sea (SCS) respectively. The estimation of divergence time between two lineages was 2.1–3.8 Ma, corresponding to a period of the early Pleistocene to late Pliocene. The vicariance of the two lineages was connected to the historical isolation of marginal seas and sea surface temperature (SST) gradient, pointing that SST might play an important role in maintaining phylogeographical patterns of M. petachialis. Significant overlaps between two lineages were observed in 23° to 29° N, located at the adjacent area of the ECS and SCS, which might be promoted by the connectivity of China Coast Current. However, the influence of ocean currents on mixings between two lineages was limited. In comparison, significant relationships were found between genetic distances and geographic distances if the North and South populations were analyzed separately, result of which might be due to some small reciprocal, rotating flows along coastal areas and special geographical conditions.


Introduction
Historical climatic changes in sea level, water temperature and ocean currents can have a powerful impact on the distribution and genetic structure of the marine biota [1,2]. During the Pleistocene, glacial-interglacial climate fluctuations resulted in falling and rising sea levels that alternately exposed and submerged shallow seafloors [3], which in the end might have led to genetic divergence and demographic expansion, respectively [4,5]. As an important marine centre of origin, the Northwest Pacific [6], especially China Seas, have gained much attention 29.3‰ in pelagic larval stage and 9.2 to 36.0‰ in larval period, and was more successful in lower salinity environments [38]. Culture of M. petechialis has been developed rapidly to meet the huge demand of market, and consequently this clam had become one of the most commercially maricultured bivalves in China. In order to ensure sustainable exploitation, management and conservation must rely on population genetic structure and gene flow among populations of the species [39]. M. petechialis was often misidentified as M. meretrix. A previous study reported that M. meretrix was only distributed in the South China Sea, while M. petechialis was more widely distributed throughout the coasts of China [40]. Genetic surveys of M. petechialis populations have been performed in China Seas [41,42], however, sampling locations were limited in these previous studies and no accurate factors that affected the genetic population structure were elucidated.
Mitochondrial DNA is well established to study population genetics and phylogeography due to its comparatively fast rate of evolution and lack of recombination [43]. Because of ease  Table 1. Two lineages are color-coded separately: Red, North lineage; Yellow, South lineage. CCC, China Coastal Current; WKCC, West Korea Coastal Current; SCC, Subei Coastal Current. This map is rendered with ODV v4.7.6 [35]. https://doi.org/10.1371/journal.pone.0183221.g001 Phylogeography of Meretrix petechialis PLOS ONE | https://doi.org/10.1371/journal.pone.0183221 August 16, 2017 of amplification using the polymerase chain reaction method and universal primers [44], COI (the first subunit of cytochrome coxidase) is one of the most frequently used mitochondrial genes and has been widely applied to a number of phylogeographical studies for marine species including Ruditapes philippinarum [32], Cellana toreuma [34], the scyphozoan jellyfish Aurelia spp. [45] and the Rock Shell Thais clavigera [26]. MtDNA markers are maternally inherited, representing only one evolutionary history [46]. Combining nuclear and mitochondrial DNA markers can improve the power of molecular data to test phylogeographic hypotheses and provide complementary information about population structure [25,47]. ITS1 gene (the internal transcribed spacer region 1) is a common nuclear marker used with COI gene in phylogeographic studies [30,31,33]. In the present study, we used COI and ITS1 markers to determine phylogeographical pattern of M. petechialis populations in the northwestern Pacific.
This work aims to (1) understand the contemporary genetic structure of M. petechialis; (2) calculate approximate divergence time between major lineages; (3) elucidate the factors that may affect the genetic population structure.  (Table 1), and were preserved in 95% ethanol. Specific permit was not required to collect M. petechialis from these locations/activities, because M. petechialis is not an endangered or protected species and were only collected from public access areas.

Sample collection
DNA extraction, PCR amplification, cloning, and sequencing Genomic DNA was extracted from the ethanol-preserved adductor muscle using an improved phenol-chloroform procedure, which was described by Li et al [48]. A fragment of COI gene was amplified for all samples, while a fragment of ITS gene was only amplified for partial samples to verify the results of COI (Table 1). Polymerase chain reaction (PCR) was carried out in a 50 μL volume containing 2 U Taq DNA polymerase, 100 ng template DNA, 1 μM each primer, 0.2 mM of each dNTP, 1×PCR buffer, 2 mM MgCl 2 , and 4% DMSO. For a fragment of COI gene amplifications, we used the primers: LCO1490 5'-ATTATTCAGAACCAATCA TAAAGATATTGG-3' and HCO2198 5'-TGTAGGAATAGCAATAATAAAAGTTAC-3' [44]. For a fragment of ITS gene amplifications, we used the primers: PEF-10 5'-TAGAGGAAG GAGAAGTCGTAACAAGG-3' and 5.8R 5'-CAAKRTGCGTTCRARRTGTCGATGWTCA-3' [49]. PCR was performed in a GeneAmp 1 9700 PCR System (Applied Biosystems, Foster City, CA, USA) under the following conditions: an initial denaturation for 3 min at 94˚C, followed by 35 cycles of 40 s at 94˚C, 40 s at 48˚C in COI or 40˚C in ITS, 40s at 72˚C, and completed with a final extension for 5 min at 72˚C.
For COI, PCR products were directly sequenced using respective primers on an ABI 3730 automated sequencer. For ITS, PCR products were purified by DNA gel extraction kit UNIQ-10 following the recommended protocol, then were cloned into cloning vector pEASY-T1 (TaKaRa). For each individual, we selected two to five clones to sequence by M13 primers.

Sequence analysis
Sequences of COI and ITS were assembled and checked by SeqMan in DNASTAR software. The consensus sequences were then aligned with BioEdit [50] using ClustalW [51] under default settings. Sequences were cut into the same length using MEGA 6.0 [52]. Software RDP4 [53] was used to detect genetic recombination by the default parameter settings within ITS sequences. Haplotypes of sequences were determined with the software program DnaSP v5 [54], and sites with gaps/missing in the ITS-1 sequences were considered. All sequences of haplotypes were deposited in GenBank database (Accession numbers: KY318078-KY318174, KY318176-KY318461). The program Arlequin 3.5 [55] was used to analyze molecular diversity indices including haplotype diversity (h), nucleotide diversity (P), and mean number of pairwise differences (k). For the net average genetic distance between major lineages, the Kimura 2-parameter (K2P) model [56], pairwise deletion of gaps data, and "Transitions + Transversions" options were used to calculate in MEGA 6.0 [52].

Phylogenetic analyses
The neighbor-joining (NJ) tree with bootstrap analysis (1000 pseudoreplicates) was constructed using MEGA 6.0, and parameters setting was same as distance analyses. Bayesian inference (BI) was performed in MrBayes 3.2.0 [57], in which the Markov-chain Monte Carlo (MCMC) search was run with four chains for 50 million generations with a sampling frequency of 1/1000 trees. Substitution models were inferred using jModeltest 0.1.1 [58] using the Akaike information criterion (AIC). M. lamarckii (KY318175) was used as outgroup.

Genetic diversity and population structure
Hierarchical analyses of molecular variance (AMOVA) [59] were performed using Arlequin 3.5 on both COI and ITS data set to evaluate population structure. The variance components, sum of squares and F statistics were calculated between lineages, among populations within lineages and within populations, respectively. The significance of F-statistic analogs was evaluated with 10000 random permutations, and the Tamura-Nei model was selected to correct for multiple substitutions because the HKY+I model is not implemented in Arlequin. Haplotype network showing the genetic relationships was generated with PopART 1.7 [60] at the provincial level based on COI gene. Pairwise F ST values among populations within lineages were evaluated using COI gene by Arlequin 3.5, and the statistical significance was tested by 10000 random permutations. Mantel tests to test the relationship between genetic distances and geographic distances were performed on the Isolation by Distance Web Service (IBDWS at http://ibdws.sdsu.edu/~ibdws/) [61] with 10,000 permutations using COI sequences.
To infer if there had been a past population expansion or selective sweep of M. petechialis, Fu's Fs [62] and Tajima's D [63] statistics along with their statistical significance (1000 permutations) were carried out using DnaSP v5. For the purposes of eliminating the possibility of a selective sweep, the McDonald-Kreitman test [64] was conducted in DnaSP v5 using M. lusoria as outgroup species. Because ITS-1 was not a protein-coding gene, only COI sequences were analyzed for the McDonald-Kreitman test. In cases where expansion was evident based on neutrality statistics, statistics in mismatch distributions were performed in DnaSPv5 to test whether a sudden population expansion occurred in the history. Sum of squared deviation (SSD) [65] and Harpending's raggedness index (RI) [66] were calculated in Arlequin 3.5.

Divergence time estimation
Divergence time between the major distinct lineages was estimated in Beast v1.7.5 software [67] only on the more stable COI marker with a calibrated molecular clock and HKY+I models. We supposed an uncorrelated lognormal distribution for the clock and speciation: Yule Process as the tree prior. Analyses were run for 10 million iterations and sampled every 1000 iterations with a final burn-in of 10%. The lack of a fossil record for M. petechialis precludes an estimate of a species-specific molecular clock, which will influence any estimates of divergence time. Molecular clocks have been estimated for the COI gene in several bivalves such as Cyclina sinensis (0.7-2.4%/Myr) [31], Ruditapes philippinarum (0.9-3.35%/Myr) [32], Tegillarca granosa (2-2.4%/Myr) [33] and Atrina pectinata (2.4%/Myr) [68], but the rates vary widely. In our study, we used molluscan specific mtCOI divergence rates calibrated for bivalve (1%/Myr) [19,69] and gastropoda (2%/Myr) [70,71] to estimate the divergence time between two lineages. Although large variance might exist in the divergence rate of 1%/Myr to 2%/Myr for M. petechialis, the rate was used to estimate a rough time for elucidating phylogeographic hypotheses. The mean substitution rates defined in our analyses were obtained by dividing the mean sequence divergence rates by two [72]. Convergence diagnostics were conducted in Tracer v1.5 [73] and the effective sample size (ESS) for each parameter exceeded 200. TreeAnnotater v1.7.5 was used to generate the maximum credibility tree.

Genetic diversity
The aligned COI and ITS-1 gene sequences were obtained with lengths of 747 bp and 860 bp, respectively. The obtained 305 COI gene sequences without indels contained 112 variable sites and 68 parsimony informative sites, yielding 97 haplotypes, of which 72 (74.2%) haplotypes were singletons, being represented by a single sequence in the sample. Wenzhou population showed the highest values of nucleotide diversity (0.01339 ± 0.00701) and mean number of pairwise differences (9.99900 ± 4.70302) because lineage overlap was observed in this population. Overall values of COI gene for haplotype and nucleotide diversities were 0.9483 ± 0.0054 and 0.03364 ± 0.01638, respectively. The obtained 345 ITS-1 gene sequences contained 194 variable sites and 65 parsimony informative sites yielding 286 haplotypes, in which 260 (90.9%) haplotypes were singletons. Overall values of ITS-1 gene for haplotype and nucleotide diversities were 0.9968 ± 0.0009 and 0.00734 ± 0.00386. No genetic recombination was detected within ITS-1 sequences. Overall, values of haplotype diversity, nucleotide diversity and mean number of pairwise differences for ITS-1 gene were higher than those of COI gene (Table 1).

Phylogenetic analyses
Phylogenetic reconstruction with the COI gene based on NJ and BI methods separated M. petechialis into two lineages with high nodal support (Fig 2), which were identified as North and South lineages, each corresponding to one of the biogeographic areas. North lineage included samples from Bohai Sea, Yellow Sea and East China Sea, while samples of South lineage were from East China Sea and South China Sea. The overlap zone was observed in Zhejiang and Fujian provinces. But no obvious difference was found between two lineages in phylogenetic reconstruction with the ITS-1 gene. The net average genetic distance (± SD) between lineages was 5.85 ± 0.95% for COI and 0.08 ± 0.03% for ITS-1.

Genetic diversity and population structure
When populations were grouped by sea basins, for COI gene, hierarchical analyses of AMOVA (Table 2) indicated that the most of the total genetic variation (93.22%) was from differences between two lineages (P < 0.01). A smaller (5.98%, P<0.01) amount of genetic diversity occurred within populations, and the smallest (0.80%), yet significant (P < 0.01), of the total variation was present among populations within lineages. For ITS-1 gene, the result showed that the most of the total genetic variation (84.20%, P<0.01) was apportioned within populations, followed by a smaller (9.00%, P<0.01) amount occurred between groups, and only a small part (6.81%, P<0.01) among populations within groups.
Based on the COI, two haplogroups were formed in the network (Fig 3), which was consistent with the BI and NJ trees. The topology of haplotype network was characterized by multiple star-like type. Significant overlaps between two lineages were observed in 23˚to 29˚N, where is located at the adjacent area of the ECS and SCS.
Pairwise F ST values based on COI sequences among populations within lineages were shown in Tables 3 and 4. For North lineages, Changyi and Weihai populations were found significantly different from others except for Busan. For South lineages, Xuwen, Qishui, Qinzhou and Haiphong populations located in Beibu Gulf, were remarkably different from Wenzhou, Xiamen, Shantou and Zhanjiang populations. The IBD analysis displayed a significant   Table 5). The observed mismatch distributions of two lineages based on COI and ITS data (except for Southern lineage on COI) uniformly displayed unimodal distributions, which fit well with the model of sudden expansion (Fig 5). The RI and SSD for two lineages all matched the null hypothesis of sudden expansion model with nonsignificant values (Table 5). Via an inspection of the network of lineage south (Fig 3), two haplogroups could lead to the bimodal distribution.

Divergence time estimation
Results of the COI divergence time estimations for M. petechialis are shown in Fig 6. The divergence age estimate between North and South lineages was approximately 2.1-3.8 Ma,

Phylogeographical patterns
Based on the network relationships and phylogenetic reconstruction for COI markers, significant population differentiation were observed between northern (ECS) and southern groups (SCS), and overlapped at the adjacent area of the ECS and SCS (Fig 1). In addition, genetic differentiation between two lineages accounted for a large proportion according to the AMOVA analysis. Therefore, we concluded that phylogeographical pattern of M. petechialis was connected to the historical isolation of marginal seas, which is consistent with the observations in  Phylogeography of Meretrix petechialis mitten crab Eriocheir sensu stricto [74], marine fish (Chelon haematocheilus [75] and Bostrychus sinensis [76]) and bivalves (Cyclina sinensis [31] and Tegillarca granosa [33]). The divergence time between two lineages dated back to 2.1-3.8 Ma, generally near to early glaciations intensified time and the start time that sea level major drops 60 to 120 m (2.5 Ma) [8,77]. Because of the sea level falling, ECS and SCS were isolated, so a physical barrier was created for M. petechialis, giving rise to divergent lineages. The divergence time in our study is similar with two subspecies of oyster (2.0-3.6 Ma) [43] and two species in the mitten crab (2.24 Ma) [78]. Above all, in the northwestern Pacific region, historical glaciation during Pleistocene played an effective role in generating intraspecific genetic divergences [25]. In postglacial period, sea level rising and land bridge disappearing contributed to rapid population expansion of intertidal species and secondary contact at the adjacent areas [16,79,80]. The signal of the demographic expansion of M. petechialis was also detected through star-like type of network topology, mismatch distribution analysis and neutrality tests. One major cooling period was occurred in the 2.2-1.0 Ma, which are reflected in the increasing SST gradient in the subtropical area, where SST cooled significantly (decreased by 3.3-5.4˚C and 1.0-2.1˚C in winter and summer respectively), while SST showed little or no cooling (decreased by 0.9˚C and 0.6˚C in winter and summer respectively) in the tropical western Pacific [81]. SST gradient between the ECS located in the subtropical area and the SCS located in the tropical area might lead to genetic divergence combined with glacial isolation. Furthermore, SST possibly plays an important role in maintaining phylogeographical pattern of M. petechialis. Firstly, the boundary of two M. petechialis groups is approximately near to Taiwan Strait. In the northwestern Pacific, the Tropic of Cancer is the boundary between the warm-temperate and tropical regions, and 20˚C isotherm of annual mean temperature running through the Taiwan Strait, which is a major biogeographic boundary for many species [82]. Secondly, analogous to two lineages of Mugil cephalus [17], the distribution range of northern group is consistent with gradients of the cold waters from the Subei Coastal Current (SCC) flowing southward along the ECS coast, while the southern group appeared be restricted to the warm water of the China Coastal Current (CCC) flowing from the SCS into the ECS. Thirdly, because of the temperature difference from the north to south, northern populations have a later reproduction period. These phenomena indicated that SST might limit the  [83,84], which will potentially promote species to migrate northward along Chinese coastline [29].
Similarly, Tegillarca granosa and Cyclina sinensis, plankton larval duration of M. petechialis is about 5-8 days, however, M. petechialis migrates from mid-tidal region to low-tidal region along with tidal currents extending from late May to late June and mid-September to late-September. This specific feature may increase opportunity for gene flow among different populations. During the breeding season of M. petechialis (spring and summer), the China Coast Current with a common velocity of 20 cm/s [34] flowing northward, transported pelagic larvae from SCS to ECS. Overlapping phenomenon was found at the adjacent area of the ECS and SCS (23˚to 29˚N) because of secondary contact, which might be promoted by the connectivity of China Coast Current [25]. Significant relationship was found between genetic distances (F ST ) and geographic distances when the North and South populations were analyzed respectively, after fitting IBD model. In addition, some pairwise F ST values were significant between populations within each lineage. For North lineage, the genetic difference between some Bohai Sea populations (Changyi and Weihai) and East Sea populations were statistically significant. For South lineage, some populations located in Beibu Gulf (Xuwen, Qishui, Qinzhou and Haiphong) were significantly divergent from Wenzhou, Xiamen, Shantou and Zhanjiang populations. This is most likely attributed to the low-dispersal ability of the clam and the existence of small but persistent reciprocal flow and rotating flow along coastal areas, which impede the exposure of larvae to currents and their ability to be transported effectively [31].
He et al.'s [85]study indicates that the South China Sea should be considered to be a colonization origin of Periophthalmus modestus because of an older coalescence time, a putative ancestral haplotype, higher genetic diversity, more proportion of private haplotypes and higher species diversity of the genus Periophthalmus. Similarly, we proposed that South China Sea is the earlier refuge for M. petechialis. First of all, north lineage originated from refugium in the Okinawa Trough about 0.7 Ma, whereas South lineage derived from refugium in South China Sea about 0.8 Ma. Although they are approximate times, the relative coalescence times of two populations are reasonable. Moreover, the portion of private haplotypes from South lineage (75.9%) was higher than that from North lineage (72.1%), which supported the South China Sea was the colonization origin [86]. In addition, haplotype and nucleotide diversity (h: 0.8996, P: 0.00543) of South lineage had higher values compared with Northern lineage (h: 0.8939, P: 0.00274). Generally speaking, older colonized regions are expected to exhibit higher genetic diversity due to the longer evolutionary history [87]. Finally, although the bivalve genus Meretrix Lamarck, 1799 are broadly distributed in the West Pacific and the Indian Ocean, M. petechialis is the only Meretrix species in the ECS regions.

The discrepancy between mitochondrial and nuclear topologies
MtCOI gene analysis separated M. petechialis into two lineages, while nrITS data showed a slightly different topology. The net average genetic distance (± SD) between lineages was 5.85 ± 0.95% for COI, but only 0.08 ± 0.03% for ITS. Result of AMOVA also displayed that nrITS were not significantly differentiated between groups. This inconsistent pattern has also been found in black-throated tits Aegithalos concinnus [88] and the pen shell Atrina pectinata [89]. A large proportion of cases that discordance between mitochondrial and nuclear data likely appeared following geographic isolation and secondary contact, which attributed to incomplete lineage sorting and post-glacial introgression, respectively [47]. However, they are difficult to be distinguished [90]. COI gene will have a four times faster coalescence time than that of ITS gene [1], because this rate is inversely proportional to the effective population size [91] (a fourfold smaller effective population size [92]). In our study, the divergence time between two lineages dates back to a relatively young divergence (Pleistocene). It is reasonable that ITS gene keep homogeneous and will have longer time to spit, while COI gene have achieved coalescent. Furthermore, under incomplete lineage sorting, no discernible biogeographic pattern among populations is expected to be found [91]. In our study, phylogenetic reconstruction with the COI gene separated M. petechialis into two lineages, but no obvious difference and discernible geographic pattern was found between the two lineages using the ITS gene. We also found samples from different mtDNA lineages shared the same haplotype for ITS gene, for example, samples from Busan, Beihai and Caotan shared the same haplotype. Therefore, incomplete lineage sorting is more likely to be the reason of the discepancy between mitochondrial and nuclear topologies.

Conclusion
To better elucidate phylogeographic pattern of M. petechialis, we sequenced the COI and ITS1 genes from 311 individuals sampled from 22 locations along the coastal areas of China, Korea and Vietnam. The results provided evidence for strong genetic divergence between two lineages, which might be attributed to variance of sea level changes and SST gradients. Here, the combinatory use of nuclear ITS gene and the mitochondrial COI gene helped elucidate the phylogeographical pattern of M. petechialis, such as population structure and population expansion analyses. With the objectives of getting more detailed information about phylogeography of M. petechialis with higher genetic resolution, analyses with multiple nuclear genes would be conducted.