Origin and expansion of the world's most widespread pinniped: Range-wide population genomics of the harbour seal (Phoca vitulina).

The harbour seal (Phoca vitulina) is the most widely distributed pinniped, occupying a wide variety of habitats and climatic zones across the Northern Hemisphere. Intriguingly, the harbour seal is also one of the most philopatric seals, raising questions as to how it colonized its current range. To shed light on the origin, remarkable range expansion, population structure and genetic diversity of this species, we used genotyping-by-sequencing to analyse ~13,500 biallelic single nucleotide polymorphisms from 286 individuals sampled from 22 localities across the species' range. Our results point to a Northeast Pacific origin of the harbour seal, colonization of the North Atlantic via the Canadian Arctic, and subsequent stepping-stone range expansions across the North Atlantic from North America to Europe, accompanied by a successive loss of genetic diversity. Our analyses further revealed a deep divergence between modern North Pacific and North Atlantic harbour seals, with finer-scale genetic structure at regional and local scales consistent with strong philopatry. The study provides new insights into the harbour seal's remarkable ability to colonize and adapt to a wide range of habitats. Furthermore, it has implications for current harbour seal subspecies delineations and highlights the need for international and national red lists and management plans to ensure the protection of genetically and demographically isolated populations.

INTRODUCTION

The success of an organism is often evaluated according to its distributional range and the breadth of niches it can occupy (Afkhami et al., 2020; Hutchinson, 1957). Widely distributed generalist species typically have the capacity to expand or shift their range into new areas, whereas narrowly distributed and highly specialized species may struggle to track their preferred habitat in space and time. Assessing the dynamics of dispersal and range expansions of species is central to our understanding of the processes of adaptation, speciation and extinction. Also, monitoring species distributions, richness and interactions, as well as evaluating a species’ resilience to environmental change and disturbance from human activities are vital to good management practice (Liu et al., 2020; Martin et al., 2019). This is particularly relevant given the ongoing biodiversity and climate crises, in which species are experiencing range shifts or contractions at increasing rates (Singh, 2002; Urban et al., 2016).

Pinnipeds (Phocidae, Otariidae and Odobenidae) comprise a diverse, highly mobile and widely distributed group of marine mammals, which through adaptive radiation have come to occupy a wide range of breeding habitats and utilize diverse prey resources across polar, temperate and tropical environments. Similar to marine birds, sea turtles and marine iguanas, but uniquely among mammals, pinnipeds are characterized by a marked separation of core life activities; they forage in the aquatic realm and breed, moult and rest in terrestrial (or sea ice) habitats (Bartholomew, 1970; Cassini, 1999). These characteristics have implications for the distributional range of each pinniped species, which will be the product of its origin in time and space, as well as the distribution of suitable breeding and feeding habitats, the species’ ability to disperse and cross barriers, and whether it has been subject to recent environmentally or anthropogenically induced declines in abundance and distribution.

The harbour seal (Phoca vitulina) is the most successful pinniped species according to several criteria. It occupies a wide range of climatic zones and habitats of diverse geomorphology, hauling out on sea ice, beaches, rocky coasts and intertidal flats in marine, brackish and freshwater systems from the Subtropics to the Arctic (Davies, 1958). Despite being subject to population declines and displacements resulting from hunting, habitat destruction, environmental pollutants, interspecific competition and pathogens, the harbour seal has demonstrated a remarkable ability to recover and recolonize former ranges when efficient conservation measures have been implemented (Brasseur et al., 2018; Cammen et al., 2018; Härkönen et al., 2006; Olsen et al., 2018). Taken together, these observations might suggest that harbour seals, similar to other northern phocids, are characterized by a high dispersal ability and limited population structure (Andersen et al., 2009; Carr et al., 2015; Coltman et al., 2007; Davis et al., 2008; Fietz et al., 2016; Klimova et al., 2014; Thompson et al., 1996). However, this is not the case. Multiple genetic and tagging studies of local and regional populations have clearly demonstrated that harbour seals exhibit significant genetic differentiation across relatively fine geographical scales and that animals tend to be highly philopatric and rarely move more than 50–100 km from their natal terrestrial haul‐out areas (Andersen et al., 2011; Andersen & Olsen, 2010; Carroll et al., 2020; Cordes et al., 2017; Dietz et al., 2013; Goodman, 1998; Härkönen & Harding, 2001; O’Corry‐Crowe et al., 2003; Olsen et al., 2014, 2017; Rosing‐Asvid et al., 2020; Stanley et al., 1996b; Westlake & O’Corry‐Crowe, 2002).

How then did the harbour seal come to colonize such a broad, environmentally heterogeneous geographical range? Where did the species originate, which dispersal routes did it follow and did it encounter any barriers along the way? What patterns of genetic diversity and population structure did these processes create? What does this imply for current subspecies classifications? And what may this mean for species‐, subspecies‐ and population‐level responses to future stressors? To address these questions, we generated a comprehensive data set consisting of genome‐wide biallelic single nucleotide polymorphisms (SNPs) genotyped from 286 harbour seals sampled across the species range. Based on these SNPs, we used principal component analysis (PCA) and clustering analyses to assess genome‐wide patterns of genetic variation and to characterize population structure. Moreover, we inferred the historical relationship of the sampled harbour seal populations and quantified their genetic diversity, and finally we assessed the effect of distance and dispersal barriers in shaping the current diversity of the species. In addition to shedding light on the origin and expansion of harbour seal populations, our results have implications for the current subspecies debate. Furthermore, our results provide a baseline for future studies of the phenotypic and genotypic characteristics that have enabled harbour seals to prosper throughout most temperate coastal areas of the Northern Hemisphere and for understanding how the species might respond to future stressors.

MATERIALS AND METHODS

Sample collection and DNA sequencing

