Sewage treatment plant associated genetic differentiation in the blue mussel from the Baltic Sea and Swedish west coast

Sampling

Mussels were sampled during June 2012 and 2013 from five locations in the Baltic Proper (BP) and off the Swedish west coast (WC): Trosa/Nynäshamn Archipelago (east coast of Sweden) (ASK), Tvärminne in the Gulf of Finland (TVA), Karlskrona (Southern Sweden) (KAR), the Gulf of Gdańsk, Poland (GDA) and Kristineberg (the Swedish west coast, Skagerrak) (KRI) (Fig. 1, Table 1). Within each location three sites representing different pollution types were chosen based on information from reports and local authorities’ assessment schemes. These sites were located close to a sewage treatment plant (STP) effluent, the vicinity of a harbor (HAR), and a reference (REF). The sewage treatment plants studied include large plants with >100,000 persons connected to the plant (ASK_STP and GDA_STP), a medium plant >40,000 persons (KAR_STP) and two small plants each with approximately 7,000 persons connected (TVA_STP and KRI_STP), but all STP sites had similar amounts of nitrogen and phosphorus emission (see Table S1 for details). The sampled vicinities of harbors included different types of harbors; ferry ports (ASK_HAR, TVA_HAR), cargo and ferry Gdynia harbor (GDA_HAR), a refinery port at Preemraff (KRI_HAR), and a naval port (KAR_HAR) (see Table S2 for details). STP and HAR sites were paired with a reference site (REF) that was not influenced by sewage effluents or harbor associated pollutants. All REF sites were situated within 5–30 km from the STP and HAR sites (Tables 1 and 2). The geographical distances between sites (Table 2) were estimated using the most direct marine route in Google Earth version 7.1.2.2014 (Google, Mountain View, CA, USA). Temperature regimes and salinity differed between the different locations (west coast, and southern and northern Baltic Proper) but were similar among the three sampled sites within each of the locations (Table 1).

At each of 15 sites, 30 blue mussels (24–30 individuals used in the subsequent AFLP analyses, Table 1) with a mix of both sexes were sampled. The mussels were sampled using either a benthic sledge, a triangular bottom scraper or by hand, from a depth of between 1 and 13 m depending on local distribution. The age of each individual was estimated by counting growth rings (see Haskin, 1954) and all individuals included in this study had an estimated age of 2–5 years. From each individual the adductor muscle was dissected, snap frozen and immediately stored in −80 °C, prior to DNA isolation.

DNA isolation

Total genomic DNA was isolated from a small piece of the muscle tissue ∼2 × 2 mm of each individual using the E.Z.N.A. Mollusc DNA Kit (OMEGA Bio-Tek, Norcross, GA, USA), with small modification of the manufacturers protocol. In short, tissue was incubated overnight at 56 °C in 250 µl ML1 lysis buffer and 25 µl proteinase K followed by a chloroform:isoamyl (24:1, 300 µl) extraction step with 3 min centrifugation (at 14,000 rpm). DNA suspended in the aqueous phase was carefully separated from the interphase (200-300 µl) and 300 µl MBL buffer with 10 µl RNAse was added, and samples were incubated for 10 min at 70 °C. The DNA was precipitated from the solution with 99.7% ethanol, transferred to DNA HiBind spin columns and centrifuged to remove any traces of alcohol. Purified DNA was eluted with 50 µl 10 mM Tris–HCl buffer (pH 8.5) (preheated to 70 °C). DNA quality was visually estimated using agarose gel 1.5% electrophoresis (100 V, 20–30 min) stained with (5 µl/100 ml, 30 min) SYBRSafe (Life technologies).

Introgression analyses

