Genetic structure and diversity of the seagrass Zostera marina along a steep environmental gradient, with implications for genetic monitoring

Ries, Stefanie R.; Faust, Ellika; Johannesson, Kerstin; Jonsson, Per R.; Moksnes, Per-Olav; Pereyra, Ricardo T.; Jahnke, Marlene

doi:10.3389/fclim.2023.1303337

ORIGINAL RESEARCH article

Front. Clim., 18 December 2023

Sec. Climate, Ecology and People

Volume 5 - 2023 | https://doi.org/10.3389/fclim.2023.1303337

This article is part of the Research Topic Seagrasses and Climate Change: Challenges and Opportunities View all 5 articles

Genetic structure and diversity of the seagrass Zostera marina along a steep environmental gradient, with implications for genetic monitoring

$\r\nStefanie R. Ries$ Stefanie R. Ries¹

Ellika Faust¹^†

Kerstin Johannesson¹

Per R. Jonsson¹

Per-Olav Moksnes²

Ricardo T. Pereyra¹

Marlene Jahnke¹^*

¹Tjärnö Marine Laboratory, Department of Marine Sciences, University of Gothenburg, Strömstad, Sweden
²Department of Marine Science, University of Gothenburg, Gothenburg, Sweden

Zostera marina (eelgrass) is a foundation species in coastal zones in the northern hemisphere. Eelgrass is declining across its distribution, a trend likely to accelerate under climate change. In Sweden, eelgrass is a species of particular concern in management and conservation. Here, we provide information on genetic variation, an important component for the potential persistence and adaptation of any species in a changing environment. In particular, the steep salinity gradient over which eelgrass is distributed along the Swedish coast (26 psu on the west coast to 5 psu on the east coast) calls for a better understanding of genetic diversity, connectivity, and potential for local adaptation. To assess genetic variation and population genetic structure, we genotyped individuals with 2,138 single nucleotide polymorphisms (SNPs) from 15 eelgrass meadows spanning the whole Swedish distribution. We found a geographic population genetic structure from west to east parallel to the salinity gradient and with a clear genetic break at the entrance to the Baltic Sea. Meadows along the low salinity east coast consisted of a few or only one clone. Eelgrass on the west coast had higher genotypic richness, higher genetic variation, and showed population differentiation on smaller geographic scales. With their low genetic variation, the east coast meadows are especially threatened amidst global changes. Lack of sexual reproduction and the capacity to generate new genotypes is an issue that needs to be seriously considered in management and conservation. In addition, the lack of sexual reproduction renders clonal eelgrass less likely to recover and recolonize after disturbance, and more challenging to restore. The here provided information on genetic clusters, clonality, and genetic variation can be included for prioritizing meadows for conservation and for identifying meadows for restoration purposes. Most importantly, genetic monitoring is urgently needed to assess temporal genetic changes of eelgrass along the Swedish coast and elsewhere facing climate change.

1 Introduction

The seagrass Zostera marina (eelgrass) is a foundation species in coastal areas in the northern hemisphere and provides several important ecosystem services (Moore and Short, 2006; Boström et al., 2014). For instance, eelgrass meadows serve as fish nurseries (Bertelli and Unsworth, 2014; Kennedy et al., 2018), provide essential habitats for various marine organisms (Orth et al., 1984; Lilley and Unsworth, 2014), play an essential role in maintaining water quality (Moksnes et al., 2018; Reusch et al., 2021), and provide long-term storage of nitrogen and carbon (Röhr et al., 2016; Moksnes et al., 2021). Eelgrass meadows are declining worldwide, driven by various threats, mainly anthropogenic activities and climate change (Waycott et al., 2009; Boström et al., 2014). Along the Swedish west coast and in the southern Baltic Sea up to the islands of Öland and Gotland eelgrass is commonly found on shallow (<6 m) soft sediment bottoms, forming large, often monospecific meadows. Further north - up to the archipelago of Stockholm - it is less common and found in mixed meadows with limnic flowering plants (Boström et al., 2003, 2014). Eutrophication in combination with overfishing has been identified as the main driver for the historic losses in Swedish waters (Moksnes et al., 2008; Baden et al., 2010; Eriksson et al., 2023) and remaining meadows are in addition threatened by increasing small-scale coastal exploitation (Eriander et al., 2017). Future predicted drivers for eelgrass loss along the Swedish coast are direct or indirect effects of decreasing salinity and increasing temperature caused by climate change (Torn et al., 2020; Jeffery et al., 2023).

