Genetic Diversity and Structure Revealed by Genomic Microsatellite Markers of Mytilus unguiculatus in the Coast of China Sea

Simple Summary Microsatellite markers are widely used in genetic breeding and population genetic structure analysis. In this study, microsatellite markers are used to analyze the genetic diversity and population genetic structure of Mytilus unguiculatus in seven coastal areas of China. The results showed that M. unguiculatus has high genetic diversity with genetic structure differences observed between the Qingdao population and the other six populations. The genetic structure of M. unguiculatu is relatively weak. These findings provide a molecular biological basis for the rational development and protection of wild germplasm resources of Mytilus unguiculatus in China and can serve as a data reference for the formulation of reasonable breeding programs. Abstract The hard-shelled mussel Mytilus unguiculatus plays an important role in mussel aquaculture in China due to its characteristic and nutritive value. In this study, 10 microsatellite loci are used to study the genetic diversity and genetic structure of seven location populations of M. unguiculatus in coastal areas of China. The results of amplification and genotyping show that the observed heterozygosity (Ho) and the expected heterozygosity (He) are 0.61~0.71 and 0.72~0.83, respectively. M. unguiculatus has high genetic diversity. The inbreeding index (FIS) of M. unguiculatus is significantly positive (FIS: 0.14~0.19), indicating that inbreeding might exist within populations. The genetic structure of M. unguiculatus is weak within populations from the East China Sea All results showed that genetic differences existed between the Qingdao population from the Yellow Sea and other populations from the East China Sea. It does not detect a population bottleneck event or expansion event in the populations. The results from this study can be used to provide important insights in genetic management units and sustainable utilization of M. unguiculatus resources and provide a better understand of genetic structure of marine bivalve with similar planktonic larval stage in the China Sea.


Introduction
Mytilus unguiculatus (Mollusca; Bivalvia; Mytiloida; Mytilidae; Mytilus) is an offshore warm-temperate benthic shellfish. The shell of M. unguiculatus is tip and thin, and the shell surface is often brown. In East Asia, it is widely distributed in the littoral of the Japan Sea, Bohai Sea, Yellow Sea, East China Sea, and South China Sea [1,2]. In China, it has a long breeding history, and the main breeding area is in Zhoushan City, Zhejiang Province, where the number of breeding currently exceeds 6 billion grains per year [3,4]. Due to its high nutritional value and delicious taste, M. unguiculatus is widely consumed by people [5]. M. unguiculatus is one of the important aquaculture shellfish species in China. With the increasing demand of people, M. unguiculatus have become one of the key fishing targets. However, the generation replacement rate of M. unguiculatus is slow, and large-scale and high-intensity harvesting is bound to reduce mussel resources [6,7]. Therefore, studying the genetic structure of wild mussels can provide valuable information for improving the quality of provenances and cultured mussels, which can effectively protect resources of genetic variation in wild populations and improve the yield of M. unguiculatus.
Microsatellite DNA (Simple Sequence Repeats, SSR) is a molecular marker that is widely distributed in the genome of eucaryon. It consists of a repetitive DNA sequence with a certain number of base pairs (1 to 6 bp) [8,9]. Microsatellite DNA has a higher mutation rate than other regions of DNA [10]. There are different repeats of DNA motifs between individuals within species because of the influence of the environment or other factors. Therefore, this characteristic makes significant variation between individuals and populations, which means the microsatellite DNA is one kind of high polymorphism molecular marker [11]. Additionally, microsatellite DNA follows the Mendelian law and is inherited in a codominant manner, which makes it more powerful than previous markers in studying the relationship among alleles in an individual. It can also distinguish between homozygous and heterozygous alleles [12]. It is easy to obtain a microsatellite as it is ubiquitous in the genome. PCR can be used to obtain sequence details, which can then be analyzed with genetic statistical methods for further analysis. Therefore, microsatellite markers are used in many fields, including genetic breeding [13], genetic resource status assessment [14,15], genetic linkage analysis [16,17], and genetic mapping [18,19].
Currently, the population genetics of M. unguiculatus are mainly studied using mitochondrial molecular markers [20][21][22]. There are few reports of using SSR to analyze M. unguiculatus. The objectives of this study were to investigate the genetic diversity and population genetic structure of M. unguiculatus using microsatellite markers and compare the results with those obtained using mitochondrial molecular markers. This study aims to provide information on the geographical distribution patterns and evolutionary processes of this species, as well as to fill knowledge gaps in related research fields. Additionally, the new data obtained will aid fishery resource managers in assessing present resources and developing conservation strategies.