The study was initially based on 311 tissue samples collected from 22 locations spanning the geographical range of the harbour seal (Figure 1; Table S1). The samples were stored frozen under various conditions and preservation media (e.g., DMSO, ethanol, salt‐dried or in no medium). For most of the samples, DNA was extracted using the KingFisher Cell and Tissue DNA Kit following the manufacturer's protocol with a few modifications, including incubation at 56°C for 12–24 h using 25 µl proteinase K and 200 µl lysis buffer mix solution. Samples from New England were extracted using a Qiagen DNEasy Kit, as described in Cammen et al. (2018). Extraction blanks were added in every extraction round to control for contamination. All extractions were verified for the presence of a high‐molecular‐weight band on a 1–2% agarose gel and concentrations were measured using a Qubit 2.0 Fluorometer (Life Technologies).

FIGURE 1

Global and regional genetic structure of harbour seals. (a) Map of the sampling localities colour coded by cluster assignment. (b) Principal component analysis (PCA) of the global population of harbour seals (N = 286) demonstrating the existence of five major genetic clusters: Northeast Pacific, Northwest Pacific, Northwest Atlantic, Atlantic Arctic and Northeast Atlantic. (c) PCA of North Atlantic harbour seals (N = 230) illustrating a marked west–east and north–south pattern of genetic structure, as well as additional substructure within the Arctic North Atlantic cluster, corresponding to samples from Greenland, Svalbard and Iceland. (d) PCA of Northeast Atlantic harbour seals (N = 109), suggesting some separation between northern UK harbour seals and the rest, as well as a north–south genetic gradient from the southern North Sea to northern Norway. Sample size of each locality is given in parentheses

DNA samples with concentrations of 10–100 ng µl–1 were processed using genotyping by sequencing (GBS) at the Institute for Genomic Diversity at Cornell University following the method of Elshire et al. (2011). This included extraction and genotyping controls to test for contamination. First, RNAse A was added to all extractions to remove RNA, and then GBS reduced library construction was performed using the methylation‐sensitive “six base cutter” restriction enzyme PstI with recognition site 5′‐CTGCA↓G‐3′. Individual samples were digested and ligated with specific indexing adaptors before libraries were pooled and sequenced on an Illumina HiSeq2000 using single‐end 100‐bp chemistry and Illumina 1.8+ technology. The samples were run on a single Illumina lane for each 96‐well plate, resulting in single‐end reads of 101 bp including the 8–10‐bp barcodes.

Mapping

The raw reads were demultiplexed using the process_radtag function of stacks version.2.53 (Catchen et al., 2013). We trimmed low‐quality bases and adaptor sequences using fastp (Chen et al., 2018) with the following setting: a minimum value of 15 for qualified base (‐q 15), a maximum of 50% unqualified bases allowed per read (‐u 50) and a sliding window size of 5 bases (‐W 5). We further trimmed all clean reads to a uniform length of 91 bp in order to remove multiple bases of low quality at the 3′ end. The final pruned reads were mapped to a female harbour seal reference nuclear genome (NCBI GenBank accession: GCA_004348235.1) using bwa mem version 0.7.17 with default settings (Li, 2013). This reference genome is assembled on the scaffold level and consists of 2,363,596,144 bp over 5,541 scaffolds (N50: 41,024,070). For mapping statistics, see Table S2.

Data quality filtering

To remove low‐quality data before performing population genetic analyses, we used an iterative approach. First, we removed samples that were misidentified during sampling or storage, or were sequenced at too low depth. We then used the raw data from all the remaining samples to compile a set of filters that allowed us to discard low‐quality regions and loci in the genome. Next, using these filters we performed additional sample quality filtering where we removed samples that were either duplicated or closely related. Finally, we updated the quality filters for genomic regions and loci, so they were based only on the samples that passed all the sample quality filters.

Sample quality filtering: Low sequencing depth and mislabelling

We first excluded samples with a mean depth across all loci below 3 (n = 15). Furthermore, we excluded two samples that were originally labelled as harbour seal, but proved to be other species based on a preliminary PCA (Figure S1). After these filters a total of 294 samples were left for further analysis (Table S3).

Quality filtering of genomic regions and loci: Mappability and repeats

We estimated the mappability of the reference genome with genmap version 1.3 (Pockrandt et al., 2020) by calculating mappability scores with K‐mers of 90 bp and allowing three mismatches (‐K 90, ‐E 3). Loci with a mappability score less than 1 were excluded. In addition, we discarded loci in lower case in the reference fasta file, which were identified as repeats and low‐complexity DNA sequences using windowmasker (Morgulis et al., 2006).

Quality filtering of genomic regions and loci: Sex‐linked scaffolds

To identify sex‐linked scaffolds, we first calculated the average depth per scaffold per individual using the output of samtools idxstats and normalized the values by the mean depth of the 10 longest scaffolds (Li et al., 2009). Subsequently, we performed a PCA based on the normalized depth using scaffolds of a minimum length of 1 Mb and identified two groups of samples (Figure S2), which potentially represented the two sexes (Nursyifa et al., 2021). For a subset of the samples, the identification of sex was validated by consistency with morphological sex identification performed during sampling. Given that harbour seals are diploid animals, we further designated sex‐linked scaffolds as those in which the average depth of the assigned male was below 60% of the assigned females. These sex‐linked scaffolds were excluded from the downstream analysis (Table S4).

Quality filtering of genomic regions and loci: Sequencing depth and data missingness

Given the relatively high stochasticity in presence/absence of data often seen in GBS approaches, we excluded loci for which less than 50% of all 294 samples were sequenced with at least one read. Moreover, since the choice of this exact cutoff value was somewhat arbitrary, we also tested a stricter cutoff of 25%, which led to almost identical results in the downstream analyses, suggesting that the choice of cutoff within this range did not affect our final conclusions (results not shown). We further excluded extremely shallowly or excessively sequenced loci where the average depth of sequenced samples (with at least one read) was below 2 or above 50, with cut‐offs based on visual inspection of the distribution. Finally, we discarded the loci with substantial variation in depths among the sequenced samples (SD >30).

