Integrating genetics, biophysical, and demographic insights identifies critical sites for seagrass conservation

,


INTRODUCTION
Eelgrass (Zostera marina) forms an important habitat along temperate to arctic coastlines throughout the northern hemisphere, but huge declines have been observed across the entire distribution over the past century , Bostr€ om et al. 2014. Eelgrass wasting disease in the 1930s caused the disappearance of 90% of eelgrass along the Atlantic coasts of North America and Europe (Muehlstein 1989). Although eelgrass recovered in many shallow areas in the 1960s, it never obtained its former distribution, particular in deeper areas (e.g., Bostr€ om et al. 2014). With increasing nutrient pollution and other anthropogenic threats, large losses of eelgrass occurred again in the 1970s and 1980s and continue today in many areas , Bostr€ om et al. 2014.
As foundation species, seagrasses provide critical habitat to a large variety of species (Bertelli andUnsworth 2014, Ellison 2019), forming a level of productivity on par with coral reefs and rain forests. Eelgrass is an annual or perennial monoecious (male and female flowers on same individual) flowering plant (Ackerman 2006). It is a facultatively sexual reproducer with underwater flowers, pollen, and seeds (Ackerman 2006). Eelgrass is outcrossing, but also self-compatible resulting in self-fertilization, but little is known about the relative contributions over time and space (Reusch 2000). Within meadows, clonal extension also plays a major role (Reusch et al. 1999). Along the Swedish Skagerrak coast, the historic decline of eelgrass has resulted in significant losses of ecosystem services, e.g., reductions of cod catches and lower capacity for carbon and nutrient sequestration (Cole and Moksnes 2016). Loss of meadows and poor connectivity can lead to a reduction of genetic diversity and fitness leading to an extinction vortex (Procaccini et al. 2007, Allendorf 2017, where ever declining population sizes and genetic diversity are predicted to reinforce decline, leading to a decrease in evolutionary potential (Leimu et al. 2006), and finally extinction (Benson et al. 2016). The resulting habitat loss also affects connectivity and genetic diversity of associated species, potentially also trapping associated communities in an extinction vortex (Hughes et al. 2009, Sorte et al. 2017.
Improved management of seagrass is increasingly urgent in order to achieve biodiversity targets (Ellison 2019). One important strategy in conservation management for achieving this goal is to design connected landscapes or seascapes (Foley et al. 2010), and such studies are an important way to identify vulnerable areas as well as particularly valuable sites for connectivity and persistence of metapopulations. Assessments of connectivity on relevant spatial scales are an important part of seascape studies, as dispersal is fundamental to ecology and evolutionary dynamics of metapopulations. Understanding dispersal in the marine environment is difficult due to the high spatial and temporal variability of ocean currents (Cowen and Sponaugle 2009), but mutually enriching genetic and biophysical approaches have been used successfully to assess patterns of seagrass dispersal on small (Sinclair et al. 2018) and large geographic scales (Hernawan et al. 2016, Jahnke et al. 2016, 2017, Triest et al. 2018. Landscape/seascape studies become particularly valuable when different approaches to assess connectivity and identify particularly vulnerable sites are combined. Important seascape approaches applicable to marine spatial management include biophysical modeling (Jonsson et al. 2016, Grech et al. 2018, seascape genetics (Riginos andLiggins 2013, Jahnke et al. 2018) and demographic modeling (Padr on and Guizien 2015, Matz et al. 2018). In particular, a combined use of these different approaches is able to give information on the past (genetics), the present (biophysical modeling) and the future (demographic modeling). The practical value of data from seascape studies is their use to manage marine habitats and to model their persistence (Bostr€ om et al. 2011, Pittman 2017. Although landscape studies are commonly employed in terrestrial ecosystem management and conservation (Turner 1989), they are still rarely applied to coastal seascapes (Li and Mander 2009), including seagrass meadows (Bostr€ om et al. 2011). Worldwide, governments are requesting research to help plan regional-scale and local-scale habitat connectivity (Brodie et al. 2016, Albert et al. 2017. There is thus a need for a general framework to support the development of tools for managers and stakeholders who, ultimately, are the users of information on seagrass seascape ecology. Along the Swedish northwest coast, over 60% of eelgrass meadows have been lost since the 1980s (Baden et al. 2003). The losses have largely been attributed to the effects of coastal eutrophication in combination with overfishing of large predatory fish, causing a trophic cascade and an increase in ephemeral macroalgae that cover the eelgrass meadows (Moksnes et al. 2008, Baden et al. 2012. Despite successful efforts to reduce nutrient loading over the past 20 yr, no recovery of eelgrass coverage has been observed (Nyqvist et al. 2009). Instead, losses of eelgrass continue, in particular in the Marstrand fjord region in the southern part of the Swedish northwest coast where less than 8% of the meadows remain today compared to the 1980s .
Here, the losses of eelgrass have destabilized the bottom and increased sediment resuspension and turbidity, resulting in local regime shifts that appear to be spreading to neighboring areas . Substantial cumulative losses also occur due to small-scale coastal exploitation for recreational boating (Eriander et al. 2017). To meet these challenges, Swedish authorities have recently initiated a national action plan for eelgrass involving both increased spatial protection and large-scale restoration of eelgrass habitats . However, these measures are presently constrained by a lack of understanding of connectivity and genetic diversity, useful for prioritizing meadows for protection and restoration. Information is also lacking regarding how the extensive losses in the Marstrand region have affected the genetic diversity of eelgrass in the area, and their ability to recover. Here, we focused on these well-studied eelgrass meadows along the Swedish west coast to address the challenge of efficient conservation at a spatial scale relevant to local management, by developing an integrative framework based on (1) genetic, (2) biophysical, and (3) demographic modeling methods. We apply the workflow and the toolbox we present here to two fjord regions along the Swedish West coast, resulting in specific management recommendations for the area, and discuss how the approach may be developed into a management tool applicable in other areas and for other species.

Study area and sampling
The study area included two large fjord regions along the Swedish northwest coast in the eastern part of the North Sea: the Marstrand fjord region (57.8-58.1°N and 11.7-11.8°E), which encompasses the Hakefjord, € Alg€ ofjord, and S€ al€ ofjord and has suffered extensive losses of eelgrass; and the Gullmarsfjord region (58.2-58.5°N and 11.3-11.8°E), which comprises Abyfjord, Brofjord, and Gullmarsfjord and has experienced little decline (~5%) since the 1980s (Baden et al. 2003; Fig. 1). For both areas, eelgrass distribution data are available from the early 1980s and 2000s (Baden et al. 2003, Nyqvist et al. 2009). These reference points were used to identify extant and historic meadows for the biophysical modeling ( Fig. 1), and to select meadows for genetic sampling.

Genetic sampling and analyses
Genetic diversity, population structure and connectivity were analyzed at 22 locations using 20 microsatellite loci and a suite of analytical approaches. Eleven sites in the Gullmarsfjord region and 11 sites in the Marstrand region were sampled (Figs. 1 and 2). At each site, 40 shoots of Z. marina were collected (with the exception of G6, for which N = 38). Shoots were collected haphazardly by snorkeling or diving, at 1.0-1.5 m intervals and depths ranging from 1.2 to 5.3 m. Sampling was carried out during July and August 2016.

DNA extraction and microsatellite amplification
DNA was extracted from 5 to 10 mg of silica-gel-dried leaf tissues and 22 microsatellite loci were multiplexed for PCR amplification. See Jahnke et al. (2018) for further details, including for clone identification and data quality checks. Some of the genetic distance metrics used for downstream analyses assume that microsatellite loci are following a stepwise mutation model (SMM), where new alleles arising by mutation have one repeat unit difference to the original allele (Ohta and Kimura 1973). To test this assumption, a likelihood approach implemented in MISAT (Nielsen 1997) was used and analyses were performed for each locus in the two fjord regions separately.

Genetic diversity measures
We calculated common genetic diversity measures using GenAlEx 6.5 (Peakall and Smouse 2012). Genotypic richness (R), a proxy of the amount of clonal growth, was calculated with the formula (MLG À 1)/ (N -1) (Dorken and Eckert 2001), which relates the number of multilocus genotypes (MLGs) to the number of ramets (N) sampled. For example, 40 ramets were collected at a site, but genotyping revealed only three unique genotypes. This would indicate extensive clonality, whereas if genotyping revealed 40 unique genotypes, then sexual reproduction was dominant. Allelic richness (Ar), a measure of the number of alleles per locus accounting for sample size, was calculated with standArich in R 2.15.3 (Alberto et al. 2006). Standardization was made for a sample size of 13 MLGs, the lowest number of different MLGs found at a site.

Assessment and modeling of reproductive parameters
As genetic diversity is slow to react to changes and might not give a realistic picture of current and recent population dynamics, we also investigated demographic parameters more directly, by using Genetic Diversity Spectrum (GDS; Rozenfeld et al. 2007) and modeling approaches to investigate rates of outcrossing, clonality, and self-fertilization in each fjord region. This approach relies on genetic distance measures, and we used the results as a base for biological assumptions in the demographic modeling with SLiM. See Appendix S1 for a detailed description of methods and scripts are provided on Zenodo (see Data Availability).

Genetic differentiation and isolation by distance
To investigate gene flow within and among the two fjord regions, we looked at population differentiation. To assess population structure, we used a Bayesian, model-based clustering method implemented in TESS 2.3 (Chen et al. 2007). See Jahnke et al. (2018) for further details on the TESS analysis. To assess Isolation by Distance (IBD) within each fjord region we used a GDS-based approach (Arnaud-Haond et al. 2014) and Mantel tests (Mantel 1967). Sea-distance, distances between sampling sites without crossing land, was calculated using the package marmap in R 3.3.0 (Pante and Simon-Bouhet 2013). Custom-made scripts for the GDS approach are provided on Zenodo (see Data Availability). Mantel tests were carried out using the package ncf in R 3.3.0 (Bjornstad 2009) with 100,000 permutations, and the pairwise genetic differentiation indices (F ST and D) were calculated with the diveRsity package (Keenan et al. 2013).

Genetic network analyses
Network analyses provide a way to assess network hubs without relying on population genetic assumptions (Greenbaum et al. 2016). Instead, networks take advantage of all information embedded in the data set to let the data define their own topology (Rozenfeld et al. 2007). Network properties can be quantitatively estimated and they can reveal interesting properties about the connectivity patterns in the studied populations. Networks were built at the percolation threshold, which corresponds to the point below which the largest part of the network becomes fragmented (Stauffer and Aharony 1991), in order to keep only the essential links, yet maintaining the connectivity between the largest components of the network. As networks are dependent on the genetic distance metric used to build them, we used several genetic distance measures: Goldstein, D PS , Reynold's D, and cGD. Goldstein and Reynold's D pairwise estimates were calculated using EDENetworks (Kivel€ a et al. 2015), D PS with MSA 4.05 (Dieringer and Schl€ otterer 2003) and cGD with the popgraph package in R 3.3.0 (Dyer 2014). For all genetic distances, the percolation threshold was assessed using EDENetworks (Kivel€ a et al. 2015), and the pairwise genetic distance matrices were imported in R 3.0.0 and plotted using the  Table 1. Genetic clusters (green, blue, and red) of the background coloration show the spatial interpolation of ancestry coefficients (Q values or proportion of individuals belonging to each cluster) based on the TESS analysis with (a) K max = 4 and (b) K max = 2; the gradient within each color indicates percentage of group membership belonging to the genetic clusters. The overlaid colored dots (yellow, green, orange, and light blue) represent release points of particles in the oceanographic modeling. The different colors indicate the different oceanographic clusters identified by a clustering method based on modeled dispersal probabilities. Dots with the same color indicate areas that have an internal connectivity above the dispersal restriction, and the transitions of colors thus indicate partial dispersal barriers. Both areas are separated into three oceanographic clusters. packages igraph and ggmap (Csardi andNepusz 2006, Kahle andWickham 2013).
To investigate network properties further, we used a range of measures to define the networks in each fjord region (see Appendix S1: Table S1). At an element (node) scale, we used betweenness centrality to identify meadows that serve as stepping stones. Highest values of betweenness centrality correspond to network hubs, which in the case of genetic networks, can be understood as those nodes acting as a source or sink, relaying and/or receiving gene flow from other nodes (Rozenfeld et al. 2008, Becheler et al. 2010. We compared networks of each fjord with a random network using the Erd} os-R enyi model in the Rpackage igraph to assess their network type. A total of 1,000 random networks were created using the same number of nodes and edges as for the fjord networks. These network property analyses were carried out with the R packages ggplot2 and ggmap (Kahle andWickham 2013, Wickham 2016) and scripts can be found on Zenodo (see Data Availability).

Biophysical model of drifting reproductive shoots
The local connectivity of eelgrass meadows in the Gullmarsfjord and the Marstrand fjord regions (Fig. 1) was estimated with a biophysical model of shoot dispersal. Two high-resolution, three-dimensional hydrodynamic models were set up in each fjord regions. The commercial software MIKE 3 FM, with an unstructured computational grid, was used to model the local hydrodynamics (DHI 2017,software available online. www.mikepoweredbydhi.com). 2 The grid-cell size ranged from about 100 m close to land contours and up to 1,000 m offshore (see Appendix S1: Fig. S1). The vertical resolution ranged from 1 m at the surface to 10 m at the deepest bottom layers. Water movement is constrained by seabed bathymetry (see Appendix S1: Fig. S1) and forced by gridded atmospheric data (wind, air temperature, and pressure), freshwater input based on measured data together with modeled water velocities, sea-surface height, and water density (temperature and salinity) at the model boundaries, produced by a large-scale regional model (MIKE3 FM, Baltic Sea and North Sea model). Water transport was modeled during August and September and repeated for 5 yr (2011)(2012)(2013)(2014)(2015), a period spanning the range of annual North Atlantic Oscillation Index, known to influence large-scale circulation patterns in the eastern North Sea (Hurrell and Deser 2010).
The dispersal of eelgrass shoots was simulated with a Lagrangian Agent Based Model (ABM) implemented in the MIKE ECO Lab module (Kuusem€ ae et al. 2018). Particles representing drifting reproductive shoots with seeds were released from a large number of sites (meadows) with present or historic occurrence of eelgrass ( Fig. 1). Release of particles occurred every 5 d during August (seven occasions) and the position of particles was recorded each day for a total drift time of 30 d. Eelgrass reproductive shoots are positively buoyant and particles in the model drifted in the surface. Period of flowering and detachment, as well as duration that shoots stay afloat were based on empirical field studies along the Swedish Skagerrak coast (K€ allstr€ om et al. 2008, Infantes and. Transport of objects floating close to the water surface is influenced both by advection from ocean currents and by wind forcing of the near-surface layer (windage). To model this, we added 3% of the wind velocity to the transport of simulated particles (Wu 1983). We assumed that the negatively buoyant seeds may be dropped from the floating spathes at any time during the 30 d of drift time, and the position in the seascape of each simulated spathe was recorded once every day representing positions of potentially dropped seeds.
Existing occurrence data of eelgrass from field surveys were digitized as polygons in a GIS software (QGIS 3.4; Quantum GIS Development Team 2014), where each polygon was considered an eelgrass meadow (Fig. 1). Occurrence of eelgrass representing present and historic conditions were collected in 2015 and 1980, respectively. The simulations of the ABM included 140 eelgrass meadows in the Gullmarsfjord region, which were assumed to have been relatively stable for the past four decades. For the Marstrand region 237 meadows were included to represent historic conditions around 1980, and 126 meadows were included in the model to represent what remains today (Baden et al. 2003, Nyqvist et al. 2009Fig. 1). Each meadow was seeded in a grid every 100 m. For each release time 130,000 and 308,400 particles were released in the Gullmarsfjord and Marstrand fjord regions, respectively. In total, across all meadows, release times, and years, the dispersal of 1.5 9 10 7 particles were simulated in the ABM resulting in 4.5 9 10 8 positions along shoot dispersal trajectories. From initial release positions and recorded dispersal positions, a connectivity matrix was calculated for each day, release time and year. The connectivity matrices, normalized to the number of released particles, represent the probability of dispersal from meadow i to meadow j. For the eigenvalue perturbation and barrier analyses the connectivity matrices were averaged to obtain one mean matrix for each of the three sets of meadows: Gullmarsfjord, and the Marstrand present and historic data sets.

SUBPOPULATIONS
We employed a clustering method to identify subpopulations and partial dispersal barriers based on modeled dispersal probabilities in the seascape (Nilsson Jacobi et al. 2012). Only dispersal between eelgrass meadows (present or historic) was considered. This theoretical framework analyzes the seascape connectivity matrix to find clusters as a signature of partially isolated subpopulations. Identification of subpopulations is formulated as a minimization problem with a tunable penalty term for merging clusters that makes it possible to generate population subdivisions with varying degrees of dispersal restrictions that are typical for both genetic divergence and demographic independence (here, intra-cluster connectivity was about 100 times higher than inter-cluster connectivity). The results can therefore be used to identify spatial units relevant for conservation and management of marine populations, as an alternative or complement to genetic approaches (Nilsson Jacobi et al. 2012).

Ranking of meadows using eigenvalue perturbation theory (EPT)
Eelgrass meadows, and locations with known lost meadows, may be viewed as a network with exchange of seeds where each meadow represents a node. We applied eigenvalue perturbation theory (EPT) to rank the importance of individual meadows for the connectivity and growth of the whole network of meadows Hanski 2003, Nilsson Jacobi andJonsson 2011). The method ranks each meadow according to their strength as both donor and recipients of propagules, which is critical for maximizing the size of the metapopulation, particularly in small and threatened populations (Nilsson Jacobi and Jonsson 2011). This may be viewed either demographically (as the contribution of each individual meadow), through dispersal (to the overall growth rate of the metapopulation), or evolutionarily (as the contribution of individual meadows to the spread of genes; Nilsson Jacobi and Jonsson 2011). The EPT analysis was applied to the three mean connectivity matrices and each meadow received a score indicating its dynamic network importance. In the Marstrand fjord region, we used the EPT procedure also to identify and rank the most important lost meadows for connectivity and persistence of the present metapopulation.

Demographic modeling based on genetic and biophysical data
The availability of both genetic data and dispersal probabilities from biophysical modeling provides a unique opportunity to explore demographic and genetic stability of the assessed meadows in the future. In a first step, we modeled the past population dynamics to see if we could obtain genetic diversity and allele frequency spectra (AFS) similar to our observed data, and then extended this model into the future using the empirical genetic and biophysical data.

Demographic modeling of the past
We used individual forward simulations in SLiM (Haller and Messer 2019) to identify a demographic scenario that could explain our observed genetic data. We used a non-Wright-Fisher model and simulated a neutral 2 9 10 8 genome as one chromosome (i.e., the approximate size of the Zostera marina genome; Olsen et al. 2016) assuming no background mutations, a microsatellite mutation rate of 0.002 per allele/per generation and a recombination rate of 1 9 10 À8 (Stapley et al. 2017). Guided by the unique AFS of the real populations with a high proportion of high-and low-frequency alleles (see Fig. 4), the focal population was simulated for 20,000 generations with a carrying capacity of 100 individuals. To introduce the high number of low-frequency alleles into the focal population, we modeled a second population with a carrying capacity of 1,000 individuals, which supplied migrants to the focal population. Clonality was allowed to fluctuate in the new recruits of each generation around 0.3, guided by the genotypic richness observed in the two fjord areas (Table 1). Clonal lineages were tracked, and only individuals produced by clonal reproduction were allowed to age beyond 10 yr. No self-fertilization was implemented, because of the low levels inferred in the GDS analysis. Age structure was inferred from the real data, by multiplying clone extent with a vertical extension rate of 0.15 m/yr (Olesen et al. 2017). Mortality rates per age class were set to produce an age distribution similar to what we observed in our data, mimicking long overlapping generations created by clonality of up to 100 yr (equivalent to the largest extent of 13.5 m of a clone found in this study). Microsatellite evolution was modeled as in Haller and Messer (2016: Modelling microsatellites). After 20,000 generation genetic diversity measures were stable. At this point, the two populations experienced a bottleneck of 90% reflecting metapopulation size reduction during colonization of the Kattegat and Skagerrak by marine species after the last glacial maximum (LGM,~8000 yr ago). Now migration at the level of 0.001 was set between the two populations. Populations were allowed to recover with a growth scaling factor of 1.1 until they reached their original carrying capacity. After a further 7,900 generations, the populations experienced another bottleneck of 75%, the minimum decline reported after the wasting disease during the 1930s in Denmark (Rasmussen 1973) and at the scale observed for the recent decline in the Marstrand fjord region since the 1980s. Custom-made bash and R scripts were used to analyze the AFS and genetic diversity measures of the modeled data and compare it to the observed data. Several other models of a single population, three populations and different numbers and degrees of bottlenecks were explored, but produced AFS that fitted the observed data less well. Each model was run with 10 successful replicates (i.e., runs where no population went extinct during or after a bottleneck), and outputs were examined over the entire model run, i.e., for 28,000 generations. After confirming that replicates gave similar estimates in the first 27,000 generations, we plotted the run with the best fit to the observed data at generation 28,000 (SLiM, bash, and R scripts can be found on Zenodo; see Data Availability).

Demographic modeling of the future
To explore the stability of populations in each fjord region in the future, we extended the model beyond 28,000 generations and included dispersal probabilities from the biophysical modeling, effective population sizes (N e ) calculated from the genetic data and decline rates recorded during eelgrass monitoring for each meadow as input parameters. At generation 28,001, we split the focal population into the 11 sampling sites in each fjord region. We allowed continued migration from the outside population at 0.001 and set migration rates between the 11 sampling sites as estimated from the biophysical model (see Appendix S1: Tables S2 and S3). Population sizes were set according to the effective population size (N e ) calculated with the ldNe function in R package strataG (Archer et al. 2017). Note that we do not necessarily believe these N e values to be "true", particularly as we are dealing with a partially clonal species, but we consider that N e ratios among sites are meaningful. We modeled the Gullmarsfjord with stable population sizes (N e was set as following: G1, 63; G2, 19; G3, 11; G4, 39; G5, 395; G6, 7; G7, 3; G8, 19; G9, 12; G10, 7; G11, 13), and modeled the Marstrand meadows with constant decline rates as recorded for each site over the last four decades (N e and decline rates were set as following: M1, 500, 5% per yr; M2, 500, 4% per yr; M3, 10, 0% per yr; M4, 10, 4% per yr; M5, 5, 0% per yr; M6, 64, 0% per yr; M7, 1, 11% per yr; M8: 17, 0% per yr; M9, 116, 11% per yr; M10, 31, 11% per yr; M11, 45, 11% per year; SLiM, bash, and R scripts can be found on Zenodo; see Data Availability).

Toolbox combining the different approaches
We employed a number of different approaches including genetic analyses, biophysical modeling comprising an EPT ranking and demographic modeling. Importantly, these different tools are able to look at different time scales (past, current, and future) and can be used complementary to each other for many questions. After confirming that genetic and oceanographic cluster analyses suggested very similar subpopulations, we divided the 266 extant meadows into the seven suggested management units (MUs) in the two fjord regions. In both regions we ranked all extant eelgrass meadows within each MU according to their EPT scoring. In the Marstrand region we also ranked all historic meadows according to their EPT scoring for the historic as well as the present metapopulations. For the 22 sites for which genetic data was available, we moreover gave site-specific values or MU-wide scores for past metapopulation dynamics, past connectivity, and the future potential. For the past metapopulation dynamics, we report allelic and genotypic richness per site or per MU. For the past connectivity, we report scoring based on the inbetweenness centrality values from the genetic network analysis and scored the entire MU as high, moderate, or low. Investigating the future potential using demographic modeling in the stable Gullmarsfjord region, we used a scoring based on predicted allelic richness values 90 generations into the future, and scored MUs into either stable or moderate. For the impacted Marstrand fjord region, we scored the 11 assessed sites according to the predicted time until extinction (highest score for sites that go extinct last) and time until recolonization after extinction from other modeled sites (highest score for sites that become recolonized earlier). MUs were scored as declining (time to extinction) or moderate (time to recolonization).

Genetic analyses
Clonality and data checks.-To get a clean genetic data set, we carried out several data checks to assess clonality and to decide which genetic markers to include in the analysis. We removed 441 individuals representing clonal replicates within sites, resulting in 1,195 genets from the 1,636 ramets sampled. The number of MLGs per sampling site ranged from 13 to 38 (Table 1). The software Microdrop revealed an allelic dropout rate of 11% in locus D2, which was therefore removed for all downstream analyses. All other loci had null allele rate estimates below 1%. Deviations from Hardy-Weinberg proportions on the remaining 21 loci per population were significant in 5.1% after correcting for multiple comparisons. Significant linkage disequilibrium was found in 0.6% of comparisons. Locus GA35 was identified as an outlier by the neutrality tests in both LOSI-TAN and BayeScan for both fjord regions (see Appendix S1: Fig. S2). It appeared as being under balancing selection due to having lower than expected genetic differentiation among populations. However, this locus has the highest mean number of alleles (see Appendix S1: Fig. S3), and has consequently lower values of F ST due to a downward bias of loci with high polymorphism. The reliance of both LOSITAN and BayeScan on F ST might cause GA35 to be identified under balancing selection only because it is highly polymorphic. As the exclusion of locus GA35 did not have a considerable effect on PCAs (not shown), it was left in the data set for consecutive analyses. The likelihood ratio test from MISAT rejected the null hypothesis (SMM) in locus G4 for both Gullmarsfjord and Marstrand regions (P values were 0.001 and 3.7 9 10 À7 , respectively) after Bonferroni corrections. This locus was removed from the data set for all further analyses.
Twenty microsatellite loci were used for all subsequent analyses.
Genetic diversity.-Allelic (basically reflecting population size and age in a biogeographic context) and genotypic richness (basically reflecting sexual reproduction vs. vegetative spread) were high and similar in the two fjord regions (Table 1; Appendix S1: Fig. S4). Genotypic richness ranged from 0.31 in meadow M6 to 0.95 in two meadows in the Gullmarsfjord (Table 1; Appendix S1: Fig. S4). Allelic richness standardized to 13 MLGs was 4.30 AE 0.32 (mean AE SD) in the Gullmarsfjord and 4.26 AE 0.22 in the Marstrand region.
Assessments of outcrossing, self-fertilization and clonality.-GDS analyses and modeling indicated that clonality, outcrossing, and self-fertilization rates are similar in the two fjord systems (see Appendix S1: Figs. S5-S8). The outcome of iterating a specific reproduction mode over many generations was consistent among the sampling sites and when simulating all reproductive modes together, the GDS mostly overlapped with the "real" data or deviated towards the left (see Appendix S1: Figs. S7 and S8). In all simulated scenarios, there was a gap between zero and larger genetic distances, as seen in the real data (see Appendix S1: Figs. S5, S7 and S8). Therefore, both fjord regions show high levels of outcrossing and vegetative propagation, but little self-fertilization.
Although not directly relevant for eelgrass management or management of other species with simpler reproduction strategies, this knowledge is an important contribution to understanding the biology of eelgrass in the region, and allows more realistic demographic modeling.
Genetic differentiation and isolation by distance.-The TESS analysis of all sites (n = 22), showed differentiation between and within fjord regions (see Appendix S1: Fig. S9). Looking at each fjord region separately, the analysis identified four clusters in the Gullmarsfjord region, and two clusters in the Marstrand region ( Fig. 2; Appendix S1: Fig. S9). In the Gullmarsfjord region, each of the three fjords ( Aby-, Bro-and Gullmarsfjord) formed separate clusters. In addition, the two meadows sampled in the inner part of the Gullmarsfjord (G6 and G7) formed a fourth cluster. In the Marstrand fjord region, the sampled meadows in the Hakefjord and € Alg€ ofjord formed one cluster, whereas the two meadows in the S€ al€ ofjord formed a separate subpopulation (Fig. 2). The two subpopulations in the Marstrand fjord region show low levels of gene flow between some site (see Appendix S1: Fig. S9). Isolation by distance was evident for the GDS-based analysis (see Appendix S1: Fig. S10) and for the genetic differentiation measure Jost's D in Mantel tests (P < 0.05), but not F ST (see Appendix S1: Tables S4-S7 for Jost's D and F ST values).
Genetic network analyses.-Networks of both fjord regions together retained several links connecting them (not shown). The networks of each fjord have very different appearances: the Gullmarsfjord network retained very few connections (links) at the percolation threshold and was rather linear (see Appendix S1: Fig. S11). The Marstrand fjord region in contrast, exhibited more connections among its meadows (Appendix S1: Fig. S12). As the networks built with different genetic distance measures agreed on higher connectivity in the Marstrand region, we restricted the comparison of network properties to Goldstein distance using the percolation threshold of each network (see Appendix S1: Table S1). The network of the Marstrand region, where eelgrass has been declining drastically, showed more connectivity with 11 nodes and 23 links, while the Gullmarsfjord region only retained 7 links at its percolation threshold (Appendix S1: Table S1). In the Marstrand fjord region the network connected all sites, while in Gullmarsfjord region, four sites fell out of the network (see Appendix S1: Figs. S11 and S12). The graphs exhibit the properties of a small-world network and the networks were considerably more clustered than random networks (Appendix S1: Table S1). In the Gullmarsfjord region, nodes had a maximum of two links (Appendix S1: Fig. S11), while in the Marstrand region, nodes had up to eight links, M4 being the most connected meadow (Appendix S1: Fig. S12). The betweenness centrality measures identified four to five meadows in each region as critical stepping stones for genetic connectivity (G3, G4, G8, G11 and M1, M2, M3, M9, M10; Table 2 and Data S1, List S1).

Biophysical analyses
Biophysical modeling of seed dispersal.-The biophysical model showed that dispersal is considerable among meadows within both fjord regions, with non-zero dispersal for 27%, 51%, and 50% of all possible connections in the Gullmarsfjord, Marstrand present, and Marstrand historic data sets. Also for the sites that were assessed for genetics, modeled connectivity was high, with the Marstrand fjord region having more connections (see Appendix S1: Tables S2 and S3). Each meadow showed 38 AE 11, 64 AE 19, and 117 AE 45 (mean AE SD) connections to other meadows for the Gullmarsfjord present, Marstrand present, and Marstrand historic data sets. Self-recruitment (propagules remaining in the meadow of origin) was frequent with 45% of meadows in the Gullmarsfjord showing some self-recruitment, and 83% and 80% for present and historic meadows in the Marstrand fjord region. The maximum dispersal distance between meadows was 30 km for the Gullmarsfjord and 45 km for the Marstrand fjord region, which represent the approximate length of the two model areas.
Identification of biophysical dispersal barriers and subpopulations.-The biophysical analysis of dispersal barriers showed clusters that were very similar to population genetic differentiation (Fig. 2). In the Gullmarsfjord region, dispersal barriers where found between the Aby-, Bro-, and the Gullmarsfjord,  Fig. 2a) of the Gullmarsfjord region based on (I) current metapopulation dynamics and connectivity, i.e., the EPT-analysis ranked on the MU scale, (II) past metapopulation dynamics based on allelic (left) and genotypic richness (right), (III) past connectivity based on in-betweenness centrality in genetic network analyses, and (VI) the future potential based on demographic modeling using genetic and biophysical empirical data, ranked on a fjord region scale. consistent with the genetic results. In contrast to the analysis of genetic differentiation, no dispersal barrier was identified in the inner Gullmarsfjord (i.e., meadow G6-7; Fig. 2a). The biophysical barrier analysis based on the present distribution of eelgrass in the Marstrand fjord region also resembled the genetic results, with a dispersal barrier between S€ al€ ofjord in the south and € Alg€ o-and Hakefjord in the north. In addition, a third cluster was formed by the small meadows located on the more exposed, western side of islands at the mouth of the Hakefjord (orange dots in Fig. 2b), where no genetic sampling was performed. For the historic distribution in the Marstrand area, no stable solution for barriers was found, likely because of high internal connectivity.
Ranking of meadows.-The EPT ranking based on oceanographic dispersal probabilities of each meadow's contribution to the whole metapopulation identified many extant and historic meadows as important for the network dynamics (Fig. 3). In the Gullmarsfjord region, highly ranked meadows were found in each fjord, including most meadows sampled for genetics. In the Gullmarsfjord itself, highly ranked meadows were mostly found in the inner fjord or the mouth (Fig. 3a).
In the Marstrand fjord region, the EPT-analysis showed a high proportion of highly ranked meadows in the 1980s (Fig. 3b). Among the present Marstrand meadows, most highly ranked meadows were found on the western side of the Hake-and € Alg€ ofjord, but only few on the mainland side, where most eelgrass has been lost ( Figs. 1 and 3c). No S€ al€ ofjord meadows were ranked highly (Fig. 3c). In contrast, highly ranked meadows were present in both areas in the 1980s (Fig. 3b). The EPT-analysis that included lost meadows, to identify those that have the largest positive impact for connectivity and persistence of the present metapopulation in the Marstrand region, identified seven lost historic meadows that would best enhance connectivity and metapopulation growth among the remaining eelgrass meadows in the Marstrand fjord region (blue squares in Fig. 3c). These sites are prime candidates for restoration.

Demographic modeling
Demographic modeling of the past.-When modeling the past population dynamics, we found that a scenario of a small focal population that has first been isolated for tens of thousands of generations, and has then received gene flow at low levels from a separate population, fits our empirical data well, both in respect to the AFS of the microsatellites and the genetic diversity measures  Fig. 4; Appendix S1: Fig. S13). A scenario of two bottlenecks (strength 75-90%) is consistent with the observed data (see Appendix S1: Fig. S14). As the last bottleneck is recent (~100 generations ago), we observe a lot of variability among model runs (see Appendix S1: Fig. S13), as we indeed also observe among our sampling sites (Table 1, Fig. 4; Appendix S1: Fig. S13), indicating that the metapopulation has not yet reached a stable state after the last bottleneck.
Demographic modeling of the future.-The demographic modeling of the Gullmarsfjord with stable population sizes, showed, unsurprisingly, stable dynamics at each meadow (Fig. 5). As the effective population size (N e ) differed among sites, genetic diversity varied among meadows more than among replicate runs (Fig. 5). In general, predicted allelic richness and heterozygosity of the future 90 generations were lower than modeled for generation 20,000-28,000 (see Fig. 4 vs. Fig. 5c). It is possible that even at stable population dynamics, the Gullmarsfjord metapopulation is starting to head towards an extinction vortex of declining genetic diversity. Calculated effective population sizes are small (N e , 5-395), and the biophysical dispersal probabilities we The modeling of future genetic diversity patterns of (a, b) allelic richness, (c, d) observed heterozygosity, and (e, f) expected heterozygosity was performed with SLiM and using effective population sizes as calculated from the actual genetic data, dispersal probabilities as estimated from the biophysical modeling and decline rates as observed from ecological monitoring as input. The 11 meadows in each fjord region are indicated by 11 different colors. The replicated model runs are shown in the same color for each meadow.
Xxxxx 2020 SEASCAPE APPROACH TO CONSERVATION Article e02121; page 11 report here (see Appendix S1: Table S2), do not seem to supply enough gene flow in the Gullmarsfjord to increase genetic diversity at each site. The model for the Marstrand fjord region with declining populations, showed much higher variability in predicted heterozygosity values (Fig. 5). Notably, H obs tended to be higher than H exp , a sign of population bottlenecks (Piry et al. 2004). Population sizes reduced quickly at all modeled meadows (see Appendix S1: Fig. S15). Many meadows became extinct after dozens of generations in replicated runs (Fig. 5), and surprisingly, this included the sites for which the decline rate was set to 0% (M3, M5, M6, M8). On a more hopeful note, all modeled sites could recover after an extinction through recolonization from other modeled sites given the high recorded oceanographic dispersal probabilities in the Marstrand fjord region (Fig. 5).

Toolbox
The ranking of sites was performed in each fjord region separately and split into MUs (four for Gullmarsfjord and three for the Marstrand fjord region; Fig. 2) as suggested by both the genetic differentiation and oceanographic dispersal barrier analyses. A complete list of all extant 266 meadows included in the study is shown in Data S1, List S1. An example of one small management unit (G2) in the Brofjord of the Gullmarsfjord region is shown in Table 2. While biophysical and genetic approaches showed very good agreement on the management unit scale, there was more variation in the ranking of sites according to the metric considered on a site level. This indicates higher levels of uncertainty of the EPT analysis on a meadow scale, but for a number of sites there was also good agreement among several methods (e.g., G5 and G11, Table 2, Data S1, List S1).
A summary of different metrics per water body in the Marstrand region is shown in Table 3. Although losses of eelgrass have been high in most water bodies and management units, the analyses show that losses have mainly affected the connectivity in the S€ al€ ofjord water body severely, which has lost all high-ranked meadows (Table 3). It is also the area where restoration of lost meadows would have the best impact on the connectivity of the metapopulation in the Marstrand region.

DISCUSSION
Understanding past and present spatial population structure and identifying areas with restricted and high connectivity are critical information for spatial management, conservation, and restoration (Foley et al. 2010, Bostr€ om et al. 2011, Pittman 2017). Here, we combined different approaches to obtain this information for both past, present, and future eelgrass populations. The combination of different methods and assessing impact of eelgrass losses over long time periods on small spatial scales allowed us to identify meadows that are valuable for connectivity within the metapopulation, and those that are vulnerable to extinction at scales relevant for management. We cross-validate multiple seascape approaches and bring them together in a toolbox to prioritize limited resources based on information on value, connectivity and vulnerability. We show how the results from the study can provide direct advice for management, using our case study of eelgrass along the Swedish west coast as an example.
In our assessment, key results directly applicable to management include the genetic analysis that shows that genotypic richness is similarly high in both areas, despite dramatic losses in the Marstrand region. Both the genetic and biophysical approaches agree spatially very well and indicate higher connectivity in the declining Marstrand region compared to the stable Gullmarsfjord region. This is likely explained by the more constrained topography in the Gullmarsfjord with long and narrow fjords. The discovery clearly points out the importance of a seascape approach, as managing meadows Note: The data are presented per waterbody, the spatial unit used for assessments within the European Water Framework Directive, which overlap with the management units (MUs) M1 and M2 identified in our analyses (green and yellow dots in Fig. 2b).
Dashes indicate that there is no data available.
independently would overlook the influence of and reliance on neighboring sites (Fullerton et al. 2016) and highlights the benefit or replicating seascape genetic studies across multiple study areas (Castillo et al. 2016). However, demographic modeling of the past was able to show that the metapopulation is still affected by a recent bottleneck of~100 yr ago, which would have been difficult to infer directly from genetic data showing high diversity. Considering the indication of a delayed genetic response of historic losses, and the negative future demographic predictions for the Marstrand region, managers must work to stop further losses in both regions, and use restoration and other measures to prevent transgressing a possible extinction threshold (Hanski and Ovaskainen 2002). Combining genetic and biophysical analyses allowed us to validate the connectivity results and to identify suitable spatial management units of the metapopulations, as well as meadows most important for the connectivity within these units that could be targeted for protection. Using biophysical dispersal modeling, network analysis in combination with data on historic losses in the Marstrand area further allowed us to identify the most vulnerable areas for extinction, where measures such as restoration are most critically needed. The same data also enabled us to identify lost meadows where restoration would best benefit the present metapopulation, and genetic analyses helped identifying suitable donor meadows for the restoration.

Local management recommendations
The multiple types of analyses performed here provide guidance for agencies tasked with restoring eelgrass, as has been suggested in other areas of the world (Thom et al. 2018).
For management on a local scale, the possibly most important result was the identification of partly isolated subpopulations of eelgrass, which was supported by both genetic and biophysical analyses (Fig. 2). These subpopulations could be used as spatial management units (MUs) to ensure high connectivity and persistence of eelgrass and genetic diversity (Palsbøll et al. 2007, Allendorf et al. 2013. Interestingly, the suggested MUs overlap well with the existing coastal water bodies used for status assessment and management according to the EU Water Framework Directive (Figs. 1 and 2; Data S2, List S2), which should facilitate implementation of MUs. To exemplify how the tools provided here can be used for management within MUs, we have compiled the results for four water bodies in the Marstrand fjord region with high-quality data of eelgrass distribution and genetic data (Table 3; Data S2, List S2). The losses of eelgrass differ strongly within the region from around 5% in the North (Asker€ ofjord) to over 98% in the southern part of S€ al€ ofjord. Although the loss has decreased biophysical connectivity in the whole region, the impact is much more severe in the southern part. While losses of eelgrass have been large in the central area (Hakefjord and € Alg€ ofjord with 60% and 78%, respectively), average connectivity among remaining meadows are still high, and many of key meadows for the connectivity in the region are still present (Fig. 2, Table 3). In contrast, connectivity among the small meadows present today has decreased with 81-93% in the S€ al€ ofjord, isolating surviving meadows. Historically, this water body had the largest eelgrass area in the region, including several highly ranked meadows for the connectivity of the metapopulations in both MU-M1 and MU-M2, but all of these key meadows have been lost (Fig. 3b, c). Thus, using a combination of data on historic losses and different connectivity metrics allowed us to identify the S€ al€ ofjord water body as the most vulnerable area for extinction, where measures such as restoration are most critically needed.
Moving to a smaller spatial scale of individual meadows, we used the EPT-analysis and the fact that we have high-quality mapping data available for the 1980s to evaluate if certain areas where meadows were lost, should be prioritized for restoration actions in the declining Marstrand fjord region. For the entire Marstrand fjord region, we identified the two highest ranked lost meadows in the southern part of the S€ al€ ofjord in MU-M2, and also two highly ranked meadows in the northern part of the S€ al€ ofjord and in the € Alg€ ofjord, both belonging to MU-M1 (Fig. 3c). All four sites historically harbored very large eelgrass meadows (40-240 ha), which now are lost, and should, based on their importance for connectivity, be targeted for restoration. The genetic analyses show that local eelgrass have sufficient genetic diversity to serve as healthy donor material for restoration, demonstrating that the historic losses have not yet affected genetic diversity and fitness. Therefore, donor material for transplantations should be sourced from within the same MU as the restoration site, to favor the inclusion of locally adapted genes and avoid a potential adaptation mismatch (Jahnke et al. 2015). However, recent studies assessing the environmental conditions for eelgrass restoration at three of the four identified sites found that the loss of large eelgrass meadows have caused a local regime shift where increased sediment resuspension and drifting algal mats presently prevent eelgrass growth and restoration . Thus, other, additional measures would be required before restoration is possible, making protection of remaining meadows even more critical.
The approach presented here is also well suited for assessing the efficiency of present Marine Protected Areas (MPAs) and the need for future extensions of the MPA network. With the oceanographic EPT analysis we ranked all existing eelgrass meadows according to their importance for the whole metapopulation based on connectivity over one generation in the two study regions (in total 266 meadows). In addition, the genetic and demographic model analyses also provided information about historic connectivity and the future vulnerability of 22 meadows (see Appendix S1: Table S1 for a complete list).
Together this data can be used by managers to select optimal areas for spatial protection, and also be considered in decisions of environmental cases where eelgrass is affected by human activities, e.g., coastal exploitation for piers and marinas (Eriander et al. 2017). In the Marstrand region, with the exception of the S€ al€ ofjord, the EPT analysis identified several high-ranked sites in all water bodies (Fig. 3c), which could be candidates for spatial protection. Many of the sites on the mainland side of the water bodies and on the islands south of Nordre River are presently included in Natura 2000 MPAs. However, few MPAs are located on the west side of Hakefjord and on the islands in S€ al€ ofjord, where the largest and healthiest eelgrass meadows are presently found, which could be critical for the survival of the metapopulation within this MU. Of the sites presently not included in an MPA, the oceanographic and genetic or demographic analyses identified the sites M3, M5, M8, and M9 as important for the connectivity in the metapopulation (Fig. 2b, see Data S1, List S1), which are good targets for spatial protection.
In the Gullmarsfjord region, highly ranked meadows were found in almost all assessed water bodies and MUs (Fig. 2a, see Data S1, List S1). With the exception of MU-G2 in Brofjorden, all MUs are either completely or partially included in Natura 2000 sites. The additional important sites identified within the Brofjord, G4 and G5, are good candidates for further spatial protection. Very little documented loss has occurred in the Gullmarsfjord region since the 1980s, with the exception of a 30% loss of eelgrass on the island G as€ o at site G11. Here, large-scale restoration is presently taking place (P.-O. Moksnes, personal communication). Interestingly, this site obtained the highest overall ranking in importance for both present and historic connectivity in the Gullmarsfjord region, making it a good choice for supportive restoration. Guided by the analysis presented here, which showed high genetic diversity and high connectivity of site G11, the restoration was carried out using local donor material from within the site.

Workflow as management tool
Our study has several new features compared to other studies with similar objectives. These include (1) combination of seascape genetics and seascape ecology at spatial and temporal scales relevant for practical management; (2) biophysical modeling along very complex coasts (fjords, archipelagos); and most importantly (3) the ranking of sites according to key genetic, demographic, and ecological metrics. The comparison among methods clearly shows that the suite used together is most informative, as the combination can give information on the past (genetics), the present (biophysical modeling) and the future (demographic modeling). Our recommendation is therefore to use a combination of all three methods wherever feasible. Where time or resources are limited, high-resolution genetic assessments could in the future provide sufficient information on connectivity, as genetics combined with demographic modeling can give information on past and future gene flow. However, if high spatial resolution is needed for management, biophysical modeling might be a better sole alternative, since extensive high-resolution genetic sampling will be a challenge in most circumstances. In addition, the biophysical modeling brings the extra benefit of providing information on current changes, as well as the option to include sites in the analysis that have already been lost. Our approach of assessing a subset of samples for genetics, and then validating and extending these results on a wider spatial scale using biophysical modeling, is a cost-effective compromise between effort and results meaningful for management. We argue that the most effective way to estimate the connectivity and identify valuable and vulnerable eelgrass meadows is through an integrated, multidisciplinary framework. The framework of all approaches that we propose here provides the most comprehensive knowledge for management actions, including restoration, to ensure persistence of healthy eelgrass meadows. Nevertheless, each approach can give particular insights as detailed in the following subsections.

Genetic approaches
The genetic approaches used here give insights on evolutionary time-scales, i.e., are integrated over many generations. Because of the many possible ways in which genetic diversity may be linked to ecosystem functioning, e.g., adaptive potential and extinction risk, conservation strategies often call for preserving areas of high genetic diversity (Selkoe et al. 2016). Despite different decline rates in the two fjord regions over several decades, genotypic richness is relatively high both in the Gullmarsfjord and the Marstrand region (on average 0.69 and 0.68, respectively) and similar to a larger scale study of the Kattegat and Skagerrak (0.77; Jahnke et al. 2018). As genotypic richness does not necessarily reflect actual clonality in the populations (see our SLiM modeling and Arnaud-Haond et al. 2019), these values cannot directly be compared to other studies that potentially used a different sampling design. In comparison with other eelgrass meadows in Europe (Olsen et al. 2004;Olsen J. L. et al., unpublished data), allelic richness was also high in both fjord regions (4.30 and 4.26, respectively), even higher than what was observed previously for the Kattegat and Skagerrak (3.93;Jahnke et al. 2018). All measures of genetic diversity were variable among meadows, probably reflecting the fact that populations have not yet reached stable dynamics after the bottleneck due to the wasting disease in the 1930s or the ongoing decline (see SliM modeling). Genetic differentiation was considerable in the data set, considering the fine spatial scale assessed (sea distances between sampling sites in Gullmarsfjord are 3-67 km and Marstrand, 5-35 km; F ST 0.09 for Gullmarsfjord and 0.05 for Marstrand) with four and two genetic clusters in the two fjord regions, respectively (Fig. 2). Significant genetic differentiation on similar spatial scales has been found in an Australian seagrass using assignment tests and F ST (Sinclair et al. 2018), and in our study may be partly explained by the geography of narrow fjords that limit genetic exchange.

Biophysical approaches
Biophysical modeling can offer a more complete spatial coverage compared to genetic methods, and gives more direct demographic information on dispersal, which may better highlight immediate conservation concerns such as identifying a site at imminent risk of extinction. In contrast to genetic methods, it can also model dispersal and connectivity with historic meadows and provide information on how the loss of meadows affect connectivity (Jahnke et al. 2018), which is critical information for e.g., selecting areas for restoration. Understanding the spatial dynamics of dispersal of marine propagules, especially in the complex near-shore environment, remains a fundamental challenge in marine ecology (Carr et al. 2003, Kritzer andSale 2004). Despite the inherent difficulty to parameterize and model dispersal in near-shore environment, the biophysical model implemented here produced results that were consistent with genetic analyses, indicating the model was able to accurately approximate real dispersal phenomena of eelgrass at small geographic scales (5-67 km). The oceanographic barrier analysis identified dispersal barriers that fit the genetic breaks very well. Jahnke et al. (2018) could show a similar good fit between genetic and oceanographic barriers on a larger scale (up to 400 km) in the Kattegat and Skagerrak, but it is surprising that there is an equally good fit at the scale of kilometers in a very complex fjord system. The good agreement confirms that migration at the spatial scale of dozens of kilometers is indeed governed by oceanographic dispersal probability of rafting flowering shoots over one generation. We are aware of only one other study from Western Australia (Sinclair et al. 2018), that could show such detailed validation of genetic and biophysical estimates at the scale of individual meadows, typical and relevant for practical management.

Demographic modeling
Demographic modeling based on genetic and biophysical data is still a new approach in applied conservation studies, but is a valuable addition, as it can provide inferences about future potential for persistence (Matz et al. 2018). The allele frequency spectrum of our data and the demographic modeling clearly indicate that the whole metapopulation (in the Kattegat and Skagerrak) must have been small over thousands of generations, but received gene flow from an outside bigger population, perhaps from the North Sea proper. The demographic modeling of the past supports the assumption that the metapopulation reached an equilibrium after colonization of the area~8,000 yr ago, and shows that the assessed populations still show an imprint of a recent bottleneck.
The demographic modeling of the future is directly valuable for conservation management as it can highlight particularly vulnerable or resilient sites. Although modeled future population dynamics in the Gullmarsfjord region were stable, genetic diversity was lower than before the bottleneck. This indicates that the eelgrass in the Gullmarsfjord region could be affected by a genetic diversity extinction debt following habitat loss (Tilman et al. 1994), where the threshold condition for survival is no longer met, but extinction of alleles has not yet occurred because of the time delay in the response to environmental change (Ashander et al. 2019). In the Marstrand region, where eelgrass is declining, most modeled meadows are predicted to become extinct over the next dozens of generations if the decline continues at today's rates. This has effects on all measures of genetic diversity ( Fig. 5; Appendix S1: Fig. S17). Worryingly, even sites with stable population dynamics are frequently predicted to go extinct in the Marstrand region. However, the high dispersal probabilities predicted by the biophysical model in the Marstrand region indicates that recolonization and recovery is possible. The demographic modeling results are limited by only including the 22 sites with genetic data, and therefore underestimated the total connectivity, which warrants some caution.
Our combined approach presented here relies on a number of assumptions, which are reasonable, but uncertain. This is especially true for the demographic modeling as it combines results and assumptions from the genetic and oceanographic results, as well as having its own assumptions. Nevertheless, we suggest that the combined approach can be used as a roadmap for studies on conservation of seagrass and other threatened habitats such as coral reefs and forests.