Genetic Population Structure of Macridiscus multifarius (Mollusca: Bivalvia) on the Basis of Mitochondrial Markers: Strong Population Structure in a Species with a Short Planktonic Larval Stage

The clam Macridiscus multifarius with a planktonic larval stage of about 10 days is an ecologically and economically important species in the coastal regions of China. In this study, 3 mt-DNA markers (COI, 12S rRNA, and ND1) were used to investigate the population structure and demography of wild M. multifarius populations in 3 coastal localities of the East China Sea (ZS and ZP populations) and Beibu Gulf in the South China Sea (BH population). Sequences of 685 bp in COI, 350 bp in 12S rRNA, and 496 bp in ND1 were determined. High level and significant F ST values were obtained among the different localities on the basis of either COI (F ST = 0.100–0.444, p < 0.05) or 12S rRNA (F ST = 0.199–0.742, p < 0.05) gene, indicating a high degree of genetic differentiation among the populations. F ST values were significant but weak for the ND1 gene because it is highly conservative. The median-joining network suggested an obvious genetic differentiation between ZS and BH populations, and the finding is consistent with the results of our demographic analyses using the unweighted pair group method with arithmetic mean. Our study unraveled the extant population genetic structure of M. multifarius and explained the strong population structure of a species with a short planktonic larval stage species; this information could be useful for fishery management measures, including artificial breeding and conservation.


Introduction
The life history of most marine organisms includes a planktonic stage during which larvae disperse to distances that range from several meters to hundreds of kilometers from their location of release [1]. In general, species with a long planktonic larval stage are capable of moving great distances and dispersing widely via ocean currents. The population genetic theory predicts that a longer planktonic larval stage would result in increased gene flow, and consequently, decreased levels of population differentiation. For example, a solitary coral species with a brief planktonic larval stage showed a stronger population genetic structure than another coral species with a longer larval stage [2]. Shanks et al. [3] found that the mollusc Haliotis rubra, with a 6-day pelagic larval phase, has a dispersal distance of less than 15 km, that Ensis directus, with a 16-day pelagic larval phase, can disperse to about 111 km, and that Perna perna, with a 15 to 20-day pelagic larval phase, can disperse to about 235 km. Shanks et al. observed that propagules that spent in the water column dispersed further [3]. However, Galarza et al. [4] showed that the pelagic larval phase often fails to achieve complete dispersal potential. Thus, the relationship between the pelagic larval phase and gene flow in planktonic developers may be quite complex and needs to be studied further [4].
The clam Macridiscus multifarius L. F. Kong, Matsukuma & Lutaenko, 2012 [5], also called the sandy clam, belongs to the family Veneridae (Mollusca; Bivalvia; Veneridae; Macridiscus) and is known to be native to the coastal waters of the Western Pacific region, which extends from the coast of Japan to Northwestern Australia [6]. In China, M. multifarius is found along the coast, and it has a shell length of 30-40 mm. It has a brief pelagic phase of about 10 days at 24-29°C and 20-25 ppt of salinity. Spawning of M. multifarius occurs in summer and peaks in July and August [6]. The adults are benthic and relatively immobile, so patterns of larval dispersal and recruitment help in understanding the population connectivity of M. multifarius. This clam is an important and commercially exploited species; however, M. multifarius populations have declined largely because of over-exploitation and habitat destruction [6]. A previous study showed that the shell length of most M. multifarius specimens collected from Nanji Island, Zhejiang Province, was less than 15 mm [6]. The decline in M. multifarius populations has caused the price to increase, for example, in Zhoushan City, Zhejiang Province, the price has increased from ¥20 to more than ¥100 per kilogram in a few years [7]. Therefore, there is a great need to conserve and manage M. multifarius populations.
Enhancement practices that involve the release of M. multifarius juveniles into natural coastal areas were initiated by the Zhejiang Marine Aquaculture Research Institute in 2011 (about 6,156,000 individuals) and 2012 (about 2,155,000 individuals). M. multifarius is also bred artificially; however, very little is known about the population dynamics. Management and conservation efforts could be improved by examining how environmental factors influence population connectivity and patterns of population genetic structure [8]. Population genetic structure reflects the evolutionary history as well as the evolutionary potential of a species. Lack of information on the population genetic structure may lead to over-exploitation and subsequent resource collapse. The population genetic structure of marine organisms has been partially explained by biophysical factors such as the biology, ecology, and behavior of the species as well as hydrographical barriers to dispersal [9,10].
Genetic markers are useful tools for measuring genetic variation and gene flow among populations [11]. In the present study, we used the gene sequences of partial mitochondrial cytochrome c oxidase subunit I (COI), ribosomal 12S subunit (12S rRNA), and dehydrogenase subunit 1 (ND1) to assess the genetic diversity and population structure of natural M. multifarius populations in China. The objective of this study was to test the influence of planktonic larval stage duration on the genetic structure of M. multifarius. Our approach could add in the design of spatial fishery management and conservation strategies for species that inhabit the East China Sea and South China Sea.