Due to the complex system of occurrence of M. edulis and M. trossulus in the Baltic Sea and the west coast, and extensive mixing of the M. edulis and M. trossulus genomes (Väinölä & Strelkov, 2011) in the Baltic Sea, we used the Glu-5′ (ME 15/16) (Inoue et al., 1995) species marker to evaluate any remaining species identity effects on the genetic structure. This marker is commonly used and distinguishes between all three Mytilus taxa found in the studied area (Kijewski et al., 2006). Individuals were genotyped at locus Glu-5′ (ME 15/16) (Inoue et al., 1995), and each individual was coded as homozygous M. edulis, M. trossulus, M. galloprovincialis or as heterozygotes between them. The AFLP method is suitable for analyses of closely related species but homoplasy, the non-homology of bands with the same electrophoretic properties, is expected to be more frequent in hybrid populations (Egger et al., 2007). Occurrence of homoplasy will have four main consequences for population genetic inferences: (1) an overestimation of the frequency of the presence allele; (2) an underestimation of the differentiation between populations; (3) an overestimation or underestimation of the heterozygosis, depending on marker frequency; and (4) reduction in the power to detect selective loci (Caballero, Quesada & Rolán-Alvarez, 2008).

The Glu-5′ PCR fragments were amplified in 10 µl reactions containing: 20–40 ng DNA, 25 mM MgCl₂, 1 µl 10X PCR-buffer, 5 U AmpliTaq (Applied Biosystems), 1 µl BSA, 10 mM dNTP (Thermo Scientific Molecular Biology), and100 µM of each primer. Forward primer (ME15) 5′-CCAGTATACAAACCTGTGAAGA-3′ 5′-end was labeled with fluorescein and reverse (ME16) 5′-TGTTGTCTTAATAGGTTTGTAAGA-3′ (MWG Biotech). The amplification cycle consisted of an initial denaturation for 5 min at 94 °C followed by 38 cycles at 94 °C for 30 s, 55 °C for 30 s and 72 °C for 45 s, followed by elongation for 10 min at 72 °C. DNA fragments were separated on an ABI-3730XL capillary electrophoreses unit at Uppsala Genome Center with separation medium POP7™ Polymer (Applied Biosystems), size standard GeneScan™ 500 ROX™ (Applied Biosystems), injection time 15 s (1.6kV), run time 1,600 s and array length 50 cm.

Generation of AFLP markers

Markers were generated as described by Vos et al. (1995) with minor modifications as described by Bensch et al. (2002). Briefly, genomic DNA (20–40 ng/ml) was digested with restriction enzymes for 1h at 37 °C. Each 20 µl reaction contained 0.5 µl of 50U EcoR1 (Thermo Scientific Molecular Biology) and 0.5 µl of 50U Tru1 (Thermo Scientific Molecular Biology), 2 µl 10x TA-buffer (Thermo Scientific Molecular Biology), 1 µg BSA and ddH₂O. Followed by ligation for 3 h at 37 °C, in; ligation buffer (10X), 100 mM E-adaptor 5′-CTCGTAGACTGCGTACC-3′, 3′-CATCTGACGCATGGTTAA-5′ (MWG Biotech), 100 mM M-adaptor 5′-GACGATGAGTCCTGAG-3′, 3′-TACTCAGGACTCAT-5′  (MWG Biotech) and 10 U/ml T4 ligase (Thermo Scientific Molecular Biology).

Pre-amplification was performed in 20 µl reactions containing 10 µl DNA template (cut, ligated and diluted), 100 µM E-primer 5′-GACTGCGTACCAATTCA-3′ and 100 µM M-primer 5′-GATGAGTCCTGAGTAAC-3′, 25 mM MgCl₂, 2 µl PCR-buffer (10x), 0.8 µl BSA, 1 mM dNTP, 5U Taq polymerase and ddH₂O. Pre-amplification PCR conditions included an initial denaturation step for 2 min at 94 °C followed by 20 cycles of 94 °C for 30 s, 56 °C for 30 s and 72 °C for 60 s, and a terminal step at 72 °C for 10 min.