Sweden is surrounded by two central water bodies, in the east the Baltic Sea, and the west the Skagerrak-Kattegat, which are connected by the Öresund (Snoeijs-Leijonmalm et al., 2017; Figure 1). A steep salinity gradient between the Baltic Sea (east coast) and the Skagerrak-Kattegat (west coast) ranges from 3 to 26 psu (Snoeijs-Leijonmalm et al., 2017; HELCOM, 2018; Figure 1B). Due to its low sea surface salinity, the east coast presents a particularly demanding habitat for marine species. As the critical salinity growth level for eelgrass is about 5 psu, eelgrass is already at the edge of the species range along the east coast (Boström et al., 2003). Predicted freshening is likely to push eelgrass over its tolerance limit, which has been modeled to result in its disappearance from large parts of the Baltic Sea (HELCOM, 2013; Törnqvist et al., 2019). Warming and increasing frequency and intensity of marine heatwaves may further affect Swedish eelgrass, given that high temperatures reduce shoot growth and may result in shoot loss after marine heatwaves (Bergmann et al., 2010; Sawall et al., 2021; Wernberg et al., 2023). Therefore, it is important to develop a foundational understanding of the contemporary status of eelgrass meadows, encompassing key aspects such as population differentiation and genetic variation. This knowledge is important for predicting resilience, prioritizing conservation strategies, and establishing a baseline for tracking temporal changes (Evans et al., 2017).

FIGURE 1

Figure 1. Map of the 15 sampling sites along the Swedish coast covering the Swedish eelgrass distribution, which is limited to ca. 60 N along the Swedish Baltic coast. (A) shows the modeled average sea surface temperature along the Swedish coast over the period May to September over three decades (1976–2005) and (B) shows the modeled average sea surface salinity over the period May to September (1976–2005) for the upper 3 m of the water column.

Like many flowering plants, eelgrass reproduces sexually through flowering and seed production, but also reproduces clonally through vegetative growth via horizontal rhizome extension (Ackerman, 2006). Many meadows exhibit high levels of sexual reproduction, but several others are composed of a single large clone (Reusch et al., 1999; Olsen et al., 2004; Becheler et al., 2010; Reusch and Boström, 2011; Jahnke et al., 2018, 2020). In the Baltic Sea, a high dominance of clonal reproduction has been reported, where large eelgrass clones were estimated to be more than 1,000 years old (Reusch et al., 1999; Yu et al., 2020). Moreover, flowering shoots are rarely seen in the northern part of the Baltic Sea (Boström, 1995; Möller and Martin, 2007). Clonality may be advantageous when colonizing a new area and during range expansions (Rafajlović et al., 2017; Pereyra et al., 2023) and might favor survival in extreme environments (Edgeloe et al., 2022), but in terms of conservation a high dominance of clonal reproduction is generally seen as a disadvantage, particularly under global change (Pipithkul et al., 2021). Nevertheless, some ancient clones have persisted for centuries to millennia, and due to their presumed high fitness in the local environment are referred to as all-purpose-genotypes. However, the extent of their adaptability and/or persistence under climate change remains uncertain. Impacts of clonality on genetic variation and species survival are currently mostly studied using conceptual approaches, but high rates of sexual reproduction are more efficient at purging deleterious mutations, and forming new genotypes by recombination, aspects which will likely aid in adaptation under a changing environment (Otto and Lenormand, 2002; Orive et al., 2017). Assessing the number of clones in a population (also referred to as genotypic richness) is therefore one important level of genetic variation that provides the building blocks for a species to persist and adapt in a changing environment (Ehlers et al., 2008; Reynolds et al., 2016). As loss of genetic variation represents a threat in itself, by making organisms more vulnerable to environmental changes (Exposito-Alonso et al., 2022), spatiotemporal genetic monitoring assessments are being initiated (Schwartz et al., 2007). The aim of such spatiotemporal monitoring programs is to detect changes in genetic variation in response to climate change, conservation actions, or anthropogenic activities (Schwartz et al., 2007; Gaitán-Espitia and Hobday, 2021). Establishing a baseline of current genetic variation in a given geographic region is a crucial part of detecting genetic changes and it serves as a reference point to evaluate conservation efforts (Sandström et al., 2019). Genetic diversity assessments not only assist in assessing changes, but also provide valuable information for prioritization in marine spatial planning for protection and restoration (Jahnke et al., 2020).

