Population Genetic Divergence among Worldwide Gene Pools of the Mediterranean Mussel Mytilus galloprovincialis

Simple Summary The smooth-shelled marine mussel Mytilus galloprovincialis maintains its genetic integrity as a species on a worldwide scale. Current population genetic analyses confirm the largest divergence of M. trossulus compared to the rest of congeneric species and place M. chilensis as an intermediate taxon between M. galloprovincialis and M. edulis. Unlike previous reports, M. galloprovincialis from the Atlantic Northeast is the most likely source of exotic settlements worldwide. As a super-adaptive species, M. galloprovincialis should not be considered invasive in a human-like supremacist manner, but rather as a flexible evolutionary species (FES). The worldwide distribution of this species suggests that it is naturally endowed with plastic adaptation. Therefore, it could counteract stressful conditions and provide intergeneric ecological opportunities in the face of climatic rarefaction of world coasts. Abstract The Mediterranean mussel Mytilus galloprovincialis is distributed in both hemispheres either natively or introduced. The updated population genetic distribution of this species provides a useful knowledge against which future distribution shifts could be assessed. This study, performed with seven microsatellite markers and three reference species (M. edulis, M. chilensis and M. trossulus), aimed to determine the scenario of genetic divergence between 15 samples of M. galloprovincialis from 10 localities in Europe, Africa, Asia, Australia, North America and South America. In agreement with previous data, M. trossulus was the most divergent taxon of the genus, but M. chilensis appeared as an intermediate taxon between M. edulis and M. galloprovincialis, though closer to this latter. M. galloprovincialis from the Atlantic Northeast appears as the most likely source of worldwide exotic settlements instead of the previously thought Mediterranean population. The successful worldwide establishment of M. galloprovincialis suggests it is a flexible evolutionary species (FES), i.e., a species or population whose genetic background allows it to rapidly adapt to changing environments. This natural endowed plastic adaptation makes it a candidate resilient species amidst the ongoing climatic change.