For selective amplification eleven primer combinations were evaluated using a subset of individuals from all the geographical regions and habitats. The two chosen combinations MCTG_EAGC and MCTT_EACT, showed polymorphic bands, with a congruent pattern in duplicated individuals. The PCR reaction for selective amplification was conducted containing; 2.5 µl of the diluted pre-PCR template, 25 mM MgCl₂, 1 µl 10x PCR-buffer, 1 µl BSA, 1 mM dNTP, 100 mM E-primer 5′-GACTGCGTACCAATTCNNN-3′ with a fluorescein 5′-end, 100 mM M-primer 5′-GATGAGTCCTGAGTAANNN-3′ and 5U Taq polymerase and ddH₂O to a volume of 10 µl. The selective amplification cycle consisted of denaturation for 2 min at 94 °C, 12 cycles at 94 °C for 30 s, 56 °C with a touch-down temperature −0.7 °C for 30 s and 72 °C for 60 s, followed by 23 cycles at 94 °C for 30 s, 56 °C for 30 s, 72 °C for 60 s and a terminal step at 72 °C for 10 min.

The DNA fragments (AFLP) were separated on an ABI-3730XL capillary electrophoresis unit at Uppsala Genome Center with separation medium POP7™ Polymer (Applied Biosystems), size standard GeneScan™ 500 ROX™ (Applied Biosystems), injection time 15 s (1.6kV), run time 1,600 s and array length 50 cm.

AFLP analyses

To generate a genotype matrix for each primer combination, AFLP data from the two primer combinations were analyzed separately in Genemapper 4.0 (Applied Biosystems). Analyses settings: primer MCTG_EAGC: 250–500 bp, (Primer MCTT_EACT: 230–500 bp), bin width 1.5 bp, peak detector 170 rfu, normalization using sum of signals from all individuals within one primer combination and no smoothing. For error rate analysis 7–14% of the samples were randomly chosen and duplicated from extraction and/or ligation step. The duplicates were manually checked for congruency and the recorded duplicated loci where used to calculate error rate i.e., the number of inconsistent loci per primer combination divided by total number of loci resulting in the percentage of error. The two genotype-matrices were merged into one dataset used in all subsequent analyses.

Introgression and species identity analyses

To score the target allele length at the Glu-5′ locus, GeneMarker 2.5.2 (SoftGenetics LCC) software with size standard GS500 and analysis type Fragment (Animal) was used. In total, 364 individuals were scored, by the allele lengths as M. galloprovincialis (allele length 124 bp), M. trossulus (165 bp) or M. edulis (177 bp) Twenty percent of the samples were duplicated to check for congruency.

Statistical analyses

AFLP

A total of 396 individuals were used in the analyses and the two primer combinations yielded 354 polymorphic loci (MCTG_EAGC: 195 loci, MCTT_EACT: 159 loci). For error rate analyses 14% of duplicated samples for MCTG_EAGC and 7% of duplicated samples for MCTT_EACT were randomly chosen and yielded a congruent pattern and an error rate of approximately 6% and 7% respectively.

Gene diversity and genetic differentiation

Gene diversity (H_E) for each site and mean gene diversity for each pollution type is presented in Table 1. There was a significant overall partitioning of genetic diversity (F_ST = 0.027, P-value <0.01). No significant difference in mean gene diversity between the pollution types was found; including the west coast sites (mean_REF = 0.131, mean_STP = 0.130, mean_HAR = 0.124, ANOVA: df = 2, F = 0.933, P-value = 0.42) and excluding the west coast sites mean_REF = 0.128, mean_STP = 0.130, mean_HAR = 0.124, ANOVA: df = 2, F = 0.387, P-value = 0.690). Significant F_ST pairwise differentiation (after FDR correction, FDR = 0.05) was found between the Baltic Proper sites and the west coast sites, and between two pairs of STP and REF sites (ASK_STP/GDA_REF and ASK_STP/KAR_REF) and between one pair of HAR sites (ASK_HAR/KAR_HAR), within the Baltic Proper (Table 2).

When sites were pooled within pollution types, significant genetic differentiation (F_ST) between the pollution types STP and REF was observed, both including and excluding the west coast sites (Table 3). No significant genetic differentiation was found between the pollution types HAR and REF either including or excluding the west coast sites (Table 3).

The comparison of genetic differentiation (F_ST) between the five locations, ASK, TVA, KAR, GDA and KRI (i.e., 10 pairwise comparisons) for each of the three pollution types showed a significant difference after FDR corrections (FDR = 0.05) between all the Baltic Proper locations and KRI and one significant comparison within the Baltic Proper i.e., between ASK and KAR for HAR (Table 3).

