Eelgrass shoots were collected from 15 locations along the Swedish west (Skagerrak-Kattegat and Öresund) and east coasts (Baltic Sea) spanning the entire Swedish distribution of eelgrass (Figure 1; Supplementary Table 1). Salinity and temperature during sampling were measured at most sites (Supplementary Table 1) and surface (top 3m) temperature and salinity mean for the summer periods 1976–2005 for each sampling site were retrieved from the regional climate model RCA4-NEMO (Gröger et al., 2019) with a horizontal resolution of 3.7 km (2 NM). Sampling was conducted between June and September 2020–2022 in shallow-water meadows (1–4 meters deep). The 15 sampling sites were selected based on their geographic position along a salinity gradient from the west to the east coast of Sweden and at approximately similar geographic distances (mean distance along the coastline between sampling sites was ca. 100 km). At each location, 19–20 shoots were collected by snorkeling in a “roughly linear swim” sampling one shoot every 1.5 m to minimize the probability of sampling the same clone (Arnaud-Haond et al., 2007; Jahnke et al., 2017, 2019). On the shore, four to five pieces of 2 cm leaf fragments from the clean inner leaves near each leaf sheath were collected, and epiphytes were removed with a scalpel. Samples were placed in bags and directly stored at −20°C and freeze-dried before DNA extraction. Ca. 2 cm of leaf tissue was ground with a mixer-mill (Retsch MM 301) in 1.5 ml Eppendorf tubes with tungsten beads (4 min at 30 Hz). The NucleoSpin ^® Plant II kit (Macherey-Nagel) was used to extract genomic DNA of 298 samples, followed by sodium-acetate precipitation.