Introduction
Smooth-shelled mussels of the genus Mytilus are among the most cosmopolitan genera inhabiting marine and estuarine coastal areas over temperate and sub-polar regions [1].Three Mytilus species from the Northern Hemisphere have been profusely studied, i.e., M. edulis Linnaeus 1758, M. trossulus Gould 1850 and M. galloprovincialis Lamarck 1819 [2,3].Those species show distinct latitudinal ranges with patching patterns and hybridization at overlapping areas [4].The temperate M. galloprovincialis evolved in the Mediterranean and later expanded along Atlantic Northeast shores as far as the British Isles and Northern Animals 2023, 13, 3754 2 of 15 Africa [5].Exotically distributed populations of this species are believed to be introduced to Australia [6], Argentina [7], Brazil [8], California [9], Central Chile [10], Southeast Asia [11] and South Africa [12].At those exotic regions, M. galloprovincialis hybridize with native taxa, e.g., in Australia, South America and California, i.e., with southern Hemisphere M. galloprovincialis, M. edulis platensis and M. trossulus, respectively.
Its widespread distribution combined with local environmental pressure on shell shape have traditionally produced a confused taxonomy within this genus [13].Despite their interspecific hybridization and their potential planktotrophic larvae dispersal, Mytilus species maintain their genetic integrity over large geographical areas [14,15].The species status was allowed to develop specific molecular tools related to commercial interests of traceability in fisheries and aquaculture [16][17][18].Also, the reliable identification of species is a prerequisite to determine the natural or Introduced distribution of the Mytilus species, to detect hybridization and introgression, as well as to gauge adaptive responses of conservation pertinence [19,20].Recent advancement of genetic technologies has provided a wide variety of specific molecular markers useful to clarify taxonomic uncertainties within Mytilus spp., e.g., in many parts of Australia and New Zealand [19], East Asia [21], Northern Africa [22] and South America [23].Those studies indicate that the main genetic divergence is between species, but a great deal also exists intraspecifically [24].
In the case of M. galloprovincialis, global concerns about sustainable fisheries and aquaculture management have raised interest in its distribution and dispersal patterns through its native range.Such distribution is generally the byproduct of natural population dynamics for many marine shellfish, i.e., the larval dispersal ability along the coasts determines the gene flow intensity that finally shapes the metapopulation gene pool scenario [25].For instance, the regional genetic distribution of this species has been described on the coasts of the Iberian Peninsula [26], along the crossroads between Southern Europe and Northern Africa [22,27,28], along the Mediterranean and the Black Sea [29,30] and at a local regional scale, such as Galicia, where it has a pivotal ecosystem role [31].Phylogeographic studies have shown that regional gene pools are hardly ever genetically homogenous in their range [32].Although selective forces prompting local adaptation cannot be excluded as causative of the regional divergence, they have rarely been experimentally shown [33].Meanwhile, the gene flow-gene drift balance usually explains satisfactorily the metapopulation patterns observed, which generally fit isolation-by-distance scenarios, excepting those at transitional barriers between biogeographic regions [22].
Understandably, the threat that the expansion of M. galloprovincialis represents to local genetic resources has fueled studies on population dynamics and conservation solutions in many regions, such as South Africa [34], California [35], Brazil [36] and Chile [37].From our personal biological perspective, this species has been too frequently demonized as one of the worst invasive species because of its rapid adaptive success to exotic locations [38,39].Would these negative arguments still be levelled against M. galloprovincialis if within a few decades it became the only intertidal mussel resilient to climate change?The updated global population genetic distribution of M. galloprovincialis provides useful knowledge against which future distribution shifts could be assessed.To date, few studies have been conducted comprehensively on its whole range.One of them is an in silico data mining study on the mtDNA COI sequence distribution, which showed a complex dispersal pattern likely involving a combination of natural and anthropogenic dispersal, coupled with local adaptation and hybridization events [40].A second global study dealt with the genetic background of this species for temperature resilience.Therein, authors reported that the adaptive genetic composition was significantly different among populations and is associated with temperature variables in the Northern Hemisphere [41].
If present knowledge on the distribution of M. galloprovincialis allows the assessment of future distribution shifts, in this study we aimed to determine the scenario of genetic divergence between 15 populations of M. galloprovincialis sampled from 10 localities in Europe, Africa, Asia, Australia, North America and South America.Given the conservation of nDNA microsatellite flanking regions between close congeneric species, [42] as occurs in Mytilus [43,44], we hypothesize the feasibility of identifying sister species, congeneric hybrids and the genetic purity of M. galloprovincialis groups from each geolocation, as well as the putative original sources of its actual exotic distribution.

Sampling
Aiming to screen the global genetic diversity of M. galloprovincialis, the sampling effort was accomplished in 2007 on intertidal areas of Southwestern Europe (Spain), North Africa (Morocco), South Africa (Cape Town), Pacific Northwest (Japan), Pacific Southwest (Australia), Pacific Northeast (California) and the Pacific Southeast (Chile) (Figure 1).congeneric hybrids and the genetic purity of M. galloprovincialis groups from each geolocation, as well as the putative original sources of its actual exotic distribution.

Sampling
Aiming to screen the global genetic diversity of M. galloprovincialis, the sampling effort was accomplished in 2007 on intertidal areas of Southwestern Europe (Spain), North Africa (Morocco), South Africa (Cape Town), Pacific Northwest (Japan), Pacific Southwest (Australia), Pacific Northeast (California) and the Pacific Southeast (Chile) (Figure 1).Three external reference samples of congeneric species were also included in the analyses, i.e., M. edulis (Denmark), M. trossulus (Canada) and M. Chilensis (Chile) (Table 1).The mantle tissue of 492 specimens were conserved in 95% ethanol until DNA extraction and purification following the FENOSALT protocol [45].Three external reference samples of congeneric species were also included in the analyses, i.e., M. edulis (Denmark), M. trossulus (Canada) and M. Chilensis (Chile) (Table 1).The mantle tissue of 492 specimens were conserved in 95% ethanol until DNA extraction and purification following the FENOSALT protocol [45].