Introgression and species identity

In total, 364 individuals were scored as M. galloprovincialis (allele length 124 bp), M. trossulus (165 bp) or M. edulis (177 bp). Twenty percent of the samples were duplicates and all of the duplicated pairs of samples had the same allele combination. Three individuals contained M. galloprovincialis alleles, all in heterozygotes, two from KRI_REF and one from KRI_HAR, but since our analyses are based on variation between and within groups, a group of three is far too small to give reliable results and these three individuals were excluded from subsequent analyses. The remaining individuals were assigned a species identity coded as homozygous M. trossulus (2), homozygous M. edulis (0) or heterozygotes between the two (1) and species identity was used in subsequent analyses to explore and account for any differences in species identity (Table 1).

To test the species identity effect on the overall genetic structure a constrained principal coordinate analysis (cPCoA) was used, using both the whole data set and a data set excluding the west coast sites. A significant effect of species identity on the genetic structure was found when including the west coast sites (df = 1, F = 3.90, P-value = 0.001), but not when analyzing only the Baltic Proper sites (cPCoA: df = 1, F = 0.90, p-value = 0.739). These results showed a general species identity effect between the west coast and the Baltic Proper, but not within the Baltic Proper.

Genetic structure

Multivariate constrained principal coordinate analyses (cPCoA) were used to assess the significance of the explanatory variables (constraining variables) location, pollution type and species identity. cPCoA analyses including order dependent, permutation based ANOVAs were conducted on different aspects on the dataset. The cPCoA including all STP and REF sites, using location, species identity and pollution type as constraining variables, was tested with the order dependent ANOVA and showed a significant structuring of the genetic variation both for location and pollution type after FDR correction but not for species identity or the interaction between species identity and pollution type (Table 4). The cPCoA, including STP and REF sites from the Baltic Proper only, using location, species identity and pollution type as constraining variables, showed a significant structuring of the genetic variation after FDR correction for pollution type but not for location, species identity or the interaction between species identity and pollution type (Table 4). The cPCoA including only sites from the west coast (i.e., KRI only) using pollution type as constraint, showed no significant effect of pollution type STP/REF (Table 4).

The conditioned cPCoA, using pollution type, STP and REF as constraint and location and species identity as conditions also showed a significant effect of pollution type on the genetic structure in the Baltic Proper and was used to visualize the effect of pollution type (Table 4, Fig. 2).

In contrast to the STP and REF comparisons, the cPCoA for the HAR and REF data showed no genetic differentiation between the pollution types. When including all sites (both west coast and Baltic Proper) a significant effect (after FDR correction) of location was found, but not for species identity, pollution type or the interaction between the two (Table 4). No effects could be found when analyzing the HAR/REF contrast using the sites from the Baltic Proper (Table 4). The analyses including only sites from the west coast (i.e., KRI only) using pollution type and species identity as constraining variables, showed no significant effects (Table 4). These results indicated that the location effect found when including all sites (both west coast and Baltic Proper) was based on the differentiation between the west coast and Baltic Proper.

To compare the genetic differentiation between the locations for each pollution type independently, three additional cPCoAs were conducted using only sites from Baltic Proper. No significant genetic differentiation between the REF sites (df = 3 F = 1.05, P-value = 0.274) or the STP sites (df = 3, F = 1.14, P-value = 0.080) could be found, while there was a significant difference between the HAR sites (df = 4, F = 1.23, P-value = 0.010).

Search for specific loci associated with differentiation between the STP and REF

To identify specific marker loci associated with the genetic differentiation between the STP and REF sites four different approaches were used, varSelRF, DFDIST, BayeScan and Cochran-Mantel- Haenszel test. DFDIST and varSelRF were used for both an overall comparison between STP and REF and for pairwise site comparisons. The overall comparison using VarSelRF identified 20 loci and DFDIST 12 significant loci (Table S3). The within location based pairwise comparisons in varSelRF identified between 3 (ASK_STP/REF) to 5 (GDA_STP/REF) loci and DFDIST between 3 (GDA_STP/REF) to 17 (ASK_STP/REF) significant loci (Table S3). In the BayeScan analysis none of the loci had a posterior probability above 0.17, indicating a very low Bayes Factor. However, many of the loci with the highest Bayes Factor (yet not significant) corresponded to loci identified using the other two methods (Table S3). In the Cochran-Mantel- Haenszel test, initially 16 loci were identified as significant, but none of them remained significant after FDR correction. Also here, many of the identified loci (but not significant) corresponded to loci identified using the other methods (Table S3).