2 Materials and methods

2.1 Sampling, DNA extractions, and library preparations

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.

2.2 Clonality

2.2.1 Multilocus lineage identification and WPGMA dendrograms

The first step in the genomic analysis of facultatively sexually reproducing species is the assessment of clonality. To distinguish clonemates from non-clonemates, we identified multilocus lineages (MLLs), i.e., clones with relatively small genetic distances which stem from somatic mutations and genotyping errors. As this is a crucial step in the analysis, we performed several complementary analyses to corroborate the results (see below). All analyses were carried out for 241 individuals and 18 technical replicates to establish error thresholds. To determine the cut-off threshold for identifying an MLL, we calculated pairwise genetic distances between samples using the function “bitwise.dist” (Hamming distance, or the number of differences between two alleles) and constructed dendrograms using the WPGMA (weighted pair group method with arithmetic mean) method with the R package “poppr 2.9.4” (Kamvar et al., 2014). To assign clones to MLLs, we calculated the maximum genetic distance cut-off using the maximum distance of the 18 technical replicates (7% of all samples) as a threshold (0.01543499; Supplementary Figure 1). Individuals with genetic distances below the set threshold were considered to belong to the same MLL and only one individual from each MLL was kept for further analyses filtering with the “poppr 2.9.4” package function “mlg.filter.”

2.2.2 Genetic similarity index (SH index)

To confirm the number of multilocus lineages established through the Hemming distance we also calculated the pairwise genetic similarity index (SH), following Yu et al. (2022). SH is calculated by dividing the number of genetic markers where two samples are identically heterozygous (N_sh) with the number of heterozygous loci (N_Het) from two samples (X₁ and X₂; Yu et al., 2022; pipeline available at https://github.com/leiyu37/Detecting-clonemates).

2.2.3 Linkage disequilibrium and inbreeding coefficient

We further corroborated the presence and level of clonality by testing the presence of linkage disequilibrium (LD) using the standardized index association r_d (Agapow and Burt, 2001) for each sampling site by calculating the association with the expectation that linkage disequilibrium is lower when excluding clonemates from the analysis (81 MLLs - the slightly higher number identified based on genetic distances) compared to assessing all sampled shoots. The R package “poppr 2.9.4” was used to measure linkage disequilibrium (LD) using the average pairwise LD of 999 randomly chosen variant nucleotide sites with one hundred replicates. We also calculated the inbreeding coefficient (F_IS) to measure the excess of heterozygotes indicated by negative values of F_IS. Clonality is expected to deviate from HWE, and including clone mates is expected to result in strongly negative F_IS values (Allendorf et al., 2022). We calculated F_IS = {1 − [mean (H_o)/mean (uH_e “unbiased expected heterozygosity”)]} per locality for the dataset only containing unique genotypes (n = 81) and the dataset including all samples (n = 241) using “gl.report.heterozygosity” embedded in the R package “dartR 2.9.7” (Gruber et al., 2018; Mijangos et al., 2022).

2.2.4 Genotypic richness (R)

Once MLLs had been established, genotypic richness (R) - the inverse of clonality-was calculated according to Dorken and Eckert (2001): R = (number of MLLs-1) / (number of samples-1) for each meadow, each identified genetic cluster (see below) and each coast (Figure 1). Genotypic richness values were then plotted as a heatmap, using the R package “wingen 1.1.0” (Bishop et al., 2023). We excluded the site Gottskär (GOT) from the interpolated heatmap because of the low number of successfully sequenced samples at this site (Table 1).

TABLE 1

Table 1. Genetic variation of eelgrass for the 15 assessed sites, the east and west coast, and the five identified genetic clusters.

2.3 Genetic differentiation

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.

2.4 Genetic diversity

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.

3 Results

3.1 Clonality

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 %).

