Seascape genetics and biophysical connectivity modelling support conservation of the seagrass Zostera marina in the Skagerrak–Kattegat region of the eastern North Sea

Abstract Maintaining and enabling evolutionary processes within meta‐populations are critical to resistance, resilience and adaptive potential. Knowledge about which populations act as sources or sinks, and the direction of gene flow, can help to focus conservation efforts more effectively and forecast how populations might respond to future anthropogenic and environmental pressures. As a foundation species and habitat provider, Zostera marina (eelgrass) is of critical importance to ecosystem functions including fisheries. Here, we estimate connectivity of Z. marina in the Skagerrak–Kattegat region of the North Sea based on genetic and biophysical modelling. Genetic diversity, population structure and migration were analysed at 23 locations using 20 microsatellite loci and a suite of analytical approaches. Oceanographic connectivity was analysed using Lagrangian dispersal simulations based on contemporary and historical distribution data dating back to the late 19th century. Population clusters, barriers and networks of connectivity were found to be very similar based on either genetic or oceanographic analyses. A single‐generation model of dispersal was not realistic, whereas multigeneration models that integrate stepping‐stone dispersal and extant and historic distribution data were able to capture and model genetic connectivity patterns well. Passive rafting of flowering shoots along oceanographic currents is the main driver of gene flow at this spatial–temporal scale, and extant genetic connectivity strongly reflects the “ghost of dispersal past“ sensu Benzie, 1999. The identification of distinct clusters, connectivity hotspots and areas where connectivity has become limited over the last century is critical information for spatial management, conservation and restoration of eelgrass.


| INTRODUCTION
Eelgrass (Zostera marina L.) is one of the most widely distributed species of seagrass in the northern hemisphere and the dominating species of the temperate North Atlantic (Short, Carruthers, Dennison, & Waycott, 2007). It is a benthic foundation species within shallow coastal areas where it provides habitat and numerous ecosystem services, such as stabilization of the coastline and improved water quality, increased fish production and uptake of carbon and nitrogen Orth et al., 2006). Large-scale losses of eelgrass have occurred worldwide (Waycott et al., 2009), including northern Europe (Boström et al., 2014), causing significant decreases of ecosystem services. For example, along the Swedish Skagerrak coast, a reduction of 120 km 2 has resulted in large losses of cod catches, and release of sequestered carbon and nitrogen, to an estimated total cost of more than 600 million US$ . In response, Z. marina has been classified as a "threatened and declining habitat" in the North East Atlantic and the Baltic Sea under regional marine conventions (HELCOM 2010;OSPAR 2017) and is also indirectly protected under several EU directives including the Habitats Directive (92/43/EEC) and its Natura 2000 network.
The largest known areal distribution of Z. marina in Europe is found in the Skagerrak-Kattegat-Belt Sea region in the eastern part of the North Sea (Boström et al., 2014). As in most parts of the North Atlantic, a dramatic loss of Z. marina occurred in the area in the 1930s, as a result of the wasting disease (Rasmussen, 1977). Eelgrass recovered in most areas in the 1960-1980s, but never obtained its historic distribution (Boström et al., 2014). In the following decades, Z. marina distribution in Denmark decreased again, probably as a result of nutrient pollution (Boström et al., 2014). It is estimated that eelgrass in Denmark today constitutes 10%-20% of its historic distribution and that the depth distribution has become more shallow by approximately 50%, resulting in a loss of most offshore populations (Boström, Baden, & Krause-Jensen, 2003;Boström et al., 2014). Along the Swedish Skagerrak coast, over 60% of meadows have been lost since the 1980s (Baden, Gullström, Lundén, Pihl, & Rosenberg, 2003;Nyqvist, André, Gullström, Baden, & Åberg, 2009). These losses have largely been attributed to coastal eutrophication and overfishing of large predatory fish, causing a trophic cascade and an increase in ephemeral macroalgae that smother Z. marina (Baden, Emanuelsson, Pihl, Svensson, & Åberg, 2012;Moksnes, Gullström, Tryman, & Baden, 2008). To mitigate the ongoing loss and assist recovery of Z. marina, a number of measures are presently being discussed, including the establishment of new networks of marine-protected areas (MPAs) and large-scale restoration SwAM 2015).
The Skagerrak-Kattegat region is particularly well suited to study the potential impact of large-scale decline on connectivity, due to the availability of detailed historic data of eelgrass as well as unusually good mapping of the current distribution (Boström et al., 2014). In addition, the oceanographic features of the area are unique, with a strong outflow of surface water from the Baltic Sea into the Kattegat, creating an asymmetric circulation along the coasts with strong effects on connectivity, also creating a barrier between the Kattegat and Skagerrak (Jonsson, Corell, André, Svedäng, & Moksnes, 2016;Jonsson, Nilsson Jacobi, & Moksnes, 2016;Leppäranta & Myrberg, 2009). Biophysical dispersal models, such as Lagrangian trajectory models (Cowen & Sponaugle, 2009;Grech et al., 2016;Selkoe et al., 2010), can model such directional dispersal based on biologically realistic assumptions, for example, time of propagule release, drift duration and depth, and may be superimposed on layers of habitat preference.
Here, we aim to establish a dynamic model of seascape population structure and connectivity for Z. marina meadows in the Skagerrak-Kattegat region. Our assessment includes a temporal comparison based on oceanographic dispersal modelling of extant and historical distribution data of Z. marina for the region, and we investigate the hypothesis that the large observed decline has resulted in decreased connectivity and lower genetic diversity. In addition, we examine how oceanographic and genetic barriers fit with present administrative borders such as countries and sea basins. We combine and crossvalidate genetic and hydrodynamic modelling approaches in order to infer the importance of dispersal in shaping population structure and to compare trade-offs and synergies offered by the integration of the two approaches as applied to management and mitigation.

