Origins of the Greenland shark (Somniosus microcephalus): Impacts of ice-olation and introgression.

Herein, we use genetic data from 277 sleeper sharks to perform coalescent-based modeling to test the hypothesis of early Quaternary emergence of the Greenland shark (Somniosus microcephalus) from ancestral sleeper sharks in the Canadian Arctic-Subarctic region. Our results show that morphologically cryptic somniosids S. microcephalus and Somniosus pacificus can be genetically distinguished using combined mitochondrial and nuclear DNA markers. Our data confirm the presence of genetically admixed individuals in the Canadian Arctic and sub-Arctic, and temperate Eastern Atlantic regions, suggesting introgressive hybridization upon secondary contact following the initial species divergence. Conservative substitution rates fitted to an Isolation with Migration (IM) model indicate a likely species divergence time of 2.34 Ma, using the mitochondrial sequence DNA, which in conjunction with the geographic distribution of admixtures and Pacific signatures likely indicates speciation associated with processes other than the closing of the Isthmus of Panama. This time span coincides with further planetary cooling in the early Quaternary period followed by the onset of oscillating glacial-interglacial cycles. We propose that the initial S. microcephalus-S. pacificus split, and subsequent hybridization events, were likely associated with the onset of Pleistocene glacial oscillations, whereby fluctuating sea levels constrained connectivity among Arctic oceanic basins, Arctic marginal seas, and the North Atlantic Ocean. Our data demonstrates support for the evolutionary consequences of oscillatory vicariance via transient oceanic isolation with subsequent secondary contact associated with fluctuating sea levels throughout the Quaternary period-which may serve as a model for the origins of Arctic marine fauna on a broad taxonomic scale.

INTRODUCTION