Genetic differentiation in blue mussels between reference sites and sewage treatment plant sites in the Baltic Proper

The major genetic features of blue mussels in the sampled area were the expected pronounced genetic differentiation between the west coast and the Baltic Sea, and the weak differentiation within the Baltic Proper caused by strong gene flow. Genetic differentiation in the blue mussel between the west coast and the Baltic Sea has previously been identified in several studies (Johannesson & Andre, 2006; Wennerström et al., 2013; Larsson et al., in press). In the present study, the pairwise F_ST analyses (Table 2) confirm this geographical pattern. The present study also identified a higher frequency of M. edulis alleles in the west coast sites (Table 1), and an overall significant effect of species identity on the genetic structure, however the effect of species identity on the genetic structure was not significant when geographical location was accounted for in the analyses, indicating that differences other than species origin might be present similar to what was found in Larsson et al. (in press). Within the Baltic Proper no effect of species identity on the genetic structure was found, also similar to what was found in Larsson et al. (in press). The current study shows a weak over-all genetic differentiation within the Baltic Proper as shown by very few significant pairwise F_ST values (Table 2). When analyzing genetic differentiation between locations within each pollution type separately, no significant differentiation between locations was found within the pollution types STP and REF as shown both by the location based F_ST analyses (Table 3) and by the cPCoAs. These results confirm the earlier findings of low differentiation between locations in the Baltic Proper, due to high oceanographic connectivity and strong gene flow (see Larsson et al., in press).

Despite the low general genetic differentiation within the Baltic Proper a significant differentiation between pollution type STP and REF was detected. This was shown by several analyses such as the cPCoA (Table 4, Fig. 2), the pollution type based F_ST analyses, and in two significant pairwise F_ST analyses (Table 2). This effect seems to be restricted to the Baltic Proper as no significant difference between STP and REF on the west coast was found (Table 4). Sewage effluents have also been shown to be associated with population genetic effects in other systems. Puritz & Toonen (2011) detected that genetic structure of the bat star (P. miniata) correlated with major sources of sewage effluents and storm water. The bat star has a similar life history to the blue mussel with large populations and strong gene flow caused by larval dispersal over large distances. Despite this strong gene flow Puritz & Toonen (2011) were able to show that differences between sites (i.e., pairwise F_ST) could be explained by sewage water causing a barrier to larval dispersal.

The significant differentiation between STP and REF sites despite strong gene flow is in line with studies finding significant genetic differentiation among populations when subjected to high gene flow (Nosil, Funk & Ortiz-Barrientos, 2009; Lamichhaney et al., 2012; Gagnaire et al., 2012; Ulrik et al., 2014; Laporte et al., 2016; Guo, Li & Merilä, 2016). Taking into account the generally low genetic differentiation within the Baltic Proper and the nested study design, evolutionary forces like restricted gene flow, directional gene flow and genetic drift are considered less likely to have caused the observed genetic differentiation between STP and REF. The most likely evolutionary force acting on these populations is selection. As we do not have any data on local recruitment over generations, we conservatively describe the observed pattern as a within-generation footprints of local directional selection, in parallel with the findings in the American eel (A. anguilla) and European eel (A. rostrata) (Gagnaire et al., 2012; Ulrik et al., 2014; Laporte et al., 2016).

