Spatial genetic structure in a crustacean herbivore highlights the need for local considerations in Baltic Sea biodiversity management

Abstract Incorporating species' eco‐evolutionary responses to human‐caused disturbances remains a challenge in marine management efforts. A prerequisite is knowledge of geographic structure and scale of genetic diversity and connectivity—the so‐called seascape genetic patterns. The Baltic Sea is an excellent model system for studies linking seascape genetics with effects of anthropogenic stress. However, seascape genetic patterns in this area are only described for a few species and are completely unknown for invertebrate herbivores, which constitute a critical part of the ecosystem. This information is crucial for sustainable management, particularly under future scenarios of rapid environmental change. Here, we investigate the population genetic structure among 31 locations throughout the Baltic Sea, of which 45% were located in marine protected areas, in one of the most important herbivores of this region, the isopod crustacean Idotea balthica, using an array of 33,774 genome‐wide SNP markers derived from 2b‐RAD sequencing. In addition, we generate a biophysical connectivity matrix for I. balthica from a combination of oceanographic current models and estimated life history traits. We find population structure on scales of hundreds of kilometers across the Baltic Sea, where genomic patterns in most cases closely match biophysical connectivity, indicating passive transport with oceanographic currents as an important mean of dispersal in this species. We also find a reduced genetic diversity in terms of heterozygosity along the main salinity gradient of the Baltic Sea, suggesting periods of low population size. Our results provide crucial information for the management of a key ecosystem species under expected changes in temperature and salinity following global climate change in a marine coastal area.

