Limited population genetic variation but pronounced seascape genetic structuring in populations of the Mediterranean mussel (Mytilus galloprovincialis) from the eastern Adriatic Sea

Abstract Population genetic analysis of variation at five neutral microsatellite loci for Mediterranean mussels (Mytilus galloprovincialis) from 18 sites along the eastern Adriatic Sea revealed little or no spatial variation. In contrast, seascape genetics analysis revealed a pronounced locus‐specific gradient in allelic and genotypic frequencies across the study region. At a sixth locus, MGE7, the frequencies of two alleles, MGE7243 and MGE7249, were strongly associated, negatively and positively, respectively, with a single environmental variable – minimum salinity (minSAL). The frequency of the MGE7243/243 homozygous genotype was strongly negatively associated with minSAL, whereas the frequencies of the MGE7246/249 and the MGE7249/249 genotypes were strongly positively correlated with minSAL. Interpretation of these pronounced gradients is confounded by the fact that minSAL and another environmental variable, maximum sea surface temperature (maxSST), are highly correlated (R = −.911) and are therefore not necessarily acting independently. BLAST searches of the MGE7 locus against M. galloprovincialis whole genome shotgun sequence returned an alignment with contig mg10_S01094 (accession UYJE01010330.1) and 7 predicted M. galloprovincialis proteins VDI82194.1 ‐ VDI82200.1. Conserved domain searches revealed a similar structure to the transcriptional regulator Msx2‐interacting protein. The BLASTp search also returned significant alignments to Msx2‐interacting proteins in Mytilus coruscus, Crassostrea virginica, and Haliotis rubra. The existence of the MGE7 gradient highlights the role that environmental variation may play in retarding gene flow among wild M. galloprovincialis populations, and also how the success of collection of young mussels (spat) from one site and their transfer to another site (the farm) may be influenced by a single factor such as minSAL or maxSST on a localized scale.