Quality filtering of genomic regions and loci: Excessive heterozygosity

We estimated preliminary genotype likelihoods separately for each locality using angsd (Korneliussen et al., 2014) with the following settings: minor allele frequency (MAF) ≥ 0.05, SNP p‐value < 2×10−6, and including only the genomic regions that passed the above‐mentioned quality filters. We subsequently calculated the site‐based inbreeding coefficient (F), ranging from −1 (all genotypes are heterozygous) to 1 (all genotypes are homozygous), and performed a Hardy–Weinberg equilibrium (HWE) likelihood ratio test in pcangsd (Meisner & Albrechtsen, 2019). Finally, loci exhibiting heterozygosity excess and significant deviations from HWE (F < 0.999 and p < .05) were excluded.

Final quality filtered sets of loci

The combination of all of the above quality filters on genomic regions and loci (mappability, repetitive regions, sex‐associated scaffolds, read depth and excessive heterozygosity) led to 8,960,582 accessible bases out of 2.39 billion bases present in the reference genome (Table S5). Based on this initial set of quality filtered loci, we performed sample quality control by filtering out duplicates and closely related samples (for details see below). With this filtered sample set, we then updated the following reference quality filters: (i) sequencing depth, (ii) missing data and (iii) excessive heterozygosity. We retained all other filters. This resulted in a final set of filtered regions and loci, passQCloci, comprising 8,969,175 bp that all the downstream population genetic analyses were based on unless otherwise mentioned. Note that many of the population genetic analyses performed in this study only rely on common (MAF >0.05) biallelic polymorphic loci. To identify these loci, we used the options ‐SNP pval 2 × 10−6 and ‐minmaf 0.05 in angsd, which led to data sets of ~13,500 common biallelic loci within the filtered regions. The exact number of SNPs varied depending on what samples were included and on how we filtered for missing data in different analyses (for details see Table S6).

Sample quality filtering: Duplicated and closely related samples

Sample duplicates and first‐degree relatives within populations were identified by calculating summary statistics based on the estimated two‐dimensional site frequency spectrum (2d‐SFS) for each pair of samples as described in Waples et al. (2019). Specifically, the summary statistics R0, R1 and KING‐robust kinship were calculated for all pairs of samples within each population, and duplicated samples and first‐degree relatives were identified based on the suggested threshold of KING‐robust kinship (Waples et al., 2019; Figure S3). Since some population genetic analyses can be affected by the presence of duplicated and first‐degree relatives, we then removed one sample from each such identified pair, yielding a total of 286 unrelated individuals for further analyses.

Input data for analyses: Genotype likelihoods and single read sampling

In order to take putative uncertainties in genotype calls due to the low sequencing depth into account, almost all of the population genetic analyses in this study relied on genotype likelihoods or single read sampling. Both were obtained using angsd limited to the quality filtered loci in the postQCloci set described above and using only reads with mapping scores of at least 30 and bases with a Phred quality score of at least 20. The genotype likelihoods were estimated using the SAMtools mode (‐gl 1) in angsd. The single read sampling was performed using angsd ‐doIBS 1, which randomly samples a single base at each position for each sample (for details about the angsd commands used to generate the input files for each analysis, see Table S6).

PCA and admixture analysis

We investigated the global population structure of all the harbour seals first by applying PCA, using pcangsd based on genotype likelihoods for all the loci in postQCloci that were polymorphic and had an MAF above 0.05 (Meisner & Albrechtsen, 2018). We set the maximum number of iterations to 1,000 to ensure convergence. We also estimated admixture proportions for different groups of populations by applying ngsadmix (Skotte et al., 2013) to the same data set with the number of assumed populations (K) ranging from 2 to 20. For each K, we performed the analysis independently 200 times, in order to assess convergence, which we defined as having the top five maximum likelihood results within 2 log‐likelihood units. We subsequently evaluated the model fit of the admixture analysis using evaladmix (Garcia‐Erill & Albrechtsen, 2020), which calculates the residuals under each model, and plotted the pairwise correlation of the residuals matrix between individuals.

To examine local population structure among the 18 sampling localities in the North Atlantic, as well as specifically for the 11 localities in the Northeast Atlantic (Table S1), we repeated the pcangsd and ngsadmix analyses for samples from these regions only with the assumed population number (K) ranging from 2 to 18, and from 2 to 10, respectively. Higher K values were not explored due to convergence issues (see Results).

Historical relationships among populations

To infer historical relationships among the sampled harbour seal localities, we built a population tree using treemix (Pickrell & Pritchard, 2012). treemix uses allele frequency data to infer the order in which different populations split from each other historically, and its output is analogous to phylogenetic trees for representing the evolution of species. To prepare the allele frequency input, we first took the single read sampling data from all the loci in the postQCloci set that were polymorphic with an MAF above 0.05 and had data missing from fewer than 30 samples. Then, within each locality we counted the number of major and minor alleles, which were determined for all samples at each of the loci, and used these counts to estimate allele frequencies for each locality. To improve the rooting of the inferred tree, we also included data from three grey seal (Halichoerus grypus) samples that were extracted, sequenced and bioinformatically processed in the same manner as the harbour seal samples. We performed the treemix analysis based on the estimated allele frequency data with a window size of 500 SNPs to account for linkage disequilibrium (‐k 500), specifying the grey seal samples as an outgroup. The harbour seal samples from Tromsø (TRO) were excluded from the treemix analyses due to insufficient sample size (n = 3).

Heterozygosity