Molecular Analyses
All mussels were genotyped with seven polymorphic microsatellites, five of which were previously described ( [43]; Mgµ1, Mgµ2, Mgµ3, Mgµ4, Mgµ5) and employed to genotype M. galloprovincialis from the Iberian Peninsula [26].Two additional markers were employed, microsatellite Mech8 [44] and an unpublished one from M. galloprovincialis, Mgµ8 (forward primer 5 -ATGTCTCCTCAATCTGG-3 and reverse primer 5 -AAATCGTT AAAAAGCAAT-3 ), annealed at 55 • C and 1.7 mM MgCl 2 .PCR amplification consisted of an initial denaturing step at 95 • C for 5 min, followed by 35 cycles at 94 • C for 1 min, 1 min at the annealing temperature, and 1 min at 72 • C for extension.A final extension step was performed at 72 • C for 15 min.The amplified fragments were electrophoresed in an ALFexpress-II automatic fragment analyzer (GE Healthcare, Barrington, IL, USA) and the allele calling was helped by molecular ladders.

Data Analyses
The number of alleles and their frequencies, the allelic richness per locus (R S ) and the fixation indexes F [46] were calculated with FSTAT 3.9.5 [47].The probability test associated to F IS was calculated with the Markov chain method implemented in GENEPOP 4.2 [48], using 20 batches of 5000 iterations each.The expected heterozygosity (H E ), the observed heterozygosity (H O ) and the Hardy-Weinberg equilibrium per sample were also calculated with GENEPOP.Correction for multiple tests was performed with the false discovery rate approach (FDR) [49].Putative frequencies of null alleles co-segregating in the allelic systems and their confidence intervals were checked per locus and sample using the EM algorithm [50] and 1000 permutations, as implemented in FreeNA [51].
Fisher's exact test and Pearson's chi-square test were used to estimate the statistical power of the sampling system at refuting the hypothesis of genetic homogeneity (combining t generations of drift and effective size values to test for a specific F ST , through one batch of 1000 replicates), as well as to estimate the proportion of false significant tests (Type I error, p < 0.05) in combined test statistics (1000 replicates with the same effective size) as implemented in POWSIM 4 [52].The interpopulation fixation index F ST was calculated with FSTAT, and the differentiation parameter D EST [53] between samples was calculated in DEMEtics 0.8-7 [54], as implemented in R-package 2.12.1., using 1000 bootstrap replicates to estimate statistical significance.
Sample relationships upon variance components on the genotype matrix were visualized in a bi-dimensional space using a principal coordinates analysis (pCoA) as available from the statistical package GenAlEx 6.5 [55].A locus-by-locus AMOVA, as implemented in ARLEQUIN 3.5 [56], was used to distribute hierarchically the genetic variance as per six major regions (Atlantic Southwestern Europe and North Africa, South Africa, Pacific Northwest, Pacific Southwest, Pacific Northeast and Pacific Southeast) and per hemisphere.Variance distribution was also computed for congeneric species as reference samples.Statistical tests for each fixation index were based on 1023 permutations.The Bayesian inference on the number of gene pools was explored with BAPS 6 [57], considering both the allele frequencies and the number of genetically divergent groups as random variables, and either an admixture analysis based on 100,000 Bayesian iterations or a mixture model [58].
The first pCoA explained 40% of the divergence between samples and separated unambiguously the four species comprised in the analysis (Figure 2).were observed in M. galloprovincialis, nineteen in M. chilensis and one in each of M. edulis and M. trossulus.Most of the microsatellites exhibited a heterozygosity higher than 70% over all samples, with average HE ± SD = 0.723 ± 0.183 in M. galloprovincialis.The values of FIS ± SD ranged from 0.193 ± 0.116 to 0.398 ± 0.213, and genotypic disequilibrium was observed in most loci, with the putative null allele frequency averaging 0.108 ± 0.082 across samples (Table S2).
The first pCoA explained 40% of the divergence between samples and separated unambiguously the four species comprised in the analysis (Figure 2).The significant admixture estimates from the Bayesian clustering inferred by BAPS showed one gene pool for each species involved in the analysis.Five gene pools were significant within M. galloprovincialis (Figure 3).The largest gene pool includedsix samples, The significant admixture estimates from the Bayesian clustering inferred by BAPS showed one gene pool for each species involved in the analysis.Five gene pools were significant within M. galloprovincialis (Figure 3).The largest gene pool includedsix samples, i.e., North Atlantic Spanish (MgRi and MgSa), South Africa (MgCt), Chile (MgDi) and California (HgtMb and HgtRf).The most heterogeneous pool composition was observed in sample MgRi from Galicia (Spain) and in the two Californian samples.The other four pools from Morocco, South Africa, Australia and Japan showed less admixture.The highest divergence between the k = 8 BAPS pools using the Nei genetic distance was observed between M. galloprovincialis and the rest of the species, i.e., M. trossulus (1.152), M. edulis (0.657) and M. chilensis (0.506).The M. galloprovincialis samples formed a clade that was the sister group of the rest of the species.The clade of M. galloprovincialis included two subclades, one grouping the East Asia samples (Japan and Australia) and the other joining all the Atlantic North samples of Iberia and Morocco and including South Africa, Chile and California (Figure 4).The significant admixture estimates from the Bayesian clustering inferred by BAPS showed one gene pool for each species involved in the analysis.Five gene pools were significant within M. galloprovincialis (Figure 3).The largest gene pool includedsix samples, i.e., North Atlantic Spanish (MgRi and MgSa), South Africa (MgCt), Chile (MgDi) and California (HgtMb and HgtRf).The most heterogeneous pool composition was observed in sample MgRi from Galicia (Spain) and in the two Californian samples.The other four pools from Morocco, South Africa, Australia and Japan showed less admixture.The highest divergence between the k = 8 BAPS pools using the Nei genetic distance was observed between M. galloprovincialis and the rest of the species, i.e., M. trossulus (1.152), M. edulis (0.657) and M. chilensis (0.506).The M. galloprovincialis samples formed a clade that was the sister group of the rest of the species.The clade of M. galloprovincialis included two subclades, one grouping the East Asia samples (Japan and Australia) and the other joining all the Atlantic North samples of Iberia and Morocco and including South Africa, Chile and California (Figure 4).The fixation indexes within M. galloprovincialis showed a high genetic variation between samples (F ST = 0.096*).Such variation was higher within large continental regions (F CT = 0.058, p = 0.0078) than among regions (F CT = 0.046, p = 0.0089), although both were significant (Table 2).No variation was observed between hemispheres (F CT = 0.007, p = 0.258).The largest interspecific variance was observed between M. galloprovincialis and M. trossulus (F CT = 0.230, p = 0.0019) as compared to other pairwise comparisons.

