Genetic diversity and connectivity within Mytilus spp. in the subarctic and Arctic

Abstract Climate changes in the Arctic are predicted to alter distributions of marine species. However, such changes are difficult to quantify because information on present species distribution and the genetic variation within species is lacking or poorly examined. Blue mussels, Mytilus spp., are ecosystem engineers in the coastal zone globally. To improve knowledge of distribution and genetic structure of the Mytilus edulis complex in the Arctic, we analyzed 81 SNPs in 534 Mytilus spp. individuals sampled at 13 sites to provide baseline data for distribution and genetic variation of Mytilus mussels in the European Arctic. Mytilus edulis was the most abundant species found with a clear genetic split between populations in Greenland and the Eastern Atlantic. Surprisingly, analyses revealed the presence of Mytilus trossulus in high Arctic NW Greenland (77°N) and Mytilus galloprovincialis or their hybrids in SW Greenland, Svalbard, and the Pechora Sea. Furthermore, a high degree of hybridization and introgression between species was observed. Our study highlights the importance of distinguishing between congener species, which can display local adaptation and suggests that information on dispersal routes and barriers is essential for accurate predictions of regional susceptibility to range expansions or invasions of boreal species in the Arctic.

can exhibit different adaptations and tolerance limits to specific environments (Limborg et al., 2012;Nielsen, Hemmer-Hansen, Larsen, & Bekkevold, 2009;Thyrring, Rysgaard, Blicher, & Sejr, 2015), which make it important to know their current distribution and the connectivity between populations and the processes governing the distribution of genetic variation. Through genetic analysis, it is possible to determine the level of genetic variability within both threatened and newly established populations, the origin of migrating individuals, direction of gene flow, and possible adaptive evolutionary changes associated with climate change (Brodersen & Seehausen, 2014;Hansen, Olivieri, Waller, & Nielsen, 2012;Laikre, Schwartz, Waples, & Ryman, 2010). All key factors needed to make predictions for the likely impact of climate change.
Bivalves of the genus Mytilus are frequently used as environmental indicators, as they are semi-sessile, have a relatively long life span, and are widely distributed in coastal regions in both Northern Hemisphere and Southern Hemisphere (Goldberg, 1986;Gosling, 2003;Rainbow, 1995;Thyrring, Juhl, Holmstrup, Blicher, & Sejr, 2015). Mytilus spp.
have been the subjects of genetic studies for decades as the different species are morphologically difficult to distinguish. Consequently, the population structure of individual Mytilus species has been difficult to establish. Mytilus edulis L. 1758, Mytilus trossulus Gould 1850 and Mytilus galloprovincialis Lmk. 1819, all belong to the M. edulis species complex and are known to coexist and hybridize with conflicting patterns on the fitness for hybrids. Some studies did not observe any depressed fitness (Doherty, Brophy, & Gosling, 2009;Koehn, 1991;Riginos & Cunningham, 2005;Toro, Thompson, & Innes, 2006), while others (Gardner & Thompson, 2001;Toro, Innes, & Thompson, 2004;Toro, Thompson, & Innes, 2002;Tremblay & Landry, 2016) found a difference in fitness between parental types and hybrids and backcrosses. These findings and numerous studies on introgression between them (Fraïsse, Belkhir, Welch, & Bierne, 2016;Roux et al., 2014) have challenged the isolation species concept (White, 1978); however, they are generally considered to be different species, as they remain ecological distinct despite semipermeable barriers for gene flow and introgression Fraïsse, Roux, Welch, & Bierne, 2014;Saarman & Pogson, 2015). Mytilus trossulus is thought to have invaded the Arctic Ocean from the Pacific Ocean around 3.5 million years ago (mya) through the Bering Strait (Rawson & Hilbish, 1995Vermeij, 1991). As the Bering Strait closed during glacial periods, allopatric speciation resulted in the evolution of M. edulis in the Atlantic. Mytilus edulis has since spread to large parts of the Atlantic and due to apparent low gene flow (at least for some loci); M. edulis populations on each side of the Atlantic are genetically distinct Riginos, Hickerson, Henzler, & Cunningham, 2004). Speciation between M. edulis and M. galloprovincialis most likely occurred through allopatric isolation approximately 2.5 mya (Quesada, Gallagher, Skibinski, & Skibinski, 1998;Rawson & Hilbish, 1995 with secondary contact and introgression occurring around 0.7 mya . Between interglacial periods 46,000 and 20,000 years ago, M. trossulus reinvaded the Arctic Ocean (Rawson & Harper, 2009). From here, it invaded both sides of the Atlantic founding M. trossulus/M. edulis hybrid zones along North American and European coasts (Riginos & Cunningham, 2005).
populations are found all along the west coast, and southern populations from Tartoq and Narsarsuaq have been shown to be genetically distinct from European M. edulis displaying higher resemblance to Canadian and North American M. edulis populations . Few genetic analyses have been performed on Mytilus spp. in Greenland, and most studies have assumed these mussels to be M. edulis without genetic verification despite observations of variations in metabolic response to low temperatures between populations from NW and SW Greenland (Thyrring, Rysgaard, et al., 2015).
Moreover, in 2004, subtidal M. edulis were discovered at the mouth of Isfjorden in Svalbard after 1,000 years of absence (Berge, Johnsen, Nilsen, Gulliksen, & Slagstad, 2005). These mussels were hypothesized to have been transported from Norway by the West Spitsbergen Current in 2002, but their origin has never been confirmed through genetic analysis. Mytilus trossulus is less common in Arctic waters. Väinölä and Strelkov (2011) found that M. trossulus had a scattered distribution in the White Sea and the Norwegian Sea, and Feder, Norton, and Geller (2003) found live M. trossulus in Arctic Alaska in the 1990s. Furthermore, Wenne, Bach, Zbawicka, Strand, and McDonald (2016) has recently reported a NW Greenlandic fjord at Maarmorilik (71°N) to be inhabited by M. edulis, M. trossulus, and their hybrids. Mytilus edulis and M. trossulus hybrid zones have also been found and studied on the European and N American Atlantic coasts. Riginos and Cunningham (2005) reviewed the literature on the subject to look at local adaptation and species segregation and found conflicting patterns of species segregation across the Atlantic. In the western Atlantic, M. trossulus was found on wave-exposed open coasts, whereas M. edulis appeared to dominate in sheltered areas of low salinity. However, European M. trossulus populations from the Baltic Sea appeared to be locally adapted to the prevailing low salinities. The latter is in line with the findings of Wenne et al. (2016), who found a higher prevalence of M. trossulus in the inner Maarmorilik fjord compared with the more saline outer fjord.
Mytilus galloprovincialis normally inhabits warmer waters, but in recent years the species and M. galloprovincialis/M. edulis hybrids have been observed along the Norwegian coast (Brooks & Farmen, 2013;. This could be related to human activities like ship traffic in rural areas enabling faster invasion of waters otherwise not directly accessible to them (Anderson, Bilodeau, Gilg, & Hilbish, 2002;Geller, Carlton, & Powers, 1994). Furthermore, it has been demonstrated that M. galloprovincialis is capable of tolerating low water temperatures (Inoue et al., 1997), highlighting the potential for this species to occur in the Arctic.
Utilizing 81 nuclear SNPs, we examined the distribution of Mytilus spp. in subarctic and Arctic regions ranging from the eastern Baffin Bay to the Pechora Sea with a special emphasis on the spatiotemporal population structure of M. edulis. We further aimed at identifying the source population or populations for the newly discovered M. edulis population in Svalbard and whether the observed differences in temperature response found in W Greenland Mytilus populations (Thyrring, Rysgaard, et al., 2015) could be caused by genetically based local adaptation.

| Study sites and sampling
Nineteen Mytilus spp. samples, consisting of 509 individuals in total, were collected from thirteen subarctic and Arctic sites (Table 1 and  and Qaanaaq (QAS and QAL). From Tromsø, Nuuk, and Qaanaaq, mussels of different size classes were collected as size can be used as a proxy for age class and hence indicate possible short-term genetic change over time; TRS, NUS, and QAS were smaller mussels in the size range 15-30 mm in length, while TRL, NUL, and QAL being mussels larger than 50 mm.
Samples were stored at −19°C (TRS, TRL, SRI, KOB, UPE, QAS, and QAL) or in 96% ethanol at 4°C (LOF, SV1, SV2, SV3, SV4, WS1, WS2, PSW, and PSE). Measurements of shell length, width, and height of frozen mussels were conducted at the laboratory facilities at DTU Aqua in Silkeborg, Denmark, whereas the measurements of ethanol preserved specimens were taken at the sampling locations.
Reference samples of M. trossulus and M. galloprovincialis were provided from the Sea of Okhotsk, Russia, and from around Galician Rías in NW Spain, respectively, to evaluate the species status of the sampled mussels and to identify potential hybrids.

| DNA extraction
A minimum of 30 mg (wet weight) of mantle tissue was dissected from each mussel, and DNA was extracted using the Omega EZNA Tissue DNA kit (Omega Bio-Tek, Norcross, GA, USA) according to the manufacturer's instructions for tissue. DNA content in the extracts was verified on a NanoDrop spectrophotometer (Thermo Scientific, Waltham, MA, USA). genetic structure (Zbawicka et al., 2012(Zbawicka et al., , 2014 and 77 new SNPs originated from RAD sequencing of genomic M. edulis DNA (EBI Sequence Read Archive (SRA) study ERP006912) at the University of Stirling (Table S1).

| Identification of hybrids
Based on the generated pairwise F ST estimates, the grouping of samples was visualized in a multidimensional scaling plot applying the cmdscale function in R (R Core Team, 2015). Additionally, a principal component analysis scatter plot was created in the R package Adegenet (Jombart, 2008;Jombart & Ahmed, 2011) to illustrate the genetic relationships among individuals across all samples. Structure v2.3.4, utilizing the Bayesian MCMC clustering approach (Pritchard, Stephens, & Donnelly, 2000) was used to visualize species integrity and identify possible hybridization among Mytilus spp. using a variable number of predefined clusters (K) for grouping individuals. This was also done to positively identify M. edulis individuals and subsequently create a reduced data set exclusively aimed at investigating population structure within this species. Considering the close genetic resemblance of Mytilus spp. and the assumed low gene flow between geographically distant samples , simulations were run for a number of predefined K values. Based on an initial analysis of K up to 18, we found the highest likelihoods for K = 3-5. Accordingly, we used this as the basis to identify the major groupings within the species complex. For all simulations, a burn-in of 10,000 iterations was used followed by 100,000 MCMC repetitions.
To evaluate the power of designating individuals as pure or hybrids, we followed the procedure described in Nielsen, Hansen, Ruzzante, Meldrup, and Grønkjaer (2003) using the program Hybridlab (Nielsen, Bach, & Kotlicki, 2006). Briefly, we simulated 1,000 individuals of each of the following classes: parentals, F1/F2, and backcrosses. This

| Population structure of Mytilus edulis
A reduced data set was used to assess the population structure in M. edulis. Based on the results from the analysis of simulated parentals and hybrids (see results section), we chose to only include individuals T A B L E 1 Summary of sample information including sampling location and year, sampling code indicating location and possibly size, estimates of expected (H e ), and observed (H o ) heterozygosities and the allelic richness the Norwegian, Svalbard, and Russian samples to infer the likely origin of Svalbard mussels. The PCA scatter plots were generated in R, using the cmdscale function and the package Adegenet. "outlier" data set and a "neutral" data set were created to test the importance of outlier loci for defining the inferred population structure of M. edulis; that is, the true connectivity among populations based on neutral processes (drift and migration) could be obscured by loci under divergent selection. For both data sets, overall and pairwise F ST estimates were generated in Genepop 4.2, while Structure v2.3.4 with K = 2 (using settings as above) was used to investigate whether the population structure found in M. edulis based on all loci could be identified from both the "outlier" and "neutral" data sets, or whether they displayed contrasting patterns.

| Summary statistics
Three loci (159069_A, 171383_A, and 170478_A) deviated significantly from HWE in ten samples or more, and they were discarded from further analyses.

| Identification of hybrids
The multidimensional scaling plot ( Fig. 2A)  proach enabled the construction of an exclusive "Mytilus edulis" data set. When using K = 3 the analysis was unable to split the samples into the three a priori defined species groups (Appendix 2B).

The inferred proportion (using the 20% criterion) of the different
Mytilus species and hybrids in each of the geographical samples ( Fig. 1) show that M. edulis is the most common species within the sampled   Table 1 3

.4 | Outlier analysis
The outlier tests identified six loci as F ST outliers, with six loci significant at the 5% level and three at the 1% level. All of these outliers are high F ST outliers (Fig. 6) indicating diversifying selection (Beaumont & Nichols, 1996), although a few of them could represent the upper tail of the neutral F ST distribution. Also, a strong genetic cline as observed here is known to sometimes overestimate the number of loci under diversifying selection (Strand, Williams, Oleksiak, & Sotka, 2012). Furthermore, introgression between Mytilus spp. has been found to cause high F ST outliers (Gosset and Bierne 2012). Pairwise F ST values ranged from 0 to 0.059 for the "neutral" data set and from 0 to 0.474 for the "outlier" data set (Tables S4 and S5). The Structure analyses for both the "neutral" and "outlier" data set also supported the initial population structure separating Greenlandic samples from the Norwegian, Svalbard, and Russian samples and with the Icelandic sample of admixed origin (Appendices 4 and 5).

| Distribution of Mytilus spp. in the Arctic
Baseline information of species distribution and their genetic composition is imperative in order to quantify the impacts of climate change on species distribution ranges, biodiversity, and the effects of hybridization between species and populations (Gardner, Zbawicka, Westfall, & Wenne, 2016). Molecular genetic knowledge is a key measure to identify the distribution of invasive congener species (Geller, Darling, & Carlton, 2010), which may cause cascading ecosystem effects. Despite congener species appearing morphologically similar, interspecific variation in ecology and physiology may impact population fitness (Fly & Hilbish, 2013;Fraïsse et al., 2016;Somero, 2005). In the Arctic, baselines studies on genetic variation and species abundance are largely absent but urgently needed Wassmann, Duarte, Agusti, & Sejr, 2011).  (Berge et al., 2005;Renaud, Sejr, Bluhm, Sirenko, & Ellingsen, 2015). For instance, the northward flowing current regimes (such as the Norwegian Current) allows non-Arctic species to extend their range into the Arctic from the Atlantic or Pacific Ocean Fetzer & Arntz, 2008;Sirenko & Gagaev, 2007). Ocean currents also explain why Mytilus spp. remain absent in NE Greenland despite the environmental resemblance of NW Greenland with regard to temperatures and ice conditions (Sejr, Blicher, & Rysgaard, 2009 is considered biogeographically different from the rest of Greenland (Piepenburg et al., 2011). The absence of Mytilus mussels in NE Greenland is likely a result of dispersal barriers due to the lack of an downstream source population, as the East Greenland Current flows from north to south, exemplifying how outflow shelves may respond slowly to climatic changes (Renaud et al., 2015). This is further supported by the presence of Mytilus mussels in SE Greenland, at Tasiilaq (Ammassalik, 65°N) (Ockelmann, 1958), which is influenced by a branch of the Irminger Current from the Atlantic Ocean.
The present study highlights the need for further genetic studies in the region as a M. trossulus population was found in the most northern  (Maggs et al., 2008;. Second, there could be a contemporary spread of M. trossulus from the Pacific Ocean. Jones et al. (2003) found that waters around NW Greenland contained high levels of phosphate indicating Pacific water being transported into this area.
Also, there are a few reports of live M. trossulus in Arctic Alaska and Canada (Feder et al., 2003), so the spread of planktonic larvae from the Canadian Arctic could be possible. A third scenario could be that M. trossulus spread to Arctic Greenland from the East coast of Canada.
However, as the West Greenland Current moves along the coast from south to north, and Mytilus mussels are expected to disperse with rather than against currents, this scenario seems unlikely (McQuaid & Phillips, 2000). Finally, Mytilus spp. are known to disperse by human activities and can survive long distances and fluctuating temperatures (Lee & Chown, 2007). Qaanaaq is situated less than 150 km from the US Thule Air Base, which receives supplies by US ships; this is providing an alternative dispersal route of M. trossulus from the north Pacific.

The invasive M. galloprovincialis appeared widespread from
Greenland to the Pechora Sea. In Norway, M. galloprovincialis appears common along the coastline (Brooks & Farmen, 2013), and the discovery of M. galloprovincialis in Svalbard suggests colonization by ocean currents as hypothesized by Berge et al. (2005) or ship traffic from the Norwegian mainland (Ware et al., 2014).

| Mytilus hybrid zones in the Arctic
Most sampling locations displayed varying degrees of hybridization and introgression between different Mytilus spp. and only four locations contained apparently pure populations (Fig. 1). Introgression can affect a population's fitness and vulnerability to climate change.
In the study region, hybrid zones were found in Norway, Svalbard, and Greenland, with the highest abundance of the invasive M. galloprovincialis found along the Norwegian coast, especially in Lofoten (68°N) further supporting the findings by Brooks and Farmen (2013) and . Additionally, a surprisingly high amount of M. galloprovincialis was found at Svalbard. We also found evidence of limited introgression of M. galloprovincialis in the Russian and Icelandic samples, and the ecological consequences of invasive mussels in these regions need to be studied further. In the White Sea, M. trossulus individuals were only recorded in one of two loca- tions. This small-scale regional variation in species composition was also observed by Väinölä and Strelkov (2011), who also found M. trossulus and M. edulis/M. trossulus hybrids but to a much lesser extent than M. edulis. It is believed that the expansion of M. trossulus in the White Sea is most likely facilitated by ships (Väinölä & Strelkov, 2011

| Population structure of Mytilus edulis
The Iceland from the east, it is perhaps more likely that Iceland would be recruiting spat from East Atlantic populations. This is also inferred by Riginos and Henzler (2009), who found postcolonization gene flow from northern Europe to Iceland.
The outlier tests identified six loci as F ST outliers at the 5% significance levels. All of these outliers are high F ST outliers (Fig. 6) indicating diversifying selection (Beaumont & Nichols, 1996). However, a strong genetic cline as observed here is known to sometimes overestimate the number of loci under diversifying selection (Strand et al., 2012).
Furthermore, introgression between Mytilus spp. have been found to cause high F ST outliers (Gosset & Bierne, 2012), and this result should be interpreted with some caution. Still, we find that the pattern of population structure is the same for the "neutral" and the "outlier" data sets (Appendices 4 and 5), suggesting that patterns of neutral population structure is correlated with adaptive evolution in response to divergent local environmental conditions. Temperature influences the large-scale geographical distribution of species (Sunday, Bates, & Dulvy, 2011); however, on a local scale other factors including predation, the presence of sea ice, suitable habitats, water current, and salinity can influence the distribution of intertidal species (Høgslund, Sejr, Wiktor, Blicher, & Wegeberg, 2014;Kroeker et al., 2016;Paine, 1974), and these conditions are very different between W Greenland and the other sampling sites (Rayner et al., 2003). Still, the high divergence between samples from the Eastern Atlantic and Greenland cannot be explained alone by loci subject to selection. F ST values for the "neutral" data set are still high (Table S4) suggesting a high degree of isolation between groups. This isolation in turn may have facilitated local adaptation at this rather large geographical scale. For more specific insights on the environmental factors responsible for local adaptation, the geographical scale, and its genomewide significance, a more elaborate sampling design is warranted including more regional samples and a higher degree of genomic coverage.

| Implications for conservation of marine species in the face of climate change
The effects of global warming increase the spread and associated threat of nonindigenous species across the globe (Gardner et al., 2016;Hellmann, Byers, Bierwagen, & Dukes, 2008;Saarman & Pogson, 2015). A study by Wisz et al. (2015)  Prior to the current investigation, multiple studies have assumed Mytilus mussels in the Arctic to be exclusively M. edulis (Berge et al., 2005;Hansen, Hanken, Nielsen, Nielsen, & Thomsen, 2011;Jensen, 1905;Strand & Asmund, 2003). The identification of three Mytilus spp.
across the Arctic has implications for ecological and ecotoxicological research in the region. Mytilus mussels are extensively used in biological monitoring programs . However, interspecific differences in physiology and responses to environmental pollutants have been reported (Brooks, Farmen, Heier, Blanco-Rayon, & Izagirre, 2015;Fly & Hilbish, 2013), and thus, the lack of genetic knowledge could seriously affect the conclusions of ongoing biological monitoring. We therefore emphasize the importance of applying genetic tools to document species status, when conducting ecological, ecotoxicological, and physiological studies.
Moreover, assuming that the distribution and genetic connectivity between regions observed in this study is to be a first approximation Krause-Jensen for collecting mussels. The project was funded by the 15th of June Foundation, the Program of Russian Academy of Sciences "Fundamental Research to the Development of Arctic", the Framcentre flagship Fjord and Coast project "Life at the Edge", and a Norwegian Research Council project (project nr 225044). This work is a contribution to the Arctic Science Partnership (ASP), asp-net.org.

SUPPORTING INFORMATION
Additional Supporting Information may be found online in the supporting information tab for this article.  Table 1.
APPENDIX 1 H e , H o and allelic richness for the "Mytilus edulis" data set. APPENDIX 4 Structure analysis for the "Mytilus edulis" "neutral" data set.
APPENDIX 3 Confidence intervals 95% confidence intervals for the inferred ancestry for simulated parentals and hybrids analyzed with Structure (K = 4