For highly mobile marine organisms with extensive distributions, the mechanisms driving genetic divergence and speciation can be difficult to discern. The lack of obvious physical barriers to gene flow in the marine environment offers the potential for genetic connectivity over broad geographic and temporal scales with some species showing little genetic differentiation even among different oceanic basins (e.g., sharks and bony fishes: Schmidt et al., 2009; Roy, Hardie, Treble, Reist, & Ruzzante, 2014; da Silva Ferrette et al., 2015). Classically, marine speciation was largely thought to occur in allopatry through vicariance (Avise, 1992; Palumbi, 1994), whereby geologic isolating barriers (i.e., mountain chains, isthmuses) or hydrologic phenomena (oceans, rivers, drainage basins, and sea ice cover) isolate formerly large species distributions (Briggs, 2006; Corrigan & Beheregaray, 2009; Gaither et al., 2010; Knowlton & Weigt, 1998; O'Corry‐Crowe et al., 2016). Numerous challenges to the notion of allopatry as the primary driver of marine divergence exist (Rapoport, 1994; Kinlan & Gaines, 2003; Palumbi, 2004; in Norris & Hull, 2012), suggesting that constraints to gene flow and/or strong selection can also promote genetic differentiation among marine populations in the absence of physical barriers (Coyne & Orr, 2004; Horne, van Herwerden, Choat, & Robertson, 2008; Miglietta, Faucci, & Santini, 2011; Palumbi, 1994). Marine examples of “soft vicariance” (e.g., isolation in oceanic currents; see Hickerson & Meyer, 2008; Cowman & Bellwood, 2013), isolation at larger regional spatial scales (Mach et al., 2011), drift effects along environmental gradients (Ingram, 2011; Roy, Hurlbut, & Ruzzante, 2012), local adaptation (Bremer, Mejuto, Greig, & Ely, 1996; Clarke, Munch, Thorrold, & Conover, 2010), and ecological diversification in sympatry (Crow, Munehara, & Bernardi, 2010; Puebla, 2009; Rocha, Robertson, Roman, & Bowen, 2005), have all been documented in marine taxa with the potential for widespread connectivity.

In Arctic marine environments, the evolutionary histories and population demographics of numerous endemic species appear to be influenced strongly by fluctuations in the extent of sea ice associated with cooling glacial and warming interglacial cycles since the onset of the geologic Quaternary period roughly 2.6 Ma ago (e.g., O'Corry‐Crowe et al., 2010). Genetic structure in beluga whales Delphinapterus leucas (O'Corry‐Crowe, Suydam, Rosenberg, Frost, & Dizon, 1997; O'Corry‐Crowe et al., 2010) and walruses Odobenus rosmarus (Andersen & Born, 2000) provides strong indicators that ice effects during the Pleistocene (0.0117–2.588 Ma ago), coupled with additional ecological factors such as large‐scale temporal variations in levels of primary productivity, can leave lasting genetic imprints in a species' evolutionary history. Range expansions and contractions have also been shown in other Arctic and sub‐Arctic species (Andersen et al., 1998; O'Corry‐Crowe, 2008; O'Corry‐Crowe et al., 1997; O'Corry‐Crowe et al., 2010; Palsbøll, Heide‐Jorgensen, & Dietz, 1997; Palo, 2003), associated with sequential cooling and warming periods, respectively. Specifically, for beluga whales, late Pleistocene oscillations (7–36 ka) in sea ice extent may explain recent population subdivisions, with genetic support for recurrent episodes of gene flow that likely coincided with warming periods (O'Corry‐Crowe et al., 2010). Species divergence, occurring with such episodic gene flow, can have important impacts on population demographic parameters such as connectivity and size (Hey, 2010b; Moore, Gow, Taylor, & Hendry, 2007; Nosil, 2008), which can then lead to substantial alteration of a species' adaptive potential (Eizaguirre & Baltazar‐Soares, 2009; Hellberg, 2009; Pelletier, Garant, & Hendry, 2009). Estimating these parameters is, therefore, an important step toward a better understanding of the evolutionary histories and potential trajectories of different species, especially considering the dynamic environmental conditions that likely govern potential gene flow among them (Seehausen et al., 2008). Disentangling the interplay among demographic parameters, environmental fluctuations, and/or the potential for gene flow, however, can be especially complicated among closely related cryptic sister species whose taxonomic differentiation remains nebulous (Griffiths et al., 2010). Such information is nevertheless urgently needed, especially for species whose main distribution includes large Arctic marine ecosystems, which are currently experiencing major climate‐related transformations (Wassmann, Duarte, Agusti, & Sejr, 2011).

The number of currently recognized Somniosus species varies according to sources, with species regarded as geographically isolated from one another despite little differences in physical appearance having occurred among sister taxa over their evolutionary histories (Yano, Stevens, & Compagno, 2004). Greenland sharks (Somniosus microcephalus) are known from the Canadian Arctic eastward, including Greenland and Svalbard through to the Laptev Sea (Chernova et al. 2015), south to Iceland and Nova Scotia on the eastern seaboard of Canada (MacNeil et al., 2012) in the western Atlantic and south to the United Kingdom and southernmost Norway (in deepwater areas). Pacific sleeper sharks (Somniosus pacificus), on the other hand, range from the Bering Sea south throughout the deep waters of the Pacific and into the Southern Ocean (Yano et al., 2004). The smaller Somniosus rostratus is reported from the Mediterranean Sea and the Northeastern temperate Atlantic Ocean, while a fourth species, Somniosus antarcticus has been described from the Southern Ocean (Yano et al., 2004). Current species delimitations, notwithstanding accurate species identification and population delineation for the larger somniosids, (S. microcephlaus, S. pacificus, and S. antarcticus) remain challenging given their overall morphological similarities and the fact that even the few “distinguishing” characteristics for each species often considerably overlap (Benz et al., 2007; MacNeil et al., 2012).

Using partial mitochondrial (mtDNA) cytochrome b (cyt b) gene sequences, Murray et al. (2008) demonstrated that Pacific sleeper sharks and Greenland sharks can be distinguished from one another, with the former exhibiting no geographic structure among samples collected across the Pacific Ocean or more southerly latitudes. More recently, genetic characterization of juvenile Greenland sharks from the Canadian Arctic revealed individuals carrying Pacific sleeper shark mtDNA haplotypes (Hussey et al., 2015). This finding, taken with a growing list of cytonuclear discordance and/or genetic admixture among closely related elasmobranchs in the literature, suggests that hybridization among Somniosus species is possible (Hinojosa‐Alvarez, Walter, Diaz‐Jaimes, Galván‐Magaña, & Paig‐Tran, 2016; Kashiwagi, Marshall, Bennett, & Ovenden, 2012; Morgan et al., 2012; Walter et al., 2014). Questions remain whether this potential hybridization is recently developed among once more distinguishable species, whether it results from long‐term gene flow between sister species, or whether the shared genetic signatures are a product of incomplete lineage sorting (retained ancestral polymorphism).

In this work, we focus on genetic variation in the Greenland shark across its current known distribution. Our objective was to compare patterns of DNA sequence polymorphism and divergence within and among Greenland sharks sampled throughout the Arctic, and sub‐Arctic‐temperate regions, with those of the closely related Pacific sleeper shark, and other deepwater somniosids in order to develop a more comprehensive picture of the origins of the Greenland shark.

METHODS

Tissue samples

Tissue samples from 247 somniosid sharks were collected across a wide distribution including the Canadian Arctic and sub‐Arctic, Western Greenland, Svalbard, Iceland, and the North Atlantic (Fig. 1) that were identified as Greenland shark (i.e., S. microcephalus) during sample collection. Pacific sleeper shark (S. pacificus; n = 6) tissue samples were collected from Alaska. Additional putative S. microcephalus individuals (e.g., identified as S. microcephalus in the field) were collected from the Mid‐Atlantic Ridge (MAR‐Azores; n = 1) and the Gulf of Mexico (GM; n = 1), and published Somniosus mtDNA sequences (n = 22) from Murray et al. (2008) were retrieved from GenBank (EF090943–EF090963) and included for comparative purposes (n = 277 Somniosus spp. samples total).

Figure 1

Sampling locations of somniosid sharks in the Arctic, North Atlantic Ocean, and North Pacific: Somniosus microcephalus (green); Somniosus pacificus (blue) and; cytonuclear discordant individuals (yellow). Numbers inside circles are sample sizes

Molecular methods

Genomic DNA was extracted from all samples using either the standard phenol–chloroform method or the Promega Wizard Extraction kit following the manufacturer's instructions (Promega Corp, Madison WI, USA). We partially sequenced the mtDNA cytochrome b (cyt b) gene, the nuclear recombination activation gene (RAG1), and the nuclear internal transcribed spacer 2 (ITS2) with the flanking 5.8 and 28 ribosomal sequences (5.8S and 28S) herein referred to as the ITS2 locus, for each sampled individual. PCR amplification and sequencing reaction conditions for the markers used are available in the Appendix S1.

Data analysis

All derived mtDNA sequences were aligned with published GenBank sequences in SEQUENCHER 5.0 (GeneCodes, Ann Arbor, MI) and trimmed to 702 bp. Nuclear sequences were also separately aligned in SEQUENCHER and trimmed to 659 bp and 1,085 bp for RAG1 and ITS2, respectively. To estimate genealogical relationships among haplotypes, a 95% statistical parsimony network was constructed in TCS 1.21 (Clement, Posada, & Crandall, 2000) for each locus. Locus‐specific nucleotide models were determined in jModelTest 2.1.10 based on AIC/BIC model comparisons (Darriba, Taboada, Doallo, & Posada, 2012). Haplotype diversity (Hd) and nucleotide diversity (π), Fu's F statistic, Tajima's D, genetic distance (D a), and tests for recent population expansion from mismatch distributions using cyt b sequences were performed in DnaSP 5.10.01 (Librado & Rozas, 2009). Estimates of genetic differentiation, F ST (frequency‐based) and ΦST (Tamura‐Nei) among S. microcephalus sampling sites were calculated in Arlequin 3.5 (Excoffier, Laval, & Schneider, 2005) using the mtDNA data. A demographic history for S. microcephalus based on the cyt b dataset (using HKY+I model derived from jModelTest results) was reconstructed using the Bayesian Skyline Plot (BSP) model (Drummond, Rambaut, Shapiro, & Pybus, 2005) as implemented in the BEAST 1.8.0 package (Drummond, Suchard, Xie, & Rambaut, 2012). Four independent runs were performed, each using 100 mol/L generations and a sampling frequency of 5,000. Each run also used a piecewise‐constant skyline model, with 20 groups, and an uncorrelated lognormal relaxed‐clock mean rate (ucld.mean) was set to 0.0007 = 7 × 10−10 Ma−1 as per Martin, Naylor, and Palumbi (1992). Replicates were combined using LogCombiner (BEAST package), applying a burn‐in of 50% for each replicate. Finally, the BSP was visualized in the program Tracer 1.5 (Rambaut & Drummond, 2007) to detect changes in population size over time.

Coalescent analysis and divergence estimation

To estimate long‐term demographic parameters, namely N e, migration rates (M; i.e., measures of gene flow) and the time since divergence (t 0) between the Greenland shark and the Pacific sleeper shark, we used Isolation with Migration (IM) models as implemented in IMa2 (Hey, 2007, 2010a). IM models search parameter space for the most likely genetically based genealogy describing the split between two taxonomic groups, using a Bayesian framework. Estimated genealogies are generated assuming random mating within groups and that neither group exchanges genetic materials with other nonsampled groups (Hey, 2010a; Hey & Nielsen, 2007). Generated genealogies specify branch lengths and group‐specific genetic variation (estimated as θ) from which demographic parameters can be computed. IM models were run using the mtDNA cyt b data, with the HKY model of sequence evolution. An initial hypothesized phylogeny between the two species was provided, and uniform priors were set regarding the distribution of each parameter of interest using preliminary runs. Three different substitution rates, each characterizing divergence among elasmobranchs: a conservative estimate (7 × 10−10 Ma−1; Martin et al., 1992), Squalus (9.90 × 10−10 Ma−1; Winchell, Martin, & Mallatt, 2004), and Squatina (1.67 × 10−9 Ma−1; Winchell et al., 2004) at cyt b were used to estimate cyt b substitution rates in this work. These substitution rates respectively range from very conservative (Martin & Palumbi, 1993; Martin et al., 1992) to those possibly more reflective of somniosids (e.g., rates derived for Squalus and Squatina), respectively. Demographic parameter searches were performed with 20 Markov chain Monte Carlo (MCMC) simulations using a geometric heating scheme to ensure adequate mixing among them. Initial searches were set to run between 1 x 106 and 12 x 106 iterations in order for the genealogies to reach a stationary distribution from which to sample. Once stationarity was reached, genealogies were permutated 10 x 106 times and sampled every 100th step for a total of 100,000 recorded genealogies from which demographic parameters were calculated (see IMa2 manual; Hey, 2010a). The different searches were repeated three times, and all derived parameters were compared among all three iterations to evaluate consistency (Hey, 2010a). All genealogies recovered from each search were used to compare the information content and thus likelihood (L‐mode) of different IM models. To check for consistency in estimated parameters, model comparison with searches omitting some parameters was also run by specifically setting migration rate priors to zero, hereafter referred to as Isolation Only (ISO). ISO models prevent any migration rate influences on estimated demographic parameters and generate parameter estimates consistent with an incomplete lineage sorting (ILS) scenario.

Due to low levels of polymorphism at nuclear loci, IM analyses were not performed for the RAG1 and ITS2 datasets. Nuclear nucleotide substitution rates for RAG1 (4.85 × 10−10 Ma−1; Eytan & Hellberg, 2010) and ITS2 (1.50 × 10−9 Ma−1; Montoya‐Burgos, 2003) were also used to estimate RAG1 and ITS2 divergence times using a simple formula whereby the number of observed mutations at a locus is divided by the estimated gene‐wide percent divergence rate per Ma.

RESULTS

A total of 61 mtDNA haplotypes were recovered across all sleeper shark samples forming two distinct and well‐differentiated mtDNA haplogroups (Fig. 2, Tables S1 and S2). The Greenland shark samples consisted of 45 haplotypes, including 30 that are novel to this study (green circles with black text Fig. 2; see Table S2). These samples showed a haplotype diversity (Hd) of 0.7100 and nucleotide diversity (π) of 0.0028, while those recovered for the Pacific sleeper sharks were 0.9700 and 0.0053, for Hd and π, respectively (Table S2). A total of 18 haplotypes were among the Pacific sleeper shark samples, including four novel to this study. Site‐by‐site statistics are reported in Table 1. DNA sequences are deposited in GenBank under the accession numbers MF537353‐MF537393 and MF555591‐MF555594.

Figure 2

Haplotype variation at cyt b, ITS2, and RAG1 loci among sampled somniosids as depicted through a 95% statistically parsimony network. For cyt b, the size of the circles is proportional to the frequency of each haplotype in the 277 individual dataset, ITS2 and RAG1 haplotypes are not scaled for frequency. Haplotype numbers correspond to Table S1, white numbers are previously reported haplotypes, black numbers are novel haplotypes to this study. Individuals identified in the field as Somniosus microcephalus possess haplotypes in green, while those identified as Somniosus pacificus possess blue haplotypes

Table 1

Description of collection and mitochondrial sequence genetics‐based data for Somniosids included in the present study. n, sample size; S, number of polymorphic sites; h, number of haplotypes; Hd, haplotype diversity; k, average number of sequence pairwise differences; π, nucleotide diversity. Sample sites are AK, Alaska; CS, Cumberland Sound; DI, Disko Island; GF, Grice Fjord; IC, Iceland; KP, Kakiak Point; MB, Maxwell Bay; NS, Nova Scotia; RB, Resolute Bay; SI, Scott Inlet; SV, Svalbard; GM, Gulf of Mexico; MA, Mid‐Atlantic Ridge; GB, GenBank

Site	n	S	h	Hd	k	π	Admixtures
AK	6	8	5	0.933	3.467	0.00494	0
CS	74	33	20	0.805	2.422	0.00345	2
DI	3	6	3	1.000	4	0.00570	0
GF	3	2	3	1.000	1.333	0.00190	0
IC	25	10	10	0.720	1.2	0.00171	0
KP	2	2	2	1.000	2	0.00285	0
MB	24	21	6	0.543	2.032	0.00290	1
NS	4	14	3	0.833	7.167	0.01021	1
RB	8	11	6	0.893	30.36	0.00432	0
SI	75	27	17	0.764	2.282	0.00883	3
SV	29	11	10	0.695	1.596	0.00227	0
GM	1	—	1	—	—	—	1
MA	1	—	1	—	—	—	—
GB	22

Only two alleles were recovered from each of the nuclear markers, with these alleles characterized by mostly species‐specific single nucleotide polymorphisms (SNPs). For the RAG1 marker, a “C” was predominant at position 607 in identified Pacific sleeper sharks, while a “T” was prevalent in Greenland sharks. Four polymorphic positions were found in 1,085 bases of the ITS2 sequences, a single indel and three SNPs. Greenland sharks possessed an “A” insertion at position 67. For the SNPs, Greenland sharks displayed an “A” at position 72, a “G” at position 142, and an “A” at position 507; while Pacific sleeper sharks possessed a “G” at position 72, a “C” at position 142, and a “C” at position 507.

Some sampled individuals reported as Greenland sharks in the field, and subsequently genotyped, showed a mixture of both Greenland shark and Pacific sleeper shark characteristic haplotypes/alleles at various loci. Six sharks sampled in the Canadian Arctic, exhibited a Pacific sleeper shark cyt b haplotype but were homozygous for the Greenland shark alleles at both the ITS2 and RAG1 loci. An individual sampled in the temperate western Atlantic (i.e., “Nova Scotia” shark) exhibited a cyt b haplotype of a Pacific sleeper shark, homozygosity at ITS2 for the Pacific sleeper shark allele, but was heterozygous at RAG1, showing alleles for both Pacific sleeper and Greenland sharks at the RAG1 locus. Similarly, an individual sampled from the Gulf of Mexico showed a Pacific sleeper cyt b haplotype and homozygosity for the Pacific ITS2 allele, but also displayed heterozygosity for the Greenland shark and Pacific sleeper shark alleles at RAG1. Thus, although genetic sequences at the sampled loci tended to separate the somniosid species, evidence of genetic admixture between the two species was clear.

Demographic expansion

Tests for population demographic changes in Greenland sharks, performed with and without hybrids all corroborated the likely recent demographic expansion of the species. Genetic distance (D a) between Greenland sharks (including hybrids) and Pacific sleeper sharks was 0.01485. Tests for population demographic changes, with and without hybrid individuals included, showed significant negative Fu's F and Tajima's D values (F hybrids = −35.903, p < .0001, F no hybrids = −45.518, p < .01, D hybrids = −1.95, p < .05, D no hybrids = −2.19, p < .01, respectively), consistent with recent population expansion. Mismatch distribution analysis, outlining changes to species nucleotide diversity over time and effective population size (N e), also indicated a likely recent population expansion (Fig. 3a), while the BSP revealed a long period of population stability followed by a recent demographic expansion beginning around 0.5 Ma ago (Fig. 3b).

Figure 3

Genetic indicators of population expansion using the mtDNA cyt b data for Somniosus microcephalus individuals (including hybrids). (a) Mismatch distribution of pairwise nucleotide differences (dark gray columns) with expected curve fit to a model of population expansion (gray line); (b) Bayesian skyline plot for S. microcephalus using the HKY+I nucleotide substitution model

There was no significant genetic differentiation of S. microcephalus among geographic sampling sites with the exception of Cumberland Sound versus Maxwell Bay (F ST) and Maxwell Bay versus Nova Scotia (ΦST, Table 2).

Table 2

Measures of genetic differentiation among Somniosus sampling sites. ΦST values (Tamura‐Nei; below diagonal), conventional F ST based on haplotype frequency (above diagonal). All values non‐significant (p > .05), except in bold

Site	CS	DI	SV	GF	KP	IC	MB	NS	RB	SI
CS		0.0306	0.0111	−0.0800	−0.0664	−0.0152	0.0380	−0.0147	−0.0357	−0.0031
DI	0.0252		0.1661	−0.1250	0.0000	0.0698	0.2944	0.0924	−0.0193	0.0805
SV	−0.0090	0.1188		−0.0099	−0.0757	−0.0046	−0.0005	−0.0357	−0.0119	−0.0055
GF	−0.1161	−0.0461	−0.0784		−0.2000	−0.1009	0.0329	−0.0909	−0.1234	−0.0411
KP	−0.1297	−0.2026	−0.0299	−0.1999		−0.0855	0.0069	−0.1915	−0.1401	−0.0704
IC	−0.0072	0.1672	0.0085	−0.0966	−0.0713		0.0098	−0.0335	−0.0404	−0.0059
MB	−0.0089	0.1165	−0.0084	−0.1077	−0.0809	0.0139		0.0164	0.0400	0.0271
NS	0.1517	−0.0578	0.2396	−0.1484	−0.1786	0.2394	0.1373		−0.0667	−0.0319
RB	−0.0290	−0.0963	−0.0114	−0.1468	−0.2019	−0.0018	−0.0102	0.0438		−0.0303
SI	−0.0077	−0.0259	−0.0050	−0.1177	−0.1382	0.0028	−0.0030	0.0983	−0.0361

Coalescent IM versus ILS analysis and divergence estimates

Across all isolation with migration (IM) models assessing coalescent‐based demographic parameter, strongest support was recovered from models with unequal migration rates between the species. The best models were characterized by coalescent migration rates equal to zero from Pacific sleeper sharks into Greenland sharks, while Greenland sharks into Pacific sleeper sharks migration rates were greater than zero. The best model also included varying long‐term effective population sizes for extant and ancestral populations, following AIC model selection tests. Following conversion of coalescent rates to “forward‐time,” results indicate asymmetrical migration whereby gene flow likely only occurs from Pacific sleeper sharks into Greenland sharks since the species diverged from one another.

mtDNA (cyt b) analysis

Simulations of IM models (m prior upper bound = 0.2) and ISO (m = 0) models produced similar maximum‐likelihood estimates (hereafter MLEs) of species‐split time (t 0) of 1.190 (Fig. 4a) and 1.130, respectively (see Table 3 for comparative split times in years). For the IM analyses, a second older and less probable peak was identified for t 0 at 7.550 followed by a flat nonzero right‐end tail (Fig. 4a). Although overall migration rate parameters (forward‐time) were not significant following the Nielsen and Wakeley (2001) log‐likelihood ratio test, the population‐specific migration rate (2Nm; also forward‐time) for Pacific sleeper into Greenland shark was significantly different from zero (p < .05). Marginal MLE of migration were asymmetrical with m PS>GS 0.01330 and m GS>PS 0.00010 (Fig. 4b). Marginal unscaled estimates of q parameters and scaled N e are presented in Table 3. Although generation times have not been reported for either S. microcephalus or S. pacificus, we applied a relatively short estimate based on a 20‐year age at first reproduction hypothesized for another squaliform species Centroselachus crepidater (Irvine, Stevens, & Laurenson, 2006; Veríssimo, McDowell, & Graves, 2011), and an approximated 150‐year generation length following the ~156 age at first reproduction hypothesized by Nielsen et al. (2016). The longer generation time scaled N e estimates were approximately 13% smaller than the short‐generation time estimates (Table 3).

Figure 4

Marginal posterior probability density distributions of species‐split time (t 0; panel a), migration rates (m; panel b) and long‐term effective population size (N e; panel c) as estimated in IMa2 using mtDNA (cyt b). Migration directions are presented forward in time. PPD, posterior probability density

Table 3

(a) Time‐splitting estimates (t 0) obtained from IMa2 for Isolation with Migration (IM) and Isolation only (ISO) models unscaled and scaled in millions years among Somniosus microcephalus and Somniosus pacificus, using published mtDNA substitution rates: mtDNA Conservative (Martin et al., 1992), mtDNA Squalus (Winchell et al., 2004), and mtDNA Squatina (Winchell et al., 2004). Smi, Somniosus microcephalus; Spa, Somniosus pacificus. (b) Comparative long‐term N e estimates based on generation times of 20 years (N e20) and 150 years (N e150)

Unscaled		Conservative7.00 × 10−10		Squalus rate9.90 × 10−10		Squatina rate1.67 × 10−9
IM	ISO	IM	ISO	IM	ISO	IM	ISO
(a) Time‐splitting t 0
1.15	1.13	2.34	2.30	1.65	1.63	0.98	0.96

Nuclear analysis

Nuclear divergence estimates between species were comparable to the cyt b coalescent analysis. An estimated split time of 3.16 Ma was found based on the single mutation difference between Greenland sharks and Pacific sleeper shark RAG1 alleles; while the ITS2 locus estimate was slightly younger at 2.52 Ma, based on the four mutation steps between both shark species.

DISCUSSION

Our genetic analyses clearly demonstrate the existence of two genetically differentiated mtDNA and nDNA lineages of Somniosus in the Northern hemisphere, largely consistent with the current recognized species ranges (MacNeil et al., 2012). However, cytonuclear discordance (initially reported among juvenile S. microcephalus in the Canadian Arctic, see Hussey et al., 2015) was present in eight individuals in this study; seven individuals from the known S. microcephalus distribution (MacNeil et al., 2012) and in an additional individual sampled well outside the known distribution (i.e., Gulf of Mexico). These data provide evidence for at least three classes of genetically admixed (i.e., hybrid) individuals that account for roughly 2% of the somniosid individuals in this study. Despite our dataset representing the largest sample size of Greenland sharks genetically surveyed to date, the elevated number of rare mtDNA haplotypes and the few individuals representing them, along with our failure to recover two previously reported Greenland shark haplotypes (Murray et al., 2008), suggests we have yet to capture the full extant genetic diversity of Greenland sharks. The hidden genetic diversity of Greenland shark notwithstanding the observed cyt b variation (i.e., high number of recovered haplotypes and the star‐like network pattern, mismatch distribution, and BSP) is consistent with a signal of demographic expansion following divergence from the ancestral somniosid. No obvious spatial genetic variation was observed (Table 2); however, the conspicuous Atlantic residency of all individuals carrying the pure‐type Greenland shark genetic signatures compared with the extant Pacific sleeper shark (and ancestral) lineages suggests speciation within, and likely adaptation to, the north Atlantic and/or Arctic environment.

The mixed genetic signatures among individuals identified in the field as Greenland shark could result from two possible processes: (1) ongoing interspecific hybridization (introgression upon secondary contact), or (2) the retention of ancestral DNA polymorphisms at the observed loci (incomplete lineage sorting—ILS). To attempt to distinguish between these two processes we applied IMa2 analyses on the mtDNA data. All IMa2 model‐fitting analysis indicated higher support for models including gene flow as the initial species split, and in all cases, this gene flow was asymmetrical with alleles originating in Pacific sleeper shark and migrating into Greenland shark. These results are consistent with a scenario of introgressive hybridization. This could explain the occurrence of Pacific sleeper shark mtDNA in Greenland sharks sampled in Scott Inlet in the Canadian Arctic (previously reported by Hussey et al., 2015) and found in other Arctic and sub‐Arctic areas such as Cumberland Sound and Maxwell Bay. The occurrence of admixed sharks in temperate and subtropical waters begs the question as to the source and origin of the Pacific sleeper shark alleles found in Greenland shark. While the Nova Scotia hybrid was caught within the distribution recognized for Greenland shark (MacNeil et al., 2012), the Gulf of Mexico shark was captured well outside this range, in the subtropical North Atlantic region. Moreover, the MAR‐Azores sample represents the first confirmed report of a complete Pacific sleeper shark genetic signature in the eastern Atlantic Ocean. These new findings suggest that DNA sequences currently reported as Pacific sleeper shark genetic signatures are not exclusive to the Pacific Ocean and that the common ancestor to both sharks likely exhibited a pan‐oceanic, deepwater distribution. Such a finding would also suggest that the divergence between the two species is, therefore, not likely linked to a single vicariant speciation event such as the rise of the Isthmus of Panama (see e.g., Bacon et al., 2015 and Marko, Eytan, & Knowlton, 2015 for a recent discussion). From the fossil record, the earliest appearance of sleeper sharks date to approximately 100 Ma ago with fossil somniosids found in Miocene deposits from Belgium and Italy (Kriwet & Klug, 2009). This suggests the Atlantic presence of somniosids prior to late Miocene cooling and subsequent onset of the Pleistocene glacial periods.

Estimates of species‐split (divergence) time for Greenland shark and Pacific sleeper shark using mtDNA did not show a consistent time line for IM and ILS models across a range of substitution rates. The observed bimodal t 0 curve, however, is consistent with a “two attractors” hypothesis described as “a recent signal of no migration, and an older signal with higher migration” (J. Hey as cited in Choleva et al., 2014). Application of a broadly based conservative substitution rate (Martin et al., 1992) and more probable substitution rates (see Table 3) identified that an IM‐based estimated species‐split time using the cyt b locus occurred during the Quaternary period, roughly 0.98 Ma ago to 2.34 Ma ago. This time frame is consistent with the onset of the Pleistocene epoch, which may have played an initiating role in the speciation trajectories of extant sleeper sharks in the Northern hemisphere. Analysis of nuclear data alone, though exhibiting only low diversity and differentiation, also supports a more recent split. Analyses not accounting for interspecific hybridization (e.g., ILS runs, and BSP) further support more recent processes (i.e., Pleistocene) as exerting demographic impacts including increased N e as suggested by the BSP (Fig. 3b) around 0.5 Ma.

Ice‐olation with migration hypothesis

We introduce the hypothesis of oscillatory vicariance with secondary contact—herein referred to as Ice‐olation with Migration. If we consider that the ancestral sleeper shark was a deepwater species that was constrained to this habitat via temperature regimes, this species likely occupied a pan‐oceanic distribution prior to the late Miocene, but only in the deep ocean. This is also consistent with the occurrence of somniosids in tropical and temperate regions where they are typically observed at depths greater than 1,000 m and reportedly down to 2,200 m (Benz et al., 2007; Campana, Fisk, & Klimley, 2015; Herdendorf & Berra, 1995), and in contrast to somniosids sampled in the Arctic and northern Quebec, where satellite telemetry data have shown sharks occurring from the surface down to approximately 1,500 m (Campana et al., 2015; Fisk, Lydersen, & Kovacs, 2012; Stokesbury, Harvey‐Clark, Gallant, Block, & Myers, 2005).

From the late Miocene, and into the Quaternary, a drastic reduction in global temperatures occurred. This decline in sea temperature resulted in sea ice formation in polar regions, including the formation of ice shelves up to 1 km in thickness (Fastook & Hughes, 2013; Grosswald & Hughes, 1999; Polyak, Edwards, Coakley, & Jakobsson, 2001). This ice thickness coupled with the submerged mountainous terrain of the Arctic, may have effectively constrained population connectivity of ancestral somniosids between the Arctic sub‐basins/marginal seas and the remaining pan‐oceanic distribution. In support of the Arctic sub‐basin isolation, glacial scouring patterns along the Lomonosov Ridge are consistent with ice shelf formations (Polyak et al., 2001; Stein et al., 2016), and it is likely that similar constraints to connectivity occurred among the Arctic and Atlantic Ocean at the East Greenland Rift Basin (Vigdorchik, 1979; Vigdorchik & Vinkovetsky, 1976). The alternating advance and retreat glaciation stages, and repeated ice shelf formation and grounding on Arctic shelves would have yielded transient periods of isolation—typically 50–100 k years in duration—long enough to initiate genetic divergence, but likely not long enough to facilitate complete reproductive isolation. Subsequent interglacial periods restored connectivity among basins. Similar processes have been proposed regarding the 100 k year scale of isolation of marine fishes across the Sunda Shelf followed by subsequent connectivity (Reece, Bowen, Smith, & Larson, 2010).

During cooling periods, polar regions may have also exhibited a more uniform thermal profile than the stratified profiles of the temperate and tropical regions. Combined, these cooler temperatures and the uniformity of vertical temperature profiles in the Arctic may have facilitated an expansion in available habitat for sleeper sharks, from deep waters to shallow coastal habitats. Habitat expansion would have concurrently provided access to a diverse array of prey, including marine mammals. The demographic expansion suggested by the mismatch analysis and the BSP, and the more than doubling of long‐term N e estimates (from IMa2 analyses) from ancestral sharks to the current species is consistent with the expansion of available habitat and prey resources for somniosid species in the Arctic region.

Population expansion and adaptation in the Arctic

Ecological specialization has been proposed as a driver for marine speciation (Corrigan & Beheregaray, 2009; Kashiwagi et al., 2012; Morin et al., 2010; Wang, Chou, & White,1999) and is a plausible explanation for the demographic expansion of sleeper sharks in the Arctic region. Greenland sharks are considered generalists and opportunistic predators feeding on a broad range of benthic and pelagic teleost and cephalopod prey (MacNeil et al., 2012). Marine mammals including ringed‐seal (Pusa hispida) have been recovered from the stomach contents of Arctic Greenland sharks (Fisk, Tittlemier, Pranschke, & Norstrom, 2002; Leclerc et al., 2012; Lydersen, Fisk, & Kovacs, 2016; McMeans, Svavarsson, Dennard, & Fisk, 2010; Nielsen, Hedeholm, Simon, & Steffensen, 2014; Ridoux, Hall, Steingrimsson, & Olafsson, 1998), as well as in the stomachs of Pacific sleeper sharks at high latitudes (Bright, 1959; Horning & Mellish, 2014; Hulbert, Sigler, & Lunsford, 2001), identifying both species as high trophic level feeders (Hussey et al., 2014). When considering the global distribution of shark and pinniped diversity, an inverse relationship occurs, with the highest shark diversity at tropical latitudes where pinniped diversity is low, and few large shark species, including the Greenland and Pacific sleeper shark, occurring at high latitudes where pinniped diversity is highest (Ferguson, Higdon, Tallman, Fisk, & Hussey, 2014). The diversity and abundance of marine mammals in the Arctic might have served as novel food sources for sleeper sharks during their colonization of the Arctic environment.

CONCLUSION

Our results show that the nearly morphologically identical sister species Greenland shark and Pacific sleeper shark can be genetically distinguished using a combination of mtDNA and nuclear markers, but interspecific hybridization (introgression) between these two species is probable and has likely been ongoing since initial divergence. Our results confirm the fourth reported instance of interspecific hybridization among elasmobranchs (Hussey et al., 2015; see also Morgan et al., 2012; Walter et al., 2014; Cruz et al., 2015), which is likely more common among this group of fishes than previously thought. Our findings, considering a range of substitution rates for the genetic markers employed here, are consistent with a more recent Greenland shark and Pacific sleeper shark speciation event occurring between 1 and 2.34 Ma, more recently than the closure of the Isthmus of Panama. Our estimates of long‐term N es also indicate slightly higher sizes for Greenland sharks compared to Pacific sleeper sharks. However, recent radiocarbon analysis of S. microcephalus eye lenses has concluded that Greenland sharks may not reach maturity until 150 years of age, and possibly live as long as 400 years (Nielsen et al., 2016). If accurate, this would mean that effective population sizes are but a fraction of our 20‐year estimates, and using the N e estimates from the 150‐year generation time of Nielsen et al. (2016) is likely more appropriate (Table 3). The lack of geographic structuring for Greenland sharks across the Canadian Arctic and beyond, but with the recurring hybridization with Pacific sleeper sharks at the edge of its distribution in warmer waters, suggests that a continuing warming Arctic may further reduce the genetic integrity of both species. Future research focused on broad‐scale population‐level connectivity across the range of S. microcephalus is required to clarify the consequences of past climatic and present‐day environmental changes on the hybridization rates and evolutionary trajectory of these species.

CONFLICT OF INTEREST

None declared.

Click here for additional data file.

Click here for additional data file.