| Study area and sampling
The study covers the eastern Skagerrak, Kattegat and Belt seas along the eastern part of the North Sea (54-59°N and 8-13°E) with a total area of 77,000 km 2 (Figure 1). For simplicity, we refer to the assessed area as Skagerrak-Kattegat. Sampling of Z. marina was guided by a previous oceanographic barrier analysis of the region (Moksnes, Jonsson, & Nilsson Jacobi, 2015;Nilsson Jacobi, André, Döös, & Jonsson, 2012) in which seven oceanographic clusters were identified (site names in Figure 1 follow the seven clusters). At least three sites were sampled from each of the seven oceanographic clusters (from here on we, will refer to these clusters as sampling areas) to ensure sampling within and across potential barriers to dispersal. At each site, 40 shoots were collected using a "roughly linear swim" (Arnaud-Haond, Duarte, Alberto, & Serrão, 2007) by snorkelling or diving. Within sites, intersample distance was maintained at 1-1.5 m (covering a distance of 40-60 m across a meadow), a standard distance for this species and an adequate compromise to capture diversity and structure, while minimizing resampling of the same genotype (Olsen et al., 2004). Sample depths ranged from 1.2 to 5.3 m. A total of 920 sampling units were collected. Among sites, pairwise distances ranged from ~10 to 400 km among the 23 sampling sites. Similar sampling scales between sites were maintained as much as possible to best detect a slow decline in allele frequency, that is, 8% of sites have a pairwise geographic distances of up to 50 km, 21% up to 100 km, 23% up to 150 km, 22% up to 200 km, 17% up to 300 km and 8% up to 400 km. Samples were collected in July and August 2016 except for two populations from Denmark (6-NH and 5-BO), which were collected during a previous sampling campaign in 2004 (Table 1). One to three 2-cm leaf pieces per shoot were selected from the clean, inner leaves near the base of the shoot and when necessary cleaned of epiphytes with a scalpel. Samples were dried and stored in silica crystals for later DNA extraction.
F I G U R E 1 Map of sampling sites for Zostera marina in the Skagerrak-Kattegat region of the North Sea (Table 1). Green dots indicate the extant mapped distribution of Z. marina; the area enclosed by the solid green line in western Kattegat shows the estimated historic distribution of Z. marina. The background heat map in (a) shows an interpolation of genotypic/clonal diversity; (b) allelic richness standardized for 21 genotypes (A 21 ) generated with the genetic diversity plugin (Vandergast, Perry, Lugo, & Hathaway, 2011) in ArcMap 10.3 (Desktop, 2014) and QGIS 2.18 (Quantum GIS Development Team 2013)

| DNA extraction and microsatellite amplification
DNA was extracted from ~20 mg of silica-gel-dried tissue in 96well plates using a silica-based Cetyl trimethylammonium bromide (CTAB) protocol (Hoarau, Coyer, Stam, & Olsen, 2007), except that samples were incubated in CTAB for 1 hr at 60°C. Twentytwo microsatellite loci were used: the original set of eight from Reusch, Boström et al. (1999),  and used in numerous genetic surveys of Z. marina, and 14 additional loci developed from expressed sequence tags (EST) libraries (Keil, 2011;Oetjen & Reusch, 2007). Primer sequences, multiplex combinations and concentrations are provided in Tables S1 and S2. Polymerase chain reactions (PCRs) were performed in 96-well microtiter plates using the Qiagen Kit Type-IT ® in a 6.2 μl reaction volume following the manufacturer's instructions. The reaction profile consisted of 95°C for 5 min followed by 30 cycles of 95°C for 30 s, 56°C for 1 min 30 s and 72°C for 30 s, with a final extension step of 60°C for 30 min.

