Salinity Driven Selection and Local Adaptation in Baltic Sea Mytilid Mussels

Knöbel, Loreen; Nascimento-Schulze, Jennifer C.; Sanders, Trystan; Zeus, Dominique; Hiebenthal, Claas; Barboza, Francisco R.; Stuckas, Heiko; Melzner, Frank

doi:10.3389/fmars.2021.692078

ORIGINAL RESEARCH article

Front. Mar. Sci., 13 August 2021

Sec. Global Change and the Future Ocean

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.692078

This article is part of the Research Topic Adaptation and Phenotypic Plasticity to Climate Change View all 10 articles

Salinity Driven Selection and Local Adaptation in Baltic Sea Mytilid Mussels

$\r\nLoreen Knbel&#x;$ Loreen Knöbel^1†

Jennifer C. Nascimento-Schulze^2†‡

Trystan Sanders^2†‡

Dominique Zeus¹

Claas Hiebenthal²

Francisco R. Barboza^2‡

Heiko Stuckas^1*§

Frank Melzner^2*§

¹Senckenberg Natural History Collections Dresden, Senckenberg Leibniz Institution for Biodiversity and Earth System Research, Dresden, Germany
²Marine Ecology, GEOMAR Helmholtz Centre for Ocean Research Kiel, Kiel, Germany

Baltic blue mussels can colonise and dominate habitats with far lower salinity (<10 psu) than other Mytilus congeners. Pervasive gene flow was observed between Western Baltic Mytilus edulis living at high salinity conditions and Eastern Baltic M. trossulus living at lower salinites, with highest admixture proportions within a genetic transition zone located at intermediate salinities (Darss Sill area). Yet, we do not understand the impacts of low salinity on larval performance, and how salinity may act as an early selective pressure during passive larval drift across salinity gradients. This study tested whether larvae originating from two different populations along the natural salinity cline in the Baltic Sea have highest fitness at their native salinities. Our results suggest that Eastern Baltic M. trossulus (Usedom, 7 psu) and Western Baltic M. edulis (Kiel, 16 psu) larvae display better performance (fitness components: growth, mortality, settlement success) when reared at their respective native salinities. This suggests that these populations are adapted to their local environment. Additionally, species diagnostic markers were used for genetic analyses of transition zone (Ahrenshoop, 11 psu) mussel larvae exposed to low salinity. This revealed that low salinity selection resulted in a shift towards allele frequencies more typical for Eastern Baltic M. trossulus. Thus, salinity acts as a selective pressure during the pre-settlement phase and can shape the genetic composition of Baltic mussel populations driving local adaptation to low salinity. Future climate change driven desalination, therefore, has the potential to shift the Baltic Sea hybrid gradient westward with consequences for benthic ecosystem structure.

Introduction

The pronounced natural salinity gradient in the Baltic Sea, ranging from 30 psu at the transition to the North Sea to less than 4 psu in the innermost sub-basins (Bothnian Bay and inner Gulf of Finland) offers a natural context in which to explore evolutionary mechanisms underlying the colonisation of extreme environmental conditions (Johannesson et al., 2011; Wennerström et al., 2013). Apart from evolutionary questions, understanding how Baltic organisms cope with extreme environmental conditions holds the key to predicting and mitigating consequences of global change on marine biodiversity (Reusch et al., 2018). In this marginal sea, climate models predict a decrease of sea surface salinities by ca. 2 psu and a westward shift of the salinity gradient by the end of the century (Gräwe et al., 2013). Additionally, the Baltic Sea has been identified as one of the fastest warming marine basins on the planet (Belkin, 2009).

Mytilid mussels have successfully colonised the Baltic Sea presumably 9000-8000 years ago after the last freshwater stage (Berglund et al., 2005; Witkowski et al., 2005; Behre, 2007; Kostecki and Janczak-Kostecka, 2011). Today, the two species Mytilus edulis and M. trossulus occur in the Baltic Sea (e.g., Koehn, 1991). Both are important foundation species which contribute significantly to benthic ecosystem function and nutrient turnover (Norling and Kautsky, 2008; Heckwolf et al., 2021). Mytilus edulis colonizes the Western Baltic (e.g., Kattegat, Belt Sea) at salinities between 12 and 25 psu and M. trossulus is found in the Eastern Baltic (e.g., Baltic Proper) at salinities between 4.5 and 8 psu (e.g., Stuckas et al., 2009, 2017; Kijewski et al., 2019). Whilst gene flow between these species is restricted on the Atlantic coastline of Northern America, a hybrid swarm is found in the Baltic (e.g., Riginos and Cunningham, 2005; Wenne et al., 2020). However, pervasive interspecific gene flow does not lead to an amalgamation of the gene pools. Instead, a pronounced genetic structure is found where Western Baltic M. edulis and Eastern Baltic M. trossulus are dominated by their species-specific alleles (e.g., Stuckas et al., 2009, 2017; Fraïsse et al., 2016; Kijewski et al., 2019; Simon et al., 2019). Species-specific allele frequencies change gradually resulting in allele frequency clines and a transition zone that is situated along longitude 12–13°E (a virtual line between Öresund and Darss Sill). This is an area of maximum genetic admixture where repeated backcrossing results in the coexistence of individuals with highly different proportions of alleles specific to M. edulis (ME-alleles) and M. trossulus (MT-alleles). Consequently, the species identity of transition zone mussels is best described by a hybrid index, e.g., continuously ranging from 0 (100% ME-alleles) via 0.5 (F1-like hybrid) to 1 (100% MT-alleles) (Stuckas et al., 2017). Despite of pervasive interspecific gene flow between Baltic mussels, substantial dissimilarities in shell morphology and physiology between Western Baltic M edulis and Eastern Baltic M. trossulus still exist (Kautsky et al., 1990; Tedengren et al., 1990; Telesca et al., 2018).

