Molecular genetic differentiation of native populations of Mediterranean blue mussels, Mytilus galloprovincialis Lamarck, 1819, and the relationship with environmental variables

Abstract Blue mussels of the genus Mytilus are important ecosystem engineers in intertidal and subtidal communities. The distribution of Mytilus mussels is influenced by a number of benthic and pelagic environmental variables (e.g., substratum type and availability, water movement, phytoplankton production, physical disturbance) as well as interactions between these variables. Because of its broad tolerance of environmental variation the Mediterranean species, Mytilus galloprovincialis, has the greatest ability of all blue mussels to colonise new geographic regions. Understanding how population genetic variation is related to, or caused by, environmental variation is important but has long been a challenge. The present study examined the genetic differentiation of native populations of M. galloprovincialis throughout its entire geographic range in the Mediterranean Sea, the Black Sea and the Sea of Azov using 53 single nucleotide polymorphisms (SNP loci). Mussels, in total 1004 individuals collected from 36 locations, were genotyped and combined with existing SNP data for mussels from 11 reference sites. Pairwise comparisons of F ST values, correspondence analysis (CA) and STRUCTURE analysis all revealed four groups of populations: the Atlantic Ocean; the western Mediterranean; the Aegean Sea; and the Azov, Black and Marmara Seas. One population – from Algeria (Oran West) – was intermediate between the two main groups of the Mediterranean Sea and Atlantic Ocean. Seascape genetic analyses using GLM and DistLM analyses were employed to test site-specific genetic variation as a function of 13 environmental variables. The GLM identified five environmental variables that explained variation in site-specific F ST values, whilst in the DistLM best-fit model only four were significant. These analyses suggest that a complex mix of environmental variables contribute to explaining genetic variation of M. galloprovincialis populations within the Mediterranean Sea, which most likely reflects the complex geological history of formation, isolation and reconnection among the regional sub-basins of the Sea.