conservation under climate change: A response to Hodgson. In addition, effects of large-scale gradual changes such as climate change are less profound in species with high dispersal capabilities, as they are better able to shift their distribution ranges through migration (Doerr et al., 2011). Species with poor dispersal abilities living in areas with physical barriers to dispersal are particularly vulnerable, so knowledge about where these barriers are located, and in which areas isolated populations are found, is two pieces of information critical to management efforts worldwide (Selkoe et al., 2016). In the marine environment, the realized gene flow resulting from various types of larval/propagule dispersal, as well as juvenile and adult movements, has for long remained unknown for a majority of species. The apparent lack of obvious dispersal barriers in the oceans led earlier biologists to conclude that marine populations show little fragmentation. However, a multitude of studies have later shown that many marine species are spatially structured on smaller scales than previously thought (reviewed in Selkoe et al., 2016). In order to better understand the location of dispersal barriers and their consequences, seascape genetic studies are performed by connecting biophysical dispersal models, ecology, and population genetic data (Galindo, Olson, & Palumbi, 2006). This information is highly relevant for managers in conservation of species, for example, in the establishment and/or evaluation of networks of marine protected areas (MPAs) (Underwood, Smith, Van, & Gilmour, 2009). Further, future species distributions will shift as a consequence of climate change, with important management implications. For example, Loarie et al. (2009) showed that due to climate change, on a global scale only 8% of protected areas designed today would still contain the habitat they were designed to protect in 100 years from now. In coastal areas, one strategy for mitigation would be to ensure that geographic areas likely to be strongly impacted by the changing environment are well connected to other suitable areas through dispersal corridors (Jonsson, Kotta, & Andersson, 2018), where one key feature would be protection of habitat-forming species.
The Baltic Sea (Figure 1) is an ideal area to study the interaction between biophysical connectivity, eco-evolutionary dynamics, and population genetic structure. This is primarily because oceanographic patterns are relatively well described in the Baltic Sea (Hordoir et al., 2019). The Baltic Sea is a multibasin marginal sea, with rather strong oceanographic barriers and sharp environmental gradients separating the different basins (Leppäranta & Myrberg, 2009).
There is also limited oceanographic connectivity between the Baltic and the adjacent North Sea (Moksnes, Corell, Tryman, Hordoir, & Jonsson, 2014), effectively isolating Baltic populations of marine organisms and reducing gene flow among coastal populations of marine species. Furthermore, the connection to the North Sea opened recently, about 8,000 years ago (Ignatius, Axberg, Niemistö, & Winterhalter, 1981), and while the majority of marine species in the Atlantic failed to colonize it completely, others have reached differently far into the Baltic Sea and along its strong salinity gradient (Ojaveer et al., 2010). These distributional differences may be due to, for example, differences in stress tolerance (Weinberger, Buchholz, Karez, & Wahl, 2008;Wrange et al., 2014), dispersal (Sjöqvist, Godhe, Jonsson, Sundqvist, & Kremp, 2015;Urho, 1999), or reproductive abilities across species (Jaspers, Møller, & Kiørboe, 2011). In addition, what seems common to many of these marine species is a strong genetic divergence between the North Sea and Baltic Sea populations, with a sharp boundary at the North Sea-Baltic Sea environmental transition zone, which is located just southeast of the Danish straits (Johannesson & André, 2006). There are multiple po- The Baltic Sea has recently been proposed as a reference and model system (a "time-machine") for other coastal areas as it is ahead of other areas with respect to issues such as impacts of climate change, habitat loss, eutrophication, pollution, and over-fishing, as well as being at the forefront when it comes to both monitoring, scientific investigations, and management measures (Reusch et al., 2018). The Baltic Sea is one of the fastest-warming regions in the world (Reusch et al., 2018), which will be accompanied by additional increased precipitation resulting in a strong reduction in salinity (Meier, Hordoir, et al., 2012a). This is expected to lead to substantial range shifts of species with loss of marine species, expansion of freshwater species, and continued rapid introduction of new species (Ojaveer & Kotta, 2015). For management of individual species, knowledge of dispersal barriers and population fragmentation is critical to improve measures. One major concern is whether northern Baltic Sea populations possess dispersal capabilities and adaptability to higher temperatures, or to salinity reductions and shift their ranges to the south to escape extinction (Johannesson, Smolarz, Grahn, & Andre, 2011;Jonsson et al., 2018).
However, in different species of flatfish this does not seem to be the case (Le Moan, Gaggiotti, et al., 2019b;Nielsen, Nielsen, Meldrup, & Hansen, 2004), highlighting that population genetic patterns often are species-specific and depend on traits such as dispersal potential and population sizes.
One missing key piece of the puzzle of seascape genetics in the Baltic Sea is the mesograzer herbivore guild-small, mostly crustacean grazers, which use macroalgae for food and shelter and in turn provide an important food source for fish. Three species of isopods of the genus Idotea are arguably the most important herbivores in the depauperate Baltic Sea ecosystem (Leidenberger, Hårding, & Jonsson, 2012), the most common one being Idotea balthica. They are generalist grazers on different algae and sea grasses, but also have the uncommon ability to survive and grow solely on a diet of brown algae of the genus Fucus (Bell & Sotka, 2012). Indeed, in some areas in the Baltic, densities of I. balthica can rise to astonishing numbers (more than 80 individuals/100 g wet weight of Fucus) and their grazing can severely decimate Fucus in a local area (Engkvist, Malm, & Tobiasson, 2000). While these invertebrates thus form an important ecological function, no previous studies have examined population genetic patterns or dispersal potentials of these species in the Baltic Sea. Idotea species brood their young and have no pelagic larval dispersal phase, which would largely prevent dispersal. On the other hand, adults are strong swimmers, and in addition, they have the potential to raft long distances and cross open-sea barriers attached to free-floating algae (Rothausler, Corell, & Jormalainen, 2015;Thiel & Gutow, 2005), which might counteract the lack of larval dispersal.
To investigate the dispersal potential of mesograzers in the Baltic Sea, we here use a seascape genetics approach applied on I. balthica.
Our objectives are to describe the population genetic structure and identify barriers to dispersal in the Baltic Sea region. We also aim to locate hotspots of genetic diversity, which could be considered to be of higher protective value, as well as isolated areas of low diversity and low connectivity which might be at higher risk of local extinction. In order to achieve these goals, we combine biophysical F I G U R E 1 Map of the study area (the Baltic Sea), with collecting locations marked with dots colored by region transport models and population genomic data of high spatial resolution and test for concordance between genetic and biophysical model data. If patterns are concordant, it will be possible to extrapolate genetic patterns geographically in order to generate maps that could be used in management efforts. Thus, the information produced here provides important input for future considerations in the Baltic Sea MPA network.
We collected isopods from outside of the entrance of the Baltic Sea to the innermost areas of the species' distribution range (the Bothnian Sea and the Gulf of Finland), with almost half of the sites being MPA sites. Using 2b-RAD sequencing (Wang, Meyer, McKay, & Matz, 2012), we identified approximately 35,000 SNP markers randomly distributed across the genome, which we used to estimate population differentiation, barriers to dispersal, and areas of high and low genetic diversity. We also estimated a connectivity matrix using a biophysical model based on oceanographic circulation patterns combined with hypothesized dispersal traits of I. balthica, which we compared to the population genetic structure indicated by the genomic data. We hypothesized that the lack of larval dispersal in this species has resulted in strong population fragmentation along the Baltic Sea coast and that the population genetic structure follows an isolation-by-distance pattern where oceanographic connectivity largely explains the genetic patterns.
We explicitly compared genetic diversity within and outside of HELCOM MPA sites. According to international conventions, placement of MPA sites should take genetic diversity into account (e.g., EU Framework Directive 2008/56/EC). If the current marine MPA network within the study area has been optimally designed in this respect, we would expect samples from MPAs to contain higher genetic diversity than non-MPA ones. Further, we hypothesize that the I. balthica population in the Baltic Sea has a reduced genetic diversity compared to outside locations, due to periodically low population sizes either before, during, or after the colonization of the Baltic, as has been observed in other species (Johannesson & André, 2006).  (Table 1). From hereon, a group of isopod individuals collected at one specific sampling location and time will be referred to as a sample. Animals were decapitated and stored in 95% ethanol at −20°C. From most locations, DNA from 20 individuals was extracted, four of which twice for technical replication, using a Qiagen Blood & Tissue kit and following the standard protocol, with the addition of RNAse during the last 30 min of the tissue lysis step in order to avoid RNA contamination. DNA quantity was measured using a QuBit dsDNA BR assay, and quality was assessed through gel electrophoresis.

