Population Genetic Analysis of the Wild Hard-shelled Mussel, Mytilus Unguiculatus (Valenciennes 1858) in South Korea Using a Microsatellite Multiplex Assay

The Korean or hard-shelled mussel, Mytilus unguiculatus, previously known as Mytilus coruscus, is one of the most economically and ecologically important bivalves in South Korea. However, the population size of this species has drastically reduced owing to overharvesting and habitat shrinkage. Because its genetic information is poorly documented, we contributed, in this study, the genetic diversity and structural analyses of 246 adult samples of M. unguiculatus from seven populations along the coastal areas of the mainland and islands of South Korea using a microsatellite multiplex assay. Genetic diversity analyzed from eleven polymorphic microsatellite loci was consistently moderate (0.50–0.57) in all populations. No recent bottleneck was found, indicating that the number of the studied populations did not decrease to an extent that resulted in a reduction of genetic diversity. Additional tests did not reveal any genetic structure across them, possibly resulting from constant gene flow, strong dispersal of planktonic larvae, and genetic admixture between wild populations. These results suggest that M. unguiculatus populations along the coastal areas of South Korea should be managed as a single unit. Our study provides crucial information for future genetic monitoring, conservation management, and population restoration plan in preparation for the rapid decline in bivalve resources.


Introduction
The Korean or hard-shelled mussel, Mytilus unguiculatus (Valenciennes 1858), previously known as Mytilus coruscus, is a species of marine bivalve mollusks widely distributed in the coastal waters of China, Japan, and Korea.This species is not only one of the most economically important species among the favored bivalves owing to its high nutritional value and large size (An and Lee 2012;Kang et al. 2013b), but also ecologically important in intertidal and subtidal areas (Li et al. 2020).Recently, there has been a drastic decrease of its populations because of overharvesting for food stock, climate change, and habitat competition with an invasive species, M. galloprovincialis (Yi et al. 2021).To sustain and recover its natural populations through effective conservation and management strategies, it is necessary to carry out a comprehensive investigation on population genetics of M. unguiculatus (i.e., genetic diversity and structure).
Among several available molecular markers, microsatellite or simple sequence repeats (SSR) DNA markers (unique sequences of flanking microsatellite regions) have been extensively used for studying population genetics due to their advantageous features such as codominant transmission, high level of mutation rate and polymorphism, abundance in various eukaryotic genomes, and ease of use and application (Ellegren 2004;Karhu 2001;Vieira et al. 2016).Commonly, several microsatellite markers are fluorescently labeled and combined with multiplex PCR, allowing multiple microsatellite loci to be simultaneously amplified in a single reaction.Subsequently, the multiplex products are analyzed using gel capillary electrophoresis and then genotyped.The multiplex PCR assay using microsatellite markers, a time and cost-effective technique, has become a powerful tool to study the genetic diversity and structure of natural populations of bivalves including oysters (An et al. 2013;Galvão and Hilsdorf 2015), mussels (Wenne et al. 2022), clams (Hargrove et al. 2015), and scallops (Morvezen et al. 2015).However, there have been limited studies on population genetics of M. unguiculatus.
In this study, we conducted comprehensive genetic diversity and structure analyses of seven M. unguiculatus populations along the coastal areas of the mainland and islands of South Korea with two sets of a multiplex PCR assay (consisting of eight polymorphic microsatellite markers each) from previous studies (Xu et al. 2010;An and Lee 2012;Kang et al. 2013a;Guo et al. 2013;Fu et al. 2018).The findings of this study will provide the crucial basis for developing effective conservation and management strategies for this economically and ecologically important species in South Korea.

Sample Collection and Genomic DNA Extraction
A total of 246 adult samples of M. unguiculatus were collected from seven areas along the coastal coasts of the mainland and islands of South Korea where its major distributions were documented, including the southeast: Hansan-do (HSA), Geomun-do (KOM), Goseong-do (KOS); the west: Sapsi-do (SAB), Socheong-do (SOC); and the east: Uljin (UJI), Dok-do (TOK) from 2021 to 2022 (Table 1, Fig. 1).These sampling sites were determined to postulate the initial hypothesis of geographic isolation by distance.The adductor muscle of each sample was excised with sterile scissors and preserved in absolute alcohol, until its total genomic DNA (gDNA) was extracted following the protocol established by Asahida et al. (1996).Its amount and quality were checked using NanoDrop One Microvolume UV-Vis Spectrophotometer (Thermo Fisher Scientific Inc., Wilmington, DE, USA).