| Microsatellite genotyping, removal of clones, data quality checks and discrimination power
PCR products were diluted 1:100 (apart from the "4-plex," which was used undiluted), and fragment analysis was performed on an Applied Biosystems 3730 DNA Analyser with a 350 ROX internal size standard added to each well. Fragments were scored automatically using GeneMapper ® (Life technologies) and re-checked by eye for each individual and locus. Samples with ambiguous or rare alleles were reamplified and re-genotyped for confirmation. We succeeded in amplifying all individuals at all loci.
Because seagrasses can spread clonally via rhizome extension, a genetic individual (genet) may consist of hundreds or thousands of shoots (ramets) covering several metres. Even though a sampling distance of 1-1.5 m is generally adequate for Z. marina (Olsen et al., 2004), it is no guarantee that the same genet might not be sampled more than once if large clones are present. Accordingly, duplicate multilocus genotypes (MLGs) were identified and removed using RClone T A B L E 1 Genetic diversity of Zostera marina at 23 locations in the Skagerrak-Kattegat region of the North Sea The 920 individuals sampled in Denmark and Sweden were assessed with 20 microsatellites. Population names are followed by the acronyms, latitude and longitude, the number of sampled ramets (N), the number of multilocus genotypes (MLG), genotypic richness (R) as MLG-1/N-1, allelic richness standardized to 21 genotypes (A 21 ) plus standard deviation (SD), not applicable (na) due to low number of MLGs, observed heterozygosity (H O ), expected heterozygosity (H E ) and the inbreeding coefficient (F), and standard error (SE). Numbers in bold indicate significant F values. (Bailleul, Stoeckel, & Arnaud-Haond, 2016)  Linkage disequilibrium (LD) and HWE were evaluated for each locus and across all loci in each population with Genepop 4.2 (Raymond & Rousset, 1995; 100 batches and 1,000 iterations per batch plus Bonferroni corrections). To test for neutrality of our loci, "outlier" analyses were performed using both Lositan (Antao, Lopes, Lopes, Beja-Pereira, & Luikart, 2008) and BayeScan (Foll & Gaggiotti, 2008). See Figure S1 for additional information.
To analyse the statistical power of our set of microsatellites to discriminate clonal replicates, we calculated the probability of identity (PI) in GenAlEx 6.5 (Peakall & Smouse, 2012) at each site and used
To gain a first impression of genetic structure of all MLGs without a priori population genetic assumptions, we used PCA as implemented in adegenet 2.0.1 (Jombart, 2008) in R 3.3.2 after using the scaleGen function. PCAs were also used to investigate potential outlier loci and for visualization of the distribution of our sampling sites in a larger geographic context (see Supporting Information). Population genetic differentiation was calculated as the proportion of shared alleles among populations, D ps ' = 1 -ps, in MSA 4.05 (Dieringer & Schlötterer, 2003).
We chose D ps , because it is free of equilibrium assumptions (Bowcock et al., 1994). We also calculated several standard variance-based measures of population differentiation (Tables S5-S7).
Structure cannot always identify clusters accurately when geographic sampling is discrete along clines and/or when isolation-by-distance (IBD) patterns or autocorrelations dominate the data (Chen et al., 2007). TESS addresses these issues using a spatially continuous prior based on the geographic coordinates of each individual. Therefore, our main analyses rely on TESS results, while Structure was carried out as an additional support. As we had geographic information only at the meadow level, we used TESS to calculate slightly adapted geographic coordinates for each individual. TESS was run using the conditional autoregressive Gaussian (CAR) admixture model, which assumes spatial autocorrelation of genetic differentiation, using the default value of 0.6 for the strength of the autocorrelation. We first ran a test with default settings for K max = 2-25 and then repeated TESS for a range of likely K max (2-7) with a burn-in of 10,000 sweeps followed by 25,000 sweeps, with 100 independent runs conducted for each K max .
The independent runs were averaged and compared to assess convergence. The average deviance information criterion (DIC) for each value of K max was used to evaluate the most likely number of genetic clusters by determining K max at which a higher number of parameters did not improve the model significantly. We used pophelper (Francis, 2016

| Directional migration
We use the term "migration" in the population genetic sense, which includes both successful movement and contribution to the local gene pool (Lowe & Allendorf, 2010). In contrast, we use "dispersal" when discussing movement based on oceanographic modelling to reflect a passive process of transport influenced by currents but that does not necessarily result in any contribution to a local gene pool.
Directional migration rates based on the microsatellite data were estimated using two different methods: DivMigrate-online (https:// popgen.shinyapps.io/divMigrate-online/) and GENECLASS2 (Piry et al., 2004). DivMigrate is an indirect approach that extends genetic differentiation to include a directional measurement by identifying migrants based on the geometric means of the allele frequencies in each population. Directional migration rates are then inferred from allele frequencies and genetic differentiation in pairwise comparisons of G ST (Sundqvist, Keenan, Zackrisson, Prodohl, & Kleinhans, 2016

| Mapping of suitable habitat
Data for present-day distribution of Z. marina were based on national inventories in Norway, Sweden and Denmark, and were obtained in geographic information system (GIS) format from the Norwegian Environment Agency, the Swedish County Administrative Board of Västra Götaland and the Danish Nature Agency. Along the Swedish Skagerrak coast, distribution was based on satellite image analyses (Envall & Lawett, 2016), whereas the distribution in other areas was based on national field surveys and monitoring sites. In the oceanographic modelling, all grid cells that intersected with eelgrass locations were used as sources in the particle tracking simulation and subsequent construction of the connectivity matrices.
In addition to mapping the present distribution of Z. marina, we also explored the effect of including the known historic distribution on multigeneration connectivity. These data were obtained from a recent analysis of historic records of Z. marina presence in the Kattegat collected around 1900 (Petersen, 1893;Rosenvinge, 1909). Tidal harmonics define the sea surface height and velocities at the boundaries, and Levitus climatology defines temperature and salinity (Levitus & Boyer, 1994). The model has a free surface, and the atmospheric forcing is a dynamic downscaling of the ERA40 data set.
Runoff is based on climatological data based on a number of different databases for the Baltic Sea and the North Sea.
Velocity fields were updated in the model domain every three hours, and the calculation of particle trajectories was performed with as well as duration that shoots stay afloat, were based on empirical field studies along the Swedish Skagerrak coast (Infantes & Moksnes, 2017;Källström et al., 2008). Different drift durations simulate that individual spathes with seeds on the same shoot mature at different times and that some negatively buoyant seeds may be dropped and sink, while the shoot continues drifting (Infantes & Moksnes, 2017).
Particle release was repeated for 8 years (1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002), representing years with a range of North Atlantic oscillation index values (NAO, Hurrell & Deser, 2010), which is known to correlate well with the variability in circulation pattern. In total, 2.5 million particle trajectories were included. Dispersal probabilities between all sampling sites, over a single generation, were calculated by summing all the trajectories starting in site i having end positions within site j, normalized by the total number of simulated trajectories from site i. We also calculated multigeneration connectivity where stepping-stone dispersal was allowed over 32 single-generation dispersal events by multiplication of the single-generation dispersal matrix with itself 32 times producing connectivity probabilities when summed over all possible dispersal routes . Stepping-stone dispersal was only allowed between grid cells that intersected with the known extant or reconstructed historical habitat distribution (see Figure 1). Stepping-stone dispersal over 32 generations was considered sufficient to span the approximate spatial scale (~500 km distance) of the model domain. In terms of the temporal scale, 32 generations may represent as little as 32 years when assuming annual sexual reproduction of Z. marina, or >1,000 years when assuming high levels of clonal reproduction and clone longevity (Reusch, Boström et al. 1999).

| Oceanographic dispersal barrier analysis
We employed a clustering method to identify partial dispersal barriers based on modelled dispersal probabilities in the seascape (Nilsson Jacobi et al., 2012). Only dispersal between areas with present or historic distribution of eelgrass was considered. This theoretical framework finds partially isolated clusters. Identification of clusters is formulated as a minimization problem with a tunable penalty term for merging clusters that makes it possible to generate population subdivisions with varying degree of dispersal restrictions. As the focus was to compare putative dispersal barriers to genetic differentiation of Z. marina, the mean connectivity between oceanographic clusters was set low (0.004).

| Isolation by "sea distance" and oceanographic distance
To test for IBD, we correlated genetic distance with "sea distance," defined here as the shortest path possible among sampling sites at sea without crossing land. We used the R package marmap (Pante & Simon-Bouhet, 2013) to calculate "sea distance." One distance value (between 5-BO and 6-NH) had to be adjusted manually to ensure that no land was crossed. Mantel tests were carried out using the R package ncf (Bjornstad, 2009) in R 3.2.2. Matrices were resampled 100,000 times and after log 10 transformation of "sea distance." To test for isolation by oceanography (IBO), we correlated genetic distance with minimum oceanographic dispersal probabilities (calculated from the model described above) defined here as oceanographic distance. We used minimum dispersal probability to generate a symmetric matrix of dispersal, because it may arguably be best correlated with (symmetric) geographic and genetic distance (Wrange et al., 2016). We also correlated directional dispersal probabilities with asymmetric (genetic) migration rates in a Mantel test adapted for asymmetric matrices (Matlab 2016a, Mathworks Inc).

| Network analyses
Network analysis is a graphic approach with many applications, one of which is to understand landscape patterns of connectivity and prioritize areas for conservation (Engelhard et al., 2016;and references therein). We used networks to examine connectivity both for genetic (D ps ) and oceanographic distance applied to modelled dispersal probability matrices for single-generation-extant, multigeneration-extant and multigeneration-historic dispersal probabilities, and to highlight sites that are central to connectivity. Networks were drawn using the R packages igraph (Csardi & Nepusz, 2006) and popgraph (Dyer, 2014), where nodes represent populations and edges the pairwise distance among populations.
Thresholds were chosen systematically following the "intermediate threshold" method of Greenbaum, Templeton, and Bar-David  and minimum multigeneration-historic dispersal probability = 1e-12. Use of these thresholds resulted in the loss of five populations from each network, as also seen in the Bayesian clustering analysis ( Figure 2). Only edges with genetic distances below the threshold and dispersal probabilities above the threshold are shown.

| Genetic data quality checks and discrimination power
We identified 14 to 40 MLGs per population ( The outlier analyses identified several loci to be potentially under balancing and positive selection ( Figure S1). As their exclusion did not alter PCAs ( Figures S3 and S4), they were retained. A PCA that included additional sites from the Baltic and the North Sea did not reveal any indication that our study area may represent a secondary contact zone between genetically differentiated Baltic and North Sea Z. marina meta-populations ( Figure S5).
The probability of identity (PI) by chance for a 20-locus MLG was low ranging from 5.9 × 10 −7 in population 3-GH to 3.5 × 10 −10 in population 2-MV. The probability for detection of sibs was higher ranging from 1 × 10 −3 to 6 × 10 −5 . Statistical power simulations in POWSIM of the 20-locus set indicated a 100% probability of detecting an F ST as low as 0.0025, and the α error (false significance) was close to the intended value of 0.05 (Table S3).

| Genetic diversity and population differentiation
Genotypic and allelic diversity was found to be high overall ( Figure 1 and Table 1) indicating a dynamic and diverse environment characterized by predominant sexual reproduction. One MLG was shared among two populations (Gottskär, 3-GS and Gåsö, G-SG) separated by 120 km. This is within our estimates for single-generation northward oceanographic dispersal and reattachment at Gåsö. This genotype has three alleles at three loci not otherwise found at Gåsö. Population differentiation using pairwise shared allele distances (D ps ; 0.09-0.36) or F ST (0.01-0.21) was significant (p < .05; Table S4 and S5) and consistent with strong overall population structure.
Further characterization of the genetic population structure started with a PCA that generated a horseshoe-shaped cloud of apparently little differentiation based on the first two axes ( Figure S2) and even less on the third axis (not shown). Such a pattern is expected under scenarios in which allele frequencies are locally correlated and thus covariance is decaying with geographic distance in an IBD pattern (Frichot, Schoville, Bouchard, & François, 2012).
Continuing on, the spatial Bayesian analysis in TESS suggested genetic population subdivision into three clusters (K max = 3, Figure 2), which was very similar to the K = 3 scenario of Structure ( Figure S6 and S7 Figure S7).
Under the two-cluster partition suggested by Structure, the two smaller TESS clusters 1 and 2 are depicted as one.

| Directional migration rates
Directional migration based on the genetic data was estimated in two ways. Based on DivMigrate (Table S8 and Figure S8a), directionality was stronger from south to north with eight sites identified as sources and five identified as sinks. The genetic assignment test based on GENECLASS2 (Table S9 and Figure S8b) identified 31 first-generation migrants, of which only seven could be assigned to other sampling sites (Table S9). Here, directionality was predominantly south to north and west to east. Both methods confirm that long-distance dispersal occurs.

| Oceanographic dispersal based on particle modelling
The biophysical particle modelling indicated dispersal up to 200 km in a single generation, consistent with other estimates for Z. marina (Harwell & Orth, 2002;Källström et al., 2008), and more than 300 km when allowing stepping-stone dispersal over multiple generations.
Few particles (modelled seeds) dispersed between sample sites during a single generation, although local retention within the same sampling occurred ( Figure S9). The most northern sites (names starting with 1 and G) received particles from many other meadows; sites from the sampling areas 6 and 7 supplied particles to most other sites.
Probabilities for multigeneration oceanographic dispersal based on the present-day distribution were much lower than those based on single generations, but there were nonzero probabilities of dispersal among all sampling sites ( Figure S10). The site 3-HH was the best source for particles, while 1-ST and 6-LM supplied few particles to other sites. Inclusion of the historic distribution of Z. marina ( Figure 1) in multigeneration dispersal modelling only slightly changed the overall picture, but with important differences at certain sites. For instance, 6-LM and G-SG acted as much stronger sources and 6-GH, 6-LM, 6-NH and 4-HO received considerably more particles in the past ( Figures S10-S12). Notably, the lowest dispersal probabilities increased by two orders of magnitude compared to multigeneration dispersal based on the present-day distribution alone.

| Oceanographic dispersal barrier analysis
At the chosen threshold, the minimization algorithm applied to the multigeneration dispersal matrix (including the historic habitat) generated six oceanographic clusters with partial barriers among them ( Figure 2b). Connectivity within oceanographic clusters was approximately 100 times greater than among clusters. Four barriers were identified ( Figure 2): (i) at 58°N, which spatially coincides with the green genetic cluster and the division between the Kattegat and Skagerrak; (ii) a barrier along the Swedish Kattegat coast encompassing the red genetic cluster; (iii) at 57°N roughly following the gradient from "pure" to high genetic admixture observed within the blue cluster; (iv) a barrier across the south-west corner of the Kattegat.
This last barrier was not reflected in the genetic cluster analysis, but the asymmetric migration analysis also indicated low gene flow in this area ( Figure S8).

| Isolation by "sea distance" and oceanographic distance
Geographic distance defined as "sea distance" (see Material and Methods) ranged from ~10 to ~400 km. A significant pattern of isolation by "sea distance" (IBD) with genetic differentiation (D ps ) is evident (Table 2). While oceanographic connectivity based on singlegeneration dispersal probability was not correlated with the genetic differentiation measure D ps , multigeneration-extant and historic dispersal were strongly correlated (three different scenarios of IBO). IBO was strongest for the dispersal probability based on the historic distribution of Z. marina, which reached a correlation coefficient as high as −0.59 (Table 2). Patterns of F ST -related genetic indices were similar (not shown). IBO was also observed for the correlation between genetic asymmetric migration rates and the directional dispersal probabilities (Table 2). For both genetic differentiation indices, correlation coefficients are much higher for multigeneration compared to singlegeneration dispersal probabilities. This indicates that stepping-stone dispersal over several generations can explain genetic differentiation better than single-generation dispersal probability, which is limited by geographic distance (Table 2). Correlations are further improved when considering the historic distribution of Z. marina (Table 2).

| Network analyses
The four network analyses in Figure 3  The network analysis also allowed visualization of populations central to connectivity (Greenbaum & Fefferman, 2017;Rozenfeld et al., 2008). The oceanographic and genetic network analyses indicated that meadows from south-western and central Kattegat (sampling areas 4, 5 and 7) are central to connectivity (Figure 3). The offshore populations from the island Laesø (7-LS and 7-480) are of greatest interest as they are located in an area of historically extensive eelgrass meadows (and still retain high allelic richness; Figure 1).

| Distribution of genetic diversity of Z. marina in the Skagerrak -Kattegat
Our initial hypothesis of reduced genotypic and allelic diversity as a consequence of the massive losses and fragmentation over the past century was not borne out. Rather, the Skagerrak-Kattegat region harbours some of the highest diversity for Z. marina in Europe (J.L. Olsen, unpublished data). Interestingly, the highest values of allelic richness are found in the centre of the Skagerrak-Kattegat region around the Laesø Islands (7-LS, 7-480), where eelgrass was historically abundant (Figure 1b). These observations are consistent with probable glacial refugia (Maggs et al., 2008) and the original postglacial colonization of the nascent North Sea basin (Hewitt, 2000;Maggs et al., 2008;Olsen et al., 2004) including the Skagerrak-Kattegat, when the current system was established ca. 8,000 years ago (Gyllencreutz, Backman, Jakobsson, Kissel, & Arnold, 2006). At that time, the Baltic was still an isolated, freshwater ice lake and colonization of both areas came most likely from the south (Ireland, Brittany, Iberian tip), although a high North refugium in northern Norway cannot be ruled out for macrophytes in general (Coyer et al., 2011;Maggs et al., 2008;Olsen et al., 2013). No evidence for a secondary contact zone, commonly observed to coincide with biogeographic regions (Gagnaire et al., 2015), was evident between genetically differentiated Baltic and North Sea populations and the Skagerrak-Kattegat (see Figure S5).

| Population genetic structure and connectivity
The 23 sampling sites form three distinct genetic clusters based on the TESS analysis ( Figure 2). The small (red) genetic cluster consisting of the two sites, South Kråkerön (K-KR) and North St. Överön (K-SO), is located in the Marstrand area, which has lost an estimated 93% of meadows since the 1980s and losses continue .
Although these two sites are genetically isolated from the other clusters, they exhibit high allelic diversity. Gene flow to these sites may be provided from small fragmented eelgrass beds that are still found in the Marstrand area (currently under investigation). Each of the three clusters is further characterized by strong population genetic structure among all sampling sites and few first-generation migrants (Table S5). This is typical for seagrasses and caused by one or more of the following factors: partial clonality (Olsen et al., 2004), a larger role for mutations over migration due to the longevity of clones , sporadic recruitment (Becheler, Diekmann, Hily, Moalic, & Arnaud-Haond, 2010), founder-takes-all recruitment (Waters, Fraser, & Hewitt, 2013) or stochasticity of dispersal (Kendrick et al., 2012).
Dispersal among populations was further explored with the biophysical particle modelling. Single-generation dispersal explains the differentiation of the small genetic cluster in the Marstrand area (red cluster in Figure 3b) and is significantly correlated with asymmetric migration rates-but not with genetic differentiation (D ps ; Table 2).
Multigeneration dispersal explains a higher proportion of genetic differentiation and asymmetric migration rates (Table 2), and longdistance connectivity increases (Figure 3c). Inclusion of the historic distribution in the multigeneration model results in an almost perfect recovery of the genetic clusters (Figure 3a,d), and both D ps and asymmetric migration rates have an improved fit with this measure (Table 2).
Overall, this assessment with presumably neutral genetic markers suggests that the processes of migration and genetic drift explain a large T A B L E 2 Results of Mantel tests between log 10 -transformed "sea distance" or dispersal probabilities (see Material and Methods) and a genetic differentiation matrix based on the proportion of shared alleles (D ps ) or asymmetric migration rates based on G ST and calculated with DivMigrate (asymm. mig.)

"sea distance"
Single-generation dispersal probability

Multigeneration dispersal probability
Historic multigeneration dispersal probability All but the correlation between single-generation dispersal probability and D ps is significant (bold). Note that a negative correlation is expected between D ps and minimum dispersal probability, because sites with a high probability of dispersal between them are expected to show low genetic differentiation, whereas asymmetric migration rates are expected to be positively correlated with dispersal probability.
part of the observed genetic population structure, but adaptation to local physical parameters could be an additional explanation.
One major genetic and oceanographic break observed is located approximately at the border of the Kattegat and Skagerrak and is confirmed by the few previous studies of population genetic structure in the Skagerrak and Kattegat, e.g., for herring (Lamichhaney et al., 2012), harbour porpoise (Lah et al., 2016) and cod (Barth et al., 2017).
A particularly relevant study of the macroalga Saccharina latissimi, which also shows exclusively passive dispersal, indicates a similar genetic break between the Kattegat and Skagerrak (Moller Nielsen et al., 2016). The other important genetic and oceanographic break we observe is located in the Marstrand area. This has not been previously reported-probably due to lack of geographically detailed sampling for genetic studies in the area. The only studies in the Skagerrak and Kattegat that used both genetic and biophysical methods found high correlations between gene flow and oceanographic connectivity for diatoms (Godhe et al., 2013), while the correlation was lower for actively moving cod (Barth et al., 2017).
The network analysis also allows visualization and identification of populations central to connectivity (Greenbaum & Fefferman, 2017;Rozenfeld et al., 2008). The oceanographic and genetic network analyses indicate that meadows from south-western and central Kattegat

| Comparison of connectivity measures and temporal scales
The best fit among Mantel tests was obtained between genetic distance based on the proportion of shared alleles (D ps ) and multigeneration-historic dispersal, explaining ~40% of genetic variability (Table 2). This metric takes into account oceanographic dispersal distance, stepping-stone dispersal over multiple generations and historic habitat continuity. IBO that incorporates oceanographic distance and habitat discontinuity, and/or stepping-stone dispersal, was able to achieve similarly high correlations for giant kelps and fucoid macroalgae that also disperse by rafting (Alberto et al., 2011;Buonomo et al., 2017). Thus, support is mounting that IBO is a better approach than IBD to explain genetic structure and gene flow. However, the correlation between genetic differentiation and single-generation oceanographic dispersal performs worse than classical IBD. Asymmetric migration, calculated with DivMigrate, correlates more strongly than D ps with both single-and multigeneration-extant dispersal probability, but not with multigeneration-historic dispersal. This indicates that this measure of asymmetric migration is capable of capturing more recent migration rates (Sundqvist et al., 2016). This measure is relatively new, and to our knowledge, we test and show here for the first time that this metric has indeed a better fit with dispersal probability on shallow time-scales, as would be expected (Sundqvist et al., 2016). The good fit of both genetic measures and oceanographic connectivity further indicates that genetic structure of eelgrass in the Skagerrak-Kattegat is mainly driven by migration and genetic drift-and not selection.
Despite the documented high loss of eelgrass meadows in the area, the effect is not visible in the levels of genetic diversity and differentiation. Thanks to the availability of historic distribution data, we were able to compare modelled connectivity of the extant eelgrass beds with the historic distribution of ~100 years ago. As would be expected, the loss of the historic meadows has resulted in some changes in the probability of dispersal and the network structure. For example, Limfjord (6-LM) and Gåso (G-SG) have become oceanographically isolated over the last century (Figure 3c,d), but this is not visible (yet) in the genetic structure (Figure 3a). This genetic memory or "ghost of dispersal past" (Benzie, 1999)

| Complementary value of genetic and biophysical models
Cross-validations of the genetic and oceanographic modelling data show good agreement and provide different insights into the structure and connectivity of populations. The genetic survey integrates over many gene-flow mechanisms and captures regional population history through deep time. In addition, diversity metrics and population differentiation, as well as inferences about demography, can only be determined with genetic data. In general, genetic methods are less well suited for inferring the spatial component and directionality of dispersal.
In contrast, biophysical models provide insights about the generational time depth of dispersal and the shaping of populations with respect to barriers and circulation patterns. Although single-generation dispersal may be a weak predictor, the ability to simulate a range of generational time depths through stepping-stone simulations is a distinct advantage. When historical distribution records are also available, as is the case here, predictions of where populations should or could persist become very powerful. Biophysical models also offer better spatial coverage than is feasible with most genetic sampling efforts. The main shortcoming of biophysical models is that they say nothing about demographic history, adaptive potential or genetic health of the species in question. In terms of resource investment, initial front-end development of suitable oceanographic models is both time-consuming and cost-intensive, and limited to the specific region of interest. In the absence of such oceanographic models, genetic surveys remain the best alternative. Microsatellite markers are available for many seagrass species, and an assessment such as the one reported here is standardized and easy to perform.
The added value of the dual approach further strengthens conservation planning and eventual monitoring of a particular management plan, because it is possible to rerun a biophysical model with different data and under different scenarios to reflect adaptive management (e.g., McCook et al., 2010). The results of the two approaches to connectivity may also be used to rank sites, for instance according to their connectedness, whether they act as sources or sinks or their level of diversity (Jonsson, Nilsson Jacobi et al., 2016). Such results could then be used by conservation managers in spatial planning programmes such as Marxan and Zonation for prioritization in large-scale conservation efforts (Delavenne et al., 2012).

| Implications for management
Understanding spatial population structure and identifying areas with restricted as well as excellent connectivity are essential in conservation management. Here, both the genetic and hydrodynamic connectivity assessments identified dispersal barriers, creating distinct clusters that could serve as management units (Palsbøll, Bérubé, & Allendorf, 2007), which should be managed separately to ensure longterm persistence and protection of genetic diversity (Allendorf et al., 2013). While a genetic and oceanographic barrier is evident between the Kattegat and Skagerrak, no such break is visible between Denmark and Sweden. Thus, it is important to assess whether existing or proposed MPAs in the study area constitute functional networks within each of the three genetically distinct clusters and their biophysical barriers also across countries. For instance, all Z. marina meadows within the large cluster covering most of the Kattegat could be managed as a single unit within, for example, an MPA network, as replenishment from one site to another can be expected. However, the two sites that have become oceanographically isolated (6-LM and G-SG) will require further local management. Likewise, it is critical to protect the small meadows off the Island Laesø (7-LS and 7-480), remnants of the large historic offshore population that appear key for connectivity ( Figure 3). Fortunately, the meadows off Laesø are presently included in a ca. 1,000 km 2 large Natura 2000 site that includes protection of shallow water soft-sediment habitats (Moksnes et al., 2014). However, dramatic improvements of the environmental conditions in, for example, water clarity would be required for recovery towards the historic distribution in offshore Kattegat where the maximum depth distribution of eelgrass has decreased by >50% (Boström et al., 2003). It is also important to ensure that MPAs provide the intended protection to habitats and biodiversity (Almany et al., 2009). Along the Swedish coasts, for example, small-scale destruction of eelgrass meadows for construction of piers and marinas is high and occurs even inside protected areas (69% of the studied cases within eelgrass meadows were approved for construction ;Eriander, Laas, Bergström, Gipperth, & Moksnes, 2017). Hence, there is an urgent need to review regulations and management of existing MPAs.
From a local management perspective, one of the key findings is the low connectivity into and out of the Marstrand area (red cluster 1) in the Swedish Kattegat, where major losses of eelgrass have occurred and continue to occur. The oceanographic isolation could be related to two large rivers that enter the Kattegat just south of this area and may create a dispersal barrier. The genetic isolation indicates that natural replenishment from outside is unlikely, making protection of the remaining Z. marina beds crucial. The high allelic diversity suggests that the losses have not yet negatively affected fitness and that local meadows constitute good donor material for restoration (e.g., Reynolds et al., 2012), while transplantation from other areas should be avoided (e.g., Kettenring, Mercer, Reinhardt Adams, & Hines, 2014).
In fact, this study shows that genetic diversity and connectivity of Z. marina in the Skagerrak-Kattegat seem to be generally in a healthy state, but assessment such as this, in addition to assessments on a more local scale, can highlight vulnerable sites or could be used as baselines for tracking future changes.
The historic loss in the Skagerrak-Kattegat is possibly the largest reported seagrass loss in the world (P.-O. Moksnes, unpubl.). Despite protection by international conventions and directives, and large MPA networks (covering approximately 15% of the Skagerrak-Kattegat; Moksnes et al., 2014), recovery has been very limited, and losses continue (Boström et al., 2014;Moksnes et al., 2016). Thus, protection from physical impacts within MPAs may not be sufficient to ensure persistence and recovery of Z. marina, but more and new measures are needed to improve environmental conditions. In addition to increased efforts to reduce nutrient input to coastal waters, management should consider measures to enhance depleted populations of large predatory fish that would restore the trophic structure of coastal ecosystems (Östman et al., 2016), measures that can break self-generating feedback mechanisms such as sediment resuspension that lock the system in a turbid state (Maxwell et al., 2016;Nyström et al., 2012), and eelgrass restoration to facilitate a natural recovery of lost meadows (van Katwijk et al., 2016), including compensatory mitigation of eelgrass lost or damaged during, for example, coastal exploitation .

| CONCLUSIONS
Our analysis supports the notion that passive rafting of flowering shoots by oceanographic currents influence patterns of gene flow of Z. marina in the Skagerrak-Kattegat and is the main driver of observed population genetic structure and meta-population dynamics.
We show that the meta-population is driven by stepping-stone dispersal over many generations and that current genetic differentiation is best explained by connectivity considering the historic Z. marina distribution. This "ghost of dispersal past" is also evident in the distribution of allelic richness, where highest diversity is found in the Laesø Island area, where major historic losses occurred. Using two complementary methods to assess connectivity enabled us to investigate and compare dispersal and migration patterns at different temporal scales. In this study, we found strong concordance among the two methods in detecting sources, sinks and connectivity patterns. This information can be used to pinpoint areas where local protection is necessary or where populations could/should be managed in a network approach for MPAs. The temporally more dynamic oceanographic modelling was also able to highlight areas where connectivity has become limited over the last decades. Such information is additionally helpful for marine spatial management to pinpoint geographic areas where there is a need to improve environmental conditions. The large geographic scale study presented here forms a framework for future detailed assessments of connectivity and genetic diversity on smaller scales within the coastal archipelagos and fjords. Such multiscale information should aid managers at the local, national and international levels in marine spatial planning, for example, for the identification of hubs, and important extant or historic source meadows that should be targeted for protection, or key areas for restoration.