To assess the genetic diversity of harbour seals, we estimated genome‐wide heterozygosity as the proportion of heterozygous loci according to the 1d‐SFS for each harbour seal sample. To estimate each of these 1d‐SFSs, we first used angsd to generate a saf (site allele frequency likelihoods) file based on genotype likelihoods for all loci in the postQCloci set. Then, we used this file as input to realSFS in angsd, which we ran with default settings. The resulting heterozygosity estimates for the individual samples were pooled and averaged for each sample locality. We did not detect any bias in the individual heterozygosity estimates caused by low sequencing depth (Figure S4).

Genetic differentiation, nucleotide diversity and isolation by distance

To quantify the levels of genetic differentiation and diversity among and within sampling localities, we estimated a number of SFSs using the function realSFS in angsd with the default settings. Specifically, we estimated the 1d‐SFS for each sampling locality and a 2d‐SFS for each pair of sampling localities. To prepare input for realSFS, we first used angsd to generate saf files using genotype likelihoods for all loci in the postQCloci set. Once the SFSs were estimated, we calculated a global estimate of Hudson's F ST (Hudson et al., 1992), which accounts for imbalanced sample sizes among populations, directly from the 2d‐SFS for each population pair using a custom python script. We also calculated the net nucleotide divergence (Da) for each population pair, following the formula: , where Dxy denotes the pairwise nucleotide differences between the populations and π denotes the nucleotide diversity within the population (Cruickshank & Hahn, 2014). The statistics of Dxy and π were directly calculated from estimated the 2d‐SFS and 1d‐SFS, respectively, with a custom Python script. We further tested for associations between the genetic distances (F ST) and geographical distances for all population pairs using Mantel permutation tests with the R package vegan (Oksanen et al., 2020). A significant association would indicate isolation‐by‐distance between harbour seal populations. The geographical distances were calculated based on the shortest waterway distance between locations, following https://rpubs.com/bearedo/AOTTP‐Seaway‐Distances.

Estimation of effective migration surfaces (EEMS)