Study area and sampling
All experiments and animal sampling were reviewed and approved by the State Oceanic Administration of China and the Ethics Committee of Zhejiang Ocean University and performed according to national laws and regulations. Wild adult specimens of M. multifarius were collected from 3 coastal localities in the East China Sea (Zhoushan [ZS], 29°35 0 N, 122°25 0 E, Zhejiang Province and Zhangpu [ZP] 23°32 0 N, 117°58 0 E, Fujian Province) and Beibu Gulf (Beihai [BH], 20°54 0 N, 109°47 0 E, Guangxi Province). These sites do not belong to a national park or a protected sea area or a relevant regulatory body concerned with wildlife protection or a private owner. We confirmed that the field studies did not involve endangered or protected species. Geographic locations and sample sizes of all the examined populations are provided in Fig 1 and Table 1. All the samples were collected in October 2014. Tissues from the adductor muscle were dissected from fresh specimens, preserved in 95% ethanol, and frozen at -20°C until DNA extraction.

DNA extraction and PCR amplification
Total genomic DNA was extracted from the adductor muscle of M. multifarius by using the salt-extraction procedure [12] with slight modifications, and the tissues were pretreated with Proteinase K. The extracted DNA was stored in TE buffer at 4°C.
For all markers, PCR was performed using the Taq Master Mix (ComWin Biotech Co., Ltd., Beijing, China). The PCR primer sequences have been listed in Table 2. Amplification was

