Assessing Evolutionary Significant Units (ESU) of the Endangered Freshwater Pearl Mussel (Margaritifera margaritifera) in Southeast Norway on the Basis of Genetic Analysis

A total of 312 specimens of freshwater pearl mussel (Margaritifera margaritifera) were sampled from 11 populations, located in four different river systems in Southeast Norway, and analyzed for 11 simple sequence repeat (SSR) (microsatellite) markers. All study populations have landlocked brown trout (Salmo trutta) as the only possible host. Several populations had experienced recruitment failure, probably due to low pH (about 6.0) and calcium concentration. STRUCTURE clustering analysis revealed two genetic clusters, of which one cluster occurred mainly in the western river systems, and totally dominated in one population (Fallselva (A-FAL)) that had higher genetic diversity than the others. Cluster 2 completely dominated in the populations of the eastern river systems, and all of them had low genetic diversity. Bottleneck events were indicated in all populations and the inbreeding coefficient FIS was significant in all populations, except for the southernmost population (Sørkedalselva (B-SØR)), which was the only population with genotypes in Hardy–Weinberg equilibrium. FIS were especially high in the populations of the eastern river systems, and maximum shell length was negatively correlated to FIS. If artificially breeding and stocking should become necessary for future preservation, it should be based on single populations; alternatively, the eastern populations should be based on cross-breeding of populations within the cluster to increase their genetic diversity.


