Causes and consequences of fine-scale population structure in a critically endangered freshwater seal

Laboratory analyses

We used tissue specimens from Saimaa ringed seals found dead and deposited in a tissue bank maintained by the University of Eastern Finland and Natural Heritage Services of Metsähallitus[22]. The majority of the sampled individuals were pups (<1 yrs; 76%), and the main cause of death of both adults and pups was entanglement in fishing gear (66%). All specimens had been stored at -20°C.

Total genomic DNA was extracted using DNeasy Blood and Tissue Kits (Qiagen) following the manufacturer’s instructions. Microsatellite loci were amplified using primers originally designed for grey seal (Halichoerus grypus)[23–26], leopard seal (Hydrurga leptonyx)[27], and harbour seal (Phoca vitulina)[28, 29] (Additional file1: Table S1). PCR reactions were set up in 10 μl volumes, each containing 25 – 50 ng genomic DNA, 1 μM each primer, 0.5 U AmpliTaq Gold DNA polymerase (Applied Biosystems), 1X PCR buffer, 1.25 – 1.5 mM MgCl₂, and 0.1 mM of each dNTP (Finnzymes). The thermal cycler profile consisted of initial denaturation at 94°C for 3 min, followed by 30 – 35 cycles of denaturation at 94°C for 30 s, annealing at 50 – 58°C for 30s, and extension at 72°C for 30s, with a final extension step at 72°C for 5 min (Additional file1: Table S1). PCR products were run on an ABI 3730 DNA Analyzer (Perkin Elmer Applied Biosystems), and allele peaks scored using GeneMapper v. 4.0 (Applied Biosystems). Micro-Checker v. 2.2.3[30] was used to test for the presence of null alleles, allelic dropout, and scoring errors due to stuttering, with Bonferroni-adjusted 95% confidence intervals. The frequency of null alleles at each locus was estimated using the software FreeNA[31].

The microsatellite dataset (Additional file2: Table S2) was supplemented by, and contrasted with, mtDNA sequence variation in 215 individuals from the same population, obtained by Valtonen et al.[17]. The 704-bp mtDNA fragment spans the 5′ domain of the mitochondrial control region.

Spatial and temporal division of specimens

The genotyped individuals originate from different parts of Lake Saimaa and represent a time span of three decades (1980 – 2008). In order to analyse the spatial genetic structure of the population, the specimens were initially divided into four regional samples following the topography of the lake (Figure 1B): Northern Saimaa (microsatellites: N_ms = 15, and mtDNA: N_mt = 19), Haukivesi area (N_ms = 99 and N_mt = 116), Pihlajavesi area (N_ms = 43 and N_mt = 61), and Southern Saimaa (N_ms = 15 and N_mt = 19).

In addition, the individuals were divided into three temporal samples based on their collection decade: the 1980s (N_ms = 59 and N_mt = 79), 1990s (N_ms = 48 and N_mt = 54), and 2000s (N_ms = 65 and N_mt = 82). A decade was considered an appropriate time span to detect temporal changes, since the estimated generation time of the ringed seal is c. 11 years (estimated by Palo et al.[16] after Smith[32]).

Genetic diversity indices were also calculated across an alternative temporal division scheme, in which individuals were grouped based on their year of birth: 1965–1979 (N_ms = 14), the 1980s (N_ms = 54), the 1990s (N_ms = 47), and the 2000s (N_ms = 54) (see[17] for details on age determination).

Estimation of genetic diversity and effective population size