| INTRODUC TI ON
Trying to understand the causative forces of patterns of population genetic variation at different spatial (and temporal) scales has long been a challenge for marine science (Schmidt et al., 2008;Selkoe et al., 2015). Classical population genetics is based on assumptions of marker neutrality because the investigator is usually trying to interpret gene flow (connectivity) in the absence of environmental variation -that is, without including the influence of selection because selective mortality may provide a picture of gene flow that is different from the actual (neutral) genetic connectivity (Holderegger et al., 2006). However, with the advent of highly polymorphic markers such as microsatellites, and more recently single nucleotide polymorphisms (SNPs), and with new analytical approaches, there is now potential to investigate how environmental variation may influence gene flow and population genetic structure (Gagnaire et al., 2015;Manel et al., 2003). The ability to generate new knowledge from markers that are themselves under selection (or are tightly linked to genes under selection) has implications for disciplines such as fisheries, aquaculture, and species conservation (Gaggiotti et al., 2009;Vaux et al., 2021;White et al., 2010), as well as providing new and broad insights into how organisms interact with their environments.
For example, such new knowledge may be important in understanding species-specific responses to environmental variation (e.g., increased sea surface temperature or changes in concentrations of aragonite/calcite and silica) in the face of climate change (e.g., Segovia et al., 2020;Wei et al., 2013a;Zeng et al., 2020).
For many benthic marine invertebrates, the main mechanism for exchange of individuals between populations is via a pelagic larval stage, the fate of which often depends on oceanographic features such as currents, fronts, straits, and shoreline configuration (Pascual et al., 2017;Riginos & Liggins, 2013). During the pelagic (dispersive) stage larvae may be transported within distinct packets of water, so environmental variables related to seawater may be relatively stable (Selkoe et al., 2006). With settlement in intertidal or shallow coastal regions where water mixing is often high, juveniles experience the sea moving about them bringing greater variation in environmental conditions (Riginos & Liggins, 2013). For sessile intertidal species, which cannot avoid the dynamic range of environmental stresses in the intertidal zone (Tomanek & Helmuth, 2002), the transition between pelagic and benthic stages may be a time of particular vulnerability (Pechenik, 1999). If extreme environmental conditions are a regular feature of the settlement location, selective mortality of juveniles can lead to genetic structuring in the adult populations.
For example, peak spawning activity of Mediterranean mussels (Mytilus galloprovincialis Lamarck, 1819) in the Adriatic Sea occurs from late winter through to spring, a period that overlaps with times of pronounced freshwater influence (Da Ros et al., 1985;Krivokapić et al., 2011;Lipizer et al., 2014), a scenario that could lead to genetic structuring as a consequence of genotype-dependent mortality in response to salinity fluctuation. This sort of environmentally driven selection has been very well documented in blue mussel (Mytilus edulis Linnaeus, 1758) populations in Long Island Sound (New York State, Atlantic coast of the USA), where regular annual selective mortality is driven by strong selection against juveniles with the Lap 94 allele in response to low salinity (Hilbish et al., 1982;Hilbish & Koehn, 1985).
A recent review of seascape genetics studies (Selkoe et al., 2016) noted that three factors -temperature, oceanography, and geography -were equally important in explaining patterns of spatial genetic variation across 100 studies. In addition, it was noted that many different factors are likely to affect connectivity at distinct spatio-temporal scales, so while some studies may identify a single key (causative) environmental variable as explaining genetic variation, others will identify two or more variables, often acting synergistically (e.g., Silva & Gardner, 2016;White et al., 2010). It is increasingly clear that understanding patterns of spatial genetic variation, and the factors that contribute to the variation, is very much context-and species-dependent. While important generalities may be drawn from reviews and meta-analyses (e.g., Cárcamo, 2021;Selkoe et al., 2016), understanding the local situation is, of course, critical.
Our study focusses on M. galloprovincialis, a native species on the eastern shoreline of the Adriatic Sea. M. galloprovincialis is a model organism with a large body of research covering areas such as its role as an ecosystem engineer, physiological adaptations to the rocky intertidal zone, aquaculture and its potential as an invasive species (Arribas et al., 2014;Braby & Somero, 2006;Gardner et al., 2016;Kovačić et al., 2017;Wenne et al., 2022). Its success as an invader and ability to outcompete native congeners is in part due to its greater tolerance of a wide variety of environmental conditions in a time of global climate change (Branch & Steffani, 2004;Evans & Somero, 2010;Saarman et al., 2017;Tomanek, 2012). For example, on the Pacific coast of North America, it has supplanted the northern blue mussel (Mytilus trossulus Gould, 1850) in many locations from southern California (USA) to British Columbia (Canada) in waters where it is more warm-tolerant than its native congener (Geller, 1999;Schneider & Helmuth, 2007;Tomanek, 2012). Ecological genetics, Population genetics Dubois et al., 2016;Zonn & Kostianoy, 2017). The Adriatic Sea is a dilution basin for the Mediterranean Sea, contributing ~30% of the freshwater inflow for the whole Mediterranean Sea (Estournel et al., 2021). The Adriatic Sea is nonetheless very saline with seasonal average salinity for the surface water of each of the three basins varying between 37.9-38.2, 37.7-38.3, and 35.5-37.4 psu for the southern, middle, and northern basins, respectively. Evaporation of water from the Adriatic Sea itself and the inflow of high-salinity waters from the northern Ionian Sea (>38.6 psu) are the main processes that maintain this high salinity (Artegiani et al., 1997). The eastern Adriatic coastal region is characterized by a limestone karst landscape with accompanying erosion, which has led to an extensive and complex underground network of flow with many freshwater wells and springs, both on land and on the sea bottom (Bonacci et al., 2013). Mussels in these areas, as well as those in estuaries, can experience wide variations in salinity, which can change very quickly in times of heavy rain (UNEP/MAP-RAC/SPA, 2015). Mussel tolerance to salinity variation may be as tied to the magnitude and speed of changes as to the absolute maximum and minimum values reached (Hamer et al., 2008).
The Mediterranean Sea is predicted to be particularly strongly influenced by global warming trends, with increasing occurrences of marine heat waves and temperature-related mass mortality events (Di Camillo & Cerrano, 2015;Lejeusne et al., 2010;Michaelidis et al., 2014;Oliver et al., 2018). Conditions are already thought to be close to the upper limit of thermal tolerance for Mediterranean mussels in some regions with several instances of mass mortality associated with marine heatwaves in Spain, France, Greece, and in the Adriatic Sea (Galli et al., 2017;Lejeusne et al., 2010;Lupo et al., 2021;Michaelidis et al., 2014). Understanding existing levels of genetic diversity, patterns of gene flow, and the relationship between genetic variation, species ecology, physiology, and environmental factors is critical to the management of this ecologically and commercially important species. Understanding species biology in the native environment can also help in the management of the species in introduced habitats. Here, we describe the use of six polymorphic microsatellite markers to study the population and seascape genetics of 18 populations of M. galloprovincialis, along the eastern Adriatic coast.