To identify potential migration corridors and barriers, we conducted an Estimation of Effective Migration Surfaces (EEMS) analysis (Petkova et al., 2015). As an input to EEMS, we calculated the matrix of average genetic distances between all pairs of harbour seals based on the single‐read sampling data from all the loci in the postQCloci set that were polymorphic with an MAF above 0.05. The geographical limits for each EEMS analysis were drawn manually (http://www.birdtheme.org/useful/v3tool.html), and EEMS was deployed using the runeems_snps program with the settings of 18 million steps of the chain including 3 million burn‐in iterations and 600 demes. To ensure convergence, we ran 10 independent Markov chain Monte Carlo (MCMC) chains and extended the chain by 15 million iterations if the initial chain had not converged. We assessed the convergence of the EEMS analysis by using the Gelman–Rubin diagnostic based on the trace plots with the R package CODA (Plummer et al., 2006). The output of EEMS was visualized using the customized make_eems_plots R script (https://github.com/dipetkov/reemsplots2).

RESULTS

Genomic data and samples

In order to assess the probable evolutionary origin and expansion of harbour seals across the Northern Hemisphere, we generated GBS data for 311 samples from 22 localities representing the North Pacific and North Atlantic range of the species (Figure 1a; Table S1; for average depth per sample see Table S3). Rigorous data quality filtering led to a final genomic data set containing ~13,500 biallelic SNP loci with a mean depth of 10.1 (range: 1.0–50.0; median: 4.8) across 286 harbour seals, excluding 25 individuals from the original sample set due to low sequencing depth (n = 15), misidentification of species (n = 2), and first‐degree relatedness or sample duplication (n = 8). Given the low‐depth nature of these data, all population genomic analyses presented below relied on genotype likelihoods or single‐read sampling to account for genotype uncertainty.

Population structure

To investigate global and local harbour seal population structure, we performed both PCA and admixture analyses. The PCA suggested a division into five highly distinct genetic clusters comprising localities in the Northwest Pacific, Northeast Pacific, Northwest Atlantic, North Atlantic Arctic and the Northeast Atlantic (Figure 1b). Strikingly, PC1 explained 41.5% of the variation in the data set, separating Pacific and Atlantic harbour seals, whereas PC2 explained only 6% of the variation, mainly separating populations within ocean basins. In a PCA conducted solely on North Atlantic localities, PC1 and PC2 split the Arctic Atlantic cluster into Greenland, Iceland and Svalbard, revealing a marked geographical pattern across the North Atlantic, with PC1 (16.6%) separating sites by longitude and PC2 (5.9%) separating sites by latitude (Figure 1c). In both the global and North Atlantic PCA, subsequent PCs tended to separate Hokkaido, Svalbard and California, respectively, from other localities (Figure S5). Finally, the high sampling density in the Northeast Atlantic allowed us to detect further substructure (Figure 1d; Figure S5), including a separation of the Scottish localities Orkney and Lismore from all other Northeast Atlantic localities along PC1 (9.2%), a longitudinal gradient from northern Norway to the southern North Sea via Limfjorden and Kattegat along PC2 (5.4%), and another longitudinal gradient from northern Norway to the Southwest Baltic along PC3 (3.2%).

The results from the admixture analyses indicated an initial split into two main genetic clusters comprising the North Pacific and North Atlantic localities, respectively, when assuming the presence of two ancestral populations (K = 2, Figure 2). At K = 3, the North Atlantic split into a northeastern and northwestern cluster with no clear geographical division between them. At K = 4, the North Pacific split into clear northwestern and northeastern clusters comprising Hokkaido and Alaska–California, respectively, and at K = 5 the North Atlantic split into northwestern, Arctic and northeastern genetic clusters. As K increased up to 13, these five overall genetic clusters were further subdivided into clusters that clearly reflected the geography of sampling localities, including California, New England, Gulf of St. Lawrence, Newfoundland, Greenland, Iceland, Svalbard, Orkney–Lismore, the southern North Sea, Norway and Kattegat. An admixture analysis including only the 18 North Atlantic localities indicated even finer‐scale genetic clustering, illustrated by the detection of additional substructuring within Norway and Greenland (Figure S6). Evaluation of the admixture model fit using evaladmix indicated that these clusters fit the data well, as strong signals of nonzero residuals were only present in a few populations (Figure S7). Given the observed fine‐scale genetic structure, we treated all sampling localities as separate populations in all the subsequent analyses.

FIGURE 2

Fine‐scale genetic structure of global harbour seal populations illustrated by ngsadmix analysis with assumed number of ancestral populations K ranging from 2 to 13 (convergence was not obtained for any higher K values). Each small vertical bar represents a harbour seal and the colouring corresponds to its genetic ancestry. Locality names and sample sizes are listed in Table S1

Historical relationship and genetic diversity

We explored the origin and expansion of harbour seals by inferring historical relationships among the sampled harbour seal localities and quantifying their genetic diversity. Historical relationships were inferred using treemix, which infers the history of splits based on allele frequencies. The results of this analysis supported a split between North Pacific and North Atlantic localities, with Alaska and California in the Pacific being the first to split after the grey seal (Figure 3b). In the North Atlantic, Newfoundland and the other Northwest Atlantic localities were at the base of the tree, which had gradual bifurcation across the North Atlantic from west to east, with localities in the southern North Sea (i.e., the Wash and Dutch Wadden Sea) splitting off last.

FIGURE 3

Historical relationships and genetic diversity of harbour seal populations. (a) Estimates of genome‐wide heterozygosity per animal for each sampling locality demonstrate a substantial loss of genetic diversity during the harbour seal's expansion from the Northeast Pacific, indicative of a serial founder event across the North Atlantic. (b) The historical relationship among harbour seal populations inferred by treemix suggests a Northeast Pacific origin for harbour seals with subsequent spread to the North Atlantic. The horizontal branch lengths (drift parameters) are proportional to the amount of genetic drift that has occurred on the branch, resulting either from long‐term isolation or population bottlenecks. Locality names and sample sizes are listed in Table S1

In accordance with the historical relationships observed in the treemix results, the North Pacific populations exhibited substantially higher levels of genetic diversity than North Atlantic populations (Figure 3a; Figure S8). The high diversity of North Pacific harbour seals was primarily driven by Kodiak Island and Endicott Arm in Alaska, but California and Hokkaido also had higher levels of diversity than observed in any of the North Atlantic localities. Intriguingly, within the North Atlantic, there was a clear trend of declining diversity along a longitudinal west–east cline from the USA and Canada to Europe. Taken together, the treemix and genetic diversity results are consistent with a Northeast Pacific origin for harbour seals and a stepping‐stone colonization of the Atlantic in a series of founder events from the northwest toward the northeast.

Genetic differentiation, isolation by distance and dispersal barriers

To quantify genetic differentiation among the different sampling localities, we computed pairwise estimates of F ST and net nucleotide divergence (Da) between all pairs of locations (Figure 4a; Figure S9; Table S7). The two estimates showed very similar patterns and we thus limit our description to the pairwise F ST estimates. These indicate a deep genetic split between North Pacific and North Atlantic harbour seals (F ST =0.38–0.69) with some substructure within both the North Pacific (F ST =0.00–0.27) and the North Atlantic (F ST =0.02–0.38). Notably, Hokkaido (F ST =0.22–0.27) and Svalbard (F ST =0.24–0.38) appear particularly distinct from other populations within their respective ocean basins, suggesting a high degree of isolation.

FIGURE 4

The estimates of genetic differentiation indicate a deep genetic split between North Pacific and North Atlantic harbour seal localities and substantial substructure within ocean basins, driven by geographical distance and dispersal barriers. (a) Pairwise F ST estimates among all sample localities; (b) isolation by distance across the harbour seal's range plotted as pairwise F ST estimates among sampling localities against their shortest waterway distance measured in kilometres; (c) dispersal barriers (shades of brown) identified by an EEMS analysis. Locality names and sample sizes are listed in Table S1

Next, we plotted F ST against the shortest waterway distance for each pair of sampling locations to investigate the role of geographical distance in the observed patterns of genetic differentiation (Figure 4b). This revealed that the F ST values to a large extent fit a model of isolation‐by‐distance, which was also supported by a Mantel test (r = .8973, p < 2 × 10−6). However, there were two exceptions to this pattern. First, the population pairs within the Pacific Ocean all lie at the bottom of the Atlantic “point cloud,” suggesting there might be less pronounced dispersal barriers driving differentiation in the Pacific. Second, pairwise F ST estimates within the Pacific and across oceans differed substantially depending on whether Hokkaido (HOK) was included in the analysis or not. A lack of samples from much of the northwestern Pacific (including the Aleutian Islands, Kamchatka Peninsula and Kuril Islands) prevents us from determining whether this apparent distinctiveness of seals from Hokkaido is an artefact caused by incomplete sampling or results from a real pattern of Hokkaido (Kuril) seals being strongly isolated from other Pacific harbour seals.

Finally, to further investigate potential deviations from isolation‐by‐distance caused by barriers to dispersal, we performed an EEMS analysis for all of the harbour seal populations. Nine out of 10 MCMC realizations in the EEMS converged based on the diagnostic posterior trace plots, where the potential scale reduction factor of Gelman and Rubin's diagnostic is less than 1.1 (point estimate: 1.04, upper confidence interval: 1.08), suggesting no convergence issues were detected. In general, the EEMS results (averaged across the nine converged chains) fit with an isolation‐by‐distance model, but also reveal multiple effective barriers to dispersal (m < −1, see Figure 4c). Not surprisingly, the North American continent separates Pacific and Atlantic harbour seals, and intermittent periods (e.g., glacials) with extensive Arctic Sea ice have limited dispersal. In addition, wide stretches of open ocean appear to act as dispersal barriers: the Labrador Sea and Denmark Strait isolating Greenlandic harbour seals; the Northeast Atlantic, Norwegian Sea, Barents Sea and Greenland Sea separating Iceland and Svalbard harbour seals from Europe; the northern North Sea separating Scottish and Scandinavian seals; and finally, the Scandinavian mainland separating seals in the inner Danish waters from northern Norway.

DISCUSSION

Northeast Pacific origin and leading‐edge colonization of the North Atlantic

Fossil finds and phylogenetic analyses of molecular and morphological data strongly indicate that the earliest phocid seals originated in the North Atlantic and Mediterranean/Paratethys region (Berta et al., 2018; Deméré et al., 2003; Fulton & Strobeck, 2010; Rule et al., 2020). However, the exact evolutionary relationships among phocid seal species are unclear, and their respective origins and radiations into subspecies and populations on a global scale are complex (Fulton & Strobeck, 2010; Rule et al., 2020). It has been suggested that the genus Phoca diverged from Pusa in the Atlantic Arctic ~1.5 million years ago and that harbour seals subsequently colonized the North Pacific during an interglacial period (Berta et al., 2018; Deméré et al., 2003; Fulton & Strobeck, 2010). However, this is at odds with several other observations; first, the harbour seal's sister species, the largha seal (Phoca largha), is only found in the North Pacific, indicating that the harbour seal and largha seal diverged in this region, not in the North Atlantic. Deméré et al. (2003) proposed that the largha seal may be the descendant of one of the earliest Phoca populations to colonize the North Pacific. However, this would imply that harbour seals and largha seals are paraphyletic species, which contradicts current mitochondrial and nuclear DNA data (Cordes et al., 2017; Mizuno et al., 2018; Nakagawa et al., 2010; O’Corry‐Crowe & Westlake, 1997; Stanley et al., 1996a). Second, among the Phoca and harbour seal fossils reported to date, the oldest—dating back to the middle to late Pliocene and Pleistocene—are from Oregon and California (Barnes & Mitchell, 1975), lending support to a North Pacific origin of the harbour seal. In agreement with these observations, we find that contemporary populations of Northeast Pacific harbour seals are basal in the graphical representation of the historical relationships among populations provided by treemix, and that the harbour seals in this region have substantially higher levels of genetic diversity than their conspecifics in the North Atlantic. Taken together, our results point towards a North Pacific origin of the genus Phoca, and perhaps more specifically in the Bering Sea–Chukchi Sea region given the probable Arctic origin of Phocinae (northern true seals) (Rule et al., 2020).

Our treemix analysis also placed present‐day Newfoundland harbour seals basally among the North Atlantic populations, suggesting that the migration from the North Pacific to the North Atlantic occurred via the Canadian Arctic. GBS data do not allow for reliable estimates of divergence times, but the genetic composition of North Atlantic harbour seals, the species’ preference for temperate climates and the occurrence of fossils in the North Atlantic dating back to the Late Pleistocene (Deméré et al., 2003) suggest that migration may have occurred as a single founding event during a recent interglacial period. Indeed, in their analysis of the spatiotemporal occurrence of fossils of extant and extinct marine mammals, Peredo and Uhen (2016) hypothesized a substantial exchange between marine mammals of the North Pacific and North Atlantic through the Canadian Arctic during the Eemian interglacial (Marine Isotope Stage 5e; 124–119 thousand years ago), which was characterized by high global sea‐levels (Hearty et al., 2007; Rohling et al., 2007). The harbour seal's subsequent expansion across the North Atlantic appears to have occurred in a series of bifurcating founding events, one southwest along the coastline to other Canadian and US localities, and another to Greenland. From Greenland, some harbour seals appear to have successfully made the jump to Svalbard and others moved southeast to Iceland, supporting findings by Andersen et al. (2011). The colonization history of Europe is more complex, but the treemix analysis suggests that northern Scotland and Norway split off first, with harbour seals in the southern North Sea and the Kattegat–Baltic Sea region splitting off next. It is unclear how widespread harbour seals were in the North Atlantic during the last glacial period. However, it seems likely that the southern North Sea and Kattegat–Baltic Sea region was first colonized as Fennoscandia started to deglaciate 22,000 years before present (bp) and Doggerland became submerged 7,000–8,000 bp (Coles, 1998; Stroeven et al., 2016). The oldest harbour seal subfossil material from Denmark dates to 6000–8000 bp (Aaris‐Sørensen, 2009).

Our estimates of genome‐wide heterozygosity vary substantially among regions and show a marked gradient from the Northeast Pacific (Alaska) to the Northwest Pacific (Hokkaido) and Northwest Atlantic (New England, Newfoundland and Gulf of St Lawrence), Atlantic Arctic (Greenland, Svalbard and Iceland) and further to the Northeast Atlantic (European mainland) populations where the lowest levels of diversity are found. This pattern largely resembles that of a leading‐edge colonization where serial founder events result in a gradual loss of diversity towards the colonization front. Although rarely described in marine species, genetic footprints of range expansion have been described for humans (Homo sapiens), where genetic diversity decreases with distance from Africa (Rosenberg et al., 2002). Similarly, Northern Hemisphere terrestrial flora and fauna have diversity gradients along a north–south axis that is typically attributed to northward range expansions during interglacials (Hewitt, 2004; Hewitt et al., 1996). Founding events during such range expansions can result in beneficial and deleterious mutations attaining high frequencies or becoming lost at the leading edge of the colonization wave (Excoffier et al., 2009; Hallatschek & Nelson, 2008). In this context, it is intriguing that the highly diverse and hypothesized source population in the Northeast Pacific has escaped major viral disease outbreaks, whereas Hokkaido and North Atlantic harbour seals have not (Fujii et al., 2007; Härkönen et al., 2006; Stokholm et al., 2019). This raises the question as to whether genetic variation putatively involved in pathogen susceptibility has been lost or fixed during the harbour seal's expansion from the Northeast Pacific?

An efficient colonizer despite philopatry and dispersal barriers

In pinnipeds and many other marine species, the main factors limiting dispersal are geographical distance and barriers such as land, extensive sea ice, or strong temperature gradients (Faria et al., 2021; Palumbi, 1992). Our range‐wide study suggests that geographical distance has been a main driver of isolation in harbour seals, supporting previous findings on a regional scale (Goodman, 1998; Olsen et al., 2014; Westlake & O’Corry‐Crowe, 2002). However, the isolation‐by‐distance and EEMS results indicate that stretches of open water might also act as barriers to dispersal in harbour seals. This observation is supported by tagging data from both North Atlantic and North Pacific harbour seal populations, which demonstrate that harbour seals perform foraging trips of up to several hundred kilometers, but almost always return to the same haul‐out site or region and very rarely cross large open stretches of deep water (Carroll et al., 2020; Dietz et al., 2013; Peterson et al., 2012; Rosing‐Asvid et al., 2020; Small et al., 2006; Womble & Gende, 2013). In contrast, other northern phocid species, such as the sympatric grey seal, as well as the Arctic harp (Pagophilus groenlandicus), hooded (Cystophora cristata) and ringed seals (Pusa hispida), regularly undertake long distance migrations and/or movements covering thousands of kilometers and cross open water (Andersen et al., 2009; Nordøy et al., 2008; Svetochev et al., 2016; Thompson et al., 1996; Yurkowski et al., 2016). The open water barriers to harbour seal dispersal are particularly prominent in the North Atlantic, where population connectivity appears to be reduced by a heterogeneous seascape consisting of haul‐out sites distributed across multiple distinct coastlines and islands. In contrast, connectivity in the North Pacific has likely been facilitated by the existence of a long and more or less continuous coastline along most of the North Pacific rim, with the exception of the Bering Sea and Sea of Okhotsk, which might have resulted in the isolation of harbour seals in Hokkaido (and perhaps Kamchatka) from Northeast Pacific populations.

The harbour seal's success in colonising new habitats is likely due to several factors, which may collectively facilitate the occupation of more predator‐exposed, marginal and dynamic environments than other seal species. The harbour seal is a dietary generalist, which preys on the fish, squid and crustaceans that are locally and seasonally available (e.g., Scharff‐Olsen et al., 2019). Thus, when colonising new areas, harbour seals are likely not often limited by prey availability, and their relatively small body size implies that their overall energy requirements are low. Additionally, like several other phocids, the harbour seal becomes sexually mature at a young age (3–5 years; Härkönen & Heide‐Jørgensen, 1990), and is characterised by high annual pregnancy rates (up to 90%), resulting in relatively high population growth rates (up to 13%) (Olsen et al., 2010). Furthermore, and in contrast to most other phocids, harbour seal pups shed their natal white lanugo fur in utero and are born with some blubber and brown fat, allowing them to effectively thermoregulate on both land and in water, and thus have some mobility during and after lactation (Oftedal et al., 1991). We hypothesise that these anatomical, physiological and behavioural characteristics make the harbour seal an effective coloniser, capable of forming and maintaining viable populations from a few pioneering dispersers across a wide breadth of habitats. Indeed, a wide distribution, but strong population structure may be a general hallmark of philopatric, but otherwise opportunistic and generalist marine species, including several species of dolphins (Andrews et al., 2010; Louis et al., 2014; Manlik et al., 2019; Möller et al., 2011), sea turtles (Baltazar‐Soares et al., 2020; Clusa et al., 2018; Simões et al., 2021) and sharks (Feldheim et al., 2014; Klein et al., 2019; Mourier & Planes, 2013). In such species, we hypothesise that the colonisation of new habitats mainly happens in the form of rare long‐range dispersal events followed by very limited subsequent gene‐flow. It is uncertain what triggers such dispersal events. In harbour seals and many other species, juveniles ‐ and in particular pups of the year ‐ are the dispersing age class (Carroll et al., 2020; Dietz et al., 2013; Hamilton et al., 2019; Härkönen & Harding, 2001; Thompson et al., 1994). Thus, in harbour seals, colonisations may have resulted from juvenile animals embarking on exploratory trips.

Implications for harbour seal taxonomy

The taxonomy of harbour seals is a matter of debate. Up to five subspecies have been recognized: P. v. stejnegeri—or Kuril seal—in the Northwest Pacific; P. v. richardii in the Northeast Pacific; P. v. concolor in the Northwest Atlantic; P. v. vitulina in the Northeast Atlantic; and the freshwater Ungava seal (P. v. mellonae) (Davies, 1958; Rice, 1998). However, based on unclear delineations and indications of paraphyly in mitochondrial DNA data (Stanley et al., 1996), these subspecies have been questioned, and only three subspecies—P. v. richardii in the Pacific, P. v. vitulina in the Atlantic and the Ungava seal (P. v. mellonae)—are recognized by some researchers (Berta & Churchill, 2012; Committee on Taxonomy, 2021). The genomic analyses presented herein serve to resolve some of these issues, while also raising new ones. First, our results indicate that the North Atlantic consists of either a single, three (i.e., Northwest Atlantic, Atlantic Arctic and Northeast Atlantic) or five subspecies (i.e., Northwest Atlantic, Greenland, Svalbard, Iceland and Northeast Atlantic), depending on the criteria for defining marine mammal subspecies (Olsen & Galatius, 2018; Taylor et al., 2017). A conservative approach would imply abandoning P. v. concolor, while retaining P. v. vitulina to describe North Atlantic harbour seals, as has also been suggested by other recent assessments (Berta & Churchill, 2012; Committee on Taxonomy, 2021). Second, although our data have sampling gaps across much of the northern Pacific rim, we did find harbour seals from Hokkaido to be markedly genetically different from seals in the Northeast Pacific, including Alaska and California. This may be a result of our samples originating from the seemingly isolated Cape Erimo population (Mizuno et al., 2020) or may reflect Hokkaido (Kuril) harbour seals generally being different from other Pacific populations. However, we also note that while Hokkaido is monophyletic in the treemix analysis, it sits within a paraphyletic Northeast Pacific group. Thus, some of our results support the delineation of P. v. stejnegeri and P. v. richardii in the North Pacific, whereas others do not, highlighting the need for further comprehensive genome‐level analyses including samples from across the North Pacific. We were not able to obtain samples from the freshwater P. v. mellonae, so its status as a subspecies remains unclear.

Intriguingly, the division between Pacific and Atlantic harbour seals is supported by remarkably high levels of genetic differentiation, with F ST ranging from 0.38 to 0.69 for pairwise population comparisons between the two oceans. In comparison, clearly defined subspecies or ecotypes in other marine megafauna, such as harbour porpoises (Phocoena phocoena) and bottlenose dolphins (Tursiops truncatus) are characterized by F ST values up to the range of 0.157–0.270 for genomic data (Chehida et al., 2021; Louis et al., 2014; Nykänen et al., 2019), which is much lower than those we estimate between Pacific and Atlantic harbour seals. The understanding of marine speciation processes is certainly not settled (Faria et al., 2021). For example, there are known instances of species hybridizations in marine mammals (Bérubé & Aguilar, 1998; Kovacs et al., 1997; Savriama et al., 2018; Skovrind et al., 2019) and there is little consensus generally on how to delineate subspecies and species (Olsen & Galatius, 2018; Taylor et al., 2017). However, it is clear that the high levels of genetic differentiation estimated between Pacific and Atlantic harbour seals are rare in studies of marine megafauna, and could potentially point to the existence of two distinct species.

Perspectives

All harbour seal populations have historically been subject to extensive commercial hunting and culling (Bowen & Lidgard, 2013; Cammen et al., 2019; Crowell, 2016; Olsen et al., 2018) and several have experienced devastating viral disease outbreaks (Anthony et al., 2012; Härkönen et al., 2006). Fortunately, the characteristics that enable harbour seals to be efficient colonizers, typically faring better than many other pinnipeds in unpredictable and unstable habitats, also make harbour seals resilient to human disturbance, disease epidemics and environmental fluctuations. Many populations that were reduced in the past have now recovered, have recolonized former haul‐out sites and some are even approaching carrying capacity (Cammen et al., 2019; Harding et al., 2018; Sigourney et al., 2021). There are notable exceptions, including populations in Scotland, Norway, Greenland, Iceland and Hokkaido, which are endangered and still under extreme pressure from hunting, fisheries bycatch and other human activities, as well as potentially interspecific competition with grey seals (Andersen et al., 2011; Granquist & Hauksson, 2019; Kobayashi et al., 2014; Rosing‐Asvid et al., 2010). Moreover, the Faroe Islands have not yet been recolonized by harbour seals following their local extinction in the mid‐19th century (Mikkelsen, 2010). A recent population viability analysis showed how a seemingly stable harbour seal population close to carrying capacity can collapse with just slightly reduced fecundity, slightly increased hunting pressure and sporadic epizootics (Silva et al., 2021). The strong site fidelity of harbour seals should be taken into account in both international and national management plans to ensure the protection of locally adapted unique populations.

To understand how harbour seal populations respond to ongoing and future human disturbances and environmental fluctuations, studies should focus on identifying putative genotypic and phenotypic characteristics of local populations. This could include interdisciplinary approaches: to shed light on pup movement and recruitment; to elucidate the extent to which the environment, phenotype and genotype influence philopatry and dispersal strategies; to explore how such strategies may differ between leading‐ and trailing‐edge populations during range expansion; and to investigate how the diversity of dispersal strategies is maintained as an adaptation to future colonization of new areas. In addition, utilizing the full power of nuclear genome data would also shed more light on the relationships within Phoca and the timing of the migration from the North Pacific to the North Atlantic. Genomic data may also resolve the issue of harbour seal subspecies (and species) with consequences for the classification of taxa on the IUCN Red List. This could include phylogenetic analyses of mitogenome data and genetic assignment tests based on nuclear data to infer the level of (sub)species and population diagnosability (Pattent & Unitt, 2002; Taylor et al., 2017), as well as systematically reviewing other lines of evidence, such as variation in skull morphology, pelage coloration, breeding phenology and other traits that may vary between geographically distinct harbour seal populations.

CONFLICT OF INTERESTS

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

X.L., I.M. and M.T.O. designed the research; R.D., K.H., T.H., I.M. and M.T.O. provided funding; S.M.G., A.R.A., R.D., J.T., A.G., K.C., G.O.C., K.H., T.H., A.H., E.L.C., J.I.H., Y.K., M.H., G.S., A.K.F., C.L., K.M.K., L.W.A., S.G. and M.T.O. provided materials; S.R.S., K.C., G.O.C., E.L.C., Y.K. and M.T.O. performed the laboratory work; X.L., S.R.S., F.G.V., R.H., I.M. and M.T.O. analysed the data; X.L., I.M. and M.T.O. drafted the paper; all authors provided editorial inputs and approved the final version of the manuscript.

Figure S1‐S9

Click here for additional data file.

Table S1‐S6

Click here for additional data file.