We used Arlequin v. 3.5.1.2[33] to estimate the number of alleles (A), observed (H_O) and expected (H_E) heterozygosities, and Wright’s inbreeding coefficients (F_IS) for each spatial and temporal sample. Departures from Hardy–Weinberg equilibrium and the presence of linkage disequilibrium between pairs of loci were tested with Arlequin and GENEPOP v. 4.1.3[34]. Estimates of allelic richness (A_R) were obtained using the rarefaction method implemented in HP-Rare[35]. A_R was estimated separately for each spatial and temporal sample, based on the smallest sample size within each division scheme. Individual observed heterozygosities (H_O) were calculated using IR macroN4 ([36];http://www.zoo.cam.ac.uk/directory/william-amos), and the relationship between heterozygosity and year of birth was tested by performing a linear regression analysis in SPSS Statistics v. 19 (IBM).

Effective sizes for the total population, as well as for different regional and temporal subsamples within the lake, were estimated using two different methods. First, we used the single-sample method based on linkage disequilibrium implemented in the software LDNe[37], while assuming a lowest allele frequency of 4% in order to prune singletons also from the smallest regional samples (= Northern and Southern Saimaa). Second, utilising the sampling from different decades, we estimated N_e with the temporal method in TempoFS[38], assuming Waples sampling scheme 1 (see[39] and references therein). For both methods, 95% confidence intervals were obtained through jackknifing. The sample size from the Haukivesi area allowed also estimating effective population sizes separately for different decades.

Assessments of population structure

Genetic differentiation among the spatial and temporal samples was investigated using analysis of molecular variance (AMOVA)[40, 41] and by estimating pairwise F_ST values between samples using Arlequin. For the two regions with the largest sample sizes (Haukivesi and Pihlajavesi), we also ran a hierarchical AMOVA with spatial samples subdivided into temporal subsamples, and vice versa. Significance levels were estimated on the basis of 10,000 permutations. Because F_ST is influenced by the level of heterozygosity (H_S) in the studied subpopulations[42–44], we used GenAlEx v. 6.501[45, 46] to estimate overall spatial differentiation also based on D_est[42] and G"_ST[44], which correct for genetic diversity as well as the number of subpopulations in the analysis; departures from zero were inferred based on 999 permuted datasets. To facilitate comparisons between the microsatellite and mtDNA datasets (see[47]), we estimated global D_est also for the mtDNA sequence data using a script[48] employing the seqinr[49] and ape[50] packages in R 3.0.1[51]; statistical significance was assessed based on 10,000 permutations of the data. The distribution of microsatellite variation across individual seals from the main regions of the lake was visualized with factorial correspondence analysis (FCA) implemented in the program GENETIX v. 4.05.2[52]. The statistical significance of differences in FCA axis 1 and 2 scores among regional groups were tested using a multi-response permutation procedure (MRPP) test in PC-ORD v. 5.33[53].

Genetic structuring of the population was further studied using Bayesian genotype-assignment approaches. Analyses in Structure v. 2.3.4[54, 55] were run without prior information on the sampling locality. We employed the admixture model with correlated allele frequencies, and implemented a burn-in period of 100,000 followed by 500,000 MCMC iterations. We set the number of clusters (K) from 1 to 8, and ran 20 replicate analyses per each K to ensure convergence. Subsequently, we determined the most likely number of clusters by estimating the logarithm of the data probability given each number of K (logP(X|K);[54]), and on the basis of ΔK values estimated using the ad hoc approach of Evanno et al.[56] implemented in Structure Harvester v. 0.6.93[57].

Because collection locations were known for all but one of the sampled 172 seals, we used TESS v. 2.3.1[58, 59] for a more detailed evaluation of spatial differentiation within the lake. In these analyses, individual sampling-site coordinates were included as prior information. Before the analyses, we removed some links from the default neighbourhood system created by the program, in order to improve its match to the geography of Lake Saimaa (Additional file3: Figure S1). Thereafter, we ran TESS using a no-admixture model with 200,000 iterations, of which the first 100,000 were excluded as a burn-in, to test the number of clusters (K) from 2 to 8, with ten replicates for each K. Plots of the deviance information criterion (DIC) against K were used to identify the most likely number of clusters, which was then used in 100 replicate runs with the admixture model, using the same number of iterations as above. Finally, we averaged results from the ten runs having the highest likelihoods using CLUMPP v. 1.1.2[60], and visualised the results using Distruct v. 1.1[61]. The initial hypothesis concerning the regional division of the Saimaa ringed seal population, which was based on the topography of the lake, was then reassessed according to the Structure and TESS results.

The presence of an isolation-by-distance (IBD) pattern among seals within the lake was tested by using SPAGeDi v. 1.3[62] to contrast the degree of relatedness among individuals with their geographic distances on a logarithmic scale[63]. Using the aforementioned location information, we estimated IBD for the microsatellite (N = 171) and mtDNA (N = 214) datasets using the kinship coefficient of Loiselle et al.[64], which is considered suitable for datasets containing rare alleles[65]. For mtDNA, we also used N_ij, a kinship analogue based on genetic distances among haplotypes[66]. Between-haplotype distances were estimated based on a TrN + I + Γ substitution model using parameter values from Valtonen et al.[17]. In all three analyses, ten spatial distance classes were created using the equal-frequency method, which generates uneven distance intervals with roughly equal numbers of pairwise comparisons. Significance levels for the mean kinship coefficient (departure from zero) for each distance class, as well as for the regression slope (b), were obtained by 10,000 permutations, and standard errors were estimated by jackknifing over loci.

Estimation of migration rates and sex-biased patterns of gene flow

Migration rates among Lake Saimaa regions were estimated using the Bayesian multilocus genotyping method implemented in BayesAss v. 1.3[67]. We performed 1 × 10⁶ burn-in iterations followed by 9 × 10⁶ iterations with a sampling frequency of 2,000, with Δ parameters of allele frequency, migration rate, and inbreeding coefficient set to 0.15, 0.17, and 0.17, respectively, to achieve the recommended acceptance rates of between 40% and 60%. The analysis was repeated three times with different initial seed numbers to ensure convergence.

In order to infer whether dispersal within the lake is sex-biased, we estimated the relative amount of male and female gene flow following the approach of González-Suárez et al.[68], which takes into account the different mode of inheritance of mitochondrial and nuclear genes, as well as their differing numbers of gene copies of within a population (see also[69]). This indirect method was chosen because the small number of subadults and adults in the data[22] prevented direct assessments of sex-specific differentiation. We first estimated the amount of male differentiation using equation 1b in González-Suárez et al.[68], and then assessed the ratio of male and female gene flow rates based on equation 2c in the same paper.

Data characterization, genetic diversity, and effective population size

Genotypes for all 17 microsatellite loci were obtained for 168 individuals, while four individuals lacked data for one locus each (Additional file2: Table S2)[22]. Significant heterozygote deficit was observed at eight loci, but the finding was not consistent through regional or temporal samples. Over all samples, Micro-Checker analyses suggested presence of null alleles at six loci (Hg1.4, Hg3.6, Hg6.1, Hg8.9, Hgdii, and Pvc78) with low frequencies (r < 0.1) in all cases, most likely due to population substructure (see below). Scoring errors suggested for the loci Hg3.6 and Hg8.9 were ruled out by independent manual rechecking of the data by two researchers. In Arlequin, consistent and significant linkage disequilibrium in all temporal samples after a sequential Bonferroni correction was found between loci Hg4.2 and Pvc78, whereas Bonferroni-corrected GENEPOP results indicated no disequilibrium. We consider the latter result more reliable, because the permutational likelihood-ratio test implemented by Arlequin assumes that all loci are in Hardy–Weinberg equilibrium[70], which is not the case here (Additional file2: Table S2). Therefore, the main analyses were performed using the full 17-locus dataset.

No obvious differences in genetic diversity were found among regions (Table 1). The data also suggested diversity decrease along the decades, but there is not enough statistical power for testing the pattern. However, there was a weak but statistically significant negative relationship between individual observed heterozygosity and year of birth (H_O = 3.762–0.002*birth year; r² = 0.030, P = 0.025; Figure 2). Temporal F_IS estimates increased over the study period from effectively zero for seals born before 1980 to a significantly positive F_IS = 0.100 (P = 0.001) in the 2000s (Table 1).

Population structure in space and time

In microsatellite analyses, significant overall differentiation was detected among the four main regions of Lake Saimaa (F_ST = 0.065, D_est = 0.043, G"_ST = 0.121; all P ≤ 0.001), and the result was confirmed by pairwise F_ST estimates among regions (Table 3). These patterns can be seen in the FCA plot, in which seals from the same region tended to cluster together, even though there was wide overlap among groups (Figure 3A). FCA axis loads of the regional groups differed significantly from each other (MRPP test, P < 0.001), and all pairwise comparisons were statistically significant (P < 0.01 in each case). Microsatellite differentiation was moderate in comparison to estimates based on mitochondrial control-region sequences, which were consistently higher for overall differentiation (F_ST = 0.390, P < 0.001[17]; D_est = 0.256, P = 0.047) as well as for pairwise F_ST values (Table 3). A highly significant negative relationship between relatedness and spatial distance between pairs of individuals was found using both microsatellites (b = -0.039, P < 0.001; Figure 4A) and mtDNA (b_Loiselle = -0.055; b_Nij = -0.067; both P < 0.001; Figure 4B).

			Northern	Haukivesi		Main Haukivesi	Pihlajavesi	Southern
	N ms	N mt	Saimaa	area	Kolovesi	area	area	Saimaa
Northern Saimaa	15	19	–	0.454	0.698	0.467	0.446	0.392
Haukivesi area	99	116	0.039**	–	–	–	0.389	0.333
Kolovesi	20	21	0.161	–	–	0.544	0.534	0.527
Main Haukivesi area	79	95	0.047	–	0.170	–	0.398	0.344
Pihlajavesi area	43	61	0.096	0.063	0.236	0.057	–	0.311
Southern Saimaa	15	19	0.153	0.071	0.209	0.075	0.054	–

N_ms = sample size in microsatellite analyses; N_mt = sample size in mtDNA analyses.

Estimates based on microsatellites are given below diagonal and mtDNA variation above diagonal. Results are shown for the initial four-region division, and for the updated five-region scheme in which Haukivesi is split into the Kolovesi and Main Haukivesi areas (in italics). Differentiation in all comparisons was statistically significant at P < 0.001 after a sequential Bonferroni correction, except for **P < 0.01.

Overall microsatellite differentiation among temporal samples based on decades (1980s, 1990s and 2000s) was very weak (F_ST = 0.002, P < 0.001), and this was true also for all pairwise comparisons (all F_ST ≤ 0.002, P ≥ 0.293). Hierarchical AMOVA of the two regions with the largest sample sizes (= Haukivesi and Pihlajavesi) likewise showed that, regardless of the order of the hierarchy, most of the differentiation could be attributed to variation among regions and individuals, while the effect of decades was weak and statistically non-significant (results not shown). In contrast, mtDNA-based analyses showed highly significant differentiation among decades (F_ST = 0.384, P < 0.001; see discussion in[17]).

In Bayesian assignment analyses in Structure, mean log-likelihood was maximized at K = 4 (Additional file4: Figure S2A), but the Evanno et al.[56] approach favoured a two-cluster model (Additional file4: Figure S2B). In the two-cluster results, 22 individuals were assigned to cluster 1 and 146 to cluster 2 (red and light blue in Figure 5A, respectively), while the assignments of four individuals remained ambiguous (membership coefficient Q < 0.7, following, e.g., Kopatz et al.[71]). Notably, all individuals assigned to cluster 1 originated from the Haukivesi area (Figure 5A) and, upon closer examination, 86% of them were found to be from the Kolovesi basin in northeastern Haukivesi (Figure 1C). The Kolovesi individuals were similarly identifiable in the four-cluster plot, but the geographic pattern of the other three clusters was less apparent (Figure 5B).

In the TESS assignment analyses, DIC reached a plateau at K = 4 (Additional file4: Figure S2C). However, the spatial distribution of individuals in these clusters did not entirely follow our predefined structure (Figures 1D and5C). Instead, Northern Saimaa was the only predetermined region that was more or less identifiable in its original form in the TESS chart (Cluster 3, blue), while the distribution of Cluster 4 (yellow) individuals roughly corresponds to Pihlajavesi and Southern Saimaa. Cluster 1 (red) individuals were largely restricted to the aforementioned Kolovesi basin (Figure 1D), whereas the rest of the Haukivesi area seemed to represent an admixture zone: although 35% of individuals belonged to Cluster 2 (green), 44% remained weakly assigned.

Based on these results, we reassessed our initial, topography-based division of Lake Saimaa by splitting the Haukivesi area (Figure 1B) into two subregions, Kolovesi and Main Haukivesi (Figure 1C,D). For the rest of the lake, we retained the original division scheme, on grounds of the significantly non-zero pairwise F_ST values (Table 3). As expected, using the five-region division elevated global estimates of differentiation, based on both microsatellites (F_ST = 0.107 vs. 0.065; D_est = 0.074 vs. 0.043; G"_ST = 0.196 vs. 0.121) and mtDNA (F_ST = 0.444 vs. 0.390; D_est = 0.298, P < 0.001, vs. 0.256). Pairwise F_ST estimates showed that Kolovesi is highly differentiated from the rest of the lake (Table 3), which can also be seen in the five-region FCA plot (Figure 3B; MRPP test, overall P < 0.001; pairwise contrasts, all P < 0.001). Furthermore, the significantly positive F_IS = 0.074, P = 0.0002 of the undivided Haukivesi area disappeared when Main Haukivesi and Kolovesi were analyzed separately (F_IS = 0.024, P = 0.183; and F_IS = -0.031, P = 0.741, respectively) (Table 1). Separating the Kolovesi and Haukivesi individuals elevated LD-based N_e estimates in the new Main Haukivesi area, while estimates based on the temporal method decreased (Table 2).

Migration rates and sex-biased patterns of gene flow

Genetic diversity and effective population size

The Saimaa ringed seal has very low genetic diversity (H_E = 0.36) in comparison to marine ringed seal populations, in which microsatellite-based estimates of expected heterozygosities range from H_E = 0.80 to H_E = 0.89[72–74]. Indeed, microsatellite variation within the subspecies is, to our knowledge, the lowest that has thus far been found within the order Pinnipedia in studies in which monomorphic loci have been excluded (cf. Mediterranean monk seal (Monachus monachus), H_E = 0.40[24]; Hawaiian monk seal (Monachus schauinslandi), H_E = 0.49[75]; spotted seal (Phoca largha), H_E = 0.51[76]; northern elephant seal (Mirounga angustirostris), H_E = 0.40[77]), and heterozygosity is low even when contrasted with the mean for placental mammals that have experienced a recent demographic threat (H_E = 0.50 ± 0.03;[78], see also[79]). The reduced diversity of the Saimaa population is not plausibly explained by ascertainment bias[78, 80], which could stem from the use of loci developed for other seal species, since it should similarly affect the levels of diversity in marine ringed seals. On the contrary, heterozygosity estimates in marine ringed seals[72–74] are in most cases higher than in the species for which the loci were originally designed[25, 27, 29]. Coalescent simulations indicate that the main loss of genetic diversity in the Saimaa ringed seal predates the 20^th-century anthropogenic bottleneck[16, 17], but the question remains as to whether the current population – given its size and spatial structure – will be able to maintain the remaining variability.

Unfortunately, it seems that the answer to this question is "No": Estimates of total and regional effective population sizes were very low, ranging roughly between 5 and 113 (Table 2). Comparable values have been obtained for other isolated and/or bottlenecked populations of large mammals, e.g., Finnish wolves (Canis lupus), N_e = 37.8 – 43.0[81], Western Canadian mountain goats (Oreamnos americanus), N_e = 23.2 – 65.9[82], and Iberian lynxes (Lynx pardinus), N_e = 8.5 – 23.1[79]. It should be noted that historical, long-term N_es of several hundred, estimated by Palo et al.[16] and Valtonen et al.[17], are far higher than our contemporary estimates, suggesting that the population has been substantially larger in the past. Notably, estimates based on linkage disequilibrium were clearly lower than those based on temporal changes in allele frequencies (Table 2). This is probably due to population structure lowering LD-based estimates (see[83]), as is indicated by the rise in Main Haukivesi N_e after exclusion of Kolovesi individuals (Table 2). We therefore consider the higher estimates based on the temporal method to be more applicable for the Saimaa ringed seal. When reflected to the current census size (310 in the year 2012[15]), the temporal N_e estimates yield a ratio of N_e/N_c = 0.17 – 0.36. The lower end of this range is consistent with the average ratio observed in wild populations (0.16), and the higher end would comply with the higher-than-average N_e/N_c ratios observed in small threatened populations[84].

Nevertheless, even the higher temporal N_es apparently are not enough for safeguarding the existing diversity: our detailed analysis based on the birth years of the sampled individuals detected slowly declining individual heterozygosity in microsatellite loci during the past few decades (Figure 2), which is in line with the post-1960s decrease in mtDNA variation observed by Valtonen et al.[17]. The fact that the overall inbreeding coefficient has risen at the same time (Table 1) suggests that the increasing autozygosity is due to a progressive increase in differentiation among and/or within regions rather than to loss of alleles. Examples of similar trends are accumulating from other species that have undergone recent anthropogenic collapses[85] and range fragmentations[79, 86].

Population structure and gene flow

Gene flow among populations is often restricted or even prevented by geographic barriers, but assessing what constitutes an obstacle for a given species is not always straightforward[11, 87, 88]. In analyses of spatial genetic variation, a clear benefit of the Saimaa ringed seal population is that the topography of the lake unambiguously determines possible routes of dispersal among subpopulations, and the level of differentiation then depends on the migration rates through these routes. Telemetry studies have shown that adult home ranges may span sites nearly 40 km apart[19] and that even first-year pups can travel up to 15 km per day after being weaned[21]. However, it has hitherto remained unclear whether such occasional long-distance movements lead to gene flow among the main breeding areas.

Our results demonstrate that ringed seal subpopulations inhabiting the main regions of Lake Saimaa are differentiated with respect to both autosomal and mitochondrial variation (see also[17]). However, Bayesian assignment analyses showed that genetic structuring within the lake is more subtle than we originally envisioned, by uncovering semi-isolation of the Kolovesi basin, which essentially constitutes a labyrinthine cul-de-sac separated from the main parts of Haukivesi by a narrow strait (Figure 1C,D). All in all, the level of genetic differentiation among different parts of Lake Saimaa (F_ST = 0.065 – 0.107) is remarkable considering the short distances among the main basins (Figure 1). As a comparison, Palo et al.[72] showed that ringed seals of the Baltic Sea are essentially panmictic (F_ST = 0.000), and geographic differentiation is weak also in grey seals, in which Graves et al.[89] found Saimaa-like F_ST values only between Baltic and North Sea breeding colonies (F_ST = 0.068 – 0.097), while differentiation within the Baltic was clearly lower (F_ST ≤ 0.023). In general, interpopulation genetic differentiation in marine seals is often negligible – or at least lower than that observed within Lake Saimaa – even across thousands of kilometres[72–74, 90–92]. A potential caveat here is that the dependence of F_ST on overall heterozygosity[42, 44] means that direct comparisons are complicated by the very disparate levels of genetic diversity in the Saimaa ringed seal and its marine counterparts. We therefore used the approach of Heller & Siegismund[43] and Meirmans & Hedrick[44] to convert overall differentiation indices among marine ringed seal populations estimated by Davis et al.[73] (F_ST = 0.005) and Martinez-Bakker et al.[74] (F_ST = 0.0086) into G"_ST values, based on the reported F_ST s, mean values of H_S, and numbers of sampled populations. The converted estimates (G"_ST = 0.047, and G"_ST = 0.056, respectively) resulting from these studies, both of which include populations on different continents, are still below the overall differentiation among regions within Lake Saimaa (G"_ST = 0.121 – 0.196).

As can be expected based on the high spatial differentiation, interregional migration estimates (Table 4) were very low in comparison to values estimated for other large mammals inhabiting fragmented landscapes (e.g.,[2, 11]). Across comparable geographic distances, equally low rates have been detected only for the endangered Ethiopian wolf (Canis simensis), in which remaining populations are restricted to mountaintops separated by unsuitable lowlands[3]. Although the migration estimates were generally near zero, it appears that Pihlajavesi, the most important breeding area producing some 40 – 50% of the pups born annually[15], serves as a source for the small Southern Saimaa subpopulation. Assessing whether the estimated migration rates are high enough to maintain demographic connectivity is difficult (cf.[10]), but our results suggest at least partial demographic independence of subpopulations, which could significantly increase the risk of stochastic extinction by lowering the demographic effective population size (N_d) well below the census size[93].

Microsatellites exhibited lower levels of interregional differentiation (measured as F_ST and D_est) and a smoother isolation-by-distance pattern than did mtDNA (Figure 4A,B), suggesting that gene flow is mainly mediated by males, as is often the case in pinnipeds[5, 94, 95] and mammals in general[96]. The roughly sevenfold male-to-female gene flow ratio is comparable to values reported for, e.g., harbour seals[97] and California sea lions (Zalophus californianus)[68].

Availability of supporting data

The microsatellite dataset supporting the results of this article is available in the Dryad repository (http://datadryad.org),http://dx.doi.org/10.5061/dryad.5j754. The mtDNA sequence dataset is available from Genbank (http://www.ncbi.nlm.nih.gov/genbank), accession numbers: JX109584–JX109833.