| Sample collection
In total, 843 individuals of M. galloprovincialis were collected from 18 sites (hereafter populations) along the eastern coastline of the Adriatic Sea, a region spanning four countries (Croatia, Bosnia-Herzegovina, Montenegro, and Albania) and a variety of coastal environments ( Figure 1; Table 1). The 18 populations were divided between the three regions of the Adriatic Sea (defined by bathymetry and corresponding currents, see Figure 1) with six, seven, and five populations in the northern, central, and southern regions, respectively. Seven were farmed and 11 were wild populations. The number of individuals sampled per site varied between 31 and 50 (mean ± SD of 46.8 ± 5.6). A variety of different sizes (shell length) and therefore putative ages of mussels was collected per site, with shell length in the range of 8-85 mm. We collected all wild (native) mussels from the intertidal zone and the farmed mussels were collected from the shallow subtidal region (top 2 m of the water column F I G U R E 1 Sampling sites of mussels,  at the farm). Samples of mantle tissue were preserved in 99% ethanol for DNA extraction.

| Laboratory protocols
DNA was extracted using Geneaid genomic DNA Mini Kits following manufacturer's instructions and DNA concentrations were quantified using a NanoDrop™ ND-1000 (Thermo Scientific). Individuals were genotyped using seven microsatellite loci described in three published papers as follows: Mgμ1 (Presa et al., 2002), MGE3, MGE4, MGE5, MGE7 (Yu & Li, 2007), MT203, MT282 (Gardeström et al., 2008). The choice of microsatellites was based on amplification success and the presence of polymorphism; details, including amplification conditions that followed Westfall (2011), are summarized in Table A1. Amplified fragments were scanned using an automated analyzer ABI 3730XL by two different service providers (Macrogen and Sangon Biotech). Different versions of the data set have been used for analyses (Table A2). This is because locus Mgμ1 had data missing for ~56% of all mussels and because locus MGE7 was potentially under selection (see outlier analysis section in Results). We used seven loci for all analyses, but repeated analyses with six loci (excluding Mgμ1).

| Genetic data sets
We present results for six loci only because results using seven loci and six loci (with and without Mgμ1) were similar. For the AMOVA, Structure and AWclust analyses, we also report results for five loci (i.e., excluding Mgμ1 and MGE7) so that results using neutral loci only may be compared directly to results that included MGE7. To calculate the F-statistics used in the seascape analyses, we used a data set of 577 individuals containing complete data for six loci (i.e., excluding Mgμ1 because of missing data).
To assess whether sample sizes were sufficiently large to characterize the allelic variation of each population, the rarefaction procedure in the software package PopGenKit (Rioux Paquette, 2011) in R 3.5.1 (R Core Team, 2014) was used with 1000 jackknife replicates to estimate the total allelic diversity for each marker, each population and over all populations.
Analyses of Hardy-Weinberg equilibrium (HWE) and Linkage Disequilibrium (LD) were performed using Genepop 4.7.0 (Rousset, 2008). Departures from HWE were assessed by two methods, the probability test and the U-test with the alternate hypothesis of heterozygote deficiency, using the default parameters. False discovery rate control for multiple statistical testing was applied to p-values (Verhoeven et al., 2005).
Rare alleles were excluded from the calculation as recommended by Waples and Do (2010).