Microsatellite Marker Selection and Genotyping
A total of 65 polymorphic microsatellite markers that had previously been characterized in M. unguiculatus were tested on three random samples from each population.The optimization of PCR settings in multiplex panels was applied to the markers with consistent PCR amplification   S1).The PCR amplification of sixteen microsatellite loci was carried out in a 20-µl volume using AccuPower® PCR PreMix (Bioneer Inc., Daejeon, South Korea) including 1 µl gDNA and each 1 µl forward and reverse primers (the forward primer was fluorescently labeled with FAM, HEX, TAMRA, and ATTO).Amplification was performed in a ProFlex PCR System (Thermo Fisher Scientific Inc.) under the following reaction conditions: initial denaturation at 94 °C for 5 min, 35 cycles of denaturation at 95 °C for 20 s, annealing at 55 °C (Set A) and 50 °C (set B) for 20 s, extension at 72 °C for 30 s, and final extension at 72 °C for 5 min.The fluorescent PCR products were sized and separated on an ABI 3130xl DNA Analyzer (Applied Biosystems, Foster City, CA, USA).Allele sizes were scored against GeneScan 500LIZ Size Standard using GeneMapper v.4.0 (Applied Biosystems).

Genetic Diversity Analysis
Genetic diversity indices at locus and population levels, i.e., number of alleles per locus (Na), effective number of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), polymorphic information content (PIC), and Wright's inbreeding coefficient (F is = 1 -Ho / He) were calculated using Cervus v.3.0.7 (Marshall et al. 1998) and GenAlEx v.6.5 (Peakall andSmouse 2006, 2012).Null allele frequencies (F null ) at each locus and in each population were estimated using the Expectation Maximization (EM) algorithm (Dempster et al. 1977) with 1,000 bootstrap replicates using FreeNA (Chapuis and Estoup 2007).Micro-Checker v.2.2.3 (Oosterhout et al. 2004) was used to evaluate the microsatellite data by identifying genotyping errors, allelic dropouts, and null allele presence.Deviation from Hardy-Weinberg equilibrium (HWE) of each locus across populations was computed using an exact HW test (Guo and Thompson 1992) in Genepop v.4.7.5 (Raymond and Rousset 1995;Rousset 2008).This software was also used for testing the heterozygote excess and deficiency of each locus per population and the linkage disequilibrium (LD) between locus pairs.