Genetic Diversity of M. galloprovincialis
The degree of interspecific conservation of microsatellite primers and their polymorphism were proportional to the genetic distance between species [44].Herein, all primer pairs from M. galloprovincialis were not only amplified in the sample of M. chilensis but this latter bore a significantly higher allelic richness and number of specific alleles than the former species.This phenomenon points to both the similarity of these genomes and their own specific status [59].The average polymorphism of seven microsatellites of M. galloprovincialis from all the continents (mean H E ± SD = 0.723 ± 0.183) was congruent with previous observations in this species, e.g., H E = 0.772 ± 0.154 from six microsatellites on Iberian samples [26], but slightly higher than in seven microsatellites from the Moroccan samples (H E = 0.552 ± 0.127) [22].This result was expected because of both the sampling amplitude and the population size, i.e., the genetic diversity is maximal in Galician Estuaries (Rías Gallegas) inhabited by the largest world population of this species [31].The global deficit of heterozygotes, especially at loci Mgµ1 and Mgµ2, indicated a significant deviation from the Hardy-Weinberg equilibrium.This phenomenon was reported earlier in microsatellites of M. galloprovincialis [22,26] but was also apparent with allozymes and nuclear DNA markers [2,60] and is common in marine bivalves [61,62].The underlying causes of that deficit range from functional to technical.Some functional hypotheses are selection [63,64], stock admixture [65], inbreeding and genotype-independent spawning [66].More likely, technical causes are stochastic genotyping errors [67] such as allelic dropout or false alleles [68], as well as systematic errors, i.e., null alleles due to primer site sequence variation [69].The high putative null allele frequency (~15%) inferred in this study for two loci across most samples is a reasonable explanation for their heterozygote deficit.The other loci showed a low or moderate null allele frequency (below 10%), but the analytical exclusion of loci Mgµ1 and Mgµ2 did not produce a different outcome regarding genetic diversity, as observed in most studies [70].