Introduction
The Mediterranean Sea has a complex history of formation. Its origin dates to the Mesozoic and the Tethys Sea, when it acquired a shape similar to the present day in the Miocene (Picotti et al. 2014). At the end of the Miocene, the Mediterranean Sea was cut off from the Atlantic Ocean by a land barrier between Iberia and North Africa, which caused the water level to drop due to evaporation, which in turn gave rise to the Messinian crisis 5.75-5.32 milion years before present (M ybp) (Rouchy & Caruso 2006;Janssen & Peijnenburg 2014). The connection to the Atlantic Ocean was re-established 5.3 M ybp when seawater filled the Mediterranean trough (Krijgsman et al. 1999(Krijgsman et al. , 2018. Subsequently, during the Pliocene, a monsoon climate period was followed by cyclical sea-level changes related to periods of glaciation in the Northern hemisphere (Barsotti & Meluzzi 1968;Picotti et al. 2014). Other major regions, sometimes considered to be a basin or sub-basin of the Mediterranean Sea, also have complex and sometimes completely independent geological histories. For example, the Black Sea, prior to the establishment of a contemporary connection to the Mediterranean Sea through the Bosphorus Strait, the Sea of Marmara (including the Dardanelle Strait) that was established 10.5 K ybp, was a freshwater lake receiving water from melting glaciers. Earlier, the Black Sea was only periodically influenced by Mediterranean Sea waters (Ryan et al. 2014). These complex geological histories and the present-day multi-basin structure of the Mediterranean Sea present an opportunity to better understand how periods of isolation and then reconnection, plus the different ages of the basins and their very different environmental conditions may influence genetic connectivity and population genetic diversity.
Given the complex history of formation of the Mediterranean Sea and its sub-basins it is not surprising that the Mediterranean Sea is a biodiversity hotspot (Coll et al. 2010;Danovaro et al. 2010;Pascual et al. 2017). Two main regions of biodiversity of the Mediterranean terrestrial fauna and flora have evolved, with different geological histories: the West Iberian and the North African separated by the Strait of Gibraltar to the west, and the Balkan and Anatolian to the east. These regions differ strongly in taxonomic composition and population divergencies of land organisms as a result of isolation in older (a few M ybp) and more recent glacial refugia (Sanmartín 2003;Froufe et al. 2016). In contrast to land and freshwater organisms, marine species, including benthic sedentary organisms such as mussels, have the ability to disperse via their pelagic stages of development (larvae) and/or colonise new areas as adults on drifting items (e.g., kelp). Such dispersal potential, when realised, may mean that their populations are characterised by high connectivity and low spatial genetic differentiation (e.g., Lessios et al. 1998;Addison et al. 2008;Reisser et al. 2014). Consequently, such marine organisms may spread and reconnect isolated populations with the result that present-day marine populations in the Mediterranean basin may be only weakly genetically differentiated.
The native smooth-shelled blue mussel, Mytilus galloprovincialis, is a widely distributed species throughout the Mediterranean Sea. The species is an important ecosystem engineer in the intertidal and subtidal hard-bottom communities (Borthagaray & Carranza 2007;Arribas et al. 2014;Çinar et al. 2020), where the sessile adults are attached by byssus threads to hard substrata but the larvae may travel large distances in the water column. In addition, the mixing of Mediterranean populations may have been caused by a long history of shipping and trade activity associated with the sea (Blondel et al. 1999;Coll et al. 2010). For example, marine organisms including the mussel Mytilus have been used as food since the early Paleolithic (ca. 50 K ybp) in present-day southern Spain and Gibraltar (Cortés-Sánchez et al. 2019), the Neolithic in the Algarve, Portugal (Callapez & Pimentel 2018) and the Middle Bronze Age in present-day Italy (Zedda et al. 2003;Minniti 2005), as evidenced by preserved mussel shells in middens in caves and the vicinity of settlements. A simple form of mussel culture was initiated in Italy over 2 K ybp (Mattei & Pellizatto 1977;Smaal 2002). More recent intensive shipping activity, translocations of hatchery stocks and rafting on natural (e.g., kelp) and/or man-made (e.g., plastics) materials may have contributed to possible reduction of genetic differentiation among populations within basins of the Mediterranean Sea (Giantsis et al. 2014a).
The progenitor of modern smooth-shelled blue mussels gave rise to Mytilus trossulus in the Pacific Ocean, which subsequently invaded the North Atlantic Ocean after the opening of the Bering Strait 3.5 M ybp and gave rise to Atlantic Ocean M. edulis (Riginos & Cunningham 2005;Bach et al. 20192019; reviewed by Gardner et al. 2021). This newly evolved taxon spread in the Northwest Atlantic Ocean to reach Europe. The oldest fossils 756 R. Wenne et al. of M. edulis are found in deposits from the early Pliocene, its shells having been reported from middle and late Zanclean deposits in the vicinity of Lepe, on the Bay of Cadiz, Atlantic coast NW of Gibraltar (Belaústegui & Muñiz Guinea 2016). M. edulis colonised the Mediterranean Sea during periods of sea-level change before the Pleistocene (Vermeij 1992;Gardner et al. 2021). Glaciation in the Northern hemisphere took place 2.6 M ybp (Sosdian & Rosenthal 2009), likely resulting in a drop of sea level by 100 m. Divergence of M. galloprovincialis in the Mediterranean Sea from M. edulis in the North Atlantic Ocean probably occurred 2.5 M ybp and was followed by isolation lasting ~1.8 M y (Roux et al. 2014). Periods of glaciation forced the displacement of northern populations of boreal and colder-water species into southern Europe and resulted in the partial invasion of the Mediterranean Sea as a refugium, for example, by the bivalve Arctica islandica in the northern Adriatic Sea 1.8 M ybp (Crippa et al. 2016). During the Early-Middle Pleistocene Transition of 1.4 to 0.4 M ybp more regular cycles of glacial-interglacial periods became established and caused further changes to sea level. The global drop in sea level isolated the Mediterranea Sea biota, including its mussels, from the Northeast Atlantic Ocean biota. This scenario is supported by Riginos and Cunningham (2005) whose genetic analyses highlighted that isolation mechanisms (vicariance) were important for divergence of M. galloprovincialis from neighbouring M. edulis. It is possible that during such sea level changes the Atlantic Ocean mussel populations came into secondary contact with the Mediterranean Sea mussel populations many times in the period of (incomplete) separation that lasted ~1.7 M years (Barsotti & Meluzzi 1968). After colonisation by Mytilus, the influx of waters from the Atlantic Ocean to the present-day Mediterranean Sea was limited or cut off periodically due to tectonic changes, but mainly due to fluctuations in sea level to the west of the Strait of Gibraltar. In addition, strong evaporation of seawater caused the periodic lowering of the water level in the Mediterranean Sea and gave rise to the isolation of parts of it as separate basins during the Pleistocene (Bianchi et al. 2012). In the last 150,000 years, changes in sea level of ~130 m below to 6-15 m above the present sea level have been caused by glaciations and warming periods (Benjamin et al. 2017). This complex and dynamic geological history of the Mediterranean Sea has given rise to the allopatric divergence of Mytilus populations and the evolution of M. galloprovincialis. In addition, natural hybrid zones are created in the areas where modern populations of M. galloprovincialis come into secondary contact with M. edulis (Bierne et al. 2003;Simon et al. 2020Simon et al. , 2021. The ability of Mytilus mussels to spread on natural and artificial floating objects has contributed to their emergence in recent years in regions as far apart as Svalbard in the Arctic -M. edulis (Berge et al. 2005;Kotwicki et al. 2021) and the South Shetland Islands off Antarctica -M. platensis (Cárdenas et al. 2020). However, of all of the members of the genus it is M. galloprovincialis that shows the greatest ability to colonise new geographic regions, whether via human-mediated transfer (e.g., deliberately for aquaculture or accidentally via hull fouling or ballast water) or naturally (e.g., via kelp rafting) (Gardner et al. 2021). In recent years, it has been recorded on the Atlantic coasts of South America, where in Argentina it hybridises with the native taxon, M. platensis , and in Brazil where it inhabits farmed cultures of the native species, Perna perna (Birckolz et al. 2020;Lins et al. 2021). M. galloprovincialis has also been introduced to South Africa, the Pacific coast of North America, the Sea of Japan, China, Korea, Australia, New Zealand and Chile (Wilkins et al. 1983;McDonald & Koehn 1988;McDonald et al. 1991;Han et al. 2017;Larraín et al. 2018;Zbawicka et al. 2019Zbawicka et al. , 2022Popovic et al. 2020). It is considered to be one of the most successful invasive species of the global coastal marine biota, often displacing native species (Geller 1999;Lowe et al. 2000;Gardner et al. 2021).
The distribution of Mytilus species is controlled by a number of environmental processes involving both benthic and pelagic habitats (e.g., substratum type and availability, water movement, phytoplankton production, physical disturbance, seawater biogeochemistry) as well as interactions between these processes (Sandman et al. 2013;Kotta et al. 2015). For example, water movement can indirectly affect Mytilus spp. by modifying sedimentation rates or may directly affect sessile mussels by physically disturbing or detaching them (Westerbom et al. 2008). Furthermore, the sedentary benthic suspension feeding lifestyle of adult mussels is a life-history characteristic that requires water-borne food delivery, but the quantity and quality of such food delivery is highly variable in time and space (Dahlhoff & Menge 1996;Gardner 2000Gardner , 2013Saurel et al. 2007). An ability to withstand environmental variability is often key to a species' range expansions, whether they be natural or human-mediated. Environmental tolerance can be related to or controlled by such mechanisms of adaptation as differences in gene family expression, gene splicing, methylation, past gene duplications and other genomic mechanisms also related to pangenome functions (Pujolar et al. 2014;Malachowicz et al. 2015;Kijewska et al. 2016Kijewska et al. , 2018Bitter et al. 2019;Malachowicz & Wenne 2019;Clark et al. 2021;Corrochano-Fraile et al. 2022;Liu et al. 2022). For example, there is now a growing body of evidence to indicate that M. galloprovincialis is the best adapted smooth-shelled mussel to tolerate thermal variation, an adaptation that may be both facilitating its spread into non-native areas and its ability to outcompete congenerics (e.g., Braby & Somero 2006;Jones et al. 2010;Tomanek & Zuzow 2010). When seeking to understand associations between environmental and population genetic variation it is important to incorporate ecologically meaningful drivers into the models (e.g., direct and indirect environmental gradients, resource gradients). However, trying to identify environmental drivers that are relevant at different temporal and spatial scales is often difficult because the drivers may not act in isolation but may instead act synergistically, although not always to the same extent in time (e.g., different seasonal interactions) or in space (e.g., coastal headlands versus bays). Thus, quantifying interactive effects (i.e., how different environmental gradients separately or interactively modulate the resource-allele relationship) is still a major challenge (Palumbi et al. 2019;Hu et al. 2020;Zeng et al. 2020).
To date, genetic studies of native M. galloprovincialis populations in the Mediterranean Sea are local in nature, having covered only certain regions (e.g., Ahmad & Beardmore 1976;Edwards & Skibinski 1987;Karakousis & Skibinski 1992;Quesada et al. 1995aQuesada et al. , 1995bComesaña et al. 1998;Daguin & Borsa 1999;Giantsis et al. 2012Giantsis et al. , 2014aGiantsis et al. , 2014bBierne et al. 2003;Sammer et al. 2010;Vera et al. 2010;Lourenço et al. 2015;Fraïsse et al. 2016;Paterno et al. 2019). No genetic research has yet been undertaken at the pan-Mediterranean Sea scale. The aim of the present study was to test for genetic differentiation of M. galloprovincialis populations throughout its entire native range. We also tested if and how much environmental variables explained the population genetic diversity of M. galloprovincialis. By including key environmental variables in our testing we seek to explore links between patterns of environmental and genetic variability of the mussel populations. A strong match between genetic and environmental variation and the presence of clear genetic differentiation according to present day environmental gradients may indicate the importance of local adaptation processes in Mediterranean Sea M. galloprovincialis or barriers to gene flow. Alternatively, a relationship between environmental variation and subtle regional genetic differentiation may reflect the complex geological history of the region and the evolutionary history of the mussel populations.

The environment
The Mediterranean Sea is a land-locked, semienclosed marginal sea of the Atlantic Ocean. The Mediterranean Sea is a very complex marine environment and is oceanographically diverse with a number of distinct sub-seas. Sea surface temperature exhibits strong seasonality and pronounced spatial gradients from west to east and north to south. The Mediterranean Sea is microtidal with a typical tidal range of less than 50 cm. The salinity of the Mediterranean Sea is high throughout the basin from 35 PSU on the West to 39 PSU on the East, but is brackish in the Black Sea and Sea of Azov. The Mediterranean Sea has very low nutrient concentrations, especially in its eastern parts, which is reflected in the patterns of primary productivity. Locally, riverine nutrient input can be significant; however, most river systems discharging into the Mediterranean Sea are small by global standards (Uitz et al. 2012;Goffredo & Dubinsky 2014).

Sampling and isolation of genomic DNA from mussel samples
Mytilus spp. samples, consisting in total of 1004 individuals of mixed ages (not quantified) and sizes (5 to 50 mm shell length) were collected from 36 sites (nominally "populations") situated along the coasts of the Atlantic Ocean, the Mediterranean Sea and the Black Sea between 2004 and 2016 (Table I; Figure 1). DNA was isolated from mantle tissue that had been stored in 96% ethanol or at −70°C, using a modified CTAB method (Hoarau et al. 2002). Eleven reference site samples of M. galloprovincialis from the Atlantic Ocean coast of Spain and Portugal, the Mediterranean Sea and the Sea of Azov, M. edulis from France and M. trossulus from Canada, consisting of 316 specimens, were also included (Table I). Fifty three single nucleotide polymorphisms (SNPs) differentiating Mytilus taxa, and with the ability to identify hybridisation, were employed to assay genetic variation (Zbawicka et al. 2012Gardner et al. 2016). Samples were genotyped using the Sequenom MassARRAY iPLEX genotyping platform (Gabriel et al. 2009).

758
R. Wenne et al.  (Tamura et al. 2013). Two methods were used for the population structure analyses. First, correspondence analysis (CA; Benzécri 1992), implemented in GENETIX (Belkhir et al. 2003), was used to visualise genetic structure among populations. The results are presented as a scatter plot, with the axes representing the contribution of inertia of the data matrix in a way that can be considered analogous to the total variance in allelic frequencies (Benzécri 1992). Second, clustering and assignment testing were performed using the Bayesian-based method implemented in STRUCTURE v. 2.3.4. STRUCTURE was employed using the model assuming admixture, ignoring population affiliation and allowing for the correlation of allelic frequencies between clusters. This admixture model allows for individual structure with mixed ancestry, meaning that fractions of the genome could have come from different ancestors (Pritchard et al. 2000;Falush et al. 2007). The most appropriate number of genetic clusters (K) was determined by a diagram-based comparison of log-likelihoods for values of K. At least five runs were used to calculate each K value, following the method described by Evanno et al.

760
R. Wenne et al. (2005). Threshold q-values of 0.2 were used as a criterion to separate hybrid and pure mussels (Vähä & Primmer 2006). Individuals were considered residents if q > 0.8 for the area where they were sampled. Individuals with q-values from 0.2 to 0.8 were considered to be potentially admixed, as they could not be readily assigned as residents or migrants (Lecis et al. 2006). Individuals with an assignment probability of q > 0.8 were defined as belonging to the wild population (cf. Jonker et al. 2013), whilst those with q ≤ 0.8 were labelled as admixed. A Monte Carlo Markov Chain was run for 100,000 iterations following a burn-in period of 50,000 iterations.
Linking genetic and environmental variation. Regionalscale proxies for weather and local environmental conditions were obtained from online databases for 13 environmental variables (Table S1). As a proxy for the global climate and weather (air and water temperature, cloud, wind speed, waves) we used ERA5 reanalyses (for all methodological details of environmental variables and modelling see references in Table S1).
Reanalysis combines model data with observations from across the world into a globally complete and consistent dataset. ERA5 provides hourly estimates for a large number of atmospheric, ocean-wave and land-surface variables. To characterise biogeochemical patterns (concentrations of chlorophyll a, nitrate and phosphate, as well as marine primary production) we used the CMEMS global biogeochemical multi-year hindcast product GLOBAL REANALYSIS BIO 001 029. CMEMS provides daily estimates for all biogeochemical variables. Model inputs for the physical conditions (salinity, mixed layer thickness) were obtained from the CMEMS global ocean eddy-resolving reanalysis product GLOBAL REANALYSIS PHY-001-030. Data for the global distribution of photosynthetically available radiation at the sea surface and on the seafloor were obtained from Gattuso et al. (2020). Relationships between environmental variables and genetic variation of M. galloprovincialis populations from the Mediterranean Sea (i.e., excluding populations CHA, LOI, MSMA and KKAT) were tested using two complementary seascape genetics approaches. First, following Wei et al. (2013), we  Table I. used population-specific F ST values (the response variable and where negative F ST values were set to zero) and the 13 environmental variables in a GLM (the GLZ routine in Statistica v. 12). We employed a normal distribution and a log link function. To minimise Type I error that may be associated with stepwise (forward or backward) model building we employed the "all effects" model. Second, following Silva and Gardner (2016), we used population-specific allelic frequencies (the response variables) and the 13 environmental independent variables in a distance-based linear model (DistLM in the PRIMER + PERMANOVA v. 6 software package - Anderson et al. 2008). This test is a permutational equivalent to partial redundancy analysis (Legendre & Anderson 1999). DistLM was used to perform an ordination of fitted values from a given model and is constrained to find linear combinations of predictor variables (environmental data) that explain the greatest variation in the data cloud (population-specific allele frequencies). Permutation of residuals was carried out under a reduced (or partial) model and because this is a permutational test, there are no assumptions about data normality (Anderson et al. 2008). Marginal tests (one independent variable at a time) and sequential tests (all independent variables entered into the model based on their relative importance (most significant first) in the marginal tests) were employed to identify the environmental variables that explained the greatest variation in the genetic dataset. Model fit was tested using adjusted R 2 (i.e., adjusted for the number of terms in the model), the AIC value and the BIC value. Note that these two seascape genetic analyses use the same environmental dataset but test for population-specific variation in different dependent variables: for the GLM this is population-specific F ST values (a summary metric of population differentiation) whereas for the DistLM this is locus-specific allele frequencies (i.e., the raw genetic data).

SNP validation and Hardy-Weinberg equilibrium
Fifty-three SNPs were genotyped (Table II). ORF identification was not possible for only five SNPs. Forty-four SNPs (83.01%) were located in coding regions, of which only three (5.66%) were not synonymous. Four SNPs (7.54%) were located in non-coding regions. The MAF ranged from 0.000 (7 different loci) to 0.404 (1 locus) ( Table II). The vast majority of loci were in Hardy-Weinberg equilibrium (HWE) in all populations.
No SNPs with departures from HWE were observed in 10 populations (AGA, BLS, CHW, CIRP, MOM, ORAW, SBRB, TURK and reference CAS and LOI). Only one population (VAL) had three SNPs that were not in HWE (P < 0.01). In the remaining samples the fraction of SNPs showing significant departures from HWE was one or two for 20 and 14 populations, respectively (Table III).

Genetic diversity
The percentage of polymorphic loci (Po) in all M. galloprovincialis populations ranged from 45% in the Morocco (AGA) and Turkey (TURK) samples to 56.60% in the Croatia (DUB) and reference IMC (Italy) samples (Table III). Observed heterozygosity (H O ) for 53 loci across 47 populations was lower than expected (H E ), except for two samples: ORAE and BLS (Table III) (Table III),

762
R. Wenne et al.  which may indicate a relationship between individuals within a population, resulting from the collection of related individuals from a small area. Estimates of average pairwise differences withinpopulation diversity among the M. galloprovincialis populations revealed that the most diverse population was BLT from Tunisia and the least diverse was CHW from Croatia. In the reference samples the greatest diversity was observed in the Atlantic M. galloprovincialis BID population from Spain and the lowest in the M. edulis MSMA population from France. The same samples exhibited the greatest and the least gene diversity (Table III).
Allelic frequencies were calculated for 53 SNPs and minor allele frequency (MAF) was determined for 51 bi-allelic SNPs (Table S2) AMOVA was performed comparing groups of samples for five different scenarios where populations were defined a priori (details in Table IV).
The estimated values of the F-statistic were significant for all five scenarios, and the greatest variance was exhibited within individuals for all scenarios.  populations from the Adriatic Sea region formed a separate, but not well supported, sub-group. The remaining M. galloprovincialis populations fell outside these two main clusters with little or no evidence of groupings based on geography. Correspondence analysis (CA) of 43 M. galloprovincialis populations (of which 7 are reference populations) revealed clear separation of an Atlantic Ocean grouping, a Black Sea plus Sea of Azov grouping, and a third grouping consisting of all other populations, except the ORAW population that was located between the Atlantic Ocean and the Mediterranean Sea groups (Figure 3). After removal of Atlantic Ocean and African samples a CA carried out on 33 populations and including 3 reference populations (LID, ORI, AZO) resolved three groups: the M. galloprovincialis from the Black Sea and the Sea of Azov, M. galloprovincialis from the Mediterranean Sea, and M. galloprovincialis from the Aegean Sea ( Figure 4). STRUCTURE analysis of 46 Mytilus populations (i.e., excluding the reference M. trossulus) revealed that the largest increase of LnP(D) was for K = 2 and then for K = 3. The greatest subdivision was detected at K = 2 where the clusters corresponded to the separation of the M. edulis and M. galloprovincialis populations ( Figure 5). At K = 3, the clusters corresponded to M. edulis, and two groups of M. galloprovincialis, with a general division into lineages from the Mediterranean Sea and the Atlantic Ocean because individuals are potentially admixed. Division between the M. galloprovincialis groups was most apparent at the population level because the frequency of one of the clusters characteristic of M. galloprovincialis was >77% (average 80%) in Atlantic Ocean populations, whereas the frequency of the second cluster was >50% (average 61%) in M. galloprovincialis from the Mediterranean Sea populations (Fig. S3). From the STRUCTURE analysis when K = 3 most of the 1292 individuals (97.2%) were correctly assigned to their original sample taxon (one of the three clusters) with q > 0.8. Mussels collected from Morocco (AGA, Atlantic Ocean) and Algeria (ORAW, Mediterranean Sea) showed gene admixture characteristic of M. galloprovincialis reference groups from the Atlantic Ocean (BID, VIG, CAS and CAM). A second sample from Algeria (ORAE), located east of the Alboran Front, showed gene admixture characteristic of M. galloprovincialis populations from the Mediterranean Sea.

Linking genetic and environmental variation
The GLM testing of the "All effects" model was significant (P < 0.05), with Mixed Layer thickness (P < 0.0001), PAR at surface (P < 0.004), PAR on seafloor (P < 0.011), Cloud cover (P < 0.03) and SST (P < 0.036) as significant terms. The other eight environmental variables did not explain significant variation in the population-specific F ST values (P > 0.05). For the DistLM analysis the best-fit models built using adjusted R 2 , AIC and BIC gave the same results so we report here only the results for AIC. The marginal test results (Table V) indicated that seven of the 13 environmental variables explained significant variation in population-specific SNP locus allele frequencies. These seven variables explained a total of 75% of the variation in the SNP dataset. In the sequential testing the best-fit model contained all 11 environmental variables, of which five variables were statistically significant (Wave Ht, PO 4, SST, PAR at surface, Mixed layer) in a 13term model that explained 61.4% of the variation in the SNP dataset (Table V).

Discussion
The Mediterranean Sea is characterised by a high biodiversity of marine organisms and a large number of endemic species (Coll et al. 2010;Danovaro et al. 2010). Despite the connection with the Atlantic Ocean via the Strait of Gibraltar and the possibility of transport of planktonic larvae and motile adults, many species show genetic distinctiveness of their Mediterranean Sea populations (Patarnello et al. 2007;Pascual et al. 2017). These species include the cirriped Chthamalus montagui (Shemesh et al. 2009;Pannacciulli et al. 2017), the shrimp Palaemon elegans (Reuschel et al. 2010) and the blue mussel Mytilus galloprovincialis (Varvio et al. 1988;Quesada et al. 1995b;Kijewski et al. 2011;Zbawicka et al. 2012;Del Rio-Lavín et al. 2022).
The location of the Almeria-Oran Front is the area where two circular sea currents meet and flow from north to south-east into the Mediterranean Sea and then south-west, which directs water of Atlantic Ocean origin into the Alboran Sea. This is considered to be an isolating factor (Millot & Taupier-Letage 2005;Millot 2013;Pascual et al. 2017) for some marine species between the Atlantic Ocean proper and the Mediterranean Sea. Abrupt changes of allele frequencies in M. galloprovincialis populations related to this front have been reported for a range of difference genetic marker types, including allozymes (Sanjuan et al. 1994;Quesada et al. 1995a), mtDNA (Quesada et al. 1995b), microsatellites (Diz & Presa 2008;Ouagajjou et al. 2010;Ouagajjou & Presa 2015)    In the present study, significant differentiation between mussel populations of the Sea of Azov and the Black Sea from the Mediterranean Sea populations was identified with two of the 53 SNP loci. No statistically significant differences were observed among populations from the Black Sea and the Sea of Azov and the Dardanelle Straits. Most genotypes were common to the Mediterranean, Black and Azov Seas and what differences existed are attributable to genotype frequency differences among the populations. A few SNP loci with alleles rare in the Black Sea and Sea of Azov populations were more frequent in the Mediterranean Sea populations (e.g., BM102A "C", BM105A "G", BM26B "T"), or those that did not occur in the Black Sea and Sea of Azov, for example, BM106B "G". Only two genotypes were more frequent in the Sea of Azov than in Black Sea populations: BM147A "C" and BM32A "A". Significant differentiation between mussel populations from the Mediterranean and Black Seas has also been found using as many as 512 SNP loci (Paterno et al. 2019) although no significant genetic structuring was noted within the
Black Sea populations using genotyping of 998 SNP loci (Paterno et al. 2019). Genetic divergence has also been reported for M. galloprovincialis from Mediterranean Sea and Black Sea populations at the mtDNA level (Ladoukakis et al. 2002). RFLP analysis of COIII gene polymorphism in M. galloprovincialis populations revealed that almost all Black Sea haplotypes also occur in the Mediterranean Sea populations, that is, the Black Sea populations are a sub-set of the Mediterranean Sea populations. In addition, no significant differences in allele frequencies have been reported between populations of M. galloprovincialis from the southern Black Sea, Bosphorus Strait and Sea of Marmara by assessment of variation in the COIII mtDNA region and six microsatellite loci (Kalkan et al. 2011). It has been suggested that a shift from Black Sea-like to Mediterranean Sea-like genetic structure occurs at a location in the Dardanelle Strait, which is supported by our SNP data. The most common mtDNA haplotype identified by RFLP analysis of ND2-COIII in a population from the Sea of Azov was also present in a sample from Villafranche-sur-Mer, northwestern Mediterranean Sea (Śmietanka et al. 2004). However, a few rare haplotypes were unique. Genotyping-by-sequencing of the same region of mtDNA revealed haplotype frequency differences between populations from the Sea of Azov and Black Sea, but most common haplotypes were also present in the northwestern Mediterranean Sea (Gerona, Banylus sur Mer and Gulf of Oristano) (Śmietanka et al. 2009, 2014). Because the timing of onset of the flow of Mediterranean Sea waters into and possible colonisation by euryhaline bivalve populations of the Black Sea has not been defined precisely (Nikula & Väinölä 2003;Sromek et al. 2019), we have reviewed in detail here the information available. It is estimated that Mediterranean Sea water from the Aegean Sea filled the Sea of Marmara and replaced freshwater originating from the Black Sea approximately 12 K ybp (Algan et al. 2001). The outflow of lacustrine water from the present area of the Black Sea probably occurred between 10 to 8.4 K ybp, which was followed by pulses of Mediterranean Sea seawater entering the Black Sea (Hiscott et al. 2007). Connection of the Black Sea to the Mediterranean Sea may have taken place 9.4 k ybp, as indicated by strontium/oxygen ratio changes in molluscan shells (Major et al. 2006). Transition of the Black Sea from a freshwater lake to a marine environment most probably occurred between 8.9 and 8.5 K ybp (Tyuleneva et al. 2014). Outflow of freshwater or isolation lasted until 8 K ybp (Soulet et al. 2011) and was followed by establishment of today's hydrographic system with surface outflow of low salinity Black Sea water and inflow at depth of Mediterranean Sea high salinity waters (Kokkos & Sylaios 2016) via the Bosphorus Strait (Hiscott et al. 2007).
The first representatives of the Mediterranean mollusc fauna appeared in deposits in the Bosphorus region 5.3 K ybp, whilst Mytilus sp. appeared for the first time ~4.4 K ybp (Algan et al. 2001). In a core from the southwest Black Sea, M. galloprovincialis has been found in layers dated to 5.9-2.4 K ybp (Hiscott et al. 2007). Filipova-Marinova et al. (2013) observed the occurrence of a M. galloprovincialis shell layer in deposits of Varna Lake (Bulgarian Black Sea) that were radiocarbon dated to 7,776 to 6,183 ybp. According to Tyuleneva et al. (2014) timing of the arrival of M. galloprovincialis in the northwestern Black Sea (estimated without 14 C dating) was from the Bugazian at 10.5-8.4 K ybp, Kalamitian beds at 7.1 to 4 K ybp and Dzhemetinian beds at 4.1 K ybp to the present in sediments containing marine euryhaline species. A discrepancy in time estimation of the earliest appearance of Mytilus in Black Sea deposits may be related to elution and displacement of some sediment layers. It can be assumed that Mytilus populations settled the Black Sea sometime since 8 K ybp as a result of post-glacial expansion from refugia.
Time of divergence of mtDNA haplotypes in the Black Sea and Mediterranean Sea mussel populations has been estimated as a few hundred thousand years ago (Śmietanka et al. 2014), which is much earlier than the origin of the mussel population in the Black Sea. Therefore, euryhaline M. galloprovincialis populations evolved in the Mediterranean Sea and the Aegean Sea (possibly including the Levantine Sea) long before the contemporary marine phase of the Black Sea. Similar to the Black Sea situation, the Baltic Sea in northern Europe was a freshwater lake and became connected with the Atlantic Ocean and newly developing North Sea approximately 8-9 ybp (Berglund et al. 2005;Behre 2007). For comparison, Kostecki and Janczak-Kostecka (2011) have reported the onset of the marine water environment in Pomeranian Bay, southwest Baltic Sea, to be 8.9-8.3 ybp as determined using geochemical estimators. The Baltic Sea is characterised by lower marine species diversity in comparison to the present-day North Sea (Johannesson & Andre 2006;Wennerström et al. 2017). Genetic divergence of Mytilus populations in both regions (i.e., the Black Sea from the  Greek coasts of the Aegean Sea revealed that these populations were homogeneous (Ladoukakis et al. 2002;Giantsis et al. 2014b). Moreover, Giantsis et al. (2014b) reported that the population from the area closest to the Dardanelle Strait (Sea of Marmara, Turkey) differed significantly from the Aegean Sea populations, whereas differences between samples from the Aegean, Ionian and North Adriatic Seas were weak with the exception of one sample from the area of Zadar (Croatia). MtDNA sequence analyses revealed genetic homogeneity among all Greek populations and the clear differentiation of the only Turkish sample (Çanakkale, Dardanelle Strait, Sea of Marmara) from the Aegean populations (Giantsis et al. 2014b). On the other hand, microsatellite data suggested significant differentiation of Italian samples (from the North Adriatic Sea and Ligurian Sea, respectively) from Aegean samples (Giantsis et al. 2014a). Genotyping with a large number of SNPs in a more local study revealed population variation in the northern and southern Adriatic Sea different from that in the northern Ionian Sea (Paterno et al. 2019). The results of the present study are in agreement with the findings of these different studies and extend the microstalite DNA analyses (Giantsis et al. 2014a) concerning the genetic  Table I. Vertical black lines separate the samples.  (Mejri et al. 2011), the sea star Astropecten aranciacus (Zulliger et al. 2009) and the cockle Cerastoderma glaucum (Nikula & Väinölä 2003). We did not observe significant genetic differentiation among populations of M. galloprovincialis in the Western Mediterranean Basin, east of the Alboran Front to the coasts of Italy and Tunisia. Similarly, analysis of samples from Banyuls (France) and Haouaria (Tunisia) using 512 SNP loci did not identify differences (Paterno et al. 2019). Lack of genetic differentiation among M. galloprovincialis populations within the Western Basin can be partly explained by natural reasons, such as extensive gene flow and a common paleogeographic history of this Basin including the Balearic, Ligurian and Tyrrenian Seas separated from the Eastern Basin by the Siculo-Tunisian Strait. After each glacial phase since the Pleistocene (cycles of 41 K years prior to and 100 K years after the Mid Pleistocene Transition), M. galloprovincialis populations expanded into newly flooded areas as a result of sea level rising and then retreated with sea level falling (repeated interglacial colonisation and expansion, and retreat during glaciations). However, the biggest drop in sea level occurred 24-19 K ybp during the Last Glacial Maximum (Boavida et al. 2019). Surface waters in the central Mediterranean Sea (in the vicinity of the Siculo-Tunisian Strait) were most probably much colder in spring during glacial periods in comparison to interglacial periods (Rouis-Zargouni et al. 2010), which suggests much colder winter temperatures in comparison to the present. Strong meltwater discharge also probably influenced cooling (lower temperatures) and increased variability of the salinity of coastal Mediterranean Sea waters. Increase in rainfall during the past interglacial/glaciation periods resulted in a decrease of sea surface salinities e.g., in the northern Tyrrhenian Sea (Dixit et al. 2020) and increased suspended matter concentrations, and intensified the flux of nutrient and organic matter in the Mediterranean Sea causing deposition of sapropels (Toucanne et al. 2015) in both the Eastern and Western basins (Rohling et al. 2015). Thus, a combination of physico-chemical changes in sea water composition may explain genetic differentiation observed in contemporary Western and Eastern Mediterranean Sea populations of M. galloprovincialis. A separation of western and eastern Mediterranean populations by a possible barrier of periodically shallow water in the central Mediterranean Sea (the Siculo-Tunisian Strait) has been indicated by Chefaoui et al. (2017). A divergence between western and eastern Mediterranean Sea populations has also been found for species such as the sea bass Dicentrarchus labrax (Bahri-Sfar et al. 2000), the sea cucumber Holothuria polii (Valente et al. 2015), the crab Carcinus aestuarii (Ragionieri & Schubart 2013), the hermit crab Diogenes pugilator (Almón et al. 2021), and the cockle Cerastoderma glaucum (Sromek et al. 2019).
Seascape genetics seeks to identify associations between environmental variation and genetic variation with the ultimate aim of identifying key environmental factors that contribute to explanation of population genetic variation and regional differences in genetic structure (e.g., Selkoe et al. 2008Selkoe et al. , 2016Riginos & Liggins 2013;Wei et al. 2013;Silva & Gardner 2016;Zeng et al. 2020). Typically this approach, which is most often based on variation of neutral genetic markers, seeks to identify environmental factors that promote (e.g., currents) or retard (e.g., salinity variation acting as a barrier) gene flow. Whilst this new approach to understanding connectivity is both powerful and elegant, it is often limited by data availability -the low numbers of environmental variables that are available for multiple sites within a region of study. In the present study we used 13 environmental variables collated from 43 sites within the Mediterranean Sea. The GLM analysis identified five environmental variables (mixed layer thickness, PAR at surface, PAR at seafloor, SST, Cloud cover) that explained variation in population-specific F ST values. In contrast, the sequential tests for the DistLM analysis best-fit model, which explained 61% of the variation in the raw SNP data set, included all 13 environmental variables, although only four variables (Wave height, PO 4 , SST, PAR at surface) were significant in the multiterm model. Earlier literature suggests that salinity and temperature often drive the generic patterns of hard-bottom intertidal species including Mytilus spp. (Kaiser et al. 2011) and this result was observed for SST but not for salinity (PSU) in the M. galloprovincialis of the Mediterranean Sea. The absence of salinity from our significant results may reflect the fact that the gradient in the Mediterranean Sea is not great, and certainly not as great as in other regions (e.g., in the Baltic Sea -Kijewski et al. 2019). Overall, our results suggest that there is a complex mix of environmental variables that contribute to genetic variation of M. galloprovincialis populations in the Mediterranea Sea, rather than a simple (one or two variable) explanation as has been reported for several other coastal and deep-sea marine invertebrates (e.g., Wei et al. 2013;Silva & Gardner 2016;Zeng et al. 2020). The complexity of the seascape genetics results in the present study may reflect (1) the reasonably large number of environmental variables in our data set (n = 13) and the fact that many different environmental variables are likely to influence gene flow in any given system (i.e., we might expect to detect a complex result simply because we have a complex environmental data set) and/or (2) the complex geological history of the Mediterranean Sea and the formation of its waters and the associated sub-basins of the region leading to complex interactions between environmental and genetic variation that may be site-specific or regional rather than basinwide. In addition, as noted above, seascape genetic analyses usually focus on detecting associations across neutral loci. Many of the SNPs loci employed here are in coding regions even if they do not exhibit outlier status (refer to Figs. S1 and S2). Thus, the detection of an association between genes under (low) selective pressure and environmental variation cannot be ruled out. Elsewhere, Sun and Hedgecock (2017) have highlighted for high gene-flow species (M. galloprovincialis falls into this group) the need to better understand the role that temporal genetic variation may play in contributing to seascape genetics analysis results from a single (snap-shot) version to a temporal sampling series. Finally, it is worth noting that the GLM and DistLM analyses are both linear-based methods that do not include interaction terms. Other analytical approaches, such as boosted regression trees (Elith et al. , 2008Leathwick et al. 2006;Hastie et al. 2009;Kotta et al. 2017) may prove to be more informative because they are not constrained to detect only linear relationships and can examine interactions among variables, but they require larger (training) data sets than are often available to most researchers and to us for the present study.
The contemporary genetic structure of the Mediterranean Sea populations of M. galloprovincialis is the result of a combination of natural and anthropogenic factors. The dispersal ability of the species is expected to result in high gene flow and connectivity among populations and lead to genetic homogenisation on a broad spatial scale. Although this could explain the homogenisation observed within individual basins, it cannot explain the comparative homogenisation found among the different basins. For instance, based on the topography and oceanographic conditions of the Siculo-Tunisian Strait in the central Mediterranean Sea, one would expect to observe (at least a degree of) genetic differentiation between samples from the Adriatic Sea and the Ligurian Sea, which this study did not find. Anthropogenic activities, including hull fouling, transport of ballast water, the movements of exploration or drilling rigs and unrecorded humanmediated transplantation of spat for aquaculture may all have played an important role in overcoming natural barriers to mussel connectivity and consequently may have contributed to the shaping of the present patterns of genetic homogeneity among the Central-Western Mediterranean Sea mussel populations. For example, the M. galloprovincialis Atlantic form has been reported to have been introduced via ballast water and/or hull fouling to ports in an area from the north of France to Norway, where it has hybridised with local M. edulis (Simon et al. 2020). These hybrids have been termed "dock mussels" (because they are mostly found in docks associated with shipping activity) and have appeared in recent decades rather than centuries. In the opposite direction, Mediterranean ports may also experience strong pressure from invasive species, including M. edulis. For example, M. edulis specimens have been identified using molecular genotyping on a barge hull, which arrived from Middlesborough, northern England, and was moored near the wreck of the cruise vessel Costa Concordia (Isola del Giglio, north Tyrrhenian Sea) in 2012 (Casoli et al. 2016). However, these mussels experienced a high mortality rate at temperatures over 22°C and did not live long enough to interbreed with local populations of M. galloprovincialis. Prior to this, Beaumont et al. (2006) reported anecdotal evidence of the introduction of M. edulis to the Mediterranean coast of France for aquaculture.
Development of an aquaculture industry targeted on Mytilus spp. production and other commercially important shellfish in the twentieth century contributed to unintentional introductions and the spread of many species (Coll et al. 2010) and can also be considered as a vector of mussel transportation between basins in the Mediterranean Sea. An excellent example is the case of Chalastra, a culture area near Thessaloniki, Greece, which serves as the main mussel culture area in the Aegean Sea. Mussel spat from this area have been repeatedly translocated to other mussel farms throughout Greece, Italy and France for ongrowing. Given that these 776 R. Wenne et al. translocations are unrecorded it is very difficult to estimate the magnitude of their contribution to the current lack of geographical structure observed within Mediterranean Sea M. galoprovinciallis populations. An example is the similarity of the Adriatic sample (CHW) to the Black Sea population (BLS). Growing attention has been paid to the spread of Mytilus spp. by rafting on natural or artificial floating objects, including anthropogenic litter (Miller et al. 2018;Rech et al. 2018;Zbawicka et al. 2019; reviewed by Gardner et al. 2021). Plastic debris has been indicated as a transport vector for Mytilus sp./ spp. on long distances in the North Atlantic and Arctic Oceans (Kotwicki et al. 2021). The main contributor to this spread via rafting is mariculture-related gear, including detached bouys ( Figure 6). Interestingly, M. edulis that originated from market discards, aquaria discharged into the sea or introduced for aquaculture purposes in Venezia, Italy (Crocetta 2012), was not observed in our samples and most probably did not establish populations on the coast of Italy. In addition, our results do not confirm the existence of M. edulis (the so-called rock mussel grown in the Ebro Delta for aquaculture purposes, Wenne et al. (2022)

Disclosure statement
No potential conflict of interest was reported by the author(s).