The AFLP method, employed in the present study, has previously been shown to identify specific marker loci linked to adaptive variation among habitats (Williams & Oleksiak, 2008; Nosil, Funk & Ortiz-Barrientos, 2009; Lind & Grahn, 2011). This possibility was explored by using four different analyses. The combination of analyses identified a large number of loci (Table S3), however only a few loci were verified across methods and combinations of data set (i.e., pair-wise and pooled STP and REF). No significant loci were identified using BayeScan and loci identified using CMH were not significant after FDR correction. The identification of only a few specific loci associated with the differentiation may be due to small allelic changes among many co-varying loci rather than pronounced changes at a few loci of large effects (Laporte et al., 2016). To further explore such polygenic selection, a high resolution sequencing method (e.g., RAD sequencing or genome sequencing) in combination with approaches not based on the hitch-hiking model would be needed (Laporte et al., 2016).

The anonymous nature of AFLP and the suggested polygenic nature of STP effects make it difficult to imply any direct mechanisms behind the observed differentiation. Cage experiments exposing adult mussels outside sewage treatment plants have shown to include several direct (pharmaceuticals and endocrine disrupting compounds) and indirect stressors (nutrients load-induced algae blooms and/or temporal anoxia/hypoxia), likely reflecting the complexity of sewage effluents. The stressors have included immunological (Akaishi et al., 2007), genotoxic and lysosomal responses (Turja et al., 2015) but also endocrine disruptive effects such as feminization in freshwater mussels (Elliptio complanata) (Gagné et al., 2011) and up-regulation of vitellogenin in males in both freshwater mussels Elliptio complanata (Bouchard et al., 2009; Gagné et al., 2011), and in marine mussels Mytilus galloprovincialis (De los Ríos et al., 2013). The experiments described above were conducted on transplanted adult individuals (De los Ríos et al., 2013; Turja et al., 2015) and in this study the mussels examined were likely exposed to pollutants throughout their life, and it is therefore difficult to imply at which life stages of the blue mussels the putatively local selection occurs. Toxicity studies have shown that larvae are more sensitive towards many pollutants when compared to adults (Beiras & Bellas, 2008), and as Puritz & Toonen (2011) suggest that waste and storm water act as barrier to larval dispersal, it is therefore plausible that STP effluent effects occur at early life stages.

As the blue mussels in the Baltic Proper show a strong introgression pattern (Väinölä & Strelkov, 2011) the effect of species identity on the genome wide genetic structure was explored, while possible effects of species differences was controlled for when analyzing pollution type effects. In this study, species identity, based on the Glu-5′ marker, was not associated with genetic structure in any of the analyses when the effect of large-scale geographical distances were accounted for. The lack of association between species identity and genetic structure in the Baltic Proper supports the presence of strong introgression in the Baltic Proper (Larsson et al., in press). These results also show that species identity is not an important part of the pollution type effects on genetic differentiation. If species identity was an important factor in the STP effects on genetic structure, a significant interaction between species and STP would have been observed. It is, however, possible that the hybridization status of each individual was not fully covered by the Glu’-5 marker and that other markers based on other parts of the genome may give a discrepant pattern (see Kijewski et al., 2006; Kijewski et al., 2011). To study introgression patterns in detail, a more comprehensive full genome scan including pure M. edulis and M. trossulus as reference individuals is needed but this was beyond the scope of the present study.

Genetic divergence between harbor sites (HAR) in the Baltic Proper

In contrast to the genetic differentiation found between STP and REF sites, we did not find any results suggesting a congruent effect on the genetic differentiation between the HAR and REF sites. However, the cPCoA show significant differentiation among HAR sites in the Baltic Proper. This pattern may be interpreted as a signature of divergent local selection among HAR sites and is in contrast to a study of Mediterranean mussels where no significant genetic differentiation between harbor sites was found (Štambuk et al., 2013). A possible explanation may be that the harbors each have a diverse pollution regime, and the selection pressure therefore differs between them. The results indicate that there may be differences in the pollution regimes between the harbor types (naval port, oil harbor and ferry ports). Another factor that could affect differentiation between harbor populations is shipping. All studied harbors are connected to inter- and intra- Baltic shipping routes, so it is possible that harbor populations may be differentiated via introductions of individuals from different areas in- and outside the Baltic Sea.

Similar to what was found for the STP pollution type, no significant effects of species identity were found, and again this result indicates that species identity is not an important part of the pollution type effect on genetic differentiation.