Data analysis
For all sequence analyses (COI, 12S rRNA, and ND1), genetic similarities were evaluated using BLAST (http://www.ncbi.nlm.nih.gov/BLAST) to identify M. multifarius sequences. Then, the sequences were aligned using Clustal W [15], and individual consensus sequences were retrieved using both alignment and manual checks. The aligned DNA sequences were imported into MEGA version 5.0 [16] for sequence comparisons and variation analysis. Standard genetic diversity indices, such as the number of haplotypes, polymorphic sites, haplotype diversity (h), and nucleotide diversity (p), were calculated with DnaSP 4.0 [17]. Tajima's D [18] and Fu's Fs [19] neutrality tests were performed using Arlequin version 3.11 [20] and 10,000 random permutations to infer population expansion events and check for deviations from a strictly neutral model of evolution. Phylogenetic trees were constructed using the unweighted pair group method with arithmetic mean (UPGMA) [21] based on a matrix of the Kimura two-parameters (K2P) distance method [22] in MEGA 5.0 [16]. Statistical robustness in the nodes of the resulting tree was determined by 1000 bootstrap replicates [23]. To depict phylogenetic and geographical relationships of the haplotypic sequences, 3 haplotype networks of COI, 12S rRNA, and ND1 were created using the median-joining method in Network 4.6 software [24]. Population structures based on COI, 12S rRNA, and ND1 genes were investigated using analysis of molecular variance (AMOVA) [25] in Arlequin version 3.11 [20]. Genetic differentiation, genetic distance, and migration rate among the populations were estimated by calculating the F statistic (F ST ) between the populations and testing their significance with 1000 permutations.

Genetic variation
Sequences of the 685-bp COI gene were determined in 61 specimens, and 33 polymorphic sites and 29 haplotypes were detected. Twenty-five haplotypes were found in only one population, 2 (H1 and H18) were found to be shared by two populations, and 2 (H11 and H13) were found to be shared by all three populations (Table 3). A 350-bp fragment of 12S rRNA gene was sequenced from 61 samples, and 4 polymorphic sites and 7 haplotypes were detected. Two haplotypes (H1 and H2) were found to be shared by ZS and ZP populations; the rest were population specific haplotypes (Table 4). In addition, 496 bp of the ND1 gene was sequenced using 65 specimens, and 5 polymorphic sites and 6 haplotypes were detected. H1 and H2 were shared among all 3 populations; H3 was observed in ZS and BH populations; H4 and H5 were specific for the BH population (Table 5). Haplotype diversity (h) and nucleotide diversity (π) based on COI, 12S rRNA, and ND1 genes for the 3 populations are listed in Table 1. Not all of Tajima's D results for the 3 genes were significant (p > 0.05; Table 6), indicating balanced selection and/or a decrease in population size [18]. Most of Fu's F S results for the 3 genes were not significant (p > 0.05; Table 6), except F S results for the COI gene in the ZS *means that the nucleotide is the same as the nucleotide of H1. doi:10.1371/journal.pone.0146260.t003

Haplotype network analysis
The median joining network (Fig 2) *means that the nucleotide is the same as the nucleotide of H1.
doi:10.1371/journal.pone.0146260.t005 between ZS and BH populations, while the N m value was < 1. For ND1 gene data, the genetic structure was weak but significant and no differentiation was detected among the populations (F ST < 0.05). AMOVA for COI, 12S rRNA, and ND1 genes on the basis of haplotype frequencies revealed that 75.56%, 41.98%, and 101.1%, respectively, of the genetic variation occurred within the populations, whereas 24.44%, 58.02%, and -1.1%, respectively, of the genetic variation occurred among the populations (Table 8). AMOVA results for overall population genetic structure within and among the populations were highly significant (p < 0.001), except for ND1 gene results.  The phylogenetic relationship of M. multifarius among the haplotypes was determined using COI, 12S rRNA, and ND1 genes of Venerupis philippinarum (NC 003354.1), Paphia undulata (NC 016891.1), and Paphia amabilis (NC 016889.1) as the out-groups (Fig 3). The 3 UPGMA trees showed that most of the haplotypes were weakly associated (less than 50% bootstrap support) or unresolved, which was possibly due to low nucleotide differences among them. Haplotypes of the 3 genes were clustered into 2 obvious branches. Furthermore, the 3 UPGMA trees (COI, 12S rRNA, and ND1 gene) consistently showed the same results: M. multifarius clustered (bootstrap 98, 96, and 99, respectively) with V. philippinarum and then clustered with the branch comprising P. undulata and P. amabilis.

Discussion
In this study, significant genetic structure was observed in M. multifarius populations by using 3 mtDNA markers. This result is consistent with the finding that species with a short planktonic larval stage have a short dispersal distance, especially those with low mobility as adults. Thus, marine bivalves with a short planktonic larval stage could have a strong population genetic structure.
Genetic diversity within and between populations provide a potential genetic resource for future adaptation, and can be critical for the fitness of a population [14]. It is mainly explained by several historical and contemporary processes, such as genetic drift, effective migration, natural selection, fragmentation and range expansion [26]. Significant molecular differentiation in M. multifarius populations was identified on the basis of all 3 mtDNA markers. Results of the COI gene showed high levels of haplotype diversity (h = 0.822-0.889) and low levels of nucleotide diversity (π = 0.004-0.006). This pattern is consistent with that observed in previous studies on other species from the East China Sea or South China Sea, for example, Mytilus coruscus (h = 0.818-0.972 and π = 0.00405-0.00747) [27] and the scallop Chlamys nobilis (h = 0.641-0.966 and π = 0.00412-0.01306) [28]. For 12S rRNA (h = 0.257-0.678 and π = 0.001-0.002) and ND1 (h = 0.199-0.391 and π = 0.001-0.003) genes, a lower level of genetic variation and fewer haplotypes were obtained (7 haplotypes for 12S rRNA and 5 haplotypes for ND1), largely because they are in a low-variation region of the mitochondrial genome. In our study, the highest level of haplotype diversity was observed in the BH population. This result may be because the bio-resource of BH is more abundant in the southern tropical regions.
Researchers often use F ST to assess gene flow, and a higher F ST value indicates a lower level of gene flow (N m ) and higher genetic differentiation among populations [29]. F ST reflects the level of inbreeding within populations [30] or the extent to which populations are differentiated [31]. The presence of genetic structure is an outcome of limited gene flow and a high level of genetic drift within each reproductively isolated group. F ST values below 0.05 indicate negligible genetic differentiation, whereas values greater than 0.25 indicate high genetic differentiation within the analyzed population [32]. High and significant F ST values were obtained among different localities of M. multifarius populations on the basis of both COI (F ST = 0.100-0.444, p < 0.05) and 12S rRNA (F ST = 0.200-0.742, p < 0.05) genes; this, indicated a high degree of genetic differentiation among the M. multifarius populations. In addition, the highest F ST values for both COI and 12S rRNA genes were between ZS and BH populations, which endorsed a signal of isolation due to distance. Populations may be divided by major geographic barriers such as land barriers and oceanographic patterns. Our result also showed a higher level of variation when compared with those of other bivalve studies, for instance, Li et al. [27] characterized the genetic relationship of the mussel M. coruscus in the East China Sea by using the COI gene (F ST = -0.03612-0.09774), and it showed no significant phylogeographic structure. The planktonic larval stage of M. coruscus was 35 days [33], and it potentially disperses over larger distances. Similarly, the F ST value using the COI gene among populations of the mussel Mytilus galloprovincialis along the coast of Eastern and Southern China was between -0.0469 and 0.0509 [34]. In contrast, Wang [35] studied the genetic diversity of the clam Coelomactra antiquata (planktonic larval stage of only about 9 days) [36] on the basis of the COI gene, and found that the F ST value between populations from the Bohai Sea and East China Sea was 0.95899, which showed a highly significant genetic differentiation between populations. These results endorsed a highly discussed hypothesis that marine species with an intensive planktonic larval stage would exhibit a high level of genetic differentiation. Nevertheless, with respect to the ND1 data, F ST values were significant but weak (F ST < 0.05), and high N m values were detected among the populations; this may be because the ND1 gene is primitive and high conservative in this species.
Results of the median-joining network suggested an obvious genetic differentiation between ZS and BH populations with respect to the COI and 12S rRNA genes, while results for the ND1 gene showed several ambiguous connections among the 3 populations. Similarly, demographic analyses using UPGMA for all haplotypes show 2 similar branches: one branch contains most (COI) or all (12S rRNA and ND1) haplotypes from the BH population, and the other branch contains most haplotypes from the BH population (COI) or haplotypes from all 3 populations. A much higher degree of genetic differentiation was observed between BH population and the other 2 populations, mostly because BH is located in the Beibu Gulf, South China Sea, and the Leizhou Peninsula (Fig 1) limits the dispersal of planktonic larvae and acts as a barrier to the genetic connectivity of this species. It has been suggested that, in summer, the northward China Coast Current transports propagules from ZP to ZS, although the ZP population is also influenced by the Taiwan Warm Current.
In conclusion, our results showed that differentiations occurred among 3 populations of M. multifarius on the basis of COI, 12S rRNA, and ND1 gene data; in addition, the BH population was isolated from the other populations. M. multifarius has a short planktonic larval stage, but it has a high population genetic structure that may facilitate fishery management, captive breeding, and restocking programs for aquaculture and conservation.