Adaptation to local salinity conditions is presumed to be one mechanism underlying the maintenance of distinct phenotypic differences in Baltic Sea Mytilus. This hypothesis was tested in field transplantation experiments of adult mussels. These analyses documented negative impacts of low salinities (<10 psu) particularly on Western Baltic M. edulis, with reduced filtration and growth rates as well as extremely high mortalities (Kautsky et al., 1990; Kossak, 2006; Riisgård et al., 2012, 2013, 2014). However, growth, survival, and feeding rates of Baltic M. trossulus individuals were comparable to those of Baltic M. edulis populations at high salinities (25–30 psu) (Kautsky et al., 1990; Tedengren et al., 1990; Riisgård et al., 2014). These specific responses of Baltic Mytilus species to salinity do not only hint toward the existence of local adaptation to salinity but also suggests that Baltic M. trossulus adults have a broader salinity tolerance window.

However, the implications of salinity changes on pre-settlement stages of Baltic mussel populations are yet to be understood. This is an important research topic, as larvae can drift for long distances during their several-week long pelagic life phase, thereby being exposed to fluctuating salinities (Larsson et al., 2016; Stuckas et al., 2017). In addition, mussel early life stages have been shown to be much more susceptible to abiotic stressors than adults, with a particularly high sensitivity to low calcium concentrations in low saline basins such as the Baltic (Thorson, 1950; Thomsen et al., 2015, 2017; Ramesh et al., 2017).

Here we explored how salinity impacts pelagic larvae during the pre-settlement life stage addressing two main questions: (i) whether larvae of both Baltic Mytilus species as well as transition zone mussels show differential salinity dependent performance (mortality, growth, settlement success) and (ii) whether salinity stress acts as a selective pressure in these early life stages, possibly shaping the genetic composition of mussel populations along the salinity gradient.

We first conducted a population comparison experiment using a common garden design with larvae originating from Western Baltic M. edulis (Kiel – KIE -, genotypes dominated by ME alleles, Figure 1) and Eastern Baltic M. trossulus (Usedom – USE -, genotypes dominated by MT alleles, Figure 1). Larval performance was quantified at both native and reciprocal salinities. A second experiment utilised larvae raised from transition zone mussel parents (Ahrenshoop, - AHP -, genotypes with intermediate frequencies of both ME- and MT-alleles, Figure 1). Although methodologically challenging, we used laboratory larval rearing trials for the first time at salinities <10 psu. These larvae were reared at low salinities analogous to what they could be exposed to during eastward drift trajectories along the southern Baltic coast. Performance of larvae was monitored in response to lowered salinity and at two temperatures prior to genotyping with species diagnostic markers upon settlement. Under the assumption of local adaptation to distinctly different habitat salinity regimes, both experiments were expected to show better performance of larvae under native compared to altered salinity. In addition, we hypothesised that low salinity treatment under laboratory conditions causes selection on transition zone mussels and predicted a significant enrichment of hybrid genotypes dominated by M. trossulus specific alleles.

FIGURE 1