Genetic Differentiation between Species
The analysis of principal coordinates separated the four species comprised in the analysis, with M. trossulus as the most distinct taxon of the genus, as shown using DNA, allozyme and morphometrics markers, e.g., [71].M. trossulus is more distantly related to M. galloprovincialis than to M. edulis, i.e., average F ST (Mt-Mg) = 0.312 ± 0.037, p = 0.001 and F ST (Mt-Me) = 0.249, p = 0.001, respectively, and is believed to have been diverging from these species for about 3.5 million years [72].The interspecific relationships upon variance components are congruent with the Nei genetic distance computed after the k = 8 BAPS pools, where M. galloprovincialis diverged from the rest of species-M.trossulus, M. edulis and M. chilensis-by 1.152, 0.657 and 0.506, respectively.Previous studies situated M. galloprovincialis and M. edulis as the most closely related sister taxa of this genus, which coalesced some 2 million years ago [73].However, the present scenario suggests that M. chilensis is an intermediate taxon between M. galloprovincialis and M. edulis but closer to the former, as previously shown with allozymes [59], microsatellites [44] and mtDNA COI [74] (but also see mitogenomic analyses [75]).These population genetic studies on M. chilensis and its latitudinal morphological description [76] provided the first proofs of the specific genetic status of this taxon in Chile.Despite the geographic proximity and hybridization between M. chilensis and M. edulis platensis in Cape Horn [77], M. chilensis seems to be genetically closer to M. galloprovincialis.That similarity could be phylogenetic rather than introgressive because of both the assumed recent introduction of M. galloprovincialis in central Chile [78-80] and its lack of admixture with M. chilensis [77].

Genetic Differentiation within M. galloprovincialis
The amount of interpopulation genetic divergence in M. galloprovincialis observed with seven microsatellites (F ST = 0.102 ± 0.044) was slightly higher than that reported from thousands of SNP markers (F ST = 0.087, [18]) and is likely due to an overestimation of null alleles [70].This result shows that a moderate number of random microsatellites provide enough signal for global population genetic analysis.None of the five regional gene pools identified within M. galloprovincialis using Bayesian clustering was explained by divergence between hemispheres, as has also been reported upon COI gene sequences [40], but better by variation within regions and among regions.The variation within regions is the byproduct of a population connectivity pattern reported in many instances between samples of Northern Africa (Morocco) and Mediterranean coasts, e.g., F ST = 0.044 ± 0.006 [22], or between Atlantic, Alboran and Mediterranean Iberian coasts (F CT = 0.0281, p = 0.023; [26]).Notably, withinregion divergence equaled that among regions, although the latter was expected to be much higher under an isolation-by-distance pattern and independent regional evolution.Except for the Australian and New Zealand samples, the global scenario suggests that the M. galloprovincialis populations inhabiting exotic locations are relatively recently settled.The small genetic divergence (F ST ) and genetic differentiation (D EST ) between close samples, i.e., Japanese or California samples, as well as the high divergence between distant samples, e.g., Mediterranean vs. Australian samples, are observations congruent with previous studies using allozymes [81] and DNA markers [78].