FIGURE 2

Figure 2. (A) WPGMA dendrogram based on bitwise distance between 259 genotyped eelgrass samples from 15 different sampling sites along the Swedish coasts. Shades of blue indicate samples from the west coast, and shades of red represent samples from the east coast. The dotted line shows the threshold (0.01543499) to distinguish MLLs from sequencing errors. Individuals with nodes to the right of the threshold belong to the same MLL. Stars show technical replicates of each sampling site and the number of stars indicates matching pairs. The colors represent the sampling sites. (B) Detected clonemates of the 259 samples with the SH index. (C) Linkage disequilibrium (LD) and (D) inbreeding coefficient (F_IS) values for 241 samples and for 81 MLLs. Numbers above/below the bars indicate the number of samples (light gray) and the number of MLLs per site (dark gray). F_IS could not be calculated for meadows consisting of less than two MLLs and LD could not be calculated for meadows consisting of less than three MLLs.

3.2 Genetic differentiation

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).

FIGURE 3

Figure 3. The first (PC1) and second (PC2) components of the principal component analysis (PCA) of 81 eelgrass multilocus lineages (MLLs, i.e., clonal lineages) from 15 sampling sites along the Swedish coast. Each dot represents one individual and the symbol and color indicate its sampling site. Individuals from the west coast are colored in shades of blue, and individuals from the east coast are colored in shades of red. Barplots for eigenvalues of the first 50 PCAs are displayed, with PC1 and PC2 in purple.

FIGURE 4

Figure 4. Genetic cluster analysis of the 81 MLLs collected at 15 eelgrass meadows. In the admixture plots, each individual is represented by a vertical line on the x-axis where the color segmentation shows the individual's estimated ancestry into one of the genetic clusters. The numbers of ancestral populations shown in (A) are K = 2 and in (B) K = 5. K = 5 represents the most likely number of ancestral populations according to the cross-entropy criterion (Supplementary Figure 2). Pie charts summarize the admixture coefficients for each sampling site and project them onto a map (C) showing K = 2 and (D) showing K = 5. Acronyms of sampling sites are described in Figure 1.

3.3 Genetic diversity

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.

FIGURE 5

Figure 5. Interpolated heatmaps of clonality and genetic diversity of the 81 multilocus lineages of eelgrass along the Swedish coast. (A) genotypic richness calculated with R = {[number of multilocus lineages (nMLL)-1]/[number of individuals (nInd)-1]}. (B) The number of unbiased polymorphic sites in percent {nPol/nLoc*100*[2*nMLL/(2*nMLL-1)]}; (C) the number of unbiased private alleles per sampling site {nPa*[2*nMLL/(2*nMLL-1)]}; and (D) the values of nucleotide diversity. Bright blue indicates higher and darker blue indicates lower genetic variation.

4 Discussion