| Allelic and genotypic estimations and analyses
The frequency of alleles, number of private alleles (PA), observed and expected heterozygosities (H O and H E ), and the fixation index (F IS ) were calculated using GenAlEx 6.503 (Peakall & Smouse, 2012).
Allelic richness (A R ) and private allelic richness (PA R ) were calculated using HP-RARE 1.1 (Kalinowski, 2005

| Differentiation among regions and populations
To test if there were genetic differences among populations in the three regions in the Adriatic Sea ( Figure  Two different cluster analysis programs, Structure (Pritchard et al., 2000) and AWclust (Gao & Starmer, 2008), were used to assess genetic structure for the 18 populations. Structure uses Bayesian analysis, which assumes conformation to HWE and LD assumptions, while AWclust is non-parametric and therefore does not require conformity to HWE and LD assumptions. Structure analyses were performed with a burn-in length of 50,000 steps, run length of 100,000 steps, 10 iterations for K = 1 to 18 using the correlated allele frequency model, and both the Admixture and No admixture models with sampling locations as priors. The Evanno method (Evanno et al., 2005) as implemented in the program Structure Harvester (Earl & vonHoldt, 2012)

| Environmental data acquisition
Environmental data were collected from published articles and online databases (Table A3). Because there is no standardized environmental data collection between countries bordering the eastern Adriatic Sea, the data available varied widely among sites.  Table 1). For these latter areas, which were all in Croatia, the online resources at SeaDa taNet.org and the Croatian sea bathing water quality database (IOR, 2022), together with the models presented in Lipizer et al. (2014) and satellite data presented in Böhm et al. (2003), were used and also allowed for comparison among all sites. We did not use satellite data available from Copernicus because we could not retrieve data from enclosed bays and convoluted coastlines that accurately reflected fine-scale conditions (Đurović et al., 2018). More detailed information about environmental data types, spatial and temporal coverage, and sources is provided in Hamilton (2019). Overall, there was good agreement for trends and inter-site comparisons across all data sources, particularly for sea surface temperature and salinity data. The chlorophyll-a data were not so complete, with detailed information available for the 10 southernmost populations but only rather sparse data for the eight northern populations, for which we relied on satellite-derived data from Böhm et al. (2003) and the website SeaDa taNet.org.
Nine environmental variables were used in analyses, the minimum and maximum values of sea surface temperature, salinity, and chlorophyll-a concentrations (minSST, maxSST, minSAL, maxSAL, minCHL-a and maxCHL-a), together with delta values (deltaSST, del-taSAL and deltaCHL-a -the difference between the maximum and minimum readings) that reflect the full range of variability. Three geospatial variables were also employed, latitude (Lat), longitude (Long), and the total coastal distance (TotalCD). TotalCD was calculated following Wei et al. (2013a) as the total of the pairwise coastal distances (CD in km) between any one population and all others (where CD is the shortest distance by sea between each population pair). Sites with high values of TotalCD are, on average, at a greater distance from all others than sites with low values of TotalCD. The environmental and geospatial variables used are detailed in Table A4.  The importance of each independent variable in explaining genetic structure was evaluated by calculating the percentage of best-fitting models containing that independent variable ( Figure A1; Table A5).

| Seascape genetics -distance-based linear modeling
To test the effect of the six independent environmental/geospatial variables (minSAL, minSST, maxSAL, minCHL-a, Lat, TotalCD) on site-specific genetic variation, a distance-based linear model was employed (e.g., Silva & Gardner, 2016 Rarefaction curves showed that ~80% of allelic variation was captured in each population for six loci estimates of the effective population size, Ne, were large (Table A6). For the majority of populations, the upper limit of the 95% confidence interval was "infinite," where an "infinite" population size means that there is no evidence of genetic variation caused by genetic drift due to a finite number of parents. Only one population, IW, had a finite upper limit for the 95% confidence interval for both data sets. The probability that alleles from each population were drawn from the same distribution for six loci, calculated by exact G-tests, was <1 × 10 −20 across all loci, less than 1 × 10 −5 for individual loci MT282, MGE4, MGE7, and MT203, 0.021 for MGE5 and 0.187 for MGE3. The probability that genotypes from each population were drawn from the same distribution for six loci was 4.4 × 10 −7 across all loci, less than 1 × 10 −5 for MGE7, 0.018 for MT282 and not significant for the remaining loci.

| Allelic and genotypic estimations and analyses
There were no differences in A R or F IS between farmed and wild populations (Kruskal-Wallis H statistic = 0.740, p = .39 and H statistic = 2.628, p = .10 for A R and F IS , respectively). Regional allele frequencies are illustrated in Figure A2.

| Differentiation among regions and populations
AMOVA (Table 2)  For the six locus data set, all three pairwise comparisons of F RT values among regions were statistically significant (p ≤ .012) whereas for the five locus data set only the North-Central comparison was significant (Table 3).
Pairwise population F ST values for the six locus and five neutral locus data sets were variable with 64 of 153 comparisons significant F I G U R E 2 Allele frequencies for MGE7 for each population, each region and across all data.
after FDR correction for the six locus data set and only one for the neutral data set (Table A8). Two populations exhibited unusually high The Structure analysis runs using the Admixture model for the six locus data set had a maximum likelihood statistic, ΔK, at K = 5, with a local maximum at K = 2. When the No Admixture model was used, the maximum likelihood was attained at K = 2 ( Figure A3). The Structure runs for K = 5 did not show any obvious genetic structure but for K = 2, using F I G U R E 3 MGE7 243 allele frequency and structure results: (a) map of frequency of MGE7 243 allele (cyan) and all other alleles (orange) at each population, (b) correlation of frequency of MGE7 243 allele with latitude, (c) map of proportion of cluster membership for six loci from the structure analysis (cluster 1 orange and cluster 2 cyan), (d) correlation of proportion of cluster membership by population for six loci (green) and five neutral loci (blue).
both models, the bar plot suggested a north-south cline in the proportion of ancestry from each cluster, something that was confirmed by linear correlation tests of proportion of cluster membership as a function of latitude (R 2 = .459, p = .002, Figures 3 and A3). Structure runs on the neutral data set did not show this north-south cline (R 2 = .111, p = .178, Figure 3) and there was no evidence of any genetic differentiation.
AWclust analyses ( Figure A4) of the six locus data set revealed overlapping confidence intervals and local maxima at K = 2, 4, and 7, indicating weakly defined structure. For K = 2, the correlation between cluster 2 membership and latitude was significant (R 2 = .308, p = .017) but weaker than the corresponding correlation from the Structure program. For the five neutral loci, only one cluster was evident and when K = 2, the correlation between cluster 2 membership and latitude was not significant (R 2 = .073, p = .279).

| Correlation testing of environmental/ geospatial variables
Based on PCA and correlation analysis with a critical cut-off value of |R| > 0.85, we removed six of 12 environmental and geospatial variables from the data set ( Figure A5, Table A9). Subsequent analyses involved six independent environmental and geospatial variables (minSAL, minSST, maxSAL, minCHL-a, Lat, and TotalCD).

| Seascape genetics -GLZ analysis
GLZ analyses were performed on five genetic response indices -LinF ST5 , Although minSST and minCHL-a were included in the top-ranked model (together with minSAL) for f(MGE7 243/243 ), the three lowest ranked models of all possible models included only minSST and/or minCHL-a.

| Seascape genetics -distance-based linear modeling
Marginal analyses revealed that two of the six variables explained significant variation in the 6-locus data set (minSAL (p = .0016), Lat (p = .0434)) ( that the single best-fit model contained only one variable -minSAL (Table 5). This model explained 13.8% of the variation in the genetic data set. As expected, addition of new terms to the model improved the fit as judged by the R 2 value, which increased from 0.138 (one variable) to 0.432 (six variables). The variable minSAL was included in all six of the six best-fit models, indicating that it is the single variable with greatest power to explain variation in the genetic data set.

| Correlation of MGE7 alleles and genotypes with environmental variables
The most common allele over all populations, MGE7 243 , was negatively correlated with minSAL (f = 0.538, R 2 = .635, p < .0001).  Figure 4). Although the observed frequency of the MGE7 243/249 genotype was 0.009 (n = 7 mussels observed), the expected frequency was 0.13 (n = 100.6 mussels expected). The ratio of observed to expected for this genotype (0.07) was very low, even in a background of heterozygote deficiencies. The ratios of observed to expected for the three other heterozygous genotypes The variables minSST and minCHL-a were included with minSAL in the top-ranked model for the GLZ analysis based on the genetic response index f(MGE7 243/243 ) but neither was significantly correlated with any MGE7 allele or genotype.

| DISCUSS ION
In this study, focusing on 18 populations of M. galloprovincialis along the eastern coast of the Adriatic Sea, neutral multi-locus microsatellite analysis revealed low levels of genetic differentiation, very large population sizes and heterozygote deficiencies (consistent with a large body of work, e.g., Daniels & Litvaitis, 2017, Hedgecock et al., 2004, Wei et al., 2013b, Zouros & Foltz, 1984. These results add fine-scale detail to the literature and are consistent with other regional or basin-wide studies that showed either well-connected populations in the Adriatic Sea close to panmixia, or weak east-west and/or north-south differentiation (Giantsis et al., 2014;Paterno et al., 2019;Štambuk et al., 2013;Wenne et al., 2022). We also report a pronounced gradient in allelic and genotypic frequencies at one locus that is very strongly associated with environmental variation along the eastern Adriatic Sea coastline.
The potential role of selection in generating a gradient of ge- and on the homozygous genotype MGE7 243/243 were strongly associated minSAL, and the models were statistically highly significant (p < .0001). The geospatial variable Lat was not statistically significant in these multivariate models, this result being consistent with the stronger single variable correlation of MGE7 243 allele frequency with minSAL (R 2 = .635, p < .001) than with Lat (R 2 = .476, p = .015).
The association between the MGE7 243 allele and minSAL extended to the homozygous MGE7 243/243 genotype, but not to the heterozygous MGE7 243/X genotype. Analysis of regional subsets of the data (e.g., northern, central, southern Adriatic) and analysis with or without individual sites (e.g., farmed and wild samples from Boka Kotorska Bay) revealed that the overall relationship is robust and not dependent on the inclusion or exclusion of specific sites (results not shown). Variations in allele frequencies that coincide with changes in an environmental variable or geographic cline may arise from chance effects and historic events (Schmidt et al., 2008), but the distinctive distribution pattern of the homozygous and heterozygous genotypes is much less likely to arise by chance alone and suggests a selective pressure. The MGE7 locus may be a neutral hitchhiker locus situated close to a gene under selection or it may itself be part of the coding region of a gene under selection (Gagnaire et al., 2015). This would not be surprising, since the locus was derived from EST sequences from GenBank (Yu & Li, 2007); it was however unexpected since this locus has previously been used successfully as a neutral marker (Westfall, 2011), including in one study focused on the Croatian coast (Štambuk et al., 2013). The trend of decreasing frequency of the MGE7 243 allele with increasing minimum salinity is not balanced equally by the remaining alleles.  Both maxSST and minSAL are relevant to many biological processes that may affect individual survival, and have long been known to play roles in genetic selection related to environmental clines or habitat mosaics (e.g., Gardner & Palmer, 1998;Halpin et al., 2002;Hilbish & Koehn, 1985;Logan et al., 2012;Negri et al., 2013;Wenne et al., 2020). While some differences in tolerance to thermal and hypoosmotic stresses are linked to differences at the individual gene level (such as in the Lap gene of M. edulis in Long Island Sound, New York State, Hilbish et al., 1982), other differences are at the regulatory level of transcription or post-translational modification.
In an analysis of the transcriptional response of Mytilus congeners to acute heat and low-salinity stress, of 45 genes that responded to both stressors, the response was in opposite directions for 44 of the genes, with the one gene with a response in the same direction encoding a thioredoxin reductase enzyme involved in oxidative stress (Lockwood et al., 2015). Although relatively few changes in transcriptional regulation between M. galloprovincialis and M. trossulus resulted in substantial differences in thermal tolerances between the congeners (Lockwood et al., 2015), the severe impact of heat stress on highly evolved and conserved systems such as energy metabolism and detoxification argue that there is little room for further evolutionary adaptation to thermal stress (Michaelidis et al., 2014;Tomanek, 2012).
Environmental stresses, including thermal and hyposaline stress, lead to increased levels of oxidative stress, thus if selection is the driving force behind the variation in MGE7 allele frequencies, the mechanism may be related to pathways that deal with oxidative stress rather than mechanisms specifically related to low salinity or high temperature (Hamer et al., 2008;Michaelidis et al., 2014).
The effects of environmental stresses operate synergistically to increase overall levels of oxidative stress, thereby reducing tolerance to individual stresses. For example, heavy metal pollution increases sensitivity to thermal stress (Michaelidis et al., 2014). Environmental stresses increase energy demands generally in order to maintain cellular homeostasis, repair damaged proteins, and detoxify reactive oxygen species (ROS). There is a trade-off between energy metabolism which generates ROS and expending energy to detoxify ROS, which can reach a tipping point particularly during times of thermal stress when food availability may be at its lowest (Michaelidis et al., 2014).
Environmental stresses also increase susceptibility to disease which may trigger mass mortalities (Di Camillo & Cerrano, 2015;Lejeusne et al., 2010). Any adaptive mechanism that enhances the ability of mussels to respond to oxidative stress is expected to be favored in environments which experience extremes in both sea surface temperature and salinity (e.g., Boka Kotorska Bay and Butrint Lagoon, the locations of sampling sites CF, CW, IF, IW, and SBF) (Hamer et al., 2008, Michaelidis et al., 2014.  & Presa, 2008& Presa, , 2009Giantsis et al., 2014;Larraín et al., 2015), and a limited number of studies investigating population genetic structure of M. galloprovincialis using microsatellites.
Interestingly, markers have often not performed well when used on mussels from different geographic regions, with, for example, only a subset of markers developed for western Mediterranean mussels (Diz & Presa, 2008), performing satisfactorily for eastern Mediterranean mussels (Giantsis et al., 2014;Štambuk et al., 2013). We suspect that this may be because of the regional differences that appear to exist in the pan-genome of M. galloprovincialis (Gerdol et al., 2020), which also may explain the very high proportion of null alleles (i.e., via hemizygosity) that is so often reported for this species and why so many loci are out of HWE. More recently population genetic analyses have been performed using SNPs Zbawicka et al., 2014) including the Mediterranean and Black Seas (Paterno et al., 2019;Wenne et al., 2022). Our data are consistent with the SNP results but add to our understanding by being much more focussed geographically and using, where possible, local environmental data. Although SNPs are now the marker of choice (Fischer et al., 2017), microsatellites offer complementary information to SNPs, having a higher mutation rate and more recent temporal scale of inference (Fotsing et al., 2019;Waits & Storfer, 2015;Wu et al., 2021).
The findings of this study have important implications for aquaculture and conservation management in relation to climate and other anthropogenic change. Aquacultural practice along the eastern Adriatic coast is largely traditional with spat often transferred from natural spawning grounds to aquaculture sites with favorable conditions for growth (Giantsis et al., 2014;Kovačić et al., 2017). The traditional methods of mussel farming employed are extremely labor intensive so the matching of environmental conditions between the natal site of the spat and the aquaculture site is recommended to reduce selective mortality related to environmental stresses (Mandić et al., 2017;Ramón et al., 2007). This is illustrated by the MGE7 locus, with the frequency of genotypes differing between environments. Aquaculture sites are often in bays and lagoons (for example, sampling sites CF and IF in Boka Kotorska Bay and SBF in Butrint Lagoon) with high maxSST and low minSAL. Our results suggest that juveniles collected from areas of low environmental variability will have lower proportions of the favorable MGE7 243 allele than juveniles collected from areas of high environmental variability and are more likely to undergo selective mortality.
Climate change has altered the geographic range of many species, moving ranges poleward and has been implicated in the success of M. galloprovincialis as an invasive species outside of its native range (reviewed by Gardner et al., 2021). Within its native range of the Mediterranean Sea, reports of mass mortalities linked to extended periods of elevated temperatures are increasing (Galli et al., 2017;Lejeusne et al., 2010;Michaelidis et al., 2014;Verdura et al., 2019). For example, widespread disease triggered by longlasting high water temperatures and eutrophic waters in the northern Adriatic Sea led to mortality across multiple taxa, including M. galloprovincialis, which, in some areas, decreased in abundance from ~20% coverage to ~2% in 2013, with recovery to 9.5% by 2014 (Di Camillo & Cerrano, 2015). Loss of large proportions of filter feeding animals will impact ecosystem functioning and water quality and increase the risk of eutrophication and further mass mortality events. Conservation and aquaculture practices should therefore be directed at maintaining genetic diversity and monitoring the quality of the environment to reduce the synergistic effects of pollution, eutrophication, and thermal stress (Harley et al., 2006;Lejeusne et al., 2010).

ACK N OWLED G M ENTS
We thank Mario Lovrinov, Dragana Milošević, and Tea Tomljanović for help locating and collecting mussels, and Catarina Silva and XiangZhao Guo for help in the laboratory and with preliminary data acquisition. This work was supported by funding from the University of Zagreb Faculty of Agriculture, Croatia to MP, MF, LS, and by funding from the School of Biological Sciences, Victoria University of Wellington, New Zealand to JPAG. JSH was supported by a scholarship from VUW, NZ.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available