The restriction-site-associated DNA (RAD) sequencing method 2b-RAD was used to generate SNPs (single nucleotide polymorphisms; Wang et al., 2012). 2b-RAD libraries were digested, ligated, and amplified following an adapted protocol (by Wang et al., 2012; https://github.com/z0on/2bRAD_denovo). In brief, for the digest, 50 ng of DNA was fragmented with the restriction endonuclease enzyme BcgI (New England BioLabs). Fragments were then ligated with adapters and barcoded, followed by amplification. The amplified fragments were extracted from a 2% agarose gel and cleaned with the NucleoSpin ^® Gel and PCR Clean-up kit (Macherey-Nagel). Four different pools of a total of 283 samples (15 samples were discarded, because of insufficient DNA quality), including 23 technical replicates (i.e., some individuals from each meadow for which library preparations were replicated twice independently to asses sequencing errors), were sequenced on the Illumina Novaseq 6000 platform, generating 50-bp sequences at the Science for Life Laboratory (SciLifeLab) - SNP&SEQ Technology Platform at Uppsala University, Sweden.

Bioinformatic analysis followed a modified de novo pipeline from Mikhail Matz, available at https://github.com/z0on/2bRAD_denovo, and carried out on the Rackham cluster of the Swedish National Infrastructure for Computing (SNIC). The samples presented here were analyzed in conjunction with data from 34 more sampling sites in the region to extract a consistent set of SNPs across previous (Faust and Jahnke, 2023) and ongoing genomic studies and to avoid potential discrepancies in genetic diversity estimates derived from different SNPs sets. Due to the short read nature of 2b-RAD reads, only forward reads were used. PCR duplicates and restriction sites were removed with a custom Perl script, low-quality bases (Q15) were trimmed with cutadapt 4.0 (Martin, 2011), and reads shorter than 25bp were removed. Trimmed reads were mapped to the eelgrass genome assembly v.3.1. (Ma et al., 2021; https://bioinformatics.psb.ugent.be/gdb/zostera/) using local alignment in bowtie2 (Langmead and Salzberg, 2012). For SNP calling, the Genome Analysis Toolkit (GATK 3.8-0; McKenna et al., 2010) was used, implementing a minDP = 3 filter (minimum depth for a genotype call) before variant quality score recalibration (VQSR). Replicates with a missingness above 25 % and a heterozygosity discovery rate higher or equal 90 % were filtered out prior to this analysis. VQSR used SNPs identically genotyped among 72 technical replicates from the conjugated data set including the remaining 18 technical replicates presented here as well as 54 additional technical replicates from 34 additional sampling sites (2,668 SNPs; transition/transversion ratio: 2.479; truth sensitivity: 99 %) as “true” variants to generate a training model to allow estimation of the probability of a SNP to be “true” and recalibrate quality scores. Using vcftools 0.1.16 (Danecek et al., 2011), an initial filter to remove individuals with high missingness (>40%) was applied, followed by a filter selecting loci with maximum heterozygosity of 50% and at least genotyped in 90% of individuals (retaining 9,069 out of a possible 13,243 SNPs). A harsh genotyping cut-off was used to aim for the best quality and avoid null alleles that affected RAD loci (Ruocco et al., 2022). Additional filtering for individuals with a missingness of >25 % was performed. The final dataset used in most analyses below consisted of 241 successfully sequenced individuals (after removing the 18 technical replicates) genotyped with 2,138 polymorphic loci.

To detect population structure and to describe the genetic differentiation among MLLs and eelgrass meadows, we performed a principal component analysis (PCA) with the R package “ade4 1.7-22” (Chessel et al., 2004; Dray et al., 2007) and used the R package “LEA 3.8.0” (Frichot and François, 2015) using the “snmf” function. LEA analyses were run three times per K for 1-15 maximum number of ancestral populations (K). The most likely number of clusters was evaluated by plotting the cross-entropy for the different numbers of ancestral populations (Supplementary Figure 2) and visually inspecting admixture plots for several Ks. Additionally, we investigated the relationship between pairwise F_ST and sea distance, salinity distance and temperature distance between sampling sites separately using Mantel tests with the R package “ncf” (Bjornstad, 2009) with 100,000 replicates (Supplementary Tables 2, 3). Sea distances without crossing land were calculated in the R package “marmap” (Pante and Simon-Bouhet, 2013) using a minimum depth constraint of −0.5 m (Supplementary Table 2). At five sampling sites, coordinates had to be marginally adjusted as the map resolution wrongly assigned them to have a positive depth. Salinity and temperature values were retrieved for each sampling site from an ensemble of 3D ocean circulation models (Gröger et al., 2019), and distances between two sites were calculated by subtracting one salinity/temperature from the other (Supplementary Table 3). Pairwise Weir and Cockerham's F_ST was calculated among sites using the “dartR 2.9.7” package. For the F_ST analysis, we removed clones, to prevent a strong increase of F_ST with the proportion of clones (Balloux et al., 2003), but note that the low number of MLLs at some sites renders F_ST values unreliable for these sites as F_ST compares intra- with inter-site variation.

Genetic diversity is a broad term encompassing several measures that quantify genetic variability within and among populations, including nucleotide diversity, heterozygosity, and percentage of polymorphic loci (Hughes et al., 2008; Allendorf et al., 2022). To estimate genetic diversity, we calculated polymorphic sites (Pol) and private alleles (Pa) using the MLL data set (n = 81) based on all 2,138 SNPs, using the functions “isPoly,” embedded in “adegenet 2.1.10” (Jombart, 2008; Jombart and Ahmed, 2011) and “private_alleles()” in “poppr 2.9.4.” To obtain unbiased estimates, the percent of polymorphic (uPol%) and private alleles (uPa) were standardized relative to the number of MLLs (Table 1). To assess nucleotide diversity and heterozygosity, we used a dataset of 5,771 loci including the 2,138 loci as well as monomorphic sites. The rationale for including monomorphic sites was to avoid underestimating heterozygous sites and their derived estimates. Observed heterozygosity (H_o = number of heterozygous sites/total number of sites with available genotype calls, customized script “individualHET.py”) and nucleotide diversity (π) were calculated following Yu et al. (2022; https://github.com/leiyu37/populationGenomics_ZM). π was calculated after filtering for missingness (resulting in 1,513 SNPs for 241 individuals and 2,092 SNPs for the 81 MLLs; also see Supplementary Table 4 for H_o filtered for missingness). The “–max-missing” function was used where no missing data was allowed (defined by 1), embedded in “vcftools 0.1.16” (Danecek et al., 2011). To display different genetic variation parameters on maps, we used the R package “wingen 1.1.0.” Finally, the above-mentioned genetic variation parameters were also calculated for each coast and genetic cluster separately, as well as for the dataset including all 241 samples.

Out of 241 samples genotyped with 2,138 polymorphic loci, we detected 79–81 MLLs among the 15 sampling sites (Table 1; Figure 2; Supplementary Figure 3 for full names). The identity of the 79 MLLs calculated with the SH index were identical to the Hemming distance calculation. The package poppr additionally identified two samples (NYN-11 and SLI-05) as belonging to MLLs (Supplementary Table 5). Notably, no clones were shared between the sampling sites. Assessment of LD and F_IS further confirmed that our chosen clonal threshold was meaningful: LD was consistently lower at each meadow for the datasets of MLLs compared to all samples, except for the meadow Slite (SLI; Figure 2C; Supplementary Table 6). F_IS was consistently closer to zero for the data sets containing only MLLs at each meadow, except for Ålabodarna (ALA), where F_IS did not change considerably (Figure 2D; Supplementary Table 6). After removing clones, the F_IS values were close to zero, except for three sites, which still had considerably negative F_IS values [Grötvik (GRO), Slite (SLI) and Nynäshamn (NYN)]. The meadows on the east coast showed a slightly higher excess of heterozygotes (F_IS = −0.064) than the west coast (F_IS = −0.054). Basing estimates of clonality on the higher established estimate of 81 MLLs, genotypic richness was generally higher on the west coast (R = 0.58) than on the east coast (R = 0.17; Table 1). Complete clonal dominance (R = 0) were found at Hornsudde (HOR), Krampö (KRA), and Björkö (BJO), all located in the northern part of the eelgrass distribution along Sweden's east coast (Table 1; Figure 1). Only one meadow on the west coast had low genotypic richness [Grötvik (GRO); R = 0.06]. The number of SNP differences within MLLs, which may be the result of somatic mutations or technical artifacts is generally low (Supplementary Table 7). On a meadow level, we observed the highest differences within an MLL at Ålabodarna (ALA; 1.08 %), Koster (KOD; 0.7 %), and Hornsudde (HOR; 0.78%). The lowest differences were found in the MLLs of Grötvik (GRO; 0.17 %) and Furumon (FUR; 0.07 %).

The PCA (Figure 3) of unique MLLs showed a gradient of genetic differentiation among meadows following the west-to-east coast (see Supplementary Figure 4 for PCA with 241 individuals). The horseshoe-shaped genetic differentiation was also along the environmental gradient of salinity. Within and among site differentiation was higher in the two Northern-most meadows - Koster and Stenungsund - than all remaining meadows combined (Figure 3). Genetic differentiation (F_ST; Supplementary Table 2) was significantly correlated with sea distance (p < 0.05) but with a low r² of 0.18. There was no significant relationship between F_ST and salinity or temperature along the entire coastline, but there was a significant relationship along the steep salinity gradient for the six sites from the west coast (r² = 0.52, p = 0.03). F_ST values on the east coast were less meaningful, due to a low number of MLLs within sampling sites, which can artificially inflate F_ST values. The cluster analysis confirmed these patterns with a gradient of differentiation between the east and west coasts (Figure 4). The most supported number of ancestral populations (K = 5; Figure 4; Supplementary Figure 2) showed stronger structuring along the west coast by separating the west coast into three clusters (dark blue, blue, light blue), but little genetic structure on the east coast (Figures 4B, D). All but one Baltic meadow belongs to one genetic cluster (yellow), except Nynäshamn, which forms its own cluster (red; Figures 4B, D). Investigating higher Ks did not lead to ecologically or evolutionary meaningful clusters (but see Supplementary Figure 5 for K = 2 to K = 10).

Like genotypic richness (Figure 5A) the genetic diversity patterns among the 81 MLLs differed between the two coasts, with generally higher genetic variation on the west coast. East coast meadows demonstrated lower percentage of polymorphic sites (percentage of uPol; 35 %) than the west coast with 65 % (Table 1; Figure 5B). The levels of private alleles (Pa), i.e., alleles that are exclusive to a meadow, were twice as high along the west coast (Table 1; Figure 5C). Observed heterozygosity (H_o) was very low throughout (ranging from 0.01 to 0.021; Table 1). Nucleotide diversity was generally higher along the west coast (Figure 5D), but the highest values were found in the east coast meadow Nynäshamn (NYN; 0.020; Table 1; Figure 5D), followed by the west coast meadow Koster (KOD; 0.017). The five lowest values were found along the east coast.

Previous to our analysis, large and locally dominant clones of eelgrass have been found around Åland in the northern Baltic Proper and Greifswalder Bodden at the German coast in the Baltic Sea (Reusch et al., 1999; Gonciarz et al., 2014; Yu et al., 2020), and no or rare flowering had been reported from the Baltic Sea (Boström, 1995; Boström et al., 2014; Salo, 2014). Our assessment of clonality using genomic methods and investigating the species' whole distributional area along the Swedish coast largely corroborate and extend these findings to much of the species' Baltic Sea distribution. We show that along the Swedish coast of the Baltic Sea, several of the sampling sites comprised only one clone, and others were dominated by a small number of clones. However, no clones were shared among the sampled meadows. When looking at differences among multiple clonemates, we found they differed on average at nine sites out of the 2,138 SNPs (0.40%). If we interpret these differences as an indication of the age of each clone, this would for instance suggest that the oldest clones were found at the west coast meadow Ålabodarna, despite the higher clonal reproduction on the east coast.

Compared to previous studies, genotypic richness (the number of different clones) was generally low in all 15 Swedish meadows assessed. A recent genomic assessment over the entire distribution of eelgrass, in both the Pacific and Atlantic, only detected a small percentage of clonemates (9%; Yu et al., 2023), as well as previous genetic assessments, based on microsatellites, along the Kattegat, Skagerrak and Southern Baltic Sea (Jahnke et al., 2018, 2020; Martínez-García et al., 2021). A high incidence of clonality in the Baltic Sea has previously been shown in other species (e.g., macroalgae; Bergström et al., 2003; Johannesson and André, 2006; Pereyra et al., 2023). A potential explanation for the high incidence of clonality in multiple facultative sexual species in the Baltic Sea is that under extreme conditions, environmental factors may be unsuitable for one or several stages of sexual reproduction. Sexual reproduction may be lost, and simultaneously clonal growth is favored, as shown for instance in the land plant Iris hexagona (Eckert, 2000; Tobler et al., 2006). Indeed, although flowering eelgrass and seed germination have been observed at as low as 5 psu in the Baltic Sea (Boström, 1995; Salo, 2014; Salo et al., 2014), this salinity level might be too low for successful seed ripening or seedling establishment. In such a scenario, clonality might increase survival in these extreme environments (Edgeloe et al., 2022), and be especially competitive during colonization and range expansions into such environments (Rafajlović et al., 2017; Pereyra et al., 2023). Indeed, the Baltic Sea turned marine only recently (~8,000 years ago) and has been colonized since then by marine species such as eelgrass (Snoeijs-Leijonmalm et al., 2017). The observed patterns may therefore be explained by the recent colonization of the Baltic Sea and the fact that clonality increases the temporal persistence of founder events through strong dominance of the initial colonizers hindering subsequent colonization (De Meester et al., 2002; Waters et al., 2013). The west coast meadows in contrast showed in general higher levels of sexual reproduction. Nevertheless, the extant of clonality of some west coast meadows is comparably high to east coast meadows - despite the marine environment and likely more historic colonization history of the west coast. Clearly, local factors play a role in the incidence of clonality as well. For instance, the meadow Grötvik has a high incidence of clonality, which may similarly be explained by recent colonization, as this is a very small (<1 ha) and isolated meadow located within a recreational marina along an otherwise exposed stretch of coast with no recorded eelgrass presence. The meadow at Koster also has relatively low genotypic diversity - however without obvious links to a more recent colonization history, which highlights the complexity of identifying drivers of clonality on the local level. An important extra advantage of clonal reproduction in modular plants like eelgrass is that during clonal reproduction, new apical shoots, responsible for propagation, stay physically connected for some time to the “king” - the leading apical shoot of the clone (Ruocco et al., 2021) and allow resource/risk sharing under unfavorable conditions (Bricker et al., 2018; Ruocco et al., 2021; Arriesgado et al., 2023). Other widespread mechanisms for clonal dominance are hybridization and/or polyploidy (Bengtsson, 2009), which are the mechanisms for the switch to clonality in the seagrass Posidonia australis (Edgeloe et al., 2022). However, neither this study nor more detailed genomic assessments have detected signs of hybridization or polyploidy in eelgrass (Yu et al., 2020, 2023).

We observed significant genetic differentiation with increasing spatial distance and detected a clear genetic separation between the east and west coast, as well as additional structuring within each coast. The Öresund on the west coast marks a genetic barrier between west and east coast meadows, as seen for many species (Johannesson and André, 2006; Johannesson et al., 2020). The salinity gradient in the sound is very steep and has been suggested as a barrier to gene flow in many species through divergent selection and variation in salinity adaptation (Johannesson et al., 2020). Along the west coast, we identified three genetic clusters in a similar pattern of genetic differentiation as we reported previously using microsatellites (Jahnke et al., 2018). Along the Swedish east coast, all but one meadow - Nynäshamn - grouped as one genetic cluster. This very distinct meadow outside Nynäshamn was identified as genetically isolated and had high heterozygosity and nucleotide diversity. Interestingly, this location was previously identified to host genetically distinct populations of another marine species, the isopod Idothea balthica, a grazer that is commonly found on eelgrass and brown algae in the Baltic Sea (Pereyra et al., 2013; De Wit et al., 2020). A possible explanation may be high amounts of ship traffic in the area with frequent ferry and container ship connections to Estonia, Gotland, Latvia, Germany and Poland. As Nynäshamn is the most differentiated location of all meadows assessed here (Figure 2), Nynäshamn is alternatively a relic population of old founding individuals of the entire region or represents a more recent colonization.

Genetic diversity quantifies the variability within populations and is likely influenced by the prevalence of clonality in seagrass populations (Jueterbock et al., 2020; Procaccini et al., 2023). Field and mesocosm experiments designed to test for relationships between genetic diversity and ecosystem function in seagrasses have shown that not only clonality but also allelic richness and/or heterozygosity are important for positive population dynamics and ecosystem services (Reynolds et al., 2012). Here we found that heterozygosity, polymorphic sites, and nucleotide diversity were generally lower on the Swedish east coast compared to the west coast - in line with clonality patterns. Relatively high genetic diversity values were previously observed along the Swedish west coast (Jahnke et al., 2018, 2020; Martínez-García et al., 2021; Duffy et al., 2022) based on microsatellites. However, recent whole genome studies indicate a generally low genetic variation of eelgrass in the Atlantic compared to the Pacific (Yu et al., 2022) due to colonization bottleneck events (Jahnke et al., 2020; Duffy et al., 2022; Yu et al., 2023). Our values for π and H_o were in the same order of magnitude as reported in Yu et al. (2023) for Atlantic populations. Notably, one Swedish site (Torseröd) was included in Yu et al. (2023) and our closest site, Stenungsund is ca. 70 km sea distance away and has approximately twice as high H_o and π values (Supplementary Table 8). We believe that this discrepancy is mostly driven by the number of MLLs/individuals assessed overall and per site, the filtering settings and the different number of SNPs. Within our dataset, we observed generally lower genetic diversity of eelgrass along the Swedish east coast - in line with the lower genotypic richness reported for eelgrass here and lower genetic variation over a wide range of taxa in the Baltic Sea described previously (Johannesson and André, 2006).

One of our key findings is the high degree of clonality in the Baltic Sea. The resulting potential lack of successful sexual reproduction along the Swedish east coast likely has an impact on the ecology and management of eelgrass compared to populations with sexual reproduction. If there was an absolute lack of sexual reproduction it would imply: (1) No seed bank that can help the natural recovery of a meadow after damage or total loss resulting in a higher risk of long-term losses after disturbance; (2) The lack of dispersing reproductive shoots reduces the potential for recolonization after losses. (3) Lack of viable seeds also means that seedlings cannot help to speed up the growth of planted eelgrass during restoration resulting in slower recovery and a higher risk of failure during restoration. (4) The high clonality and low genetic diversity also implies that Baltic eelgrass might be less resilient and able to adapt to dominant local stressors (e.g., eutrophication affecting light and oxygen levels, etc.) in addition to those associated with climate change (see below). However, it is worth noting that certain clonal meadows in the Baltic Sea have endured for a long period. While we cannot predict the future persistence of clonal eelgrass meadows eelgrass is less exposed to Labyrinthula spp. in the low salinity of the Baltic Sea, which was linked to the eelgrass wasting disease, that decimated eelgrass in the Atlantic in the 1930s (Jakobsson-Thor et al., 2019). Plasticity has likely played an important role in how these clonal meadows have survived past and ongoing stressors. Preserving thriving eelgrass meadows, whether they exhibit high genetic diversity or are clonal must be a priority. When attempting restoration, donor material needs to be selected with care to achieve increase cover, density and persistence. In summary, since natural recolonization and restoration are inherently more challenging along the Swedish east coast, it is imperative that management focuses on protecting existing meadows.

Climate change adds another dimension to the vulnerable situation of Baltic eelgrass and the importance of mapping and monitoring the genetic and genotypic diversity. Climate change is a threat to eelgrass, particularly in the Baltic Sea where eelgrass is predicted to disappear from large areas (Törnqvist et al., 2019) as (1) temperature is increasing rapidly (Reusch et al., 2018); (2) salinity is predicted to decrease (Kniebusch et al., 2019); (3) eelgrass is at the edge of the species range (Boström et al., 2003); and (4) sexual reproduction rates are low (Möller and Martin, 2007) and clonality dominant (Reusch et al., 1999). In addition to a biodiversity loss over large spatial areas (Lambers, 2015), climate change will likely result in more common extreme environments and species range shifts (Germain and Lutz, 2020). Reproductive modes are rather dynamic systems responding to such changes, and a lack of sexual reproduction likely has deep demographic consequences, affecting population expansions, recolonization and adaptation (Eckert et al., 2016; Fouqueau and Roze, 2021). Eelgrass meadows in regions currently displaying substantial genetic variation such as the Öresund, could experience shifts toward dominance in clonal reproduction. Several field and mesocosm experiments have shown that a high incidence of clonality in eelgrass generally negatively influences resistance to disturbances (Hughes and Stachowicz, 2004; Evans et al., 2017), ecosystem recovery (Reusch et al., 2005) and productivity (Hughes et al., 2008). Moreover, sexual reproduction with drifting flowering shoots is essential for the availability of propagules (Jahnke et al., 2017) dispersed with ocean currents (Jahnke and Jonsson, 2022). Maintaining the resilience of coastal ecosystems and their ecosystem functions under climate change is among the aims of several international agreements including the Kunming-Montreal Global Biodiversity Framework (UN, 2022). Therefore, temporal genetic monitoring (Schwartz et al., 2007) of eelgrass is crucial for detecting such changes (Sandström et al., 2019). Our data provides a genetic baseline and transparent data for stakeholders (Pärli et al., 2021) that can be used for future spatiotemporal genetic monitoring (Phair et al., 2021).

The data provided here is also of relevance for eelgrass restoration actions currently steeply increasing during the UN decade of restoration (Ward and Beheshti, 2023) and supported by the Global Biodiversity Framework decided on in 2022 (UN, 2022). It is documented that genetic and in particular genotypic diversity correlate with successful seagrass restoration (Reynolds et al., 2012; Evans et al., 2017) and that genetic variability of donor sites is crucial for the success of seagrass transplantation projects (Reynolds et al., 2012). Here we provide information on meadows with high genotypic and genetic diversity and define different genetic clusters. We suggest choosing donor sites from the same genetic cluster with high genotypic and genetic diversity in the local context. Because it is challenging to foresee the exact effects of introduced genetic material on the local adaptation of the transplant sites and surrounding meadows, it is important to select sites for transplantation carefully. Due to recent progress in seascape genomics (Selkoe et al., 2016) and genomic offset (Rellstab et al., 2021), there is potential for improved identification of donor material with strong adaptation traits, which could enhance resistance to predicted environmental conditions under climate change.