| Genotyping
2b-RAD libraries (Wang et al., 2012) were prepared from the DNA using a modified version of the laboratory protocol developed by Mikhail Matz, available at: https ://github.com/DeWit P/BONUS_ BAMBI_IDOTEA. In brief, 100-200 ng of DNA template was fragmented using the type 2b endonuclease enzyme BcgI, after which adapters were ligated to the ends of the excised 36-bp fragments.
Fragments were then amplified with barcoded adapters, after which they were pooled equimolarly into 24-sample population pools (20 individuals + 4 technical replicates). All pools were sequenced in an Illumina HiSeq 2,500 machine, 50bp single-end, at the Swedish National Genomics Infrastructure's SNP & SEQ platform at Uppsala University.
All bioinformatic analyses were run on the University of Gothenburg computer cluster "Albiorix" (http://albio rix.bioenv. gu.se/). All commands used in the analyses can be found here: https ://github.com/DeWit P/BONUS_BAMBI_IDOTEA. An unpublished draft genome assembly for I. balthica was used as a reference for mapping the 2b-RAD data, using bowtie2 (Langmead & Salzberg, 2012) (information for the genome project can be found at https ://github.com/The-Bioin forma tics-Group/ Idotea_genome_project).
The standard GATK pipeline (McKenna et al., 2010) was followed for SNP calling (using the UnifiedGenotyper), including InDel realignment, quality score recalibration (BQSR), and variant score quality recalibration (VQSR), using sites identically genotyped across all technical replicates as a "True" training set for the machine learning algorithm. Poorly genotyped individuals and SNP sites genotyped at < 80% of all individuals were filtered out, as well as highly (>75%) heterozygous sites. The data set was pruned in order to keep only one SNP site for each RAD fragment. Finally, technical replicates examined for concordance in genotype estimates were discarded from further analysis.

| Among-sample summary statistics and differentiation
Mean pairwise F ST , as well as sample allele frequency differentiation p-values using Fisher's exact probability test, was calculated for all loci combined for all sample pairs using GENEPOP (Raymond & Rousset, 1995;Rousset, 2008) (Option 3, "Genic differentiation for all pairs of populations," with default Markov chain parameters). This method computes sample allele frequencies and seeks to reject the null hypothesis that alleles in both samples are drawn from the same distribution. Pairwise F ST values were transformed using a metric multidimensional scaling (MDS) approach, which summarizes the genetic distances in a multidimensional plane with as few dimensions as possible, after which the first two dimensions were plotted using R. Genetic diversity indices (F IS and expected heterozygosity H E ) were calculated using Arlequin v.3.5.2.2 (Excoffier & Lischer, 2010). Mean expected heterozygosity for each sample as a response variable depending on location-specific environmental factors (summer salinity, winter salinity, summer temperature, winter temperature, MPA status, and distance from the North Sea following the coastline (see Table 1)) were modeled using multiple regression in R, after which the significance of each explanatory variable was tested using ANOVA. To do this, environmental data for the top 6 m were extracted from the Rossby Centre Oceanographic circulation model (Meier, Döscher, & Faxén, 2003) and averaged across 1995-2004 for all locations except for W-1, which was outside the domain of the model. For the location W-1, empirical data from the ICES database (ices.dk) for depths ≤ 5 m were averaged for the years 1995-2004 for the nearest locations present in the database.
Summer and winter values were averages for June-August and December to February, respectively. Summer and winter salinity were found to be strongly correlated, and thus, only summer salinity was used as a factor. For all other factors, variance inflation factors were < 3, and thus, they were kept in the model.

| Principal components analysis
An identity-by-state (IBS) distance (1-IBS) matrix was calculated on an individual basis from the full 33,774 SNP dataset, using plink 1.9 (Purcell et al., 2007). A Canonical Analysis of Principal coordinates (CAP) was performed using the Vegan R package, assigning multidimensional coordinates to each isopod individual. This method is similar to a PCA or an MDS, but it rotates the eigenvector axes in such a way as to maximize the among-sample differences. Also, a test of the significance of sample as a factor was tested through an ANCOVA within the Adonis R package.
To identify loci potentially affected by natural selection, a global F ST outlier analysis was performed using two methods, Bayescan v 2.1 (Foll & Gaggiotti, 2008), using a cutoff value of an FDR-corrected p-value of 10 -4 to determine significance, and OutFLANK (Whitlock & Lotterhos, 2015), which is more stringent than BayeScan in outlier calling, and thus minimizes false positives. The output from the two methods was then examined for overlap. IBS matrices and CAP analyses were calculated separately for the shared outliers and nonoutliers, in order to determine whether outliers had disproportionate or distortive effects on the genetic structure analyses. As outliers and nonoutliers showed similar patterns, outliers were kept in downstream analyses.

| Geographically explicit population inference
Genotype likelihoods for all SNP sites were extracted from the vcf files using bcftools, which were input to NGSadmix (Skotte, Sand Korneliussen, & Albrechtsen, 2013) for population cluster analysis. NGSadmix uses a probabilistic framework to infer ancestry. K TA B L E 1 Collecting location information. Broad-scale geographic regions are color-coded as follows: Green-Western Baltic Sea (Skagerrak/Kattegatt/Danish Straits area) (W); Blue-Swedish Baltic Sea proper (SBA); Purple-Swedish Bothnian Sea (SBO); Red-Finnish Bothnian Sea (including the Archipelago Sea) (FBO); Black-Gulf of Finland and Estonia (GFE). Color codings are used consistently throughout the manuscript. Summer and winter values for temperature and salinity are averages 1995-2004 for June-August and December-February, respectively Finally, variability in genetic diversity was modeled and extrapolated geographically using the software package EEMS (Petkova, Novembre, & Stephens, 2015). This software uses an isolation-by-distance model with stepping stones ("demes") as a null model in order to infer geographic areas of high and low genetic diversity using the parameter q, a parameter which reflects the expected genetic differences among two individuals collected at the same site (Petkova et al., 2015) (not to confuse with the ancestry coefficients from the NGSadmix analysis above, which are also denoted "q"). A full-rank distance matrix was generated by imputing

| Biophysical model
A biophysical model was used to estimate dispersal probabilities and multigenerational connectivity calculated from stepping-stone dispersal across generations. All results below from the biophysical model are referred to as "connectivity," to be distinguished from results obtained from genetic data, which are referred to using the terms "gene flow" and "population structure." The biophysical model is based on the oceanographic circulation model NEMO-Nordic that produces water velocity fields with spatial resolution of 3.7 km in the horizontal and 3-12 m in the vertical, and a temporal resolu-

| Genetic versus. biophysical connectivity
The correlation between genetic differentiation (pairwise F ST ) and multigenerational connectivity estimated from the biophysical model was tested with Mantel tests (Mantel, 1967). Correlations with both untransformed and logarithmically transformed (log 10 (x + 1e-50)) connectivity matrices (64 generations) were tested. Because pairwise F ST data are symmetric, we extracted symmetric connectivity matrices by using either the minimum connectivity (minimum of i to j and j to i), the maximum connectivity, or the mean connectivity. In addition, distributions of pairwise F ST values within each region and between regions were plotted as boxplots in R, and a regression analysis was performed to examine potential isolation-by-distance patterns between mean connectivity and pairwise F ST .

| Accuracy of genotyping
The 2b-RAD libraries generated µ = 8.64 M reads ± SD 3.72 Mreads per individual, and mapping rates (not including multiple hits) were µ = 20.23% ± SD 2.26% (Table S1). This rather low percentage is not unexpected given the short nature of the 2b-RAD tags and the incomplete genome assembly used (NG50 = 12.9 kb based on assumed haploid genome length 1 Gb; 338 k contigs;  (Table S1) from the 31 different samples (Figure 1; Table 1). Genotyping correspondences among technical replicate individuals were high, µ = 97.70% ± SD = 2.49% (Table S2), and missing data were evenly distributed across samples, with only 9 loci not genotyped in any individual in one sample, and no loci with all missing data from more than one sample (Table S3).

| Among-sample summary statistics and differentiation
Overall, considering all loci, F IS values were close to zero in all 31 samples, indicating within-location panmixia (mean = 0.020 ± SD 0.029; n = 31; Table S4).
Pairwise F ST values among samples ranged from 0.001 to 0.068 (mean = 0.025 ± SD 0.012; n = 465) (  Figure S1). In an MDS plot ( Figure S2), axis 1 and axis 2 together explain much of the variance in pairwise F ST (61% and 24%, respectively). In this analysis, the Swedish Bothnian Sea samples (SBO-1 and SBO-2) in particular diverge strongly from all others along both primary axes.
There was significant genetic differentiation (based on sample allele frequencies) among most samples (Fisher's exact test p « .001; Table S5 below diagonal). Exceptions to this pattern were found in the Skagerrak/Kattegat/Danish straits area (W-1-9; marked blue in Table S5) and along the northern part of the Finnish Bothnian Sea coast (FBO-1-6; marked green in Table S5), where many samples were not significantly differentiated from each other, indicating that these samples were taken from the same population. Mean expected heterozygosity per sample (Table S4)

| Individual-based inference of genetic admixture and population structure
Our canonical analysis of principal coordinates ( Mean summer salinity was an equally strongly significant factor in the CAP axis 1 loading in both outliers and in nonoutliers (p < 2E-16; Table   S6). Due to the similar patterns observed in both types of loci, all loci were analyzed together for inferring genetic structure.
The combined CAP analysis clearly separated samples, or groups of samples, from each other, matching the geographic distribution of the collecting locations almost perfectly (Figure 3).
Indeed, sample was a strongly significant factor, explaining 12% of the genetic variance (ANCOVA, p = .001; Table S6). Just as in Finally, within-deme effective genetic diversity (q) (Petkova et al., 2015) was higher in the Skagerrak area than in the Baltic Sea, with particularly low values on the Swedish Bothnian Sea coast (SBO-1 and SBO-2) and in western Estonia (GFE-7) (Figure 5c).

| Biophysical model predictions of connectivity
The asymmetrical multigeneration connectivity matrix (64 generations) based on the biophysical particle model ( Figure 6, Table S7) showed that there were certain regions with high internal connectivity, which were less well-connected to each other. The western Baltic, the  Table S6).

| Concordance between biophysical connectivity and genetic structure
The population genomic data produced here shows the presence of genetically differentiated populations of Idotea balthica within the Baltic Sea. In some of the regions, there are significant population differences on scales of less than 100 km. This contrasts to patterns in some fish species such as sticklebacks (DeFaveri,  and herring (Lamichhaney et al., 2012) where populations are genetically structured at broad geographic scales in the Baltic Sea, while in other fish species, such as perch and whitefish, genetic differences are found also on relatively small scales (Olsson, Florin, Mo, Aho, & Ryman, 2012;Olsson et al., 2011). Instead, the strong genetic structure in I. balthica is similar to the pattern present in the seaweed Fucus vesiculosus (Ardehed et al., 2016), which has an overlapping distribution with I. balthica and provides food and shelter to the isopod (Kotta et al., 2019).
Further, the genetic patterns match the biophysical connectivity model closely, suggesting that passive dispersal is an important factor influencing the genetic structure in this species. Without pelagic larval dispersal, long-distance dispersal is most likely F I G U R E 4 (a) Admixture plot at K = 3. Bothnian Bay (the "North Quark"; see Figure 1). Idotea balthica is quite common on the Finnish side (FBO-1,2), and this area is shallow and may provide good stepping-stone dispersal for shallow-water organisms, so even species with limited dispersal ability should be able to cross over to Sweden from the Finnish side. However, isopods are not found on the Swedish side (Leidenberger et al., 2012). In fact, our SBO-2 site, approximately 300 km further south, is the most northern point of the species distribution on the Swedish coast. There are available habitats (Fucus belts) north of this site, although the distribution of Fucus spp. becomes increasingly patchy. Niche modeling studies (Kotta et al., 2019;Leidenberger, Giovanni, Kulawik, Williams, & Bourlat, 2015) indicate that the northern Bothnian Sea coast should be inhabitable to the isopods, rendering the reason for their absence in this area an interesting topic for future research. are gone today (Boström, Baden, & Krause-Jensen, 2003). Further sampling of Idotea along these shores will be necessary in order to illuminate the patterns of gene flow here.

| Main population units of Idotea balthica in the Baltic Sea
Along the Swedish coast, starting from Öresund (W-4), samples show a strong step-wise isolation-by-distance divergence going in to the Baltic. Such a pattern is expected following a colonization front moving into the Baltic Sea along the coast, with founder effects leading to increasing levels of differentiation (Excoffier, Foll, & Petit, 2009

| Reduced genetic diversity
The low heterozygosity in I. balthica in the Baltic Proper corresponds to a 12% drop compared to the most diverse sample on the Swedish west coast (W-1). Interestingly, this corresponds very well to the average drop in nuclear genetic diversity of Baltic Sea populations of various marine species (11%-12%) (Johannesson & André, 2006).
One possible explanation is the influence of genetic drift due to periodically low population sizes purging diversity in the Baltic Sea (Johannesson & André, 2006). However, some observed patterns could also be explained by secondary contact and introgression among historically separated lineages enhancing diversity locally.
Secondary contact zones in this area have been recently described for three species of fish (Le Moan, Gaggiotti, et al., 2019b). However, as the genomic reference used in this study was highly fragmented, we could not investigate genomic patterns of introgression in our dataset. Hopefully, future improvements in genome sequencing will allow for such studies in the future.
Within the Baltic Sea, the western part of Estonia stands out as a region of lower genetic diversity, in large part due to low diversity in the sample GFE-7 ( Figure 5c; Table S4). GFE-7 was collected on the highly exposed west coast of Saaremaa island, where also the biophysical connectivity analysis indicates a lower input of propagules compared to other sites along the Estonian coast ( Figure 6).
Along the isolated Swedish Bothnian Sea coast, diversity is also low, especially in sample SBO-1 (Gräsö island). These locations could be considered important from a MPA network perspective, as local reductions in diversity are due to either low population sizes, low migration, or a combination of the two factors.

| Management and conservation implications
The reduced within-species genetic diversity, and the need for adaptations to the rapidly changing Baltic Sea, call for management and conservation actions (Johannesson et al., 2011 (Laikre et al., 2016). Useful information of species genetic diversity is currently accumulating, but managers hesitate to use genetic diversity in dayto-day operations as they find this type of information difficult to interpret (Sandström, Lundmark, Andersson, Johannesson, & Laikre, 2019). One way to overcome this problem may be to provide easy-tointerpret maps of genetic diversity and population subdivision, which can be used to implement management efforts in several ways: First, the genetic clusters identified in Figure 5b can be considered as a good starting point for defining management units in this organism, where the Kattegat/Skagerrak/Danish Straits area, and also the west coast of Finland might be managed as single units, while in the other parts of the study area the observed population divergence indicates a need for more local management units. It is important to note, however, that this study does not include data from the island of Gotland, nor from the coasts of Poland, Latvia, and Lithuania. In these areas, there are (rare) historical records of I. balthica (Leidenberger et al., 2012). Large parts of the southern Baltic coast consist of sandy areas with limited Fucus-habitat, yet I. balthica can also be found in eelgrass (Zostera marina) beds, so it cannot be excluded that there are additional genetic clusters present there. Further sampling efforts will be necessary in order to resolve this issue.
Second, highly diverse populations of high priority for conservation can be identified from the genetic data. Genetic diversity is a proxy for the long-term adaptability of a population and can also be considered a sign of a relatively large and stable population size.
In this dataset, sites with higher than average genetic diversity (as Finally, as significant warming and decrease in salinity are expected in the Baltic Sea before the end of this century, local extinctions of isopod populations throughout the area, especially near current range margins, are to be expected (Kotta et al., 2019). These could in turn have immediate negative consequences for ecosystem functions as isopods are important as food for, for example, local populations of herring and perch. While local extinctions might be unavoidable at the most northern areas, populations further south could be supported by influx of genetic material, which might contain low-salinity tolerant traits, from the north. A potential management strategy to minimize the detrimental effects of climate change would thus be to facilitate dispersal southward along the salinity gradient by protecting suitable habitats (continuous Fucus belts), especially along the Swedish Bothnian Sea and northern Baltic Sea coasts, as it may otherwise be difficult for low-salinity adapted genotypes to track the receding salinity gradient (Kotta et al., 2019). This is especially important as our data show that the Swedish Bothnian Sea coast is genetically isolated. There is currently a lack of designated MPAs in this particular area aiming to protect Fucus spp., an overview of which seems motivated to maintain stepping-stone possibilities for north to south gene flow in I. balthica. It might also be considered to physically translocate northern populations to more southerly areas. While this type of manipulation is controversial, it might prove necessary under rapidly decreasing salinities in order to sustain more southern Baltic Sea populations of isopods. Surface salinity has been decreasing since 1980 and may reach 3 psu in large parts of the Baltic Proper before the end of the century (Meier, Andersson, et al., 2012b;Meier, Kjellström, & Graham, 2006).
However, it is still unclear exactly how the interactions of all projected environmental changes will affect species distributions within the next hundred years, so consequently it will be important to establish genetic monitoring, in order to be able to give early warnings for loss of genetic diversity in of both Fucus spp. and I. balthica.

DATA AVA I L A B I L I T Y S TAT E M E N T
All raw sequencing data were submitted to the NCBI Short Read Archive, BioProject PRJNA551577. Vcf files containing all genotype calls were submitted to the Dryad online repository (https ://doi.