Sample Collection and DNA Extraction
The seven sampling locations of M. unguiculatus crossed from the Yellow Sea to East China Sea (i.e., Qingdao, Zhoushan, Xiangshan, Yuhuan, Taishan, Pingtan, and Xiamen) ( Table 1 and Figure 1). In total, 50 individuals were collected and analyzed from each population (Table 1, Figure 1). Specimens were collected from the wild environment rather than areas of aquaculture breeding programs. Adductor muscle fragments were removed from each mussel and stored in 95% ethanol and below −20 • C until DNA extractions. The DNA extraction followed the method that was reported in Aljanabi and Martinez (1997) with slight modifications [23]. DNA quality was determined by electrophoresis in 1% agarose gel. Total DNA concentration was diluted to 30-40 ng/µL and stored at −20 • C for further analysis.

Genetic Data Analysis
After genotyping, the number of alleles (N a ), observed (H o ), and expected (H e ) heterozygosity were evaluated using GenAlEx v6.501 [25]. The FSTAT v2.9.3 software [26] was used to estimate the allele richness (A r ) by setting the minimum sample size to 46 and the genetics inbreeding index (F IS ). The MICRO-SATELLITE ANALYSER (MSA) v4.05 was used to estimate the genetics differentiation index (F ST ) [27]. The program GENPOP v4.5 [28] was used to test each marker for Hardy-Weinberg Equilibrium (HWE) and the p value. Arlequin v3.5 [29] was performed the term hierarchical analysis of molecular variance (AMOVA) with 10,000 permutation. Using program NeESTIMATOR v2.0 [30] to estimate the effective size of each population based on the linkage disequilibrium (each allele's frequency > 0.05 CI = 95%). The program BOTTLENECK v1.2 [31,32] performed the bottleneck test to analyze whether there had been a population event (expansion or bottleneck). The setting was as follows: model = Infinite allele model (I.A.M.) [33], Stepwise mutation model (S.M.M.) [34] and Two-phase model of mutation (T.P.M.) [35] with 10% I.A.M., 90% S.M.M., and 10,000 permutations [36]. The Wilcoxon test was used to check whether there was significant excess heterozygote [37]. Finally, the software MICROCHECKER v2.2.3 [38] was performed to test for the presence of null alleles.
The program Populations v1.2 [39] was used to construct the NJ phylogenetic tree among seven populations based on the Nei's standard genetics distance [40]. The software STRUCTURE v2.3 was performed the clustering analysis based on Bayesian method (K set 1 to 7, 20 runs, each run used an admixture model, MCMC = 1,000,000, burn-in = 25,000). The results were submitted to the online tool STRUCTURE HARVESTER [41] for evaluation of the best K value and the program CLUMPP [42] was used to plot the STRUCTURE analysis.

Genetic Diversity
Results from the analysis of genetic diversity showed the average numbers of alleles (N a ) ranged from 5. (F IS ) ranged from 0.14 (YH) to 0.19 (ZS, XS, and TS). All loci did not significantly deviate from HWE (Table 2).

Genetic Structure
Based on the geographic distribution of sea areas, M. unguiculatus were divided into two groups (Group 1: QD and Group 2: ZS, XS, YH, TS, PT and XM). Results of AMOVA showed that genetic differentiation among groups was over 4.1% of total variation, with 1.3% attributed to variation among populations within groups and the remaining 94.6% attributed to within populations ( Table 3). The pairwise F ST values revealed that Qingdao had a slight significant genetic divergence with the other six populations. In the six investigated populations from the East China Sea, most pairwise genetic distances based on fixation index (F ST ) were not significant. The gene flow (N m ) ranged from 3.80 to 833.08, indicating frequent gene exchange between populations at different sampling locations ( Table 4). The Bayesian assignment analysis supported the result of pairwise F ST . Under the best K value (K = 2, Figure 2), Qingdao was divergence with other six populations ( Figure 3). The phylogenetic tree also showed that the population of Qingdao was far from other populations (Figure 4).

Historical Dynamics
The Ne estimates for seven populations of M. unguiculatus were large, ranging from 694.5 on TS (CI = 207.6 − infinity) to infinity on QD (CI = 372.6 − infinity). In the bottleneck test, based on three different models, it was not significant that the heterozygote exceeded any other populations. The allele distribution conforms to the L-shaped distribution without bottleneck effect events. It indicates that no genetic bottleneck has recently occurred in the populations of M. unguiculatus along the coastline of the China Sea (Table 5).

Historical Dynamics
The Ne estimates for seven populations of M. unguiculatus were large, ranging from 694.5 on TS (CI = 207.6 − infinity) to infinity on QD (CI = 372.6 − infinity). In the bottleneck test, based on three different models, it was not significant that the heterozygote exceeded any other populations. The allele distribution conforms to the L-shaped distribution without bottleneck effect events. It indicates that no genetic bottleneck has recently occurred in the populations of M. unguiculatus along the coastline of the China Sea (Table 5).

Historical Dynamics
The N e estimates for seven populations of M. unguiculatus were large, ranging from 694.5 on TS (CI = 207.6 − infinity) to infinity on QD (CI = 372.6 − infinity). In the bottleneck test, based on three different models, it was not significant that the heterozygote exceeded any other populations. The allele distribution conforms to the L-shaped distribution without bottleneck effect events. It indicates that no genetic bottleneck has recently occurred in the populations of M. unguiculatus along the coastline of the China Sea (Table 5). Furthermore, the mode-shift test also indicated no evidence for bottleneck events in the recent evolutionary history of M. unguiculatus ( Figure 5).

Discussion
High genetic diversity generally enables species to adapt to environmental changes and avoid inbreeding. Microsatellite alleles and heterozygosity are crucial parameters to determine the genetic diversity of a population [43][44][45]. A population that exhibits heterozygosity greater than 0.5 is generally considered to have high genetic diversity [46]. In the present study, the Ho and He of M. unguiculatus were 0.61 to 0.71 and 0.72 to 0.83, respectively. This finding is consistent with previous studies by Liu et al. [20], Feng et al. [21], and Wei et al. [22] that investigated the same species in the same area using mitochondrial molecular markers. For other marine shellfish analyzed using microsatellite markers, Meretrix meretrix had an average Ho of 0.3878 and He of 0.7996 [43], Crassostrea angulata had an average Ho of 0.647 and He of 0.872 [47], and Sinonovacula constricta had an average Ho of 0.7525 and He of 0.6866 [48]. In comparison with the Ho and He values of M. unguiculatus obtained in this paper, the genetic diversity level of M. unguiculatus is slightly higher than that of Meretrix meretrix, which is close to that of Crassostrea angulata and Sinonovacula constricta, indicating a high level of genetic diversity. This may be related to the in vitro fertilization of M. unguiculatus, the ability to reproduce on a large scale, and the spread of larvae [21]. The results of genetic diversity showed that the current germplasm resources of M. unguiculatus were in good condition and had good breeding prospects.
The FIS value is an important index to measure the genetic dynamics of the population, reflecting the balance between Ho and He. When FIS is positive, it indicates that the heterozygote is missing; when FIS is negative, it indicates that the heterozygote is excessive [49,50]. The FIS values of the seven populations in this study ranged from 0.14 to 0.19, and

Discussion
High genetic diversity generally enables species to adapt to environmental changes and avoid inbreeding. Microsatellite alleles and heterozygosity are crucial parameters to determine the genetic diversity of a population [43][44][45]. A population that exhibits heterozygosity greater than 0.5 is generally considered to have high genetic diversity [46]. In the present study, the H o and H e of M. unguiculatus were 0.61 to 0.71 and 0.72 to 0.83, respectively. This finding is consistent with previous studies by Liu et al. [20], Feng et al. [21], and Wei et al. [22] that investigated the same species in the same area using mitochondrial molecular markers. For other marine shellfish analyzed using microsatellite markers, Meretrix meretrix had an average H o of 0.3878 and H e of 0.7996 [43], Crassostrea angulata had an average H o of 0.647 and H e of 0.872 [47], and Sinonovacula constricta had an average H o of 0.7525 and H e of 0.6866 [48]. In comparison with the H o and H e values of M. unguiculatus obtained in this paper, the genetic diversity level of M. unguiculatus is slightly higher than that of Meretrix meretrix, which is close to that of Crassostrea angulata and Sinonovacula constricta, indicating a high level of genetic diversity. This may be related to the in vitro fertilization of M. unguiculatus, the ability to reproduce on a large scale, and the spread of larvae [21]. The results of genetic diversity showed that the current germplasm resources of M. unguiculatus were in good condition and had good breeding prospects.
The F IS value is an important index to measure the genetic dynamics of the population, reflecting the balance between H o and H e . When F IS is positive, it indicates that the heterozygote is missing; when F IS is negative, it indicates that the heterozygote is excessive [49,50]. The F IS values of the seven populations in this study ranged from 0.14 to 0.19, and the F IS values were significantly positive. The F IS values estimated were between 0 and 1, which was consistent with the results that the H o was less than the H e of the study populations. This result indicates that there is a loss of heterozygosity in the mussel population. This indicates that there may be inbreeding in the population, and it may also be caused by rare base deletion caused by human interference and other factors [51,52].
The results of the genetic structure analysis and population comparison revealed a low level of genetic variation in Qingdao population and confirmed their distribution on areas of sea. The F ST value is an important index to reflect the degree of population genetic differentiation [53]. Wright and Maxson [54] distinguished four levels of genetic differentiation, based on F ST estimations, little differentiation (0 < F ST < 0.05), moderate differentiation (0.05 < F ST < 0.15), large differentiation (0.15 < F ST < 0.25), and very large genetic differentiation (F ST > 0.25). Therefore, the pairwise F ST (0.003~0.065) was quite small for the populations studied. The pairwise F ST values, including Qingdao population, were higher than others. The phylogenetic tree (Figure 4) showed that the genetic distance of Qingdao was far away from the other populations and the STRUCTURE analysis ( Figure 3) also supported that it was divided into two groups based on the geographic distribution. The results of N m can be confirmed by F ST results. Relevant studies suggest that when N m > 4, the effect of genetic drift can be ignored [55,56]. All populations are in a state of random mating, and gene flow is the main factor affecting population genetic differentiation. In this study, the N m values of Qingdao population and the other six populations were all around 4, while the N m values between other populations were far greater than 4. This shows that except for the Qingdao population, the other six populations have frequent gene exchanges and no obvious genetic differentiation. However, Wei et al. [22] used mitochondrial molecular markers to study M. unguiculatus in the same area. The results showed that the degree of genetic differentiation between the seven populations was not obvious, and no obvious genetic differentiation was found in the genetics of each population. It was caused by the origin and specificity of mtDNA. With genetic conservation of genes, researchers made it possible to use some specific mtDNA to gain mitochondrial DNA sequences/genome from unknown species [57]. However, the genetic convenience of mtDNA also limited the detection range and detection sensitivity.
The results of this study indicated that M. unguiculatus exhibits a weak genetic structure, possibly due to the influence of ocean currents. The China Sea has a complex system of ocean currents, including the Kuroshio current and its tributaries, the Coastal currents system (e.g., the northern Jiangsu current, the Zhejiang and Fujian current, and the South China Sea current), and local circulation, upwelling, and eddy currents (e.g., the central cold water mass in the Yellow Sea in summer and autumn, the northern cyclone vortex in the East China Sea, and the Yangtze dilute water) [58][59][60][61]. The influence of ocean currents on species formation and population genetic structure differentiation of marine organisms has long been confirmed [62]. Marine plankton have high passive dispersal ability, which can spread along the direction of ocean currents during the planktonic period, thus promoting genetic exchange between different populations [63][64][65]. The East China Sea is the main producing area of M. unguiculatus in China, and its mass breeding season is in winter every year. During this period, a large number of mussel larvae will be dispersed in the sea. M. unguiculatus with a long planktonic period (about 35d) [66] drift for a long distance under the impetus of marine currents, resulting in genetic exchanges and a weak genetic structure. The current artificial breeding technology for M. unguiculatus is relatively mature, and aquaculture practitioners typically obtain seedlings from the main producing areas and transport them to breeding sites. The aquaculture method involves using attached substrates for seedling culture in natural sea areas. This open breeding method leads to a large number of cultured sample larvae being dispersed in the sea during the breeding season and becoming mixed with wild mussel larvae. As a result, a weakening of the genetic structure of M. unguiculatus can occur due to the large number of cultured sample larvae mixing with wild mussel larvae and spreading with ocean currents.
The investigated seven populations of M. unguicultas did not show significant traces of bottleneck effects, indicating that crossbreeding between populations is likely the main factor contributing to the current level of genetic variation. Inbreeding among populations could result in the loss of genetic information, making it undetectable through molecular methods. Although an artificial seedling technique has been developed for farming M. unguiculatus, most of the seedling parents still originate from the natural producing area in the East China Sea. The increase in breeding demand has led to the embezzlement of natural habitats, resulting in a decline in natural samples and a reduction in genetic variation. At present, the breeding industry of M. unguicultas lacks natural population protection management, which can lead to the dispersal of a large number of cultured larvae into the sea during the breeding season. These larvae can mix with natural mussel larvae and weaken the genetic structure of M. unguiculatus in the wild population. Therefore, to conserve the germplasm resources of M. unguiculatus, we recommend reducing the genetic exchange between cultured and natural populations and minimizing the contamination of genetic information in natural populations.

Conclusions
Population genetics research on M. unguiculatus is critical in revealing its situation in the natural environment. In this study, seven M. unguiculatus populations in coastal areas of China were analyzed by microsatellite markers. The results showed that M. unguiculatus had high genetic diversity. Qingdao and other populations had genetic differentiation. The overall genetic structure was weak, and the presence and size of any bottleneck effect on genetic variation were difficult to assess. The genetic data generated in this study could provide insights into genetic management units and the utilization of M. unguiculatus resources.