Figure 1. Distribution of Baltic Mytilus species and sampling localities. (A) Sampling locations for the population comparison (Kiel, KIE; Usedom, USE) and the desalination (Ahrenshoop, AHP) experiments. Further references locations for genetic analyses are shown ECK, Eckernförde; FSD, Fehmarnsund; GRO, Grömitz; STB, Steinbeck; GWZ, Gollwitz; WMU, Warnemünde; BAR, Barhöft; RUD, Rügen Island (Dranske); USE, Usedom (see further information in Supplementary Figure 5). (B) Distribution of hybrid classes expressed as Q_MT-values (0 = genotypes dominated by M. edulis alleles, 1 = genotypes dominated by M. trossulus alleles). Experimental data (dark grey, animals sampled in 2016 in AHP, parental animals: AHP-P, larval offspring populations: AHP-F1 and reference data from Stuckas et al. (2017) (light grey, animals sampled in 2013, including average salinity values) are given as violin plots. Median Q_MT values are indicated by red diamonds. Compare to Supplementary Figure 3 for full dataset and STRUCTURE bar plots. Salinity is given for each locality based on hydrodynamic model data obtained from Stuckas et al. (2017).

Materials and Methods

Broodstock Sampling and Maintenance

Adult Baltic Mytilus individuals were collected from three populations along the Baltic Sea coastline at the localities (KIE, average salinity 16 psu), Ahrenshoop (AHP, average salinity 11 psu), and (USE, average salinity 7 psu) (Figure 1 and Table 1) in spring 2016. Average salinities at each locality were analysed over a 2 to 3-year period and this information was used to calculate average salinity conditions for experimental procedures described below (Sanders et al., 2021). Considering the high levels of interspecific gene flow (see introduction), we will refer to Baltic mussels as “Western Baltic M. edulis” (KIE), “transition zone mussels” (AHP), and “Eastern Baltic M. trossulus” (USE) throughout the manuscript. Mussels (15–50 mm shell length) were collected and transported to GEOMAR in continuously aerated, cooling boxes and kept at 10°C in 10 L plastic aquaria (N = 20 mussels per aquarium) with 20 μm – filtered sea water (FSW) at their native salinities (Table 1). Water (ca. 50% of tank volume) was changed daily.

TABLE 1

Table 1. Parental mussel sampling sites and abiotic conditions.

Population Comparison Experiment (Kiel vs. Usedom)

The experiment aimed at testing whether larval stages originating from Western Baltic M. edulis collected in KIE and Eastern Baltic M. trossulus collected in USE are adapted to local salinity regimes and whether they can rapidly acclimate to salinity changes. The experimental design is depicted in Supplementary Figures 1, 2. Two separate larval pools originating from parental mussels from each of the two populations were compared under high and low experimental salinity conditions. Larvae used in the experiment were generated using a group spawning approach, followed by pooling gametes within each population to obtain offspring representing all possible parental combinations. Although this method may potentially mask the intra-population variability compared to analyses of individual families, it guaranties absolute synchronous spawning of mussels and imitates the natural diversity of family representation in the larval pool.

All spawning and fertilisations were conducted within each population (KIE, USE) separately at their native salinity. Spawning was induced by increasing water temperature (10 to 18°C within 30 min) and individuals were immediately isolated in 1 L beakers containing filtered seawater (FSW). Once spawning started, water was gently, periodically stirred and gamete quality was checked (active sperm, circular eggs). After 30 min, gametes were counted under the microscope for each individual and pooled at equal ratios separately for males and females to ensure equal representation from each parent. Pooled gametes of 4 females and 6 males were mixed at a 1:100 egg to sperm ratio. Fertilisation success was determined by monitoring the presence of a polar body and/or cell cleavage within 30 min of gamete mixing. Embryo pools from each population were split into two salinity treatments (16 and 7 psu, 4 treatments in total; with 6 technical replicates). The adjustment to reciprocal salinities (Baltic M. edulis: from 16 to 7 psu; Baltic M. trossulus: from 7 to 16 psu) was performed at a rate of 3 psu per day conducting 50% water changes daily for the first three days (see further details of water changes below). Larvae were cultured in 10 L glass Duran bottles (Schott, Germany) at a density of 15 larvae⋅mL^–1 at 17°C for the first 36 h until shell formation. Afterwards shelled D-veliger larvae were transferred to custom made 14 L plastic (PVC) conical culturing vessels at a density of 10 larvae⋅mL^–1 (Supplementary Figure 2). Mild aeration (ca. 50 mL air min^–1) from below prevented sedimentation and mechanical damage to larvae.

From 3 days post fertilisation (dpf) onward, larval density and growth measurements were performed every third day on larvae contained in triplicate 5 mL water samples which were taken with a Pasteur pipette. Water changes were performed every third day, by gently draining experimental water through cylindrical 65 μm mesh filters keeping larvae submerged at all times. Subsequently, filters containing larvae were inverted and gently rinsed back into experimental vessels with freshly prepared FSW. Pilot experiments revealed that this method did not alter larval densities. Larvae were fed daily (from 3 to 16 dpf) with live cultures of Isochrysis galbana (Prymnesiophyceae) at a final concentration of 15000 cells⋅mL^–1 with a cell diameter of ca. 4 μm, which is optimal for early-stage mussel larvae. After 16 dpf, larvae were fed larger Rhodomonas sp. (Cryptophyceae; cell diameter of ca. 7 μm) for the remainder of the experiment at a final concentration of 3000 cells⋅mL^–1. Larvae were fixed in Para-formaldehyde (4% PFA, pH = 8, corresponding salinity) for quantitative analyses using a stereo microscope and camera at 100x magnification (Leica Microsystems GmbH, Wetzlar, Germany). Larval size was determined for ca. 15 larvae per replicate culture vessel by measuring shell length using ImageJ 1.x (Schneider et al., 2012) at multiple time points up until 23 dpf, despite this experiment running for 33 days in total. This was due to the loss of larger individuals from the water column due to settlement after 23 dpf. Larval mortality was expressed as the decrease in larval density (percentage of the initial number of larvae⋅mL^–1) over time and monitored until the termination of the experiment at 33 dpf. At this point of time, settlement success was quantified after removing settled larvae from culture tanks with a 5 cm × 5 cm × 5 cm synthetic sponge and rinsing onto a 65 μm mesh. Collected juvenile mussels were counted using a stereo microscope and settlement success was expressed as the ratio of settled individuals to initial absolute number of larvae.

Desalination Experiment (Ahrenshoop)

This experiment exposed larvae from transition zone parents collected in AHP to native habitat salinity (11 psu) and two lower salinities (9, 7 psu). These reduced salinities reflect the range that is expected during eastward passive drift and approximates the range of predicted desalination of their native habitat over the next century (see introduction). By using a parental population with intermediate frequencies of both ME- and MT-alleles, we tested for significant enrichment of species-specific alleles in the group of surviving F1 offspring reared at lower salinities.

The experimental design is depicted in Supplementary Figures 3, 4. Larval rearing, maintenance, feeding and measurements of fitness proxies were conducted as described in the population comparison experiment with the following modifications: mussels were placed individually in 100 mL plastic beakers containing 50 mL of 11 psu FSW in a 10°C water bath. Spawning was induced by increasing the water temperature to 18°C over 30 min. Gametes of 6 females (F1–F6) and 6 males (M1–M6) were mixed in a pairwise fashion in equal numbers to generate six families. After successful fertilisation, 40,000 embryos from each of the six families were pooled in 3 separate new vessels (1 L) in equal ratios (240,000 embryos per vessel). The volume of each of the 3 vessels was then divided equally into 10 × 2 L glass Duran bottles (a total of 30 experimental vessels). Each bottle contained an equal ratio of the 6 families (24,000 embryos in total; 4,000 from each family with a larval density of 12 larvae mL^–1). Salinity in 10 experimental bottles was kept at 11 psu. In the remaining 20 bottles, salinity was reduced to 9 psu (10 bottles) and 7 psu (10 bottles), each at a rate of 1 psu every 12 h. Within each salinity treatment, half of the bottles (5) were maintained at 12 ± 0.01°C while the other half (5) were maintained at 15 ± 0.04°C, following a two-day temperature adjustment period (1.5°C day^–1). In summary, the experiment consisted of 6 salinity-temperature treatment combinations with 5 technical replicates containing 12 embryos mL^–1 (Supplementary Figure 3). Tanks were aerated via 20 mL plastic Pasteur pipettes at a rate of ca. 50 mL air min^–1 to provide fully air saturated conditions without causing mechanical damage to larvae. Shell growth (relative increase in shell length over time) and mortality (relative decrease in larval density over time) were recorded until 45 dpf despite the experimental period lasting for 67–69 days. This was because of the onset of settlement. Finally, settlement success (relative number of larvae settled) was estimated at the end of the experiment at 67–69 dpf and settled juvenile mussels were collected for genotyping at this time point.

Statistics

Statistical analyses were performed with Rstudio v0.1.0.44 (RStudio Team, 2015) using R v0.3.3.2 (R Core Team, 2014). Growth rates (increases in shell length at widest point over time) were described using the von Bertalanffy growth model (VBGM) modified by Beverton and Holt (1957). In this model, growth proceeds toward a maximum asymptotic size limit, decreasing as the size approaches the maximum limit (L_∞). This was done using the equation ⋅L(t) = L_∞(1−e^−k(t−t0)), where L(t) is the expected or average length at time (age) t, L_∞ is the asymptotic average length, K is the body growth rate coefficient, and t0 is a modelling artefact to represent time/age when average length was zero. For both experiments, parameters for the general VBGM model (L_∞ and K and t0) were estimated using the FSA R package. For comparison between treatments, subset models were created using a common t0 value (one for each experiment) and L_∞ and K separately estimated for each the treatment level of both the simulated desalination (factors: salinity and temperature) and population comparison (factors: salinity, population) experiments. Designating a common t0 is justified, as larvae used in different treatments originated from the same cohort. For each individual experiment, the Akaike Information Criterion (AIC) was used to define the most parsimonious model representing the shell length data. Mortality rates (changes in larval densities over time) were analysed using two-way ANCOVA (covariable: time (dpf), fixed factors: population and salinity or, salinity and temperature), and pairwise comparison among means of the different groups were tested with the Tukey’s HSD post hoc test. Mortality/density data related to the population comparison experiment were visually inspected for homogeneity of variances (data distribution in the residual vs. fitted plot) and deviations from normality (density plot and QQ plot). Data obtained from the desalination experiment were checked for homogeneity of variances using a Fligner-Killeen test and checks for normality were performed using a Shapiro-Wilk test. Larval settlement was analysed by comparing the settlement probability between treatment using a Generalised Linear Model (GLM) with binomial distribution and logit link function for both the simulated desalination (factors: temperature, salinity) and population comparison (factors: salinity, population) experiments.

Genetic Analyses Associated to the Desalination Experiment

Genotyping

Genetic analyses of specimens aimed to identify ME- and MT-alleles. Genotyping included parental mussels sampled at AHP in spring 2016 (AHP-P, six females and six males, Table 1 and Figure 1) and successfully settled individuals (juvenile mussels) sampled at the end of the desalination experiment from different treatments (AHP-F1, experimental populations). DNA extraction from AHP-P mussels was performed as described in Stuckas et al. (2017). Juvenile mussels from each replicate were rinsed with PBS and stored at −20°C until DNA extraction in 20 μL LoTEPA buffer according to Zhan et al. (2008). Genotyping of AHP-P and AHP-F1 was performed at four species diagnostic single copy nuclear markers that targeted non-coding (Glu5′, EFbis, mac-1) and coding (M7 lysin) genomic loci. These markers were intensively used for Mytilus genotyping, species identification and hybrid zone analyses and, therefore, used to characterize the distribution of Western Baltic M. edulis, transition zone mussels, and Eastern Baltic M. trossulus (e.g., Riginos and Cunningham, 2005; Stuckas et al., 2009, 2017). Protocols were adapted to juvenile mussel genotyping (Supplementary Tables 1, 2).

Simulated Reference Populations

The analyses of whether salinity change has an effect on the genetic composition of experimental larval populations would ideally be based on genetic profiling before and after laboratory treatments. However, sampling of individuals shortly after fertilisation (before treatment) was not feasible. This problem was circumvented by predicting the before treatment genetic profile of larval populations using simulations (simulated reference populations). This was done using the genetic data from AHP-P and mimicking the breeding scheme used for the desalination experiment (see above, Supplementary Figure 3). In a first step, 4,000 offspring genotypes of each of six parent couples were simulated using Hybrid Lab v0.1.0 (Nielsen et al., 2006) and assuming random association of parental gametes and alleles. Second, simulated offspring genotypes of these six couples were pooled (24,000 simulated genotypes in total), randomly mixed, and split into 10 virtual populations (2,400 simulated individuals each). Finally, random mortality was simulated by randomly drawing 20 virtual individuals (“survivors”) from each simulated population. This procedure was repeated three times (accounting for the three salinity treatments, Supplementary Figure 3) resulting in three sets of simulated reference populations (AHP-F1, simulation 1–3). They reflect the genetic composition of larval cohorts before salinity selection under genetic drift expectations.

Hybrid Indices Using Bayesian Inference

Bayesian inference (STRUCTURE v0.2.3.3; Falush et al., 2003) was used to estimate the proportion of ME-alleles and MT-alleles in each individual. This method used genetic data from the genotyping step (see above) and assigns individuals to genetic clusters fulfilling Hardy-Weinberg and linkage equilibrium expectations. Assignment probabilities to a given genetic cluster is expressed as Q-value ranging from 0 (no assignment) to 1 (full assignment) while values between 0 and 1 indicate admixture between clusters. Thus, Q-values can be interpreted as hybrid indices expressing to what proportion the genome of each individual is composed of alleles from one or the other genetic cluster. Bayesian inference was performed for individuals of (i) AHP-P, (ii) AHP-F1 (experimental populations), and (iii) AHP-F1 (simulations 1–3). This was done with reference to genetically pure populations of M. edulis (Helgoland, North Sea, coordinates: 54.18, 07.89), M. trossulus (Penn Cove, North America, coordinates: 48.21, −122.70), and 19 populations of Baltic mussels analysed in a previous investigation (Stuckas et al., 2017). This included previously analysed samples from KIE, AHP, and USE that were collected in 2013. A complete list of samples included in Bayesian inference is shown in Supplementary Table 3.

Analyses were performed using the admixture model under the assumption of correlated allele frequencies (Falush et al., 2003). All other parameters were set to default. The possibility for K 1–10 genetic clusters was tested performing 20 independent runs for each K and a Markov Chain length of 10⁶ including a 25% burn-in period. The best K was chosen following Evanno et al. (2005). Further analyses of STRUCTURE Q-value distribution were performed with Rstudio v0.1.0.44 (RStudio Team, 2015) using R v0.3.3.2 (R Core Team, 2014). As settled larvae could only be recovered from 9 out of 15 replicate units at 12°C and from 14 out of 15 replicates at 15°C, we decided to pool temperature replicates for the analysis of the influence of salinity treatment on Q-value distributions. This was done as a Welch t-test suggested that temperature (across salinity treatments) did not significantly impact Q-value distribution. The effect of salinity treatment on Q-value distribution of the three F1 animal groups at salinities of 7, 9, and 11 psu in comparison to the three simulated F1 populations was then analysed using one-way ANOVA followed by Tukey’s HSD tests after visually verifying normal distribution and homoscedasticity (see above). Data were displayed as violin plots (using geom_violin in ggplot2).

Ethical Approval

All applicable international, national and/or institutional guidelines for the care and use of animals were followed.

Results

Population Comparison Experiment

With our improved larval culturing system, we successfully reared Baltic Sea mussel larvae at low salinities (<10 psu) to settlement for the first time in Western Baltic M. edulis originated from KIE (16 psu) and Eastern Baltic M. trossulus from USE (7 psu). Low salinity had an effect on fitness-related parameters in both Baltic M. trossulus and Baltic M. edulis larvae. Developmental delay was observed in both populations at 7 psu compared to 16 psu with a ca. 12 h delay in formation of the complete D-veliger shell, accompanied by significantly slower growth rates at 7 psu in both populations (Figures 2A,C and Supplementary Tables 4, 5). Final shell lengths were approximately 17% lower at 7 psu than at 16 psu at 23 dpf, i.e., the end of the period analysed for shell growth (Figure 2). Larval mortality estimated from larval density changes was significantly higher at 7 psu than at 16 psu (ANCOVA, F_(1,259) = 6.6, p < 0.001; Supplementary Table 6 and Figure 3A) only in Baltic M. edulis (ANCOVA, Tukey HSD post hoc test, p < 0.001; Supplementary Tables 6, 7 and Figure 3A). Mortality rates of Baltic M. trossulus did not differ significantly between both salinities (ANCOVA p = 0.41; Supplementary Table 7), indicating that Baltic M. trossulus larvae exhibit better performance at low salinity than Baltic M. edulis larvae. Settlement began in all treatments at ca. 25 dpf and the experiment was terminated at 33 dpf when no larvae remained in the water column. Low salinity decreased settlement success in both populations (GLM, p < 0.001; Supplementary Table 8 and Figure 3B), with settlement probability of Baltic M. edulis at 7 psu being particularly poor (0.3% compared with 1.89% at 16 psu). Yet, Baltic M. edulis larvae had a higher settlement success at 16 psu compared to Baltic M. trossulus (Tukey HSD, p < 0.001; Supplementary Table 9 and Figure 3B), and Baltic M. trossulus exhibited higher settlement success than Baltic M. edulis at 7 psu (Tukey HSD, p < 0.001; Supplementary Table 9 and Figure 3B).

FIGURE 2

Figure 2. Growth rates. Shell length (μm) increase over time of all treatments throughout the course of populations comparison (A) and desalination experiments (B) fitted to von Bertalanffy growth models (VBGM). Mean asymptotic average length (L_∞) and mean Brody growth coefficient (K) estimated using the VBGM for the different treatments of the population comparison and simulated desalination experiments are represented in panels (C,D), respectively. Symbols represent means and 95% confidence interval.

FIGURE 3

Figure 3. Larval density in cultures and probability of settlement success for all treatments throughout the course of the population comparison (A,B) and simulated desalination experiments (C,D). Plots (A) and (C) represent changes in larval density over time which indicates mortality and (B) and (D) the probability of settlement success. Symbols represent means and SE.

Desalination Experiment: Larval Performance

Exposing larvae from transition zone mussels (AHP, 11 psu) to desalination and warming revealed strong impacts on larval fitness components. Exposure to desalination reduced shell length growth rates of transition zone larvae (Figures 2B,D and Supplementary Tables 4, 5) with slower rates of growth at both 9 and 7 psu. Mortality rate expressed as larval density change was also higher at low salinities (ANCOVA, F_(2,558) = 8.9, p < 0.001; Supplementary Table 6 and Figure 3C) with mortality being highest in the 7 psu treatment (ANCOVA post hoc test, p < 0.001; Supplementary Table 7) and not significantly different between 9 and 11 psu (ANCOVA post hoc test, p = 1.0, Supplementary Table 7). Settlement probability was negatively impacted by both desalination and temperature (GLM, p < 0.001; Supplementary Tables 8, 9 and Figure 3D) with larval settlement being significantly different between all three salinities only at 15°C (Tukey HSD, p < 0.001; Supplementary Table 9). At 12°C, settlement success was significantly higher only at 11 psu compared to both 9 and 7 psu treatments (Tukey HSD, p < 0.001; Supplementary Table 9). At all salinities, settlement success was higher at the higher temperature (15°C) (Figure 3D). A significant interaction between salinity and temperature was also observed for mortality and shell length growth, i.e., shell growth and mortality were more strongly reduced by low salinity at 12°C, compared to 15°C (Tukey HSD, Supplementary Table 7 and Figures 2B,D,3C). This demonstrates the interactive role of both temperature and salinity in influencing fitness related parameters in Baltic Mytilus larvae.

Desalination Experiment: Genetic Analyses

A total of 296 settled juveniles (<1 mm shell length; AHP-F1, experimental populations) from the six treatments of the desalination experiment were genotyped along with parental animals (AHP-P, 6 males, 6 females). Subsequently, genotypes of parental animals were used to predict the genetic profile of their offspring (simulated reference populations; AHP-F1, simulation 1–3) reflecting genetic drift expectations. These genetic data were used for Bayesian inference along with reference animals from the North Sea (genetically pure M. edulis), and the western Pacific Ocean (genetically pure M. trossulus), and additional Baltic reference populations from a previous survey in 2013 (Figure 1A, Supplementary Figure 5 and Supplementary Table 3; Stuckas et al., 2017). These individuals were assigned with different probabilities to two genetic clusters, an M. edulis-cluster (expressed by the Q_ME-value) and an M. trossulus-cluster (expressed by an Q_MT-value). As both Q-values add to 1, we only use the Q_MT-values to describe assignment results and to express the hybrid index, i.e., the relative proportion of MT-alleles. An excerpt of the complete results (Supplementary Figure 4A; STRUCTURE bar plots) is presented Figures 1B, 4. In Figure 1B, we show the distribution of Q_MT-values of AHP-P and experimental AHP-F1, with all animals pooled for each salinity group. Animals were pooled for this visualisation, as we observed no significant differences in Q_MT-values between temperatures (Welch’s t-test, t = −0.76, df = 20.9, p > 0.45). We also added Q_MT-distributions for animals collected in 2013 for our previous study (Stuckas et al., 2017). From this depiction it becomes apparent, that AHP-P includes roughly equal proportions of animals with very low (M. edulis – dominated genotypes) and high (M. trossulus – dominated genotypes) Q_MT-values, resulting in a median Q_MT - value around 0.6. Interestingly, this median value and distribution is more similar to that of populations further west of AHP sampled in 2013 (e.g., Warnemünde – WMU -, Gollwitz - GWZ -, Figure 1). However, experimental offspring (AHP-F1) from all three salinity treatments closely resemble Q_MT distribution of the 2013 AHP animals and populations further east (e.g., Barhöft – BAR -, Dranske at Rügen Isand - RUD -, USE), with a Q_MT-median in all three salinity groups >0.75, indicating a shift to more M. trossulus - dominated genotypes. This is also illustrated by the shift in violin shape in Figure 1B. Simulated offspring populations (AHP-F1, simulation 1–3) are characterised by a similar distribution of Q_MT-values as their parents (AHP-P) with a low median Q_MT-value of <0.25 (Figure 4A). Experimental AHP-F1 larvae raised to settlement at all three salinities were characterised by significantly higher Q_MT-values than the AHP-F1-simulation groups, with substantial enrichment of animals with Q_MT-values > 0.75 (ANOVA; F_(5,₄₇₎ = 12.44, p < 0.0001, Tukey HSD p < 0.01; Supplementary Tables 10, 11, Figure 4B). There were no significant differences in Q_MT - values between the different F1 salinity groups (Tukey HSD p > 0.05).

FIGURE 4

Figure 4. Distribution of hybrid classes expressed as Q_MT-values among parental animals (AHP-P), pooled simulated reference populations (AHP-F1 sim 1–3), and pooled treatment data at different salinities (AHP-F1_7psu, AHP-F1_9psu, AHP-F1_11psu). (A) Violin plots represent the distribution of individual Q_MT-values and medians (red diamonds). (B) Comparison between experimental (AHP-F1) and simulated offspring (AHP-F1-sim) populations depicting culture (replicate) bottle mean Q_MT-values (grey dots), with temperatures pooled for each salinity step, overall treatment mean values (red dot) and SD.

Discussion

Salinity Dependent Local Adaptation and Differential Survival in Baltic Mytilus Larvae

To our best knowledge, this is the first study in which Baltic Mytilus larvae were spawned and reared at extremely low salinities (<10 psu) under laboratory conditions. We explored whether larvae of Baltic Mytilus species from three Baltic model populations are adapted to local salinity regimes, and to test whether individuals with a high proportion of M. trossulus alleles have a higher tolerance of low salinity.

Our findings on pre-settlement stages of mussels representing different Baltic Mytilus species point towards local adaptation to comparatively high (>10 psu; Baltic M. edulis, KIE) and low salinities (<10 psu; Baltic M. trossulus, USE). This is primarily shown in the population comparison experiment where larvae originating from both localities show a better performance (mortality, growth, settlement success) under ambient compared to non-ambient conditions. Additional support comes from the desalination experiment, where transition zone mussels from AHP generally perform better at their ambient salinity (11 psu) compared to lowered salinities (9, 7 psu). Transition zone mussels show differential performance under low salinity conditions (7–11 psu) associated with a generally higher settlement success of larvae with a high proportion of M. trossulus specific alleles, as indicated by high Q_MT-values. In other words, hybrids with M. trossulus – dominated genotypes perform better under experimental low salinities. Our findings are in line with the existing knowledge on adult stages which suggests that salinity driven selection during the entire mussel life cycle can maintain genotypic and phenotypic differences between Baltic M. edulis and Baltic M. trossulus, despite of pervasive gene flow.

In a strict sense, the local adaptation observed in our experiments must be interpreted as a combination of genetically fixed phenotypic differences and acclimation and (transgenerational) epigenetic mechanisms (Palumbi et al., 2014; Donelson et al., 2018). Accordingly, future experimental approaches following multigenerational laboratory breeding of populations under common garden conditions (e.g., Sanford and Kelly, 2011) are required to help quantify the contribution of these factors to the observed performance trends, and to investigate the underlying mechanisms (e.g., Dupont et al., 2013; Shama et al., 2016).

Salinity driven selection in the laboratory shifted the genetic composition of experimental AHP-F1 animals toward the median hybrid index observed at low salinity (<11 psu) in eastern Baltic habitats, i.e., from median Q_MT – values of approx. 0.6 in AHP-P to values > 0.75 for experimental AHP-F1 resembling distributions observed for AHP, BAR, RUD, and USE (Figure 1B). The Q_MT – distribution of parental animals (AHP-P) sampled for our experiment in 2016 did not resemble that of the AHP reference population sampled in 2013 (Figure 1B; Stuckas et al., 2017). In fact, if we had started our experiment with a similar Q_MT distribution as in the reference sample from 2013, we most likely would not have detected a selection signal toward M. trossulus - dominated genotypes. The sampling bias in the AHP-P population toward a high number of M. edulis - dominated genotypes could be either due to random effects associated with the low sample size (12 parental animals) (Figure 1B). On the other hand, a shift in the position of the hybrid zone further eastward between 2013 and 2016 seems equally likely: the winter of 2014 experienced the third strongest wind-driven influx event of high saline North Sea water into the Baltic Sea in the last century, elevating salinity in the Darss Sill area (Mohrholz et al., 2015). This could have led to an increase in the fraction of M. edulis – dominated genotypes in the transition zone near AHP. AHP salinity monitoring data records indicate that salinity fluctuates strongly at this site and that salinity was relatively high in 2016 (Sanders et al., 2021).

Despite of the genetic shift toward high Q_MT - values observed in our experiment, differences between experimental AHP-F1 and naturally occurring populations still exist (Figure 1B), with a substantial number of M. edulis-like genotypes still present. This was observed in all experimental cohorts and the outcome was similar irrespective of the different salinity and temperature conditions. These observations allow forwarding two testable hypotheses: We suggests that below a salinity threshold of ca. 11 psu M. trossulus-like hybrids have a generally better performance than other hybrid classes. Presuming a threshold salinity as a trigger for differential performances explains why different experimental salinities basically resulted in similar genetic composition of settled larval cohorts. This explanation also fits to previous observations on naturally occurring populations east of Darss Sill (Figure 1B). Furthermore, given the effect of salinity driven selection on the pre-settlement phase, genetic differences between experimental AHP-F1 and natural occurring populations are best explained by the action of additional, post-settlement selection, as observed in a US East coast Mytilus population (Koehn et al., 1980).

Positive effects of increased temperature on growth and settlement success were observed in the simulated warming and desalination experiment - but also negative effects such as increased mortality. This may suggest that higher metabolic rates induced by increased temperatures may not fully buffer the negative effects of low salinity (Sprung, 1984; Sanders et al., 2018). Previous studies support the positive effect of higher temperatures on larval growth and developmental pace in Mytilus (Sprung, 1984; Beaumont et al., 2004; Sánchez-Lazo and Martínez-Pita, 2012). The warming treatment used in this experiment was within the lower part of the range of current interannual temperature variability in the natural habitat during the reproductive season (Pansch et al., 2018). Whether future warming will influence larval mussel survival and developmental energetics largely depends on whether spawning windows will shift toward earlier time points in the year. It is not unlikely that such phenological shifts will enable mussels to reproduce in thermal environments that are comparable to the ones they experience now. Strong phenological shifts with warming have been observed for Baltic plankton communities (Sommer and Lewandowska, 2011). More research on potential phenological shifts of mussel spawning windows under simulated warming is necessary to understand the combined impacts of desalination and warming on future larval mussel ecology in the Baltic Sea.

Overall, in combination with previous knowledge from post-settlement stages, this study has significance for predicting consequences for Baltic Mytilus abundance and distribution in the context of climate change. Given our experiments, predicted desalination and a westward shift in low salinities across the Baltic Sea have the potential to shift the M. edulis-M. trossulus distribution range and genetic transition zone westward (Meier et al., 2006; Gräwe et al., 2013). Even though elevated temperatures may benefit larval development by increasing developmental pace and reducing the length of the planktonic larvae phase, negative synergistic effects of high temperatures and extremely low salinities may be observed when temperatures exceed the species’ tolerance threshold.

What Mechanisms Underlie Salinity Driven Selection in Baltic Mytilus Species?

Our results suggest that M. trossulus – dominated hybrids are the outperforming genotypes under predicted future salinity conditions in the Baltic Sea. However, the putative mechanisms underlying these physiological differences between genotypes in the Baltic Sea remains elusive. Marine euryhaline molluscs are osmoconformers, meaning their extra- and intracellular environments remain iso-osmotic with surrounding seawater. Intracellular osmolarity is maintained primarily through adjustments in intracellular concentrations of compatible organic osmolytes (amino acids and amino acid derivatives), but also changes in inorganic ion concentrations (Berger and Kharazova, 1997; Yancey, 2005).

Under low salinity conditions, the depletion of intracellular ions, organic osmolytes, or both, to critical concentrations may impede cellular function (Podbielski et al., 2016). Thus, genotypic differences in salinity tolerance may depend on species-specific abilities to maintain physiological function at low intracellular osmotic pressures. Physiological costs associated with overcoming extremely low intracellular osmotic pressures such as increased active transport of inorganic ions (Maar et al., 2015) or reduced energetic efficiency of protein turnover (Tedengren and Kautsky, 1986) may underly the poorer performance of Baltic M. edulis genotypes at salinities <11 compared to Baltic M. trossulus. Biomineralisation processes may also play a key role in low salinity adaptation. In fact, reciprocal transplantation experiments demonstrated the reduced ability of adult Baltic M. edulis to calcify under low salinity conditions (Riisgård et al., 2014). Accordingly, salinities below 11 psu were shown to invoke higher costs of calcification in juvenile Baltic Mytilus with these costs potentially being related to limited availability of [Ca²⁺] (<4 mmol⋅kg^–1) below 11 psu (Sanders et al., 2018; Thomsen et al., 2018; 2021). Furthermore, a higher tolerance to lowered calcium concentrations has been observed in Baltic M. trossulus larvae when compared with Baltic M. edulis (Thomsen et al., 2018). This supports our observations of local adaptation to low salinity/[Ca²⁺] in this study and suggests that calcification substrate transport pathways may be under selection during adaptation to low salinity habitats. Interestingly, high genetic divergence in genomic regions associated with ion and osmoregulation have been identified in Baltic cod (6 psu) compared to North Sea cod (35 psu) suggesting these processes are targets of selection in low salinity environments, at least in vertebrates (Berg et al., 2015). Overall, it is likely that alleles at genetic loci associated with key physiological processes (e.g., osmoregulation, or biomineralisation) are targets of salinity driven selection and cause an indirect genome wide allele shift of species diagnostic markers. In fact, adaptive phenotypes related to physiological processes such as biomineralisation are polygenic traits (Hüning et al., 2016; Arivalagan et al., 2017) and complex selection regimes are expected to affect the entire genome.

In summary, the above cited example studies list putative mechanisms underlying salinity driven selection and show that there already is some evidence for species-specific responses to Baltic Sea salinity variation. Identifying genomic loci that are associated with salinity adaptation will enable us to better characterize the mechanistic basis for salinity adaptation in Baltic Mytilus species.

Conclusion and Future Perspectives

This study provides the first experimental evidence suggesting that (i) pre-settlement stages of Baltic Mytilus populations exhibit local adaptation to their native salinities (Baltic M. edulis, KIE; Baltic M. trossulus, USE) and that (ii) larvae with a high proportion of M. trossulus alleles (transition zone mussels, AHP) show a higher tolerance to low salinity conditions (11 psu and lower). Future work utilising complementary approaches and (multi-generational) laboratory experiments are key to elucidate ongoing selection and adaptation in this system. One task is to unravel the physiological mechanisms which modulate mussel performance under different environmental conditions. This is important to further test our prediction of a westward shift of the M. edulis - M. trossulus distribution range. Given the morphological and physiological differences between both species, these range shifts may have impacts on ecosystem functioning (nutrient and organic matter turnover) and biodiversity. A further task is to identify genomic loci causing salinity adaptation. In fact, while our genetic assay had a taxonomic focus (i.e., distinguishing between Baltic M. edulis, Baltic M. trossulus, and their hybrids), the molecular targets of salinity driven selection (genomic loci and their alleles) have not yet been identified. Overall, this study highlights the value of using the Baltic Sea environmental gradient as a model system to investigate the role of environmental factors in shaping physiological and genetic differences between populations and species (Berg et al., 2015; Guo et al., 2015; Reusch et al., 2018).

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below as: https://doi.org/10.1594/PANGAEA.892875.

Author Contributions

HS and FM conceived and designed the study. LK, JN-S, TS, DZ, and CH performed laboratory and field work. LK, JN-S, TS, DZ, FB, HS, and FM analyzed the data. HS and FM wrote the manuscript with support by LK, JN-S, and TS. All authors contributed to revisions.

Funding

Funding was provided by the Marie Curie ITN network “CACHE” (EU 7’th Framework Programme grant 605051) and the Paul Ungerer Stiftung. Technical funding and support was provided by the Kiel Marine Organism Culture Centre (KIMOCC) of the Kiel Cluster of Excellence “Future Ocean.” FB received financial support from the German Academic Exchange Service (DAAD) through the Doctoral Programmes 2015/16 (57129429).

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.

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.

Acknowledgments

The authors thank Ulrike Panknin, Julia Haegyeong Lee, Jörn Thomsen, and Anja Rauh for technical support. Earlier versions of this manuscript are part of two Ph.D. theses by TS (Sanders, 2018) and LK (Knöbel, 2019). The authors also thank three reviewers for valuable suggestions to improve the manuscript.

Supplementary Material

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

References

Arivalagan, J., Yarra, T., Marie, B., Sleight, V. A., Duvernois-Berthet, E., Clark, M. S., et al. (2017). Insights from the shell proteome: biomineralization to adaptation. Mol. Biol. Evol. 34, 66–77. doi: 10.1093/molbev/msw219

PubMed Abstract | CrossRef Full Text | Google Scholar

Beaumont, A. R., Turner, G., Wood, A. R., and Skibinski, D. O. F. (2004). Hybridisations between Mytilus edulis and Mytilus galloprovincialis and performance of pure species and hybrid veliger larvae at different temperatures. J. Exp. Mar. Biol. Ecol. 302, 177–188. doi: 10.1016/j.jembe.2003.10.009