Introduction
Freshwater bivalvia is a group with several endangered species, and the freshwater pearl mussel (Margaritifera margaritifera) is one of the most threatened [1][2][3][4]. It is consequently on the European Red List [5]. The freshwater pearl mussel is distributed from the arctic and temperate regions of western Russia, through Europe to the northeast North America [6], and their preferred habitat is oligotrophic streams where they play an ecological role by filtering water [7][8][9][10]. In populated areas, mussel habitats are often subject to physical interventions, and in some areas they are affected by eutrophication and sedimentation [11], in other areas by acidification [12]. Drastic declines and extinctions have taken place in numerous locations [3,[12][13][14], and in addition to habitat deterioration, exploitation by humans for harvesting pearls has put the species at special risk [15].
In the management of biological resources, conservation of threatened species is an important issue, and extinction due to habitat loss may be a scenario with survival in captivity as the worst case alternative [16]. Several concepts are used in management and research, such as conservation of fish in southeast Norway [41], groups of populations may have different ancestry. Conservation measures such as semiartificial infections of host fish and supportive breeding must be based on local strains [28], and this study explored genetic structure to reveal the populations' ancestry, degree of isolation, and their genetic diversity to assess inbreeding and the populations' evolutionary potential. This was further meant to identify the local evolutionary significant unit (ESU) [19,28] of the species. Potential inbreeding effects on recruitment and shell length of populations were also considered.

Study Area
The sampling locations were chosen in known freshwater pearl mussel locations in four main river systems (the western A and B, and the eastern C and D) in Southeast Norway (Figure 1), representing the border of the species distribution in the inland district. In total, 312 specimens of freshwater pearl mussels were sampled from 11 different locations. The upper marine limit in the area is approximately 200 m above sea level (m a.s.l.), and the sampling sites are situated at 150 to 304 m a.s.l. pH ranged from 5.3 to 7.6 (the low values in spring) and [Ca] from 1.4 to 16.0 mg/L. With the exception of one location (River Hunnselva, C-HUN), all the streams exhibited low conductivity and low to moderate alkalinity (Table 1). Due to acidification, lakes upstream the following sampling locations: Lomsdalselva (A-LOM), Fallselva (A-FAL), Leira (C-LEI), Kampåa (C-KAM), Bråtåa (D-BRÅ), and Løvhaugsåa (D-LØV), were lime treated, i.e., limestone powder was added to increase pH, from the early1990s, but all treatments were terminated before 2013 [42].  Table 1. Physical and chemical variables at the sampling sites, minimum-maximum scale length (L) in our samples, and status score 1 of the mussels of the sampling streams. Labels A-D in the abbreviation indicate the main river system of the locations. Larsen and Magerøy [14] reviewed freshwater pearl mussel surveys from Norway, covering more than 400 populations including the populations in this study, and categorized the populations' statuses. This was based on six criteria of importance for the long-term survival of populations: (1) population size, (2) average density, (3) prevalence, (4) smallest shell, (5) proportion of shell lengths <20 mm (assumed to represent age groups <10 years), and (6) proportion of shells <50 mm (assumed to represent age groups <20 years). More details are given in Table S1. Populations are given 0-6 points of each criteria, and the points of each criterion are added for the population. Living populations may attain a total of 4 to 36 points, i.e., the former means that old specimens are present (0 means extinct), and the latter means six points of each criterion. The populations of the present study scored from 6 to 22 points, with a mean of 11.3 (Table 1), and four scored ≤ 7 and were characterized as threatened [14]. The study samples were not representative of the smallest specimens as they were avoided for DNA sampling.
The streams Lomsdalselva (A-LOM), Etna (A-ETN), and Fallselva (A-FAL) drain to the Lake Randsfjorden and further to the Drammenselva River (A) and the Oslo Fjord. There is one migration obstacle to fish downstream (130 m a.s.l.) of the Lake Randsfjorden and one obstacle (250-275 m a.s.l.) below the sampling site of A-FAL. All the A populations (the westernmost populations hereafter) had recruitment failure, and measures were conducted to improve recruitment by breeding larvae in artificial streams for stocking after four years [14]. The population in the Sørkedalselva River (B-SØR, the southernmost population), with an outlet to the northern part of the Oslofjorden, has fairly good recruitment [44], and there is a migration obstacle (140 m a.s.l.) downstream the sampling site.
In the large Glomma river system (C), relatively few locations harbor freshwater pearl mussels, and most of them are sparse with poor recruitment [44]. Samples were collected from four streams in the river system, mentioned from the west: Hunnselva (C-HUN), Leira (C-LEI), Kampåa (C-KAM), and Gjerda (C-GJE). There is a migration obstacle to fish downstream of the sampling site in C-LEI at 225 m a.s.l., and one obstacle in C-HUN at 150 m a.s.l., close to Lake Mjøsa (120 m a.s.l.). Natural recruitment was lacking in C-HUN, whereas there was weak recruitment in C-LEI and C-GJE, and fairly good recruitment in C-KAM, with >25% of the specimens <60 mm (<20 years), coinciding with a period of liming [45]. There are several migration obstacles to fish in the Glomma River downstream of the stream inlets of Leira and Kampåa, both natural waterfalls and power plant dams, and salmon are only present in the lowermost 17 km of the Glomma River, downstream of the waterfall Sarpsfossen.
The streams Finnsrudåa (D-FIN), Bråtåa (D-BRÅ), and Løvhaugsåa (D-LØV) (hereafter the easternmost populations) draining into the Lake Vänern, and the recruitment was characterized as good in D-FIN, whereas recruitment was weaker in D-BRÅ and D-LØV, and was probably positively affected by liming of upstream lakes in the 1990s [14,46,47].

Sampling and Analysis
From 2009 to 2011, haemolymph samples (0.1-0.3 mL) were collected from the foot with 1 mL syringes attached to 1.00 × 40 mm 21Gx2" sterican needles and stored in Eppendorf tubes with 400 µL RNALater (Ambion, Austin, TX, USA). We attained permission to sample mussels from the local County Environmental Administration (Fylkesmannen in Hedmark, reference 2016/311). Bivalvia is not included in the Norwegian Act of animal welfare, but the collected mussels were treated with care and kept in fresh water, except for 1 minute or so during sample collection, and put back at the sampling site. The use of a needle nevertheless represents a risk of injury, and during 2011-2016, samples were taken with the non-invasive mucus sampling with Q-tips (overlapping the first method in C-GJE) gently rolled on the foot of the mussel, and stored dry and cool in Eppendorf tubes until freezing (−20 • C) as soon as possible. Nuclear DNA was isolated with chelex [48], and 11 SSR markers were PCR-amplified in two separate multiplexes: (multiplex-1, Table S2) Mm2209, Mm2230, Mm2233, Mm2235, Mm2236, Mm2238 [49] and (multiplex-2) MarMa3116, MarMa3621, MarMa4277, MarMa4315, and MarMa5167 [50]. Data were analyzed and evaluated with the GeneMapper (ver. 4.0, Thermo Fisher Scientific, Waltham, MA, USA) software. Randomly chosen DNA samples (≈10%) were subjected to a second round of PCR and electrophoresis to assess the consistency of the assay. The mean scoring success across samples was 96.9% (± SD = 1.70), ranging from 94.5 to 99.7%. The primary results are given in Table S3.

Statistical Analysis
The software STRUCTURE 2.3 [51] was used to infer the most likely number of population clusters (K) constituting each sample, and the number of clusters was assumed to correspond to the number of ancestor populations [51]. Each individual was assigned a membership coefficient (Q) for each inferred cluster. Ten different runs were performed for each K (1-12, i.e., 1-n + 2) simulated, assuming an admixture model. A burn-in period of 50,000 iterations and a Monte Carlo Markov Chain (MCMC) of 50,000 iterations were used. The optimum number of clusters K was determined by means of the STRUCTURE-HARVESTER software [52], as described by Evanno et al. [53]. The estimated cluster membership coefficient matrices for the best fitted K was permuted so that all replicates have as close a match as possible using the CLUMPP software version 1.1.2. [54].
Genetic diversity indices including allele frequency, expected (H E ) and observed heterozygosity (H O ), number of alleles per locus (A L ), inbreeding coefficient (F IS ), and deviations from Hardy-Weinberg equilibrium (HWE) were explored by means of the web based GenePop software. Linkage disequilibrium and indication of loci under selection (with Markov chain parameters set at the maximum dememorization number and maximum number of iterations per batch (10,000) for 1000 batches) were explored by means of the ARLEQUIN 3.1 software [55]. Allele richness (A R ) was calculated (for sample size 20) by means of the HP-Rare 1.1 software [56], and relatedness was explored by means of the ML-Relate software [57]. The GenAIEx 6.5 Excel ad in [58] was used to calculate pairwise genetic differentiation F ST [59] and Nei's differentiation index (NeiD) [60]. A dendrogram-based hierarchical clustering of allele frequencies was constructed by means of the R-software package ape [61].
The BOTTLENECK 1.2.02 software [62] was executed using an infinite allele mutation model (IAM; the least appropriate for SSR data [63]), a stepwise mutation model (SMM), and a two-phase mutation model (TPM, assumed to be the most realistic model for SSR [64]). Populations exhibiting a significant number of loci with heterozygote excess by means of a Wilcoxon sign-rank test have likely undergone a recent population bottleneck event. A second method to reveal bottleneck events is calculation of the Garza-Williamson modified index [65] across loci by means of the ARLEQUIN 3.5.1 software. This index is a ratio M = the number of alleles (k) divided by the range in allele size (r), on the basis of the assumption that the number of alleles declines faster than the range in allele size during a bottleneck. Any dataset with seven loci or more with M < 0.68 can be assumed to have gone through a recent reduction in size, i.e., a bottleneck event [65].

Results
In total, 312 specimens of freshwater pearl mussels from 11 different streams were analyzed, and the number of alleles recorded for 11 SSR markers ranged from 15 to 39. No significant linkage disequilibrium nor indication of selection were observed at any locus after Bonferroni correction.

Genetic Structure
Pairwise differentiation, expressed as F ST and NeiD, was significant in all pairs, except for the C-LEI/C-KAM pair ( Table 2). The STRUCTURE-HARVESTER analysis suggested two main clusters (Figure 2), of which cluster 1 completely dominated in the isolated A-FAL population, and comprised 43-45% of the individuals of the other westernmost populations and 30% of the southernmost population (B-SØR) ( Figure 3). Cluster 2 comprised 96-99% of the individuals of the eastern populations. From the node-less phylogenetic tree based on allele frequencies, the A-FAL population stands out as an exclusive group, with the rest of the western populations being relatively closely interrelated (Figure 4). The eastern populations, appearing as the lower right ranch of the phylogenetic tree, appear as admixed between the river systems. The two samples C-LEI and C-KAM are closely related, as demonstrated by the non-significant F ST and NeiD index. C-HUN and D-BRÅ are also relatively close, although the sampling sites were the geographically most distant among the eastern populations. Likewise, the C-GJE, D-FIN, and D-LØV were relatively close. The genetic differentiations showed no clear relationship to geographic distance, except for the "neighbors" A-LOM and A-ETN. The A-FAL and the C-HUN, sampled from two streams with neighboring catchments but of different main river systems, were found to be substantially differentiated. Table 2. Pairwise genetic differentiation of the 11 study populations expressed as F ST (below diagonal) and unbiased Nei's differentiation index (NeiD, above diagonal)

Genetic Diversity
The genetic diversity was generally higher in the western compared to the eastern populations (Table 3, Figure 5). Of the western populations, 10 to 11 loci were polymorphic, while only four to seven loci were polymorphic in the eastern populations. The number of alleles per locus (A L ) and allele richness (A R ) ranged from 1.4 to 4.1 and 1.2 to 3.3, respectively, and A L was ≥2.4 in the western and ≤1.8 in the eastern populations, and correspondingly A R was ≥2.7 in the western and ≤1.7 in the eastern populations. Expected heterozygosity (H E ) was ≥0.586 in the western and ≤0.422 in the eastern populations, and the homozygote excess was significant in all samples, except for in B-SØR, which was in Hardy-Weinberg equilibrium (HWE). A-LOM and A-FAL had the highest numbers of private alleles, 11 and 9, respectively, whereas the B-SØR sample had five private alleles. In the samples from the eastern river systems, there were null to four private alleles across samples.   The inbreeding coefficient F IS estimates were significant in all samples, except for the B-SØR population, and ranged from 0.101 to 0.170 in the samples from the western river systems. F IS was considerably higher in the samples from the eastern river systems, ranging from 0.465 to 0.825, and the maximum shell length within samples was negatively correlated to F IS ( Figure 6, F 1,9 = 8.73, p < 0.05). There was no correlation between F IS and status category (Tables 1 and 3) (p > 0.05). The proportion of probable siblings was also considerably lower among the western than among the eastern populations (Table 4). Due to heterozygote deficiency in all populations, the method employed to detect recent bottleneck events as heterozygote excess using the BOTTLENECK software was not suited. Nevertheless, the mean Garza-Williamson modified index varied from 0.017 to 0.23 and was < 0.50 at all loci, except at two markers (Mm 2233 and MarMa 3621) in the A-LOM, A-FAL, and B-SØR samples (Table S3), i.e., recent bottleneck events were indicated in all populations.

Genetic Structure and Probable Origin of the Study Populations
STRUCTURE analysis of 312 individuals of freshwater pearl mussels from four main river systems suggested two clusters, of which cluster 1 totally dominated in one of the westernmost populations, the isolated and elevated A-FAL population (in mean 98% of the individual genomes). Cluster 1 also made up 30 to 45% of individuals of the other western populations, whereas the eastern populations were completely dominated by cluster 2, comprising 96 to 99% of the individual genomes. This corresponded to the node-less phylogenetic tree, with the A-FAL population as an outlier, and the other westernmost populations, located near the Lake Randsfjorden, being more closely related to each other and relatively closely related to the southernmost (B-SØR) population. The eastern populations were differentiated from the western populations, whereas they were closely related to each other, and were admixed between the eastern river systems. The C-LEI and C-KAM populations were especially closely related, and the genetic differentiation, expressed both as the F ST and NeiD indices, was non-significant, which was different from all the other pairs tested. This leads to a suspicion that at least one of these populations was stocked by humans in relatively recent time, as pronounced structuring is commonly found among populations of landlocked freshwater pearl mussels [27].
Natural occurrence of freshwater pearl mussels presupposes a past immigration of salmon or brown trout as hosts, and in the western populations (A and B), the freshwater pearl mussel was probably hosted by anadromous salmonids in locations connected to the sea. The invasion of freshwater fish to the Glomma river system may have followed two different routes. The upper part of the Glomma River once drained to Lake Vänern in Sweden, via the Vrangselva River, and was colonized by freshwater organisms originating from the ancient freshwater Lake Ancylus [41], which once covered the Baltic Sea and parts of Sweden [66]. Brown trout, potentially hosting mussel larvae, invaded the Glomma River through this route, until the river broke through land barriers to the west and closed the connection to Lake Vänern. Since then, the Glomma River has run to the eastern side of the Oslo Fjord, and immigrants originated from the freshwater Lake Ancylus/Vänern were spread to lower parts of the Glomma River, and could also enter the Lake Mjøsa and its tributaries. The distribution of several fish species supports this invasion theory, with one example being the Arctic grayling (Thymallus thymallus L.), which has a low tolerance to salinity [67] that prevents it from entering river systems by sea migration, being present only in the Glomma river system and in systems farther east. Other species with similar distribution patterns are alpine bullhead (Cottus poecilopus) and burbot (Lota lota).
The differentiation between the western and the eastern populations suggests that the populations originated from different sources, at least in part. The mussels in the western river systems were probably introduced by anadromous brown trout or salmon at some stage after the glaciation, when the locations were available for ascending anadromous fish [41]. In the eastern river systems, the species may have been introduced by landlocked brown trout and salmon (still present in Lake Vänern), originating from the ancient freshwater Lake Ancylus; however, at some stage of the land rising, all locations below 200 m a.s.l. were physically available for anadromous fish. This leaves a possibility for the anadromous origin of salmonids hosting mussels in all river systems. Nevertheless, it remains an unanswered question as to when the invasion of mussels took place. The presence of cluster 2 in the populations A-LOM, A-ETN, and B-SØR suggested a partly common origin of three of the western populations and the eastern populations, whereas the A-FAL population, as well as cluster 1, likely originated from an earlier invasion.

Genetic Diversity
The genetic diversity of the study populations was substantially higher in the western (A R = 2.59-3.23 and H E = 0.321-0.524) than in the eastern (A R = 1.21-1. 64 and H E = 0.037-0.139) populations, and only one population, B-SØR, had genotypes in HWE and non-significant F IS . On the basis of comparative studies in central and southern Europe, researchers generally characterized A R < 2 and H E < 0.1 as low [27,28], and except for two populations with heterozygosity > 0.1, all the eastern populations fall into the category of low diversity. Karlsson and Larsen [68] found A R = 2.5 and H E = 0.377 in Fallåa (A-FAL) (on the basis of different SSR markers), and a substantially lower A R = 1.2 and H E = 0.065 in Hunnselva (C-HUN), confirming a low genetic diversity of mussels in the Glomma drainage.
Mussel populations depending on landlocked hosts have low genetic diversity and are found to be strongly structured [27], as was the case in this study. Bottleneck events, suggested for all the populations, result in loss of low frequency alleles, and subsequent genetic drift leads to differentiation between populations, and in addition, low population density leads to more frequent self-fertilizing and thereby reduced heterozygosity [27]. This may explain the low diversity of the landlocked eastern populations.
In a study on freshwater pearl mussels in river systems with anadromous brown trout and salmon as host species in northern Sweden, A R was > 3 and H e > 0.5 in 8 of 14 samples, and the genetic diversity increased downstream in the rivers [69]. If the mussel populations (C and D) of the Glomma and the Vrangselva river systems have a common immigration history with origin in tributaries to Lake Ancylus or Lake Vänern, D-FIN would be geographically closest to the source population, potentially causing the higher number of polymorphic loci in this population. Further, there were three private alleles and A R = 1.64, compared to two or less private alleles and A R ranging from 1.21 to 1.44 in the other eastern populations.
The isolated westernmost A-FAL population was special, being characterized with low abundance and poor recruitment [14] but with exceptionally high genetic diversity, suggesting a former large and well-reproducing parental population. The population also had low F IS and a low proportion of siblings as compared with the eastern populations.
The negative correlation of maximum shell length and F IS suggests an effect of inbreeding depression in the eastern populations, corresponding to what Zheng et al. [31] experienced with American bay scallops in breeding experiments. Nevertheless, the inbred C-LEI and C-KAM populations (F IS ≥ 0.448, H O ≤ 0.067) were characterized as abundant and well-reproducing, with the second and third highest status score (15 and 18) of all the eastern populations [14,45]. The C-KAM population included thousands of individuals [47], and recruitment had taken place during the last 20 years, coinciding with lime treatment of an acidified upstream lake. Likewise, the D-LØV population, with low genetic diversity, had successful reproduction in the 1990s when upstream lakes were lime treated [47], suggesting that reproduction was hampered by acidic water.
Suboptimal water quality in several mussel locations, and apparently positive effects of limestone powder on mussel recruitment in some of them, suggested that liming may be necessary to secure the freshwater pearl mussels' existence in some of the study locations [45]. When pH is about 6.0 and episodically lower and with [Ca] < 2.0 mg/L, it probably affects freshwater pearl mussel recruitment negatively. This interpretation is supported by the results of a liming project in a stream of similar quality, bringing pH above 6.2 and [Ca] above 2.5 mg/L, resulting in recruitment recovery of a freshwater pearl mussel population in Southwest Norway [70].
No negative effects of inbreeding on recruitment or viability were indicated when comparing genetics and recruitment of the study populations. It may be hypothesized that the general low genetic diversity of the eastern populations is due to an origin from settlers hosted by landlocked salmonids, different from the western populations descending from mussels hosted by anadromous host fish. There was no indication of selection on the analyzed markers, but SSRs are supposed to be neutral and not subject to selection, although there are exceptions [71]. Inbreeding in populations exposed to selection pressure may result in purging of harmful genes, favoring those conferring beneficial traits [72]. The present genotypes of the eastern populations may therefore in part be a result of selection caused by environmental stress, resulting in purging of genes in former generations, and efforts to preserve the remaining genetic diversity are urgent. Each population, with a possible exception of the C-LEI and C-KAM, should be regarded as unique ESUs, until eventually inbreeding should appear to affect the population viability. Then the eastern populations could be pooled to increase genetic diversity by experimentally admixing the populations in hatcheries. The viability of the admixed offspring should be tested by stocking in the eastern river systems in locations with suitable physical and chemical environments, but at present are not inhabited by the mussel. If this succeeds, an admixed population may be kept in captivity as a reserve stock.

Conclusions and Recommendations
Recruitment failure, low genetic diversity, and a high portion of siblings in the majority of the study populations suggested that they are prone to further loss of genetic variation and may be threatened by extinction in the future, i.e., within 100 years, except for the B-SØR, which appeared to be a healthy population. Habitat conservation is necessary, and in several study locations, the water quality, expressed as pH and [Ca 2+ ], is suboptimal for pH-sensitive organisms such as freshwater pearl mussels. Water quality and mussel reproduction in the study streams should be monitored, and lime treatment should be considered.
Excess of homozygotes and significant inbreeding coefficients suggested inbreeding, with a possible risk of inbreeding depression and future negative effects on reproductive performance [32]. Restoration measures must be based on the local populations of freshwater pearl mussels and host brown trout as the infection rate of mussel larvae depends on genetic predisposal [73,74]. Some of the eastern populations, all representing cluster 2, could be pooled if the low genetic variation and inbreeding appears to hamper recruitment and viability.

Funding:
This research received external funding from Hedmark and Oppland Counties' Environmental Administrations.