Genetic Differentiation and Structure Analysis
Determination of appropriate statistics (allele identity-based F ST or its analog, allele size-identity R ST ) for assessing the genetic differentiation of microsatellite loci was performed using the allele size randomization test (Hardy et al. 2003) implemented in SPAGeDi v. 1.5 (Hardy and Vekemans   (Pritchard et al. 2000) was used to identify genetic population structure (K genetic clusters) and assign individuals to those clusters using a Bayesian model-based cluster method.Determination of K (from 2 to 10) was carried out with 10 independent replicates for each K in 100,000 Markov Chain Monte Carlo steps and a burn-in period of 1,000,000.The results of Structure v. 2.3.4 with the most likely K value were verified by an online program, STRU CTU RE HARVESTER (Earl and VonHoldt 2012) using the delta K method established by Evanno et al. (2005).
Principal component analysis (PCA) was conducted to identify the major patterns of microsatellite loci across populations using the packages FactoMineR v. 1.3.4(Le et al. 2008) and Factoextra v.1.0.7 (Kassambara and Mundt 2020) in R software v. 4.2.1 (R Core Team 2022).In addition, GenAlEx v.6.5 was used to calculate the genetic distance matrix according to Nei (1972) and it was subsequently imported to MEGA v. 10.2.2 (Tamura et al. 2021) for constructing an unweighted pair group method with arithmetic mean (UPGMA) dendrogram.

Isolation by Distance
Isolation by distance (IBD) among populations was assessed by regressing the pairwise estimates of F ST / (1 -F ST ) against the logarithm of the pairwise geographic distance (km) (Rousset 1997) using the nonparametric Mantel test with 9,999 permutations implemented in GenAlEx v.6.5.

Bottleneck Analysis
To observe the recent bottleneck events among populations, Wilcoxon signed-rank test established by Luikart (1997) was used with 1,000 iterations in the software package BOTTLE-NECK v. 1.2.02 (Piry et al. 1999) under three different mutation models: infinite allele model (IAM), two-phased model of mutation (TPM), and stepwise-mutation model (SMM).According to Piry et al. (1999), the recommended parameters for our microsatellite dataset are TPM = 95% SMM and variance among multiple steps of 12. Mode-shift indicators demonstrating the allele frequency distribution could detect any population under mutation-drift equilibrium (approximately L-shape) or recent bottleneck (mode shift or others).

Polymorphic Microsatellite Verification
Sixteen microsatellite markers available in previous studies (Xu et al. 2010;An et al. 2013;Kang et al. 2013a;Fu et al. 2018) tested in this study were characterized across the seven M. unguiculatus populations and summarized in Table S1.All of them successfully produced PCR products for all samples (size ranging from 64 to 376 bp) of the seven populations with no evidence of scoring errors or large allelic dropout indicated by Micro-Checker v.2.2.3.The microsatellite markers were polymorphic, ranging from 3 (P31) to 42 (P05 and P39) in the number of alleles.The greatest Ho and He were found in the locus P05 (0.742 and 0.963, respectively), while these values were lowest in P32 (0.153 and 0.186, respectively).The PIC results showed that eleven loci (P26, P36, P05, P30, P25, P60, P39, P50, P52, P21, and P45) were highly informative (PIC > 0.5), three (P57, P07, and P34) moderately informative (0.25 < PIC < 0.5), and two (P31 and P32) slightly informative (PIC < 0.25).Among them, two loci had negative F is values (P31 and P34), most likely resulting from their heterozygote deficiency.Seven loci showed a significant excess of homozygosity and were flagged with null allele presence.Among them, two loci (P05 and P30) with null allele frequency (F null ) lower than 0.2 were retained due to their little influence on the average genetic diversity and structure, whereas five loci (P26, P25, P60, P39, and P21) with F null values greater than 0.2 were excluded in the further analysis.There was no genotypic LD for any locus pair across all populations, indicating that each could be treated as an independent genetic marker (data not shown).Nine loci (P26, P05, P57, P30, P25, P60, P39, P21, and P45) significantly deviated from HWE (p < 0.01).

Genetic Diversity
After the polymorphic microsatellite verification, eleven loci (P36, P05, P57, P30, P31, P07, P50, P52, P34, P45, and P32) were finally selected to study the genetic diversity of 246 individuals from the seven M. unguiculatus populations in South Korea.Table 2 summarizes the genetic variation for each locus in each population.The mean number of alleles of all loci (Na) at the population level varied from 7.45 (SOC) to 10.27 (KOM).The mean effective number of loci (Ne) was consistently lower than the mean number of alleles in all populations ranging from 4.12 (SAB) to 4.62 (KOM).There was little discrepancy between the mean Ho   and He values among the seven populations.In addition, Ho (0.50-0.57) was found to be constantly lower than He (0.56-0.59) in all populations.As a result, the inbreeding coefficient F is positively approximated to 0.00 in all populations.The mean value of Ne (5.10) was also lower than that of Na (13.27) in total (grouping all populations into one).
The mean values of Ho, He, and F is in total was 0.52, 0.58, and 0.08, respectively.Five loci (P05, P57, P30, P52, and P45) out of eleven significantly deviated from HWE across the seven populations.Among them, P05, P30, and P45 showed significant HWE deviation in all populations (except TOK and UJI at locus P30), while the loci P57 and P52 were only significant in two populations (KOM and SAB) and one population (KOM), respectively.Estimation of significant HWE deviation in total was detected at the seven loci (P05, P57, P30, P50, P52, P45, and P32).Most deviations were associated with the statistically significant result of the heterozygote deficit test (p < 0.05) (except for locus P45).

Genetic Differentiation and Structure
The pairwise F ST values (Table 3) of the finally selected eleven microsatellite loci showed a very low degree of genetic differentiation among the seven M. unguiculatus populations (ranging from 0.000 to 0.009).Most of them were not statistically significant except for the pair of SOC and HSA (F ST = 0.001) and those of TOK and HSA (F ST = 0.000).
The AMOVA results performed on the three putative groups of localities, including the southeast: HSA, KOM, KOS; the west: SAB, SOC; and the east: UJI, TOK, were displayed in Table 4.The majority of genetic variation (99.86%) was distributed within populations, whereas genetic variation among the three geographic groups was asymptotic to 0.000.All fixation indices (F ST , F SC , and F CT ) were not significant and low, indicating that none of the three groups were structured.This result was supported by those from the PCA analysis (Fig. 2).
The Evanno method from STRU CTU RE HARVESTER based on the finally selected eleven microsatellite loci across the seven populations suggested that the most likely number of groups for the dataset was K = 3 corresponding to a delta K value of 7.790 (Figure S1).However, the bar plot of K = 3, 4, and 5 similarly showed that all individuals clustered as a mixture, regardless of populations (Figure S2).The Bayesian clustering analysis found no significant structures among the populations.
The genetic distance among the seven populations according to Nei's method was substantially low (ranging from 0.012 to 0.019) (data not shown).The UPGMA dendrogram based on the Nei's standard genetic distance of the seven populations (Fig. 3) showed three distinct clusters according to three sampling areas including the southeast: HSA, KOM, KOS; the west: SAB, SOC; and the east: UJI, TOK.

IBD
The regression between pairwise population estimates of F ST / (1 -F ST ) and the logarithms of geographic distance in the seven M. unguiculatus populations (Fig. 4) slightly showed a positive correlation (R 2 = 0.0132).However, no statistical significance between them (p = 0.349) was detected from the Mantel test indicating that there was no IBD among them.

Bottleneck Analysis
Under three mutation models (IAM, TPM, and SMM), no significant heterozygosity excess was detected in any of the seven M. unguiculatus populations and in total (p > 0.05) based on the finally selected eleven microsatellite loci.Moreover, the graphical descriptor of the allele distribution exhibited an L-shape, as expected under mutation-drift equilibrium, in all populations and in total (Table 5).These results suggest that there had not been a recent genetic bottleneck in these populations.

Discussion
Genetic variation is of considerable interest for population genetic analysis in natural populations.This estimation, measured by Ho and He, serves as an indicator not only for monitoring but also predicting the maximum amount of genetic variability for population breeding, restoration, and management plans (Allendorf 1986).A high degree of genetic diversity increases evolutionary potential and adaptability to environmental changes, which is crucial for the long-term survival of populations (Zheng et al. 2022).Typically, populations experiencing a recent drastic reduction of effective population size due to sudden climate changes, habitat shrinkage, inbreeding, repeated harvesting, and marine pollution are generally characterized by low genetic diversity (Yi et al. 2021).In this study, the average genetic diversity estimated from the finally selected eleven microsatellite loci was similarly moderate in the seven M. unguiculatus populations in South Korea.Moreover, no recent genetic bottleneck was detected based on three allele mutation models (IAM, TPM, and SMM) from the Wilcoxon signed-rank test, indicating that sufficient genetic variability persisted across the examined populations under mutation-drift equilibrium, as expected.This result is consistent with a previous report by Kang et al. (2013b).Although a good indicator of genetic diversity was revealed in the studied populations, for reparation of the rapid decline in M. unguiculatus populations, continuous genetic research in these populations is necessary for planning an optimal conservation and management strategies.
Five out of eleven microsatellite loci finally selected in this study showed significant HWE deviations due to heterozygote deficiency commonly detected in aquatic organisms (An et al. 2013;Galvão and Hilsdorf 2015;Gonzalez et al. 2012;Valles-Jimenez et al. 2004;Yu and Li 2007).Many factors could cause these deficit deviations, such as biological factors (inbreeding, Wahlund effect, and gene flow), technical factors (null alleles, technical artifacts, and small sampling size) (Valles-Jimenez et al. 2004).Among them, non-random mating with close relatives (inbreeding) could be ruled out from this study, as there were no genotypic LD among microsatellite loci and a relative frequency of heterozygotes inferred from the quantitative calculation, the confirmation of a lack of bottleneck effect, and mean inbreeding coefficient F is of 0.000.The presence of null alleles could have caused heterozygote deficits.These deficits are supposedly caused by polymorphic mutations at binding sites of PCR primers resulting from the amplification failure of particular alleles, misscoring due to stutter bands, and slipped-strand mispairing of di-nucleotide repeat loci during PCR (Buttler 2005;Carlsson 2008;Oosterhout et al. 2006;Valles-Jimenez et al. 2004).Although most of these five microsatellite loci carry di-nucleotide repeat motifs, except for P57 (tri-nucleotide motif) and P30 (tetra-nucleotide motif), scoring errors due to stutter bands was not a major cause of heterozygote deficiency, confirmed by our strict  Absence of geographic barriers in the ocean, large population sizes, and high dispersal potential of the planktonic larvae phase in the life cycle generally leads to a poorly geographical differentiation among populations of marine species, especially bivalves (Launey et al. 2002).This characteristic has been observed in many marine bivalves, such as mussels (Han et al. 2017), pen shells (An et al. 2012), and oysters (An et al. 2014).The pattern of M. unguiculatus population structure varies greatly depending upon hydrological, ecological, and physical barriers, allowing gene flow across populations, which is similar to other bivalves (Palumbi et al. 1997), although its biology and ecology is poorly documented.In this study, the geographic distance between HSA and KOS was short (25.9 km) compared with that between SOC and TOK (630.4 km) among the seven M. unguiculatus populations.However, regardless of distance, low genetic differentiation and weak geographic isolation among them were revealed by the results of the pairwise F ST values, PCA, AMOVA, and IBD analyses.These results may be due to persistent gene flow either between wild populations or the passive dispersal of M. unguiculatus planktonic larvae along marine currents, detected through a strong genetic admixture.Such findings are supported by previous studies that used either mitochondrial cytochrome c oxidase I marker (Shen et al. 2009;Yi et al. 2021) or microsatellite markers (Kang et al. 2013b) on M. unguiculatus populations in South Korea and China.Absence of genetic structure detected in this study suggests that different M. unguiculatus populations across South Korea should be managed as a single unit.
In order to investigate genetic differentiation among M. unguiculatus populations, Nei's standard genetic distance, closely related to IAM, a size-independent mutation model where each mutation is likely to produce an allele of any size, is usually used (Kimura and Crow 1964;Wang et al. 2001).Many previous studies have used this parameter to infer population differentiation based on microsatellite data (Alam and Islam 2005;Ball and Chapman 2003;Barat et al. 2015;Chen et al. 2005;Gonzalez et al. 2012;Kim et al. 2012).However, calculation of genetic differentiation from microsatellite data based on IAM instead of using a size-dependence mutation model such as SMM is supposedly not appropriate (Slatkin 1995).Hence, selection of a suitable mutation model for a microsatellite dataset remains controversial.In this study, Nei's standard genetic distance was significantly low, but showed that the seven M. unguiculatus populations were clustered into the three groups according to the three geographic areas, fitting the initial hypothesis of geographic separation.This conclusion should be considered carefully because of a significantly low genetic variation among populations or unfit selection of a mutation model.However, our results may be a valuable reference for a future sampling plan to explore the genetic structure of M. unguiculatus populations in South Korea.

Conclusions
The seven M. unguiculatus populations from the three geographic groups, i.e., the southeast, the west, and the east of South Korea consistently maintained a moderate genetic diversity.Moreover, the low genetic differentiation and weak geographic separation were detected across them, suggesting that different populations along the coastal areas of the mainland and islands of South Korea should be managed as a single unit.Such findings contribute essential information for future resource management and population restoration, thus aiding optimizing inbreeding programs for this economically and ecologically important marine bivalve species.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http:// creat iveco mmons.org/ licen ses/ by/4.0/.
and obvious allelic size variation.Finally, two multiplex PCR panels (set A and set B) were obtained, each with eight microsatellite markers, generated by Xu et al. (2010), An and Lee (2012), Kang et al. (2013a), Guo et al. (2013), and Fu et al. (2018) (Table

Fig. 2 Fig. 3
Fig. 2 Principal Component Analysis (PCA) of the seven Mytilus unguiculatus populations based on the finally selected eleven microsatellite loci.Each population is marked by different colors.Each axis corresponds to the principal component (Dim1) factor 1 with

Fig. 4
Fig. 4 Correlation analysis between pairwise population estimates of F ST / (1 -F ST ) and the logarithms of geographic distance in the seven Mytilus unguiculatus populations based on the finally selected eleven microsatellite loci.The probability of the Mantel test is presented in the graph as the p value

Table 1
Sampling information of Mytilus unguiculatus in South Korea

Table 2
Genetic diversity indices based on finally selected eleven microsatellite loci of 246 Mytilus unguiculatus samples from the seven populations in South Korea

Table 3
The average pairwise F ST values in the lower diagonal and their associated geographic distance (km) in the upper diagonal of the seven Mytilus unguiculatus populations in South Korea HSA Hansan-do, KOM Geomun-do, KOS Goseong-do, SAB Sapsi-do, SOC Socheong-do, TOK Dok-do,

Table 5
Wilcoxon signed-rank tests for heterozygosity excess for the seven Mytilus unguiculatus populations IAM Infinite Allele Model, TPM Two-phase Model, and SMM Stepwisemutation Model for detection of a recent population bottleneck event within each population, HSA Hansan-do, KOM Geomun-do, KOS Goseong-do, SAB Sapsi-do, SOC Socheong-do, TOK Dok-do, UJI Uljin