Patterns of Divergence in Parapatry
Although the majority of pairwise comparisons were significant for both D EST and F ST , the higher conservativism of F ST identified a lack of divergence between very distant samples such as North Spain (MgRi) and both South Africa (MgCt) and California (MgtRf) or between this latter and Chile (MgDi).That unexpected genetic similarity was patent in two major subclades within M. galloprovincialis, i.e., one grouping being East Asia samples (Japan and Australia) and the other grouping being all Atlantic North samples including South Africa, Chile and California in both a single tree subclade and a single gene pool.These results suggests that the global genetic scenario of M. galloprovincialis is composed by two major patterns.One pattern comprises those regions diverging evolutionarily in parapatry, e.g., Iberian MgRi (Atlantic) vs. MgOr (Mediterranean) separated by the Almería-Oran Oceanographic Front [82], where strong congruence exists between genetic markers, e.g., allozymes, mtDNA and microsatellites [22,27,81,[83][84][85].A second pattern is shown between distant populations that did not diverge in parapatry but had been recently segregated either from a donor population or from one of its exotic introduction sites.Such recent settlements are understood as accidental through intercontinental trading, e.g., in ballast water or hull foiling [86][87][88], or by aquaculture interests, whether aware or not of the consequences of biological translocations.For instance, there has been much investigation on the invasive capacity of M. galloprovincialis in South Africa, California and Chile after the colonization event of this species in the last century [89][90][91][92].Current data suggest that Atlantic Southwestern Europe is the direct or indirect source of present-day populations in California, South Africa and Chile.Supporting this suggested Atlantic origin is that the largest distance observed between the Mediterranean sample MgOr and the Japanese sample MgYo was even higher (F ST = 0.186) than that observed between species, e.g., M. galloprovincialis-M. chilensis (0.116 ± 0.043).This hypothesis of an Atlantic Northeast origin of South Africa mussels is also supported by haplotype networks and F ST data from mtDNA analyses, although not recognized as such therein [40].Nevertheless, the Atlantic origin hypothesis disagrees with most previous works supporting a single Mediterranean origin of exotic M. galloprovincialis, as claimed using SNPs markers on the Chilean [18,93] or Brazilian samples [36] of M. galloprovincialis.

The Pacific Northeast
The genetic status of M. galloprovincialis from the hybrid zone of California (HgtMb and HgtRf) is congruent with its historical lack of introgression with native M. trossulus [94].This genetic status of Californian M. galloprovincialis agrees with previous studies in the Pacific Northeast from Puget Sound to the central California hybrid zone [95], as well as with the interspecific polarized distribution in the latter region, e.g., Morro Bay [35].Japan and southern Europe have been suggested as putative donors of the multiple introductions suspected to have occurred in the Pacific Northeast [89].The European origin was Mediterranean from analyses of allozyme data [9,81,96] and genomic DNA [78].However, as indicated above for other exotic locations, the Atlantic North European population is the most likely origin of those Pacific settlements, according to current microsatellites.

The Pacific Southeast
The Mytilus species inhabiting the Pacific Southeast is M. chilensis [97], which has been shown to be a genetically distinct taxon in the last decade [44,74].M. galloprovincialis is also present in the Chilean coast [10,79,98], and, to date, no evidence exists on either its expansion beyond the Gulf of Arauco in the Bío-Bío Region or its hybridization with the native M. chilensis [77].Nonetheless, this latter naturally hybridizes with the neighboring species M. edulis platensis in the Southern Cone tip [77].The natural Pacific east occurrence of Mytilus-like fossils in South America [99], as well as those in North America [9,89,94,100], does not help clarify a putative trans-equatorial historical migration of Mytilus and other taxa between these two subcontinents [101].Present knowledge allows for thesuggestion that the Pacific coasts were originally occupied by distinct species, i.e., M. trossulus or its predecessor in the Pacific north, M. californianus [102] and M. chilensis or a putative predecessor in the Pacific south andM.galloprovincialis in the Southern Hemisphere [3].Also, M. galloprovincialis seems to have been introduced multiple times to the Pacific north, via Japan or Europe [89].However, the assumed source origin of Chilean M. galloprovincialis in the Mediterranean [78], is not supported by current microsatellite data.Because of the genetic similarity after Bayesian analyses (see Figure 4) between Northern Hemisphere M. galloprovincialis and the Chilean mussel from Dichato, the origin of this latter appears to be California, Atlantic Europe or Cape town.For instance, given that phylogeographic evidence exists on the accidental introduction of M. galloprovincialis to South Africa [103], the Pacific Southeast population could have its origin in the native distribution area and/or in one of its exotic settlements.

Australia and Japan
A historical circum-Arctic migration from Atlantic European coasts has been reported to explain the existence of mussels in Australia [2,3,104].A trans-equatorial migration of mussels between the North and the South Pacific has also been proposed to explain their occurrence in those regions [81].These hypotheses are not mutually exclusive because a first circum-artic migration could have reached the Pacific Northwest and been followed by a trans-equatorial migration to Australia.In addition, knowledge on copepod parasitism of Japanese mussels [105] suggests a relatively recent human introduction of European mussels into Japanese coasts [106], likely during the Edo period of Japanese history.Whatever hypothesis is correct, the native range of M. galloprovincialis about 1 My ago would include Australia and New Zealand and possibly Chile [1,[107][108][109].Advancing in time, local evolution in parapatry and/or new introductions would have produced actual representatives of this genus, such as Southern-Hemisphere M. galloprovincialis [110] and M. chilensis [97], respectively.
Multiple introductions of M. galloprovincialis into Australia and New Zealand from its Atlantic and Mediterranean natural range have been suggested based on genetic and demographic data [111].Those introductions have led to admixtures with the native Mytilus planulatus over a large amount of the Australian coastline [112], which can explain previous scenarios of genetic heterogeneity of Mytilus samples from Australia [108,109].Current analyses indicate that the two mussel samples from Australia and Japan belong to Northern-Hemisphere M. galloprovincialis, yet they are highly divergent from each other, as well as from the rest of the gene pools of this species (see Figure 4), as has also been reported in COI sequence data [40].The inter-cluster divergence within M. galloprovincialis suggests a common origin of those two samples, while their intra-cluster divergence suggests a younger divergence between them.The above hypotheses on population sourcing from Europe to Japan or to Australia, and then from Japan to Australia or vice versa, can reasonably explain the current parapatric scenario observed with microsatellites.

Conclusions
Despite the high genetic variation exhibited by M. galloprovincialis, it maintains its genetic integrity on a global scale.Microsatellite variation confirms the higher divergence of M. trossulus from its congeneric species and places M. chilensis as an intermediate taxon between M. galloprovincialis and M. edulis.Also, microsatellite variation identifies M. galloprovincialis from the Atlantic Northeast as the most likely source of exotic settlements worldwide.The adaptive potential of M. galloprovincialis allows it to be considered a flexible evolutionary species (FES), i.e., a species or population whose genetic background allows it to adapt rapidly to changing environments.The plastic adaptation of this species makes it a resilient candidate to counteract stressful conditions and provide ecological opportunities to many intertidal genera facing global coastal rarefaction.

Supplementary Materials:
The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/ani13243754/s1,Table S1: Allelic frequencies of seven microsatellites analyzed in 15 samples of Mytilus spp.(sample codes are given in Table 1); Table S2: Genetic diversity of seven microsatellites in twelve samples of M. galloprovincialis, three samples of congeneric species (M.trossulus (MtVc), M. chilensis (MchCA) and M. edulis (MeKa) and two samples from a M. galloprovincialis-M. trossulus hybrid zone (HgtMb and HgtRf); Table S3: Pairwise estimates of genetic differentiation (D EST , above diagonal) and gene distance (F ST , below diagonal) between samples of Mytilus spp.All values except those bolded were significantly different from zero.The significance of D EST was drawn from its 95% confidence interval (Table S4).The significance threshold for F ST was generated after 100 MC batches of 5000 iterations each for alpha = 0.01; Table S4: Confidence intervals (95%) for D EST (above diagonal) and F ST (below diagonal).
Author Contributions: Conceptualization, experimentation, data analyses, and first draft, Y.O. and P.P.; data curation and interpretation of results, A.A.; review and editing, P.P., Y.O. and A.A.; funding acquisition and projects administration, P.P.All authors have read and agreed to the published version of the manuscript.

Figure 2 .
Figure 2. Principal coordinates analysis on the genetic distance FST between samples of M. galloprovincialis (Mg) relative to the control species (Mch, M. chilensis; Me, M. edulis; Mt, M. trossulus).Figure 2. Principal coordinates analysis on the genetic distance F ST between samples of M. galloprovincialis (Mg) relative to the control species (Mch, M. chilensis; Me, M. edulis; Mt, M. trossulus).

Figure 2 .
Figure 2. Principal coordinates analysis on the genetic distance FST between samples of M. galloprovincialis (Mg) relative to the control species (Mch, M. chilensis; Me, M. edulis; Mt, M. trossulus).Figure 2. Principal coordinates analysis on the genetic distance F ST between samples of M. galloprovincialis (Mg) relative to the control species (Mch, M. chilensis; Me, M. edulis; Mt, M. trossulus).
., North Atlantic Spanish (MgRi and MgSa), South Africa (MgCt), Chile (MgDi) and California (HgtMb and HgtRf).The most heterogeneous pool composition was observed in sample MgRi from Galicia (Spain) and in the two Californian samples.The other four pools from Morocco, South Africa, Australia and Japan showed less admixture.

Figure 3 .
Figure 3. BAPS posterior probability of samples of belonging to one of the five gene pools within M. galloprovincialis (Mg), to M. trossulus (Mt, dark green), to M. chilensis (Mch, dark blue) and to M. edulis (Me, yellow).Only the significant (alpha = 0.05) admixture estimates are shown.The highest divergence between the k = 8 BAPS pools using the Nei genetic distance was observed between M. galloprovincialis and the rest of the species, i.e., M. trossulus (1.152), M. edulis (0.657) and M. chilensis (0.506).The M. galloprovincialis samples formed a clade that was the sister group of the rest of the species.The clade of M. galloprovincialis included two subclades, one grouping the East Asia samples (Japan and Australia) and the other joining all the Atlantic North samples of Iberia and Morocco and including South Africa, Chile and California (Figure4).

Figure 3 .
Figure 3. BAPS posterior probability of samples of belonging to one of the five gene pools within M. galloprovincialis (Mg), to M. trossulus (Mt, dark green), to M. chilensis (Mch, dark blue) and to M. edulis (Me, yellow).Only the significant (alpha = 0.05) admixture estimates are shown.

Figure 3 .
Figure 3. BAPS posterior probability of samples of belonging to one of the five gene pools within M. galloprovincialis (Mg), to M. trossulus (Mt, dark green), to M. chilensis (Mch, dark blue) and to M. edulis (Me, yellow).Only the significant (alpha = 0.05) admixture estimates are shown.

Figure 4 .
Figure 4. BAPS reconstruction relating eight significant sample clusters upon a neighbor-joining dendrogram based on the Nei's distance (averaged over loci).

Table 2 .
Hierarchical analysis of molecular variance (AMOVA) scored at different geographic and taxonomic levels in Mytilus spp.Asterisks indicate the probability, based on 1023 permutations, that the observed values were equal to or smaller than that expected by random is p ≤ 0.01; ns: not significant.