Eelgrass provides many ecosystem services and represents a priority species in conservation in Sweden (Posledovich et al., 2021), yet there are large gaps in our knowledge of the species' clonality and genetic variation. Here we provide a first population genomic assessment of eelgrass along the entire Swedish distribution, focusing on genetic variation, genetic differentiation, and clonality. We show that there is gradually increasing population differentiation from west to east and along the steep salinity gradient along the Swedish west coast. We show that along the Swedish east coast, several of the sampling sites comprised only one clone, and others were dominated by a small number of clones. Surviving in a rapidly changing environment can be challenging for organisms (Harley et al., 2006; Cossellu and Nordberg, 2010; Phair et al., 2020), especially for sessile organisms like seagrass with a relatively long generation time and low effective population size due to clonality. The potential to adapt rapidly to environmental changes is usually strongly correlated to the amount of available genetic variation (Barrett and Schluter, 2008), making eelgrass meadows with lower genetic variation along the Swedish east coast particularly vulnerable to global change. When aiming to prevent further losses of eelgrass in the short and long term, it is important to carry out spatiotemporal assessments of genetic variation, as declining genetic variation likely reduces the potential to adapt to a changing environment. We here provide a baseline for future temporal monitoring of the genetic variation of Swedish populations of eelgrass and include a description of genetic clusters and clonal meadows. This information is also instrumental for immediate management actions, such as identifying meadows that need conservation measures, and meadows that are good candidates for targeting as donor material or recipient sites in restoration efforts.

4.1 Clonality

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).

4.2 Genetic differentiation

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.

4.3 Genetic diversity

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).

4.4 Implications for management under climate change

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.

Data availability statement

The 2b-RAD raw sequences (fastq) are openly available on NCBI's sequence read archive (SAR) (BioProject PRJNA1043091, accession numbers: SAMN38331216 - SAMN38331474). The protocol for the 2b-RAD library preparations and the bioinformatic pipeline can be found at https://github.com/stefanieries/2bRAD_Swedish_eelgrass. R scripts for genetic diversity analysis, differentiation analysis, and IBD analysis, and vcf files (2,138 SNPs: 241 individuals and 81 MLLs; 5,771 SNPs: 241 individuals and 81 MLLs) can be found at 10.17044/scilifelab.24610356.

Author contributions

SR: Formal analysis, Investigation, Methodology, Validation, Visualization, Data curation, Writing – original draft. EF: Validation, Investigation, Data curation, Writing – review & editing. KJ: Conceptualization, Investigation, Funding acquisition, Writing – review & editing. PJ: Conceptualization, Investigation, Supervision, Funding acquisition, Visualization, Writing – review & editing. P-OM: Conceptualization, Investigation, Supervision, Funding acquisition, Writing – review & editing. RP: Conceptualization, Investigation, Methodology, Validation, Writing – review & editing. MJ: Conceptualization, Supervision, Formal analysis, Funding acquisition, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study was funded by Swedish Research Council FORMAS (Grant No: 2020-008) and from funding for genetic monitoring from the Swedish Agency for Marine and Water Management (Grant No: 1717-22). The work was performed within the Linnaeus Centre for Marine Evolutionary Biology (https://www.gu.se/en/cemeb-marineevolutionary-biology) and the multidisciplinary research program Zorro (https://www.gu.se/en/research/zorro). MJ also wishes to acknowledge support through the Carl XVI Gustafs 50-årsfond för vetenskap, teknik och miljö, BalticWaters2030, The Royal Swedish Academy of Sciences and Stiftelsen Birgit och Birger Wåhlströms Minnesfond för den bohuslänska havs- och insjömiljön.

Acknowledgments

We want to thank the following people who contributed to sampling and laboratory work: Matilda Rasmussen, Henrik Möller, Swantje Enge, Bo Johannesson, Lydia Williams, and Maja Billman. Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala. The facility was part of the National Genomics Infrastructure (NGI) Sweden and Science for Life Laboratory. The SNP&SEQ Platform was also supported by the Swedish Research Council and the Knut and Alice Wallenberg Foundation.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

KJ declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fclim.2023.1303337/full#supplementary-material

References

Ackerman, J. D. (2006). “Sexual reproduction of seagrasses: pollination in the marine context,” in Seagrasses: Biology, Ecology and Conservation, eds A. W. D. Larkum, R. J. Orth, and C. M. Duarte (Dordrecht: Springer Netherlands), 89–109. doi: 10.1007/1-4020-2983-7_4