Migratory culture, population structure and stock identity in North Pacific beluga whales (Delphinapterus leucas)
Migratory culture, population structure and stock identity in North Pacific beluga whales (Delphinapterus leucas)
Greg O'Corry-Crowe 0 1
Robert Suydam 1
Lori Quakenbush 1
Brooke Potgieter 0 1
Lois Harwood 1
Dennis Litovka 1
Tatiana Ferrer 0 1
John Citta 1
Vladimir Burkanov 1
Kathy Frost 1
Barbara Mahoney 1
0 Harbor Branch Oceanographic Institute, Florida Atlantic University , Fort Pierce , Florida, United States of America, 2 North Slope Borough Department of Wildlife Management, Barrow, Alaska, United States of America, 3 Alaska Department of Fish and Game, Fairbanks, Alaska, United States of America, 4 Fisheries and Oceans Canada , Yellowknife, Northwest Territories , Canada , 5 Marine Mammal Laboratory, ChukotTINRO, Anadyr, Chukotka, Russia, 6 North Pacific Wildlife Consulting , Marine Mammal Laboratory, Seattle, Washington, United States of America, 7 University of Alaska, School of Fisheries and Ocean Science, Kailua Kona, Hawaii, United States of America, 8 National Marine Fisheries Service , Anchorage, Alaska , United States of America
1 Editor: Songhai Li, Sanya Institute of Deep-sea Science and Engineering Chinese Academy of Sciences , CHINA
The annual return of beluga whales, Delphinapterus leucas, to traditional seasonal locations across the Arctic may involve migratory culture, while the convergence of discrete summering aggregations on common wintering grounds may facilitate outbreeding. Natal philopatry and cultural inheritance, however, has been difficult to assess as earlier studies were of too short a duration, while genetic analyses of breeding patterns, especially across the beluga's Pacific range, have been hampered by inadequate sampling and sparse information on wintering areas. Using a much expanded sample and genetic marker set comprising 1,647 whales, spanning more than two decades and encompassing all major coastal summering aggregations in the Pacific Ocean, we found evolutionary-level divergence among three geographic regions: the Gulf of Alaska, the BeringChukchi-Beaufort Seas, and the Sea of Okhotsk (Φst = 0.11±0.32, Rst = 0.09±0.13), and likely demographic independence of (Fst-mtDNA = 0.02±0.66), and in many cases limited gene flow (Fst-nDNA = 0.0±0.02; K = 5±6) among, summering groups within regions. Assignment tests identified few immigrants within summering aggregations, linked migrating groups to specific summering areas, and found that some migratory corridors comprise whales from multiple subpopulations (PBAYES = 0.31:0.69). Further, dispersal is male-biased and substantial numbers of closely related whales congregate together at coastal summering areas. Stable patterns of heterogeneity between areas and consistently high proportions (~20%) of close kin (including parent-offspring) sampled up to 20 years apart within areas (G = 0.2±2.9, p>0.5) is the first direct evidence of natal philopatry to migration destinations in belugas. Using recent satellite telemetry findings on belugas we found that the spatial proximity of winter ranges has a greater influence on the degree of both individual and genetic exchange than summer ranges (rwinter-Fst-mtDNA = 0.9, rsummer-Fst-nDNA = 0.1). These findings indicate widespread natal philopatry to summering
available on our website: http://pbbegenetics.
Funding: A research grant was provided by the
Alaska Beluga Whale Committee and administered
by the North Slope Borough Department of Wildlife
departments/wildlife-management/comanagement-organizations/alaskabeluga-whalecommittee). Funding was also provided by the
National Fish and Wildlife Foundation,
grant#20080105-000, with additional support from Harbor
Branch Oceanographic Institute (HBOI) at Florida
Atlantic University and the HBOI Foundation. In
kind support was provided by the National Marine
aggregation and entire migratory circuits, and provide compelling evidence that migratory
culture and kinship helps maintain demographically discrete beluga stocks that can
overlap in time and space.
In migratory species social learning, seasonal movements and the use of geographically
separate habitats during the annual cycle can foster the cultural inheritance of migration routes
and complex patterns of dispersal, gene flow and population structure [
1, 2, 3
]. These, in turn,
have implications for gene-culture coevolution  and create novel challenges for
management and conservation, including the identification of management units, the assignment of
migrating animals to population of origin, and the assessment of risk at the population level.
The primary factors influencing population subdivision in migratory species include; (1)
Non-uniform patterns of seasonal resource distribution which can foster philopatry to
geographically discrete migration destinations (e.g., feeding or breeding grounds) that can
promote population divergence over time [
]. (2) Such population units, however, may overlap
during migration or co-occur at another migration destination or staging area for part of the
year, facilitating dispersal and interbreeding. (3) In species with extended periods of parental
care natal philopatry and thus population divergence may involve parentally directed learning
of migration routes termed migratory culture [1, 5±7]. (4) Structuring in migratory species
also depends on which sex disperses [
] with sex-biased dispersal potentially leading to the
demographic isolation of population units even in the face of gene flow and seasonal overlap
]. (5) Patterns of dispersal and heterogeneity in migratory species can be shaped by
geography and other environmental factors (e.g., water temperature, weather patterns, sea-ice cover)
that directly dictate the path and timing of migration, for example [
], and thus the level and
duration of population overlap. And finally (6), glacial history, especially at higher latitudes,
where the historical location of refugial populations and the sequence in which new habitat
emerged following deglaciation can have a lasting effect on contemporary patterns of
distribution and migration, and thus dispersal and population structure .Disentangling the
contributions of such factors to migration, dispersal and structuring in migratory species is
challenging. Genetic investigations offer huge potential in resolving the interplay between
behaviors, the environment, and demographic history in shaping seasonal movements,
dispersal and population structure in these species.
The beluga whale, Delphinapterus leucas, is a mid-sized cetacean of Arctic and sub-Arctic
waters where all these factors likely interplay to shape contemporary migration, dispersal,
breeding, and population structure patterns. In this study we use spatial and temporal patterns
of genetic variation and individual relatedness to test a series of hypotheses about these aspects
of belugas in the Pacific.
Across much of their range belugas migrate between wintering areas at or near the southern
margin of the sea-ice and summering grounds in seasonally more open water farther north
where they feed, molt and raise their young [
]. Many Arctic populations are migratory
with some completing annual circuits in excess of 6,000km while subarctic populations are less
so, often exhibiting substantial overlap between winter and summer ranges [14±19]. Highly
gregarious, beluga whales congregate by the thousands at several geographically discrete nearshore
locations following ice-breakup in summer [
]. Breeding is believed to occur primarily in
winter and early spring [
], and adult males may select different habitats than females and
2 / 32
younger animals at certain times of year [16, 23±25]. These attributes, combined with the
potential for intersecting migration routes [
] and for separate summering groups sharing the same
wintering ground during the period of peak mating [
], contribute to complex patterns of
dispersal and population structure in beluga whales, raise the possibility of stable migratory
cultures, and create unique challenges for species management, including the identification and
monitoring of stocks.
In the Pacific, beluga whales inhabit three geographically distinct regions; (1) the Gulf of
Alaska, (2) the Sea of Okhotsk, and (3) the Bering, Chukchi and Beaufort Seas, termed the
BCB region (Fig 1). Within the BCB region, there are six major summer aggregation areas
near shore that are labeled in Fig 1 as follows: Bristol Bay (#3) and Norton Sound (#4) in the
eastern Bering Sea, Kotzebue Sound (#5) and off Kasegaluk Lagoon (#6) in the eastern Chukchi
Sea, Mackenzie Delta-Amundsen Gulf region (#7) of the eastern Beaufort Sea, and Anadyr Bay
(#8) in the western Bering Sea (Figs 1 and 2). There are also three discrete summering locations
in the Okhotsk Sea (#9, #10, and #11), and a small resident population in Cook Inlet (#1) and
isolated group in Yakutat Bay (#2) in the Gulf of Alaska (Fig 1).
The available evidence indicates that beluga whales in the BCB region overwinter in the
northern and central Bering Sea at or near the southern extent of the sea ice [18, 28±31] (Fig 2).
The hour-glass shape of the Bering-Chukchi basin constricts the migration of Chukchi and
Beaufort summering groups through the narrow (85km wide) Bering Strait en route to and
from wintering areas in the Bering Sea (Fig 2). This geographic feature, in concert with the
timing and pattern of ice breakup and formation, likely exerts a substantial influence on when
whales migrate, what routes they use, where they overwinter and to what extent separate
summering groups overlap in winter. A recent study summarizing satellite telemetry data from five
of the six BCB summering groups (i.e., excluding Kotzebue Sound) revealed that whales from
discrete summering groups may also have traditional wintering areas with varying degrees of
]. Considering the peak breeding season for beluga whales is during winter and early
], there is potential for extensive genetic and individual exchange among BCB
summering groups both on migration and during the winter.
A series of genetic studies have been conducted on beluga whales in the western Nearctic
(Alaska and northwestern Canada) and in the eastern Palearctic (Russian far east). A number
revealed substantial mtDNA differentiation among geographically isolated populations and
among some summering groups that likely reflects long-established patterns of
female-mediated philopatry and demographic isolation [32±36] leading to their identification as
demographically distinct management stocks [
]. A few studies have documented lower levels
of nDNA (microsatellite) heterogeneity among geographic strata compared to mtDNA which
has been taken as evidence of extensive male-mediated gene flow among summering groups,
possibly on shared wintering grounds [
These earlier studies, however, have a number of limitations with respect to resolving
patterns of migration, dispersal and gene flow of Pacific beluga whales basin-wide. All were of too
short a duration to directly test for natal philopatry and migratory culture in this long-lived
mammal, and to assess the stability of migration and dispersal patterns over ecological time
frames. All had a regional focus that risked incomplete assessments of the molecular and
behavioral ecology of this highly vagile species in an environment with few geographic
barriers. The failure to include all populations, for example, can lead to inaccurate cluster analyses
and erroneous assignment tests. Low marker number and sample size in many cases limited
interpretations of observed patterns of differentiation. Furthermore, differing laboratory
conditions prevented the comparison of nDNA (microsatellite) data sets among studies and
populations, while few whales were sampled on migration or outside traditional coastal summering
3 / 32
Fig 1. Distribution (light blue) of beluga whales, Delphinapterus leucas, in the North Pacific Ocean. The ten major nearshore concentration areas during the summer
months are highlighted (dark blue). These areas along with a small resident group of beluga whales in the Gulf of Alaska are numbered according to Table 1.
In this study, we investigate migration behavior, dispersal, population structure and stock
identity in Pacific beluga whales. We include whales from all major coastal concentration areas
in the north Pacific for the first time. We analyze 1,444 samples for both mtDNA and eight
microsatellite loci. We conduct analyses on a further 203 Russian Far East whales from the
] for a total of 1,647 whales from across their Pacific range. We test hypotheses on
how the spatial proximity of summering and wintering areas may influence levels of dispersal
and interbreeding among discrete summering aggregations, some of the findings of which
have recently been referred to in the tracking paper by Citta et al. [
]. We include previously
unsampled and undersampled locations, greatly expand the analysis of migrating whales, and
extend the sampling time frame in many areas to encompass more than three decades from
1978±2010. With this extended time frame we assess patterns of kinship within groups and
among years, determine the stability of migration patterns and structure over ecological time
frames, and further investigate the population origins of beluga whales in one area (Kotzebue
Sound, Alaska) that has witnessed dramatic changes in the pattern of annual return of beluga
whales in summer. Finally, we explore how model-based cluster analysis of population genetic
structure and Bayesian inference of recent rates of dispersal may be influenced by sample size,
uneven sample numbers, locus number and differing degrees of subdivision.
4 / 32
Fig 2. Migration routes and sampling locations of beluga whales from the major summer concentration areas in the Bering, Chukchi and
Beaufort Seas and from the Gulf of Alaska. Summering and wintering areas, and migration routes were inferred from a combination of satellite
telemetry, aerial and shore based sightings, and Traditional Ecological Knowledge (TEK). Sampling sites are indicated by yellow circles and in the case
of migration by triangles.
Materials and methods
Sample collection and laboratory analysis
All samples were collected under the authority of U.S. Marine Mammal Protection Act permits
issued by the National Marine Fisheries Service, Russian Marine Mammal permits, and the
Dept. of Fisheries and Oceans, Canada scientific licenses. Tissue samples were collected from
1,444 beluga whales at 45 locations across 15 geographic strata in in the Gulf of Alaska, the
Bering, Chukchi and Beaufort Seas, and the Sea of Okhotsk between 1978±2010 (Table 1, Fig
1). A more detailed summary of sample numbers across the four distinct sampling periods of
this study are provided in S1 Appendix. Tissues were preserved in 20% dimethyl sulfoxide
(DMSO) saturated with NaCl. Total DNA was extracted from each sample using standard cell
lysis-high-salt extraction protocols, and each sample was screened for variation within 410bp
of the mtDNA control region according to previously published methods [
]. Total DNA
was also extracted from a number of dried tissue and teeth samples using standard silica-based
5 / 32
E. Bering Sea
E. Chukchi Sea E. Beaufort Sea
ancient DNA methods [
]. The gender of each sample was determined by PCR-based
methods  and all samples were analyzed for variation within eight microsatellite loci (S1
Table) typed on beluga whales [
] or other cetacean species [
]. Earlier studies including
an analysis of 288 belugas in the current study (including 9 known mother-calf pairs),
demonstrated that all eight loci were inherited in a Mendelian fashion and were not sex-linked [
]. Comparisons of observed to expected genotypic frequencies using the Micro-Checker
program (version 2.2.3)  for the four best sampled areas found only one of 32 tests
consistent with null alleles (i.e., homozygote excess), and no evidence of large allelic dropout (0/32
tests) or scoring errors due to stuttering (0/32 tests).
The amount and nature of mtDNA variation within the sequenced region were assessed by
determining the number of variable sites and the number of unique haplotypes using MEGA
] software and by estimating both haplotypic (H) [
] and nucleotide (π) [
using ARLEQUIN 3.5 software [
]. Estimation of genetic diversity (He, Ho and number of
alleles) and probabilities of identity (I) for nuclear loci was performed using CERVUS 3.0 [
Tests for Hardy-Weinberg (H-W) and linkage equilibria were performed in GENEPOP 4.1 [
with P-values estimated from 500,000 iterations of the data using the Mark chain method of
Guo and Thompson [
]. An analysis of both mtDNA and nDNA data revealed that the
populations were not in mutation-drift equilibrium (see [
We used both frequency-based (Fst) and distance-based statistics (Fst, Rst) to calculate the
degree of genetic differentiation among spatial and temporal strata. Fixation indices were
6 / 32
estimated using an analysis of variance framework [
], and homogeneity tests based on
50,000 permutations of the data were performed in ARLEQUIN. The model used for estimating
the evolutionary distances between pairs of mtDNA sequences was simple pairwise differences.
In the case of nuclear loci, genetic distances were based on stepwise mutation models.
Confidence intervals for Fst and related parameters were estimated via 20,000 bootstraps of the
multilocus nDNA data sets for each pairwise comparison in ARLEQUIN. Mantel tests, based on
10,000 permutations of distance data, were conducted in ARLEQUIN to assess the potential role
of geographic distance in shaping the extent of genetic differentiation among geographic
strata. Geographic distances were swim distances between the centroids of summering and
wintering ranges estimated from satellite telemetry data [
15, 16, 25, 31
The Bayesian model-based clustering method, STRUCTURE 2.3.4 [
], was used to infer
population structure and the likely number of populations (K), and to assign individuals to
population of origin. We conducted 28 distinct analyses with unique parameter settings: analyses
were run both with and without admixture, with and without sampling location included as
prior information, and they contained sample sets of varying size that were randomly sampled
from the complete sample set using R (n = 15, 30, 50, 100, and all) that included individuals
successfully genotyped at a range of loci (n 6 loci, n 7 loci, n = all 8 loci, see S1 Supporting
information). For each parameter set we used a burn-in period of 50,000 iterations followed by
2x105 iterations to collect data. Multiple runs (n = 5) were conducted for each value of K (n = 8)
to ensure convergence for a total of 1,120 separate runs. The models of Hubisz et al.  use
sampling location as prior information to reveal further underlying structure without risking
detecting structure that is not present. To identify which individuals were likely immigrants (or
descendants of recent immigrants) to their populations, we defined prior probabilities that each
individual was an immigrant (v = 0.01±0.10) and incorporated information on the geographic
origin of each sample before running the clustering analysis. We used CLUMPAK [
] to integrate
results from multiple runs for different values of K and to implement methods for choosing K.
To further explore the origins of individual whales, we assigned individuals to populations
based on estimated likelihoods of their genotype and/or mtDNA haplotype arising in each of
the sampled populations under assumptions of random assortment of alleles and
independence of loci using the programs DOH and WHICHRUN 4.1 [60±62]. A log10 ratio (LOD) score of
the most likely allocation to the second most likely of 1 (i.e., ratio of 10) denoted
confidence in the assignment to a single population [
Differences in allelic diversity among populations can result in genotypes with higher
calculated likelihoods in populations other than their source (i.e., nominal) population, even though
the genotype may still be relatively common in the source population. In such cases an
exclusion test is also required to ascertain likely population of origin. We thus examined the
estimated likelihoods of genotypes relative to those of other genotypes within each baseline
population. Using the ASSIGNMENT CALCULATOR in the DOH program  and assuming
HardyWeinberg Expectations (HWE), 1,000 new individuals were generated from the gene pools of
each population, and the natural log of the estimated probability of an individual's genotype
was compared to the likelihood distribution of the generated genotypes.
Dispersal patterns over ecological time scales were also investigated using the Bayesian
approach of Wilson and Rannala  in the program BAYESASS 1.3 (see S1 Supporting
information for details). To investigate sex-bias in genetic dispersal (i.e., male-mediated gene flow
on common wintering grounds or migration routes) we used the long-standing demographic
and reproductive isolation between Cook Inlet and Arctic populations ([
], this study) to
distinguish between the influences of Ne, locus choice and gene flow on both marker types by
comparing ratios of mean mtDNA to nDNA differentiation. Sex-bias in actual dispersal (i.e.,
individual transfer) was investigated by comparing levels of differentiation (Fst) at mtDNA
7 / 32
and levels of differentiation and average relatedness at nuclear markers (Fst, Fis, r, mAIc, vAIc)
for males and females. The tests rely on sampling individuals post-dispersal in order to exclude
the homogenizing effects of past immigration and interbreeding by either sex through the
sampling of their offspring [
]. We assumed dispersal occurs primarily at the juvenile and
sub-adult stage and so only adult cohorts were used. The randomization approach of Goudet
et al. [
] was used to test whether there was significant contemporary sex-biased dispersal.
Using FSTAT 2.9.3 we compared observed differences to a null distribution based on 10,000
randomizations of the data.
We applied the Bayesian stock-mixture method of Pella and Masuda [
] to assess the
population composition of migrating herds. This method incorporates uncertainty in the
estimation of stock proportions of the `mixture', or in this case migrant samples, and unlike
conditional maximum likelihood methods improves stock determination by using
information in the mixture sample to improve estimates of baseline relative frequencies. Using the
program BAYES [
], we ran multiple chains with different starting parameters (including different
prior stock proportions) and used the univariate statistic, the shrink factor, to monitor
convergence of chains to the posterior probability. Analyses were conducted separately on mtDNA
and nDNA data, and on mtDNA-and nDNA data combined.
Finally, because philopatry should promote high levels of relatedness among individuals
sampled across years and generations within an area we tested this by comparing estimates of
average relatedness and proportions of first and second-order relationships (based on nDNA)
within and between years using the programs COANCESTRY  and ML-RELATE [
] (See S1
MtDNA and nDNA diversity
A total of 1,444 individual beluga whales from throughout the North Pacific were screened for
sequence variation within 410bp of the mtDNA control region and adjacent proline tRNA
gene (Table 1, Figs 1 and 2). Twenty four variable sites were found, all of which were
substitutions (23 transitions and 1 transversion). A total of 36 unique haplotypes were identified, 12 of
which were represented by a single individual. To reduce sampling bias from non-random
sampling of close kin, only one member from cow-calf pairs where both members were
sampled together was included in subsequent analyses. Overall haplotypic diversity for the
remaining 1,435 whales was high (H = 0.84) due to the large number of rare haplotypes while overall
nucleotide diversity was moderate (π = 0.49%), indicating that the majority of haplotypes were
phyletically closely related. A median joining network of the unique mtDNA haplotypes
reflected this and was characterized by a series of star-like phylogenies with several rarer
haplotypes radiating from a more common central haplotype, a pattern widely interpreted as
indicative of ancient population expansions (see S1 Supporting information, S1 Fig).
A total of 1,116 individual whales were successfully screened for polymorphism at six or
more microsatellite loci (Table 1). Fisher exact tests for HWE (i.e. 8 loci x 9 geographic strata)
found 11 of 72 to be significant at α = 0.05. These deviations were partly due to combining
three areas with low sample size into one Okhotsk area and were not consistent across areas or
loci; involving 7 different loci and 8 separate locations, and no evidence of linkage was found
among the 8 loci when tested across the major concentration areas (see S1 Supporting
information). Average expected heterozygosity (He) for the 8 loci ranged from 0.511 to 0.799 across
the geographic strata (S1 Table). The number of alleles found per locus within an area ranged
from 2 for locus CS415 in Yakutat Bay, Alaska to 16 for locus DlrFCB17 in the Mackenzie
8 / 32
Delta-Amundsen Gulf, Canada, and the probabilities of identity (P(ID)) for each area were on
the order of 10−5 to 10−10.
Genetic differentiation and model-based cluster analysis
Substantial levels of mtDNA differentiation were observed among the three primary regions
that beluga whales inhabit in the North Pacific for both the distance-based (Fst = 0.105±0.321)
and frequency-based (Fst = 0.079±0.285) statistics (Table 2). All six major summer coastal
concentration areas within the BCB region: Bristol Bay, Norton Sound, Kotzebue Sound, Kasegaluk
Lagoon, Mackenzie-Amundsen, and Anadyr Bay were also substantially differentiated from
each other (Fst = 0.106±0.507) (Table 3A), the only exception being the low level of mtDNA
heterogeneity observed between Kotzebue Sound and Mackenzie-Amundsen (Fst = 0.022). The
three primary summer concentration areas within the Sea of Okhotsk also displayed substantial
mtDNA heterogeneity (Fst = 0.104±0.293). Comparing summer concentration areas and
resident populations across the entire north Pacific, Cook Inlet and Bristol Bay were found to be
the most distinct for mtDNA with average Fst values in excess of 0.39 and 0.44, respectively.
Overall, nDNA differentiation was lower than that of mtDNA but it too differed
substantially among regions (Table 2B) and among most summer concentration areas and resident
populations (Table 3B). As with mtDNA, Cook Inlet was the most distinct of all geographic
strata for microsatellite variation (F st = 0.056). No or very low levels of nDNA differentiation
were observed between Norton Sound and Anadyr Bay, between Kotzebue Sound and
Mackenzie-Amundsen and between Anadyr Bay and Mackenzie-Amundsen (Table 3B). Low
sample size precluded analyses of nDNA heterogeneity within the Sea of Okhotsk. In a hierarchical
AMOVA the proportion of total variation due to differentiation among regional groupings
compared to among populations within groups was higher for nDNA than for mtDNA (see S1
There was widespread agreement across the various parameter settings used in the
modelbased cluster analysis (S2 Table). Using no prior information on the geographic origin of
samples and no admixture among strata, the cluster analysis of genotypic data revealed that K = 2
populations was the most consistent with the data. The assignment of individuals revealed that
these two genetically discrete population clusters corresponded to a combined Cook Inlet-Sea
Gulf of Alaska
9 / 32
.222 .220 .312 .043 .314 .581 .980 .031 .101 .170 .020
0 0 0 0 0 0 0 0 0 -0 0
la l- t f
on fo in so
ia th P re
d y dn tta
e b a p
t ed ie l
o n en o
e s k p
b e c m
r a e
re p t
a re Mn
,t e o
s e h s
F ra ts i
,itc ,sn i tea
ittsa itao taum red
u tsr o
re on eh i
e d t
th se d
r a n d
b ,a n
fo , a
s t e A
e s t
lu te ro leb
a y p a
V it e T
. e r
s n e e
le e r e
a g a S
h om 01 .10
a o 0
g h n 2
u r f
o o 88
e ize 19
if u s
c a le om
a v p r
P - f
p m d
in g sa i
) in a r
(B d e
c on it e
p w d
te re taa cea
ll o tr d
te C s o
a .l ly w
s a n t
ro n O ly
c o . h
i g 5 g
m ia .0 u
t d 0 ro
gh e > e
ie th p th
s e :
so vo ed isr
A D e
N (n p o
D ts < re
l R 01 t
ia . e
r d 0 fh
d n :
n a y o
o ) re se
h A g t
o N t
it tD igh itm
m l s
m , e
n ( 1
i t 0 d
h s . e
it F 0 tr
ca ab ah p m
e s m m
d r n o
n a u c
)a ,) , y
A A 0 i
( N 0 n
n a m
S lu n
B sa aL
ay y 7 57 37 63 70 87 71 51 29 91 50 89 03 31
0 1 3 1 2 6 4 h
.000 .007 .-000 .-000 .000 .000 .-000 ,gnC
of Okhotsk stratum and the combined BCB stratum. Allowing for admixture, one population
(K = 1) was the most likely, although the same geographic clustering of Cook-Okhotsk and BCB
was still evident, especially for lower sample sizes (i.e., n = 30, n = 50), confirming that the
primary genetic divergence in the data exists at a larger than regional scale. A re-analysis of the
data incorporating prior information on sampling location detected additional structure. Most
analyses, involving a range of sample sizes (n = 15, 30, 50, 100, all), loci ( 6,
7, 8) and levels of
missing data (i.e., 0% - 4.9%), determined that five discrete population (K = 5) were most likely
given the data (Pr (K|X
Inlet, Bristol Bay, Kasegaluk Lagoon, Mackenzie-Amundsen and the Sea of Okhotsk. Norton
Sound did not form a distinct population cluster, with all individuals assigned to the two
populations corresponding to Mackenzie-Amundsen and Bristol Bay (Fig 3). Likewise, Anadyr Bay
did not form a discrete population cluster in these analyses with individuals assigned to
Fig 3. Summary plots generated in CLUMPAK of model-based cluster analysis of population structure in Pacific beluga whales using
STRUCTURE 2.3.4. The major modes for K = 4 to 6 (based on five separate runs for each value of K) are presented for the analysis using prior
sample group information and no admixture which revealed K = 5 clusters as the most likely (see panel 2). However, in a number of analyses
K = 6 was the most or second-most likely resulting in the separation of Anadyr into a discrete cluster (see panel 3). Each genotyped individual is
represented by a vertical line with estimated membership, Q, in each cluster denoted by different colors. The analysis was based on using all
individuals (n = 1032) scored at 6 or more loci (nloci 6).
12 / 32
Mackenzie-Amundsen and the Sea of Okhotsk. However, Anadyr Bay did form a separate
population in the few analyses where six populations (K = 6) were the most likely (n = 2/28) or in
analyses where K = 6 was the second most likely (S2 Table). This is clearly evident in the
CLUMPAK analyses where summary plots across multiple runs for each value of K reveal clear
clustering of Anadyr in contrast to an absence of a distinct Norton Sound cluster (Fig 3, S2 Fig). In all
analyses Kotzebue Sound was part of the same cluster as Mackenzie-Amundsen.
Mantel tests revealed positive correlations between genetic differentiation and geographic
distance within the BCB region (Fig 4). Stronger correlations with genetic heterogeneity were
found among wintering areas (rmtDNA-Fst = 0.9, p = 0.007) compared to summering areas
(rmtDNA-Fst = 0.1, p = 0.351) and for mtDNA (r = 0.1±0.9) compared to nDNA (r = 0.12±0.27)
Maximum likelihood assignment tests were performed on four geographic groups of whales
sampled on northbound spring migration (near Point Hope, eastern Chukotka, and Little
Diomede Island) or during summer and early fall outside traditional concentration areas
(near Barrow and Kaktovik) to determine population of origin (Table 1, Fig 2). Of the likely
migration destinations, WHICHRUN assigned the majority of Point Hope whales (95.4%) and
Fig 4. Mantel tests of the correlation between genetic differentiation (Fst) and geographic distance among summering and wintering grounds of
beluga whales in the Bering, Chukchi and Beaufort Seas, for both (A) mitochondrial DNA and (B) microsatellite markers.Test p values are based on
10,000 permutations of the distance data.
13 / 32
whales from the other three areas (85.7%) to the Mackenzie-Amundsen concentration area
(i.e. eastern Beaufort Sea) for mtDNA (Table 4). Individual likelihoods of these whales
based on nuclear genotypes were also higher for this concentration area over the Kasegaluk
Lagoon concentration area (i.e., eastern Chukchi Sea) with many having combined
mtDNA-nDNA LOD scores > 1 for Mackenzie-Amundsen (n = 23/45) (Table 4). Breaking
assignments out by marker type revealed that mtDNA had a greater influence than any
other locus on cumulative likelihoods (Fig 5A).
Assessing population origins of groups of migrating whales by using sample group
information in STRUCTURE, the model-based clustering method (nDNA only) assigned all the Chukotka,
Point Hope and Barrow-Kaktovik whales to the Mackenzie±Amundsen summering
concentration in the eastern Beaufort (Table 4, Fig 5B). By contrast, the majority of the whales
migrating past Diomede had higher inferred membership in, or ancestry (Q) from, Kasegaluk
Lagoon in the eastern Chukchi Sea.
The stock-mixture method, BAYES, estimated population proportions and thus likelihoods
of migrating groups being of mixed origin, and assigned individual migrants probabilistically
to baseline populations. For all runs, estimated shrink factors were close to one (1.00±1.01)
indicating convergence of multiple chains generated from different starting proportions to the
posterior density. Viewing the mtDNA results, groups of migrating whales from all four strata
had very high likelihoods (median = 0.80±0.97) of originating from the eastern Beaufort and
correspondingly low likelihoods of coming from the eastern Chukchi Sea (Fig 6, S4 Table).
Only belugas sampled near Barrow had moderate likelihoods of possible mixed composition.
The microsatellite results also found high population proportions from the eastern Beaufort
for Chukotka, Point Hope and Barrow-Kaktovik, although the density distributions exhibited
greater overlap compared to mtDNA. In contrast to the mtDNA findings, the Diomede whales
had the highest likelihood of coming from the eastern Chukchi Sea. Individual assignments for
each marker type affirmed the group findings that herds migrating past Chukotka and Point
Hope in spring were unlikely to be a mix of Beaufort and Chukchi Sea whales, with almost all
assigned to the Beaufort at P>0.9 for mtDNA and P>0.6 for microsatellites (Fig 6, S4 Table).
Assignments differed dramatically among marker types for Diomede with BAYES finding high
microsatellite likelihoods that most of the Diomede whales (8/10) were from the Chukchi but
high mtDNA likelihoods that the Diomede whales hailed from the Beaufort (Fig 5C). Finally,
while both nDNA and mtDNA assigned most of the Barrow-Kaktovik whales to the Beaufort,
some individuals were assigned to the Chukchi Sea based on mtDNA data (Fig 5C).
While the assignment of separate groupings into the same `population cluster' is not in itself
proof of random mixing, these findings and the absence of mtDNA and nDNA differentiation,
in combination with the timing of the northbound migration past Point Hope and the arrival
of beluga whales at the Mackenzie Delta and Amundsen Gulf, indicate that the whales from
both locations are part of the same population, the eastern Beaufort Sea population. It is for
this reason that analyses of dispersal and temporal patterns (see below) were conducted with
the Point Hope and Mackenzie-Amundsen whales combined in a single eastern Beaufort Sea
Genetic estimates of recent dispersal rates using BAYESASS indicated negligible immigration
(m 0.004) into Cook Inlet from other populations (S5 Table). By contrast, uncertainty over
estimates among summering concentrations within BCB suggests that, unlike the comparisons
between Cook Inlet and BCB, there was not enough information in the data to estimate rates
14 / 32
²LOD: Log10 of the ratio of the second most likely population for mtDNA and microsatellites combined. Positive values indicate assignments to the Beaufort
(Mackenzie), negative values to the eastern Chukchi
Individual whose haplotype was not in the baseline samples.
¶ Individual with incompletely scored genotype or haplotype.
of dispersal within BCB over such a short time frame (a few generations) using this method
(see S1 Supporting information).
In the assessment of sex-biased gene flow through comparison of relative ratios of nDNA to
mtDNA differentiation, estimates of microsatellite heterogeneity (Fst) for those pair-wise
comparisons involving Cook Inlet were, on average, 7.1 times lower than corresponding estimates
for mtDNA (Table 3). By contrast, mean microsatellite differentiation among the four BCB
summering areas was 22.9 times lower than corresponding mtDNA estimators. Similar results
were found with distance-based estimators. Two exceptions to this general pattern are
noteworthy. Firstly, the ratio for Bristol Bay-Norton Sound is low (Fst ratio = 7.3) due to the high
frequency of mtDNA Hap#5 in both areas. Secondly, the ratio for Norton Sound-Mackenzie is
very high (Fst ratio = 55.6) due to the very low levels of microsatellite differentiation observed
between these two areas (Table 3), a finding also reflected in the clustering analysis (Fig 3).
Assessment of contemporary sex-biased dispersal through comparisons of genetic
heterogeneity and relatedness for post-dispersal age cohorts was possible for the eastern Beaufort and
Chukchi Seas for the decade 1988±1997. MtDNA differentiation among adult ( 350 cm
standard length) males, though substantial, was lower than differentiation among adult ( 300cm
standard length) females (Fst males = 0.218 v. Fst females = 0.235, Table 5A). Similarly, using FSTAT
16 / 32
allelic frequency differences at nuclear markers between the eastern Beaufort and eastern
Chukchi were significantly lower for adult males compared to adult females (Fst males = 0.010 v.
Fst females = 0.022, p = 0.043, Table 5A), while estimated average relatedness within each area
was significantly lower for males compared to females (rmales = 0.0204 v. rfemales = 0.0419,
p = 0.042). These sex differences were more pronounced in large (i.e., older) adults (Table 5B).
Patterns of recent dispersal were also characterized using individual-assignment methods.
Using four likelihood and Bayesian inference methods very few individuals were identified as
possible migrants or of migrant ancestry (S6 Table, see S1 Supporting information). Of all the
individual whales identified as possible migrants or of having recent migrant ancestry, the
majority (74%) were adult males (S6 Table). While the number of males sampled in each
population tends to exceed that of females [
], comparison of the sex ratio for the most thoroughly
sampled population, the eastern Chukchi Sea (M:F± 1.56:1) [
], to that for likely migrants or
migrant descendants identified here (M:F± 5.67:1) revealed a significant male bias (χ2 = 4.77,
p = 0.03).
Temporal and kinship analysis
Temporal comparisons on decadal scales revealed almost no change in the geographic patterns
of mtDNA and nDNA variation for almost all geographic strata over ecological time frames.
Mean mtDNA differentiation among time periods (i.e., 1980s, 1990s, 2000s, etc.) within each
area did not differ significantly from zero (F st = 0.017 ± 0.1, p >0.05; Table A in S7 Table), and
was significantly lower than corresponding mean values among areas (F st = 0.308 ± 0.15;
ttest: p = 0.0003). Results were similar for nDNA (Table B in S7 Table), while STRUCTURE
assigned whales from each summering concentration area sampled in different time periods
into the same population clusters (S4 Fig). As with the original combined dataset, Norton
Sound exhibited low differentiation (Fst) from, and clustered with, both Bristol Bay and the
eastern Beaufort Sea. The one exception was Kotzebue Sound, where substantial mtDNA
differentiation (Fst = 0.12±0.31) was observed between beluga whales sampled in the 1980s prior
to the population decline and whales sampled in the 1990s and 2000s (Table A in S7 Table).
No nDNA data were available for the 1980s as the DNA was recovered by ancient DNA
methods from teeth and dried tissue.
The relatedness analyses found substantial numbers of closely related individuals within
groups of whales sampled together at the Kasegaluk Lagoon summer aggregation area (Fig 7).
The analysis also found that substantial numbers of closely related belugas, including likely
parent and offspring pairings, were sampled across years within summering concentration
areas. The COANCESTRY analysis revealed that the mean estimate of pairwise relatedness, r ,
within years did not differ significantly from r across years within the same summering
concentration for both moment and maximum likelihood (ML) estimators. This was the case for
sample years up to 20 years apart (Fig 8). Inferring pairwise relationships from ML estimates
of r in ML-RELATE revealed that the proportions of close relatives (i.e., parent-offspring, full-sib,
half-sib or equivalent) to distant or unrelated individuals within a given year was consistently
on the order of 20%:80% (G = 0.2±1.8, p>0.5), and this within-year ratio was not found to
differ to relationship ratios across years, even for comparisons separated by up to 20 years
(G = 0.2±2.9, p>0.5, Fig 7).
The study provided compelling evidence for widespread natal philopatric behavior of beluga
whales not just to summering concentration areas but to entire migratory circuits within
regions where there are few geographic barriers to dispersal enabling distinct subpopulations
17 / 32
Fig 5. The likely population of origin of beluga whales on spring migration sampled at four locations in the Bering, Chukchi and Beaufort
Seas. Individual assignments are represented as the relative height of stacked bars to either of two baseline populations, the eastern Beaufort Sea
(blue) or eastern Chukchi Sea (red) for mtDNA (dark shading) and nDNA (light shading). A: Maximun Likelihood assignments in Whichrun. B:
Baysian assignments using prior sampling information and admixture models in Structure 2.3.4. C: Mixed-stock assignments in Bayes. Only
individuals with complete mtDNA-nDNA profiles are shown. See Table 4 for more details.
to overlap in space and time. Closely related whales were found to aggregate together at coastal
summering areas each year, and close kin were documented at the same summering sites up to
twenty years apart. We found clearer evidence of sex-biased dispersal than previous studies,
18 / 32
Fig 6. The probability distribution of population proportions of groups of beluga whales sampled on northbound migration in
spring. Stock-mixture analysis was conducted in Bayes with the eastern Beaufort Sea (blue) and the eastern Chukchi Sea (red) as baseline
populations and the migrating groups as the potential `mixtures'.The ordinate axis indicates the number of runs.
and documented stability in migration and dispersal behavior over ecological (i.e., decadal)
time frames with notable exceptions. A much expanded nDNA dataset revealed a pattern of
limited interbreeding among many distinct subpopulations within the BCB region where the
spatial proximity of winter ranges and spring migration routes appears to have a greater
influence on the degree of both individual and genetic exchange than summer ranges. This study
supports recent satellite telemetry evidence of a series of mostly distinct, if overlapping,
wintering areas within the BCB region [
], and contrasts with earlier genetic studies [
suggested a single common wintering area and panmixia across all subpopulations within the
Bering, Chukchi and Beaufort Seas.
Extensive spatial and temporal sampling and thorough individual-based data analysis
provided a more detailed insight into the migration patterns, areas of mixing, dispersal behavior
19 / 32
415 cm (Beaufort) and/or
and gene flow of this long-lived, highly vagile species than previously possible. The
modelbased cluster analysis and assignment tests were greatly improved by larger sample sizes and
locus number, informed run assumptions, and the inclusion of all potential population
clusters. The limited information content in the microsatellite data within regions severely limited
the utility of the program BayesAss to estimate recent rates of dispersal, something recognized
by a number of other recent studies [
]. Below we discuss the investigation's findings and
their implications for beluga whale behavioral ecology and management in more detail.
Phylogeographic patterns of differentiation at nDNA and mtDNA markers among beluga
whale populations from the Gulf of Alaska, BCB, and Sea of Okhotsk regions in the Pacific
likely reflect long standing patterns of restricted gene flow and dispersal. This was especially so
for the Cook Inlet population in the Gulf of Alaska confirming that this small, geographically
isolated population is both reproductively and demographically isolated from all other beluga
whale populations, and likely has been over evolutionary time frames (see S1 Supporting
information). These basin-wide findings expand on earlier regional-level studies [32±35, 39] and
are supported by recent satellite telemetry studies that have mapped winter as well as summer
]. Together, these investigations indicate that the Alaska and Kamchatka
Peninsulas have been effective barriers to dispersal and gene flow over time. Projections for
continued climate warming across the Arctic  in concert with inherent natal homing
tendencies of beluga whales (see below) will likely act to maintain distinct regional populations in
The sequence by which discrete geographic strata emerged as distinct population clusters in
the cluster analysis may also reveal something about their origins. Bristol Bay and Kasegaluk
20 / 32
Lagoon emerged as discrete population clusters early on in the various cluster analyses, that is,
at low values of K (see Fig 3, S2 Table) while most other strata within the BCB region as well as
the small Okhotsk cluster did not differentiate out until later (i.e., higher K). Multiple factors
may contribute to such findings (see S1 Supporting information) but this may reflect an early
divergence of these populations within the BCB region and possible separate glacial refugia for
these two populations.
Philopatry, group structure and migratory culture
The substantial mtDNA differentiation we detected for both sexes among summering
concentration areas within the Bering, Chukchi, Beaufort, and Okhotsk Seas likely indicates limited
dispersal among these areas and a probable long established pattern of female-mediated
philopatry to summer migration destinations. Although this has become an almost universal
interpretation for such patterns of mtDNA variation for migratory species [7, 32, 35, 39, 71±77],
and see , it has rarely been explicitly tested. It relies on the assumptions of a simple
drift-dispersal model to explain predominantly frequency-based (Fst) differentiation. However,
multiple factors can complicate interpretations of mtDNA subdivision in terms of contemporary
dispersal patterns: local populations may not be in equilibrium, sampling may be inadequate,
and movement patterns in migratory species may change across years and generations (see
below). Using large sample sizes collected across many years we identified very low numbers
of females and males as likely migrants or of migrant ancestry and recorded consistent patterns
of mtDNA heterogeneity among most areas over decades of sampling. With age of first
reproduction estimated at 8.3yr for Nearctic belugas [
] this time series likely spanned multiple
generations and supports philopatry to natal population or summering group as a dominant
trait in beluga whales. The discovery that substantial numbers of closely related whales,
including parent-offspring pairings, were sampled across years, and even decades, within summering
areas is direct proof of natal homing.
A recent study summarizing findings from multiple satellite telemetry investigations
revealed that beluga whales from separate summering grounds across the BCB region return
to somewhat geographically discrete wintering areas in the northern Bering Sea [
] (Fig 2).
Our genetic findings viewed alongside these winter movement and habitat use data indicate
that philopatry in beluga whales may extend to entire migratory circuits including migration
routes, staging areas, molting sites, summering locations and wintering areas, which in turn
promotes the emergence and maintenance of demographically distinct subpopulations
regardless of the extent of seasonal overlap (Fig 2) or level of interbreeding (see below).
Natal philopatry has been documented or inferred in a wide range of migratory species and
long-distance travelers, including birds [
], fish , reptiles [
], and mammals [
species where there is an extended period of postnatal care or association with parents it has
been postulated that this homing behavior can entail migratory culture where migratory routes
and destinations are learned from parents [
1, 6, 7
]. Whitehead  recently declared that the
evidence for culture in non-humans is rare and that documenting cultural stability (in
nonhumans) is difficult. Furthermore, while genetic differentiation may reflect the faithful
transmission of migratory culture , it is not direct evidence of it.
Colbeck et al. [
] recently found that groups of closely related beluga whales often migrate
together in the Hudson Bay region of Canada and posited that this could facilitate the learning
of migration routes. Matthews and Ferguson [
] confirmed that beluga whale calves are
dependent on their mothers for two or more years, and thus multiple migratory cycles. Such
findings when viewed in conjunction with our findings of limited dispersal among areas, and
substantial numbers of close kin not only within herds and years at discrete summering areas
21 / 32
Fig 7. The proportion of pairwise genealogical relationships estimated for beluga whales sampled within and between years across two decades near Kasegaluk
Lagoon, Alaska. Maximum likelihood estimates of four relationship categories were estimated from genotypic data using the program ML-RELATE. The stacked bars
represent the proportions of distantly/unrelated individuals to closely related individuals (i.e., parent-offspring, full-sib and half-sib or equivalent) for a subset of the
20-year data set comprising the first three years (1988, 1993, 1994) and the last three years (2005, 2006, 2007).
but also across decades at those locations indicate the potential for lifelong associations of
closely related individuals and is compelling evidence for culturally inherited fidelity to
migration destinations in beluga whales that may extend to entire migratory circuits.
The annual return of beluga whales to traditional locations at specific times of year is likely
an adaptation to the non-uniform distribution of resources across the north Pacific and the
seasonal predictability of those resources, including food, molting sites, and calving and
breeding areas. The configuration of the migratory circuit of each beluga whale population or stock
(Fig 2), and thus their degree of overlap, has likely also been shaped by physical factors that
directly shape migration pathways (e.g., sea ice, geography) and by historical factors such as
the routes of postglacial colonization. Any shift in environmental parameters could
presumably disrupt migration routes and seasonal habitat use patterns (see below), and/or the level of
philopatry to those demographically distinct circuits. A recent study revealed that beluga
whale movements are resilient to inter-annual variation in sea ice conditions [
] which may
be related to cultural inheritance. However, the study did observe that anomalies in beluga
whale migration patterns were correlated with anomalous ice years.
Dispersal and gene flow
The increased statistical power of the expanded data set that we used allowed a more
quantitative assessment of sex-biased dispersal than earlier studies. While both sexes tend to be
philopatric, we found that when dispersal does occur it is predominantly by adult males.
Malebiased dispersal is typical in mammals and has been interpreted to be a strategy to reduce
competition for mates and avoid inbreeding [
8, 9, 84
]. Large population and group sizes in Pacific
beluga whales and opportunities to interbreed with other populations at certain times of year
(see below) may largely reduce the need for male belugas to disperse. There is growing
evidence for seasonal age and sex segregation in beluga whales and differences in habitat use
between males and females [16, 22±25, 81]. In light of the present study's findings it appears
22 / 32
Fig 8. Test of differences in mean relatedness (r ) among beluga whales within a single year compared to r between whales from two different
years for the same summering ground using COANCESTRY. The graph depicts results for Kasegaluk Lagoon 1988 compared to 1988±2007 using the
ML estimator TrioML of Wang (2007). If the observed difference (black line) falls outside the 90% (dotted lines), 95% (dashed lines), and 99% (green
solid lines) confidence intervals from the bootstrap analysis distribution the difference is adjudged to be significant.
that both sexes may remain philopatric to the same migratory network but may not always
travel in association.
While multiple factors confound comparisons of heterogeneity at markers with different
modes of inheritance we were able to attribute the much lower nDNA differentiation
compared to mtDNA differentiation among the subpopulations within the BCB region to greater
male-mediated gene flow (see S1 Supporting information). We also, however, rejected
panmixia on a common wintering ground. This differs from an earlier microsatellite study in this
region by Brown Gladden et al. [
] based on a much smaller dataset (i.e., less locations,
samples, years and loci). It also differs from recent findings on beluga whales in Canada and the
Russian Far East. Turgeon et al. [
] found no convincing evidence of nDNA differentiation
among demographically discrete summering assemblages in the Hudson Bay region
concluding this likely reflected near random mating on a common breeding area. Similarly,
Meschersky et al. [
] found no evidence of nDNA differentiation among two of the three
demographically distinct summering aggregations in the Sea of Okhotsk (also see mtDNA
findings in Table 3A this study) and concluded they share a common gene pool.
We found that genetic distance was more strongly related to distances between wintering
rather than summer grounds suggesting gene flow and dispersal most likely occurred on or
near the wintering grounds. Winter and early spring is the peak period of breeding [
the time of year when these discrete subpopulations are in closest proximity [
] (Fig 2). For
23 / 32
example, while patterns of mtDNA variation in Norton Sound were consistent with a
basinwide pattern of female-mediated philopatry, patterns of nDNA differentiation and the
outcomes of the model-based cluster analysis suggested extensive gene flow between Norton
Sound and whales in the eastern Beaufort Sea and in the Bristol Bay subpopulations (Fig 3).
While the summering grounds of these two subpopulations are far from Norton Sound
(1,000km and 2,000km, respectively) a recent telemetry study revealed that the wintering
range of the Norton Sound subpopulation may overlap to some degree with that of the eastern
Beaufort Sea and Bristol Bay whales [
]. Similarly, low nDNA differentiation and clustering
might suggest the same for the Anadyr Bay subpopulation and beluga whales in both the
eastern Beaufort and the Sea of Okhotsk. However, low sample size from the latter precludes a
well-supported inference at this time.
The assignment tests and mixed stock analysis of migrating whales revealed that while both
the wintering and summering areas of the eastern Chukchi Sea and eastern Beaufort Sea
subpopulations may overlap somewhat [
15, 16, 26, 31
], the timing of spring migration differs such
that the whales hunted at coastal sites in Chukotka, the Bering Strait (i.e., Diomede), and
northwest Alaska (i.e., Point Hope) in the spring and off Alaska's Beaufort Sea coast in
summer were predominantly from the eastern Beaufort Sea population (Fig 5). This agrees with
earlier genetic investigations  and recent telemetry studies [
] where the spring migration
of eastern Beaufort whales occurs earlier and through denser sea ice than eastern Chukchi Sea
belugas. The discovery that a few individual whales at some of these spring locations had
higher likelihoods of eastern Chukchi Sea ancestry or of mixed-ancestry, however, indicates
that the Bering Strait region is also an area of mixing in spring. This is also evident in recent
movement data where Citta et al. [
] observed that for tagged eastern Beaufort Sea whales to
migrate north in spring through the Bering Strait earlier than the eastern Chukchi belugas they
had to pass through the latter's primary wintering area.
The pattern of annual return to one traditional summering area in the BCB region, Kotzebue
Sound, stands out from all others. The number of whales returning to Kotzebue Sound in
summer dropped precipitously in the early-1980s [
]. This was followed by decades of
inconsistent returns [
]. Earlier genetic studies found evidence that the pre-decline whales were likely
a demographically distinct subpopulation and that more recent occurrences of belugas in the
Sound in the 1990s and 2000s most likely involved whales from the Beaufort Sea [
These studies, however, did not include all potential source subpopulations in the BCB region
and so the origins of Kotzebue whales remained provisional. The current study included all
major summer aggregations in the BCB region and confirmed that the pre-decline summer
assemblage in Kotzebue Sound was likely a demographically discrete subpopulation and that
beluga whales entering the Sound post-decline in the 1990s and 2000s (at least for those years
where we had large sample sizes) were different from those we originally sampled and were
most likely comprised of whales primarily from the eastern Beaufort Sea sub population.
Evolution, ecology and management
The macro-geographic patterns of genetic differentiation in beluga whales across the North
Pacific suggest divergence of the evolutionary trajectories of regional populations on
millennial timescales. Regional adaptation of beluga populations is likely over such time frames while
rescue of small and/or declining populations like that in Cook Inlet by whales from other
24 / 32
regions is probably remote. The stability of heterogeneity among spatial strata on decadal time
scales within regions and the evidence for generational-scale kinship within strata may
indicate the kind of cultural stability in beluga whale populations that is necessary for gene-culture
coevolution  where distinct cultures drive the evolution of neutral and functional genes
over much faster time scales.
The study revealed that all spatially discrete summering aggregations within regions are
demographically independent subpopulations and as such should be defined as separate stocks
under management objectives to maintain healthy, functioning populations across the species
range . The greatly increased dataset over previous studies now means that five of the six
summering aggregations within the BCB region can be more definitively defined as the
following stocks: Bristol Bay, Eastern Bering Sea (Norton Sound), eastern Chukchi Sea (Kasegaluk
Lagoon), eastern Beaufort Sea (Mackenzie-Amundsen), and Gulf of Anadyr (Anadyr Bay)
(Table 1). We found that the sixth traditional summering area, Kotzebue Sound, was most
likely a distinct subpopulation before declines in the 1980s but requires further investigation
regarding its current stock composition.
Until now, limited knowledge on levels of interbreeding combined with a paucity of
information on movements and distribution in late fall, winter and early spring shaped perceptions
of beluga whale populations, and thus management units, across much of the Arctic that are
dominated by summer distributions [
20, 37, 38, 87
]. We define stocks, for example, by their
primary coastal summering location.
This study's findings, in concert with new satellite telemetry data [
], force us to alter
this perception. The demographically distinct summering aggregations within the BCB region,
for example, in most cases: (1) do not appear to interbreed extensively, (2) return to discrete
wintering areas, and (3) disperse and interbreed over limited distances. Thus, beluga whales in
the Bering, Chukchi and Beaufort Seas comprise a series of populations or subpopulations
with discrete, traditional migratory circuits and destinations that overlap to varying degrees at
certain times of year. As well as challenging current concepts on connectivity between
subpopulations within regions, these findings will also alter our approach to estimating effective
population sizes (Ne) and how to model risk, recovery and population viability. Furthermore, when
the migratory circuits of these subpopulations are considered in light of the heterogeneity of
Pacific Arctic and subarctic marine ecosystems  the individual subpopulations likely also
have quite distinct ecologies and perhaps other unique aspects to their behavior and
population biology. It follows then that different migratory populations are likely exposed to different
Finally, the study confirms the critical importance of the Northern Bering Sea and Bering
Strait region to beluga whales. This relatively small area (Fig 2) is home to tens of thousands of
beluga whales from multiple subpopulations, some among the largest in the world, for example
], for much of the year. This region is also important for other marine mammal species
(Citta et al. in review.) and is currently undergoing dramatic environmental and ecosystem
]. Minor shifts in this region's environment and ecosystem, including sea-ice
and productivity, could have major impacts on beluga whale dispersal, breeding behavior,
population status and structure across much of their Pacific range.
S1 Supporting information. This is a file that contains supporting information text.
25 / 32
S1 Fig. Median Joining network of beluga whale mtDNA haplotypes across their North
S2 Fig. Summary plot of model-based cluster analyses of population structure in Pacific
beluga whales with Nmax = 100 (A) compared to whales with N = all (B).
S3 Fig. Estimation of the likely number of population clusters, K, using the rate of change
in the log probability of the data, the ΔK statistic of Evanno.
S4 Fig. Temporal cluster analysis of nDNA data from beluga whales in the Pacific using
STRUCTURE. For those summering concentrations where we had genotypes from more than
one decade individuals from each summering location were assigned to the same population
cluster. Each of 973 individuals is represented by a vertical line with estimated membership, Q,
in each cluster denoted by different colors. The analysis was based on eight microsatellite loci,
used prior sample group information (LOCPRIOR), and yielded similar results for both the
admixture and no admixture (shown) models of ancestry. For each geographic stratum the
decade 1988±1997 is denoted by a grey bar across the top of the figure, and the decade 1998±
2007 is denoted by a white bar.
S1 Table. Estimates of genetic diversity and probabilities of identity for eight
microsatellite loci in beluga whales. Values are given for the major beluga whale concentration areas in
the western Nearctic, and for the combined spring migration (i.e., Point Hope) and
summering concentration (I.e. Mackenzie-Amundsen) areas for the Beaufort Sea.
S2 Table. Summary of results from STRUCTURE analysis of microsatellite data from
beluga whales in the North Pacific. The likelihoods of K given the data is reported for 28
separate analyses, each with unique parameter settings. Analyses were run both with and without
admixture, with and without sampling location included as prior information, and they
contained sample sets of varying size randomly sampled using R (n = 15, 30, 50, 100, and all) that
included individuals genotyped at a range of loci (n 6 loci, n 7 loci, n = all 8 loci). For
each parameter set we used a burn-in period of 50,000 iterations followed by 2x105 iterations
to collect data. Multiple runs (n = 5) were conducted for each value of K (n = 8) to ensure
convergence for a total of 1,120 separate runs. The most likely K for each analysis is shaded and
details on the clustering is provided in the comments. Analyses where K = 6 was well
supported are shaded in blue.
S3 Table. Correlation between geographic distance and genetic distance among seasonal
ranges of beluga whales in the Bering, Chukchi and Beaufort Seas for both mtDNA and
microsatellites. Mantel tests, based on 10,000 permutations were conducted in Arlequin 3.5.
S4 Table. Point estimates (mean and median) and Bayesian posterior probability bounds
(95%) of population composition estimates, p, of beluga whales migrating past four
locations in the Bering, Chukchi and Beaufort Seas in spring and early summer using the
stock-mixture program BAYES. Baseline samples are from likely migration destinations, the
eastern Beaufort Sea and the eastern Chukchi Sea. Multiple chains (1,800±12,000 reps) were
26 / 32
run and shrink factors did not exceed 1.01.
S5 Table. Genetic estimates of recent immigration into beluga whale populations. Means
(± 95% confidence intervals) of the posterior distribution of the proportion of individuals that
are migrants (m) based on the analysis of multi-locus genotypic data in BayesAss. To insure
convergence 10 separate runs, each with 20 x 106 iterations, the first 1 x 106 iterations of
which were burn-in, were conducted. Receiving populations are in rows, delta = 0.15 and n
denotes sample size. Similar results were found when the Beaufort Sea was represented by
whales sampled only at the Mackenzie Delta-Amundsen Gulf.
S6 Table. Identification of candidate migrants and individuals with likely migrant ancestry
using Bayesian and maximum-likelihood assignment tests. Bayesian inference was based on
multilocus genotype probabilities using the model-based clustering method, STRUCTURE
(MCMCs involved a burn-in period of 30,000 rep.s followed by 1 X 106 rep.s.) and the
BayesAss program (MCMC runs with 20 X 106 rep.s, the first 1 X 106 discarded as burn-in).
Likelihoods were estimated in WHICHRUN and are reported as log10 ratios of the likelihood that
an individual's haplo-genotype was more likely in a population other than where it was
sampled assuming random assortment of microsatellite alleles and that observed mtDNA
haplotypic frequencies and microsatellite allele frequencies represent population frequencies. Relative
likelihoods were estimated with the Assignment Calculator program in DOH where P
represents the proportion of 1,000 new individuals randomly drawn from each population which
had equal or smaller likelihood values. Only results for the sampled population (i.e., baseline)
and the population the STRUCTURE and WHICHRUN analysis assigned the individual to
are presented. Only individuals that had low probabilities of non-immigrant ancestry (at either
ν = 0.05 or 0.1) and/or had high mtDNA-nDNA LOD scores ( 1) are reported.
S7 Table. Temporal patterns of genetic differentiation (Fst) within mitochondrial DNA
(A) and within microsatellite loci (B) in Pacific beluga whales. The temporal strata
comprised the following four 10-year sample windows: 1980s = 1978±1987); 1990s = 1988±
1997); 2000s = 1998±2007); and 2010s ( 2008). P-values from homogeneity tests, based
on 50,000 permutations, are represented by the following shading patterns: dark grey:
p 0.01, light grey: 0.01<p 0.05, unshaded: p>0.05. Only strata with a sample size of
n 10 are reported. The boxes bordered in blue highlight pairwise comparisons across
decades within geographic strata. 50,175 permutations.
S1 Appendix. A more detailed summary of sample sizes of beluga whales from the same
fifteen geographic strata in the North Pacific Ocean shown in Table 1. Here, the samples are
broken out by decadal time periods. The last two columns summarize the total sample sizes for
the primary sampling period from 1988±2010 and includes samples from Meschersky et al.
This study was a collaborative effort with the subsistence hunters of the Alaskan, western
Canadian Arctic, and Chukotkan coasts. We greatly appreciate their generous provision of
samples and insight, and for their continued support on all aspects of this study. Invaluable lab
27 / 32
assistance was provided by C. LeDuc, M. DeAngelis, C. Bonin, A. Nestler, and S. Rodgers.
Thanks also to M. Adams, A. Burdin, K. Burek, J. Dau, V. Eunichevun, C. Goertz, L. Lowry, R.
Hobbs, V. Melnikov, V. Nikulin, J. Orr, L. Pierce, T. Romano, G. Seaman, R. Small, H. Smith,
and A. Whiting for their advice, support and assistance with sampling. We'd especially like to
thank Willie Goodwin, Chairman of the Alaska Beluga Whale Committee, for his continued
support of beluga genetics research, his knowledge, and his assistance in obtaining samples.
Thanks also to Lloyd Lowry for helpful reviews of earlier drafts of this manuscript. All samples
were collected under the authority of U.S. Marine Mammal Protection Act permits issued by
the National Marine Fisheries Service, Russian Marine Mammal permits, and the Dept. of
Fisheries and Oceans, Canada scientific licenses.
Conceptualization: Greg O'Corry-Crowe, Robert Suydam.
Data curation: Lori Quakenbush, Brooke Potgieter, Lois Harwood, Dennis Litovka, John
Citta, Vladimir Burkanov, Kathy Frost, Barbara Mahoney.
Formal analysis: Greg O'Corry-Crowe, Tatiana Ferrer.
Project administration: Greg O'Corry-Crowe.
Resources: Robert Suydam, Lori Quakenbush.
Writing ± original draft: Greg O'Corry-Crowe.
Writing ± review & editing: Robert Suydam, Lori Quakenbush, Brooke Potgieter, Lois
Harwood, Dennis Litovka, Tatiana Ferrer, John Citta, Vladimir Burkanov, Kathy Frost.
28 / 32
29 / 32
30 / 32
62. Brzustowski J. Doh assignment test calculator. http://www2.biology.ualberta.Ca/jbrzusto/Doh.php.
Wilson GA, Rannala B. Bayesian inference of recent migration rates using multilocus genotypes.
Genetics. 2003; 163: 1177±1191. PMID: 12663554
31 / 32
Whitehead H , Rendell L . The cultural lives of whales and dolphins . 2015 ; The University of Chicago Press, Chicago.
2. Robertson GJ , Cooke F . Winter philopatry in migratory waterfowl . The Auk . 1999 ; 116 : 20 ± 34 .
3. Jorgensen SJ , Reeb CA , Chapple TK , Anderson S , Perle C , Van Sommeran SR , et al. Philopatry and migration of Pacific white sharks . Proc R Soc B . 2010 ; 277 : 679 ± 688 . https://doi.org/10.1098/rspb. 2009 .1155 PMID: 19889703
Whitehead H . Gene-culture coevolution in whales and dolphins . Proc Nat Acad Sci . 2017 ; 114 : 7814 ± 7821 .
5. Valenzuela L , Sironi M , Rowntree V , Seger J . Isotopic and genetic evidence for culturally inherited site fidelity to feeding grounds in southern right whales (Eubalaena australis) . Mol Ecol . 2009 ; 18 : 782 ± 791 . https://doi.org/10.1111/j. 1365 - 294X . 2008 . 04069 . x PMID : 19207250
6. Harrison XA , Tregenza T , Inger R , Colhoun K , Dawson DA , Gudmundsson GA et al. Cultural inheritance drives fidelity and migratory connectivity in a long-distance migrant . Mol Ecol . 2010 ; 19 : 5484 ± 5496 . https://doi.org/10.1111/j. 1365 - 294X . 2010 . 04852 . x PMID : 21083633
7. Carroll EL , Baker CS , Watson M , Alderman R , Bannister J , Gaggiotti OE et al. Cultural traditions across a migratory network shape the genetic structure of southern right whales around Australia and New Zealand . Sci Rep . 2015 ; 5 : 16182| https://doi.org/10.1038/srep16182 PMID: 26548756
8. Greenwood PJ . Mating Systems, philopatry and dispersal in birds and mammals . Anim. Behav . 1980 ; 28 : 1140 ± 1162 .
9. Dobson FS . The enduring question of sex-biased dispersal: Paul J. Greenwood's (1980) seminal contribution . Anim Behav . 2013 ; 85 : 299 ± 304 .
10. Prugnolle F , de Meeus Y. Inferring sex-biased dispersal from population genetic tools: a review . Heredity . 2002 ; 88 : 161 ± 165 . https://doi.org/10.1038/sj.hdy. 6800060 PMID: 11920116
11. Zhu X , Wu K , Wu L . Influence of the physical environment on the migration and distribution of Nibea albiflora in the Yellow Sea . J. Ocean Univ. China . 2017 ; 16 : 87 ± 92 .
Weider LJ , Hobk A . Phylogeography and arctic biodiversity: a review . Ann. Zool. Fennici . 2000 ; 37 : 217 ± 231 .
13. O 'Corry-Crowe G. Beluga Whale (Delphinapterus leucas) . In Encyclopedia of Marine Mammals. Third Edition . Ed.Wursig s B. , Thewissen J.G. M. , and Kovacs K. Elsevier . 2017 .
14. Heidi-Jørgensen MP , Sinding M-HS , Nielsen NH , Rosing-Avid A , Hansen RG . Large numbers of marine mammals winter in the North Water polynya . Polar Biol . 2016 ; 39 : 1605 ± 1614 .
15. Suydam RS , Lowry LF , Frost KJ , O'Corry-Crowe GM , Pikok D Jr . Satellite tracking of eastern Chukchi Sea beluga whales into the Arctic Ocean . Arctic. 2001 ; 54 : 237 ± 243 .
16. Richard PR , Martin AR , Orr JR . Summer and autumn movements of belugas of the eastern Beaufort Sea stock . Arctic . 2001 ; 54 : 223 ± 236 .
17. Hobbs RC , Laidre KL , Vos DJ , Mahoney BA , Eagleton M. Movements and area use of belugas, Delphinapterus leucas, in a subarctic Alaskan estuary . Arctic . 2005 ; 58 : 331 ± 340 .
18. Citta JJ , Quakenbush LT , Frost KJ , Lowry LF , Hobbs RC , Aderman H . Movements of beluga whales (Delphinapterus leucas) in Bristol Bay, Alaska . Mar Mamm Sci . 2016 ; 32 : 1271 ± 1298 .
19. Lydersen C , Martin AR , Kovacs KM , Gjertz I . Summer and autumn movements of white whales Delphinapterus leucas in Svalbard, Norway . Mar Ecol Progr Ser . 2001 ; 219 : 265 ± 274 .
20. Frost KJ , Lowry LF . Distribution, abundance and movements of beluga whales, Delphinapterus leucas, in coastal waters of western Alaska . Can J Fish Aquat Sci . 1990 ; 224 : 39 ± 57 .
21. Solovyev BA , Shpak OV , Glazov DM , Rozhnov VV , Kuznetsova DM . Summer distribution of beluga whales (Delphinapterus leucas) in the Sea of Okhotsk . Russ J Theriol . 2015 ; 14 : 201 ± 215
22. Suydam RS . Age, growth, reproduction and movements of beluga whales (Delphinapterus leucas) from the eastern Chukchi Sea . Doctoral dissertation . University of Washington, Seattle. 2009 ; 169pp.
23. Loseto LL , Richard P , Stern GA . Segregation of Beaufort Sea beluga whales during the open-water season . Can J Zool . 2006 ; 84 : 1743 ± 1751 .
24. Vv Krasnova , Chernetsky AD Kirillova OI , Bel'kovich VM. The dynamics of the abundance, age, and sex structure of the Solovetsky reproductive gathering of the beluga whale Delphinapterus leucas (Onega Bay, White Sea) . Russ J Marine Biol . 2012 ; 38 : 218 ± 225 .
25. Hauser DDW , Laidre KL , Stern HL , Moore SE , Suydam RS , Richard PR . Habitat selection by two beluga whale populations in the Chukchi and Beaufort seas . PLoS One . 2017 ; 12 : e0172755. https:// doi.org/10.1371/journal.pone. 0172755 PMID: 28235041
26. Hauser DDW , Laidre KL , Suydam RS , Richard PR . Population-specific home ranges and migration timing of Pacific Arctic beluga whales (Delphinapterus leucas) . Polar Biol . 2014 ; 37 : 1171 ± 1183 .
27. Brown Gladden JG , Ferguson MM , Freisen MK , Clayton JW . Population structure of North American beluga whales (Delphinapterus leucas) based on nuclear DNA microsatellite variation and contrasted with population structure revealed by mitochondrial DNA . Mol Ecol . 1999 ; 8 : 347 ± 363 . PMID: 10199005
28. Luque SP , Ferguson SH . Age structure, growth, mortality, and density of belugas (Delphinapterus leucas) in the Canadian Arctic: responses to environment? Polar Biol . 2010 ; 33 : 163 ± 178 .
29. Litovka DI , Hobbs RC , Laidre KL , O'Corry-Crowe GM , Orr JR , Richard PR , et al. Research of belugas whales Delphinapterus leucas in Anadyr Gulf (Chukotka) using satellite telemetry. Marine Mammals of the Holarctic . Abstracts of the conference presentations . Moscow, 294pp. Belkovich VM ed. 2002 .
30. Litovka DI , Andronov PYu , Batanov RL . Analysis of seasonal distribution of beluga whales Delphinapterus leucas and its prey objects in coastal waters of north-western Bering Sea . Bulletin of KanchatNIRO: Invest Biol Resour Kamchat NW Pac . 2013 ; 27 : 50 ± 71 . (In Russian).
31. Citta JJ , Richard P , Lowry LF , O'Corry-Crowe G , Marcoux M , Suydam R . et al. Satellite telemetry reveals population specific winter ranges of beluga whales in the Bering Sea . Mar Mamm Sci . 2016 ; 33 : 236 ± 250 .
32. O' Corry-Crowe GM , Suydam RS , Rosenberg A , Frost KJ , Dizon AE . Phylogeography, population structure and dispersal patterns of the beluga whale Delphinapterus leucas in the western Neartic revealed by mitochondrial DNA . Mol Ecol . 1997 ; 6 : 955 ± 970 .
33. O 'Corry-Crowe GM , Dizon AE , Suydam RS , Lowry LF . Molecular genetic studies of population structure and movement patterns in a migratory species: the beluga whale (Delphinapterus leucas) in the western Nearctic . In Molecular and Cell Biology of Marine Mammals. Ed. Pfeiffer C.J. Krieger Publishing Co. Malabar , Florida. 2002 : 53 ± 64 .
34. O 'Corry-Crowe GM , Lydersen C , Heide-Jørgensen MP , Hansen L , Mukhametov LM , Dove O , et al. Population genetic structure and evolutionary history of North Atlantic beluga whales (Delphinapterus leucas) from West Greenland, Svalbard and the White Sea . Polar Biol . 2010 ; 33 : 1179 ± 1194 .
35. Brown Gladden JG , Ferguson MM , Clayton JW . Matriarchal genetic population structure of North American beluga whales (Delphinapterus leucas , Cetacea: Monodontidae). Mol Ecol . 1997 ; 6 : 1033 ± 1046 . PMID: 9394462
36. Meschersky IG , Kholodova MV , Zvychaynaya EY . Molecular genetic ecology of the beluga (Delphinapterus leucas: Cetacea, Monodontidae) summering in the southern Sea of Okhotsk as compared to North American Populations . Russ J Genet . 2008 ; 44 : 1105 ± 1110 .
37. Muto MM , Helker VT , Angliss RP , Allen BA , Boveng PL , Breiwick JM , et al. Alaska marine mammal stock assessments , 2016 . U.S. Dep. Commer., NOAA Tech. Memo. NMFS-AFSC-355 2017 ; 355 : 366 p. http://dx.doi.org/10.7289/V5/TM-AFSC- 355 .
38. COSEWIC. COSEWIC assessment and update status report on the beluga whale Delphinapterus leucas in Canada. Committee on the Status of Endangered Wildlife in Canada . Ottawa. ix + 70 pp. (www. sararegistry .gc.ca/status/status_e. cfm) . 2004 .
39. Meschersky IG , Shpak OV , Litovka DI , Glazov DM , Borisova EA , Rozhnov VV . A genetic analysis of the beluga whale Delphinapterus leucas (Cetacea: Monodontidae) from summer aggregations in the Russian Far East . Russ J. Mar Biol . 2013 ; 39 : 125 ± 135 .
40. HoÈss M , PaÈaÈbo S. DNA extraction from Pleistocene bones by a silica-based purification method . Nucleic Acids Res . 1993 ; 21 : 3913 ± 3914 . PMID: 8396242
41. Rohland N , Hofreiter M. Ancient DNA extraction from bones and teeth . Nature Prot . 2007 ; 2 : 1756 ± 1762 .
42. Fain SR , LeMay JP . Gender identification of humans and mammalian wildlife species from PCR amplified sex linked genes . Proc Am Acad Forensic Sci . 1995 ; 1 : 34 .
43. Buchanan FC , Freisen MK , Littlejohn RP , Clayton JW . Microsatellites from the beluga whale Delphinapterus leucas . Mol Ecol . 1996 ; 5 : 571 ± 575 . PMID: 8794563
44. Schltterer C , Amos B , Tautz D. Conservation of polymorphic simple sequence loci in cetacean species . Nature . 1991 ; 354 : 63 ± 65 . https://doi.org/10.1038/354063a0 PMID: 1944571
45. Valsecchi E , Amos W. Microsatellite markers for the study of cetacean populations . Mol Ecol . 1996 ; 5 : 151 ± 156 . PMID: 9147690
46. Van Oosterhout C , Hutchinson WF , Wills DPM , Shipley P . MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data . Mol. Ecol. Notes . 2004 ; 4 : 535 ± 538 .
47. Tamura K , Stecher G , Peterson D , Filipski A , Kumar S. MEGA6: molecular genetics analysis version 6.0 . Mol Biol Evol . 2013 ; 30 : 2725 ± 2729 . https://doi.org/10.1093/molbev/mst197 PMID: 24132122
48. Nei M , Tajima F. DNA polymorphism detectable by restriction endonucleases . Genetics . 1981 ; 97 : 145 ± 163 . PMID: 6266912
49. Nei M. Molecular evolutionary genetics . Columbia University Press, New York 1987 .
50. Excoffier L Lischer HEL . Arlequin suite ver 3.5. A new series of programs to perform population genetics analyses under Linux and Windows . Mol Ecol Resour . 2010 ; 10 : 564 ± 567 . https://doi.org/10.1111/j. 1755- 0998 . 2010 . 02847 . x PMID : 21565059
51. Kalinowski ST , Taper ML , Marshall TC . Revising how the computer program CERVUS accommodates genotyping error increases success in paternity assignment . Mol Ecol . 2007 ; 16 : 1099 ± 1006 . https:// doi.org/10.1111/j. 1365 - 294X . 2007 . 03089 . x PMID : 17305863
52. Raymond M , Rousset F. GENEPOP Version 1.2: population genetics software for exact tests and ecumenicism . J Hered . 1995 ; 86 : 248 ± 249 .
53. Guo SW , Thompson EA . Performing the exact test of Hardy-Weinberg proportion for multiple alleles . Biometrics . 1992 ; 48 : 361 ± 372 . PMID: 1637966
54. Weir BS , Cockerham CC . Estimating F-statistics for the analysis of population structure . Evolution . 1984 ; 38 : 1358 ± 1370 . https://doi.org/10.1111/j.1558- 5646 . 1984 .tb05657. x PMID: 28563791
55. Excoffier L , Smouse P , Quattro J . Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data . Genetics . 1992 ; 131 : 479 ± 491 . PMID: 1644282
56. Pritchard JK , Stephen M , Donnelly P . Inference of population structure using multilocus genotype data . Genetics . 2000 ; 155 : 945 ± 959 . PMID: 10835412
57. Pritchard JK , Wen X , Falush D . Documentation for structure software: version 2 .3. University of Chicago. 2010 .
58. Hubisz MJ , Falush D , Stephens M , Pritchard J . Inferring weak population structure with the assistance of sample group information . Mol Ecol. Res . 2009 ; 9 : 1322 ± 1332 .
59. Kopelman NM , Mayzel J , Jakobsson M , Rosenberg NA , Mayrose I. CLUMPAK : a program for indetifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour . 2015 ; 15 : 1179 ± 1191 . https://doi.org/10.1111/ 1755 - 0998 .12387 PMID: 25684545
60. Paetkau D , Calvert W , Sterling I , Strobeck C . Microsatellite analysis of population structure in Canadian polar bears . Mol Ecol . 1995 : 4 : 347 ± 354 . PMID: 7663752
61. Banks MA , Eichert W. WHICHRUN (version 3.2): a computer program for population assignment of individuals based on multilocus genotype data . J. Hered . 2000 ; 91 : 87 ± 89 . PMID: 10739137
64. Goudet J , Perrin N , Waser P . Tests for sex-biased dispersal using bi-parentally inherited genetic markers . Mol Ecol . 2002 ; 11 : 1103 ± 1114 . PMID: 12030985
65. Pella J , Masuda M. Bayesian methods for analysis of stock mixtures from genetic characters . Fish. Bull . 2001 ; 99 : 151 ± 167 .
66. Masuda M. User's Manual for BAYES: Bayesian Stock-Mixture Analysis Program , 28 pp. National Marine Fisheries Service , Juneau, Alaska. 2002 .
Wang J. COANCESTRY : a program for simulating, estimating and analyzing relatedness and inbreeding coefficients . Mol Ecol Resour . 2011 ; 11 : 141 ± 145 . https://doi.org/10.1111/j.1755- 0998 . 2010 . 02885 .
x PMID : 21429111
68. Kalinowski ST , Wagner AP , Taper ML . ML-RELATE: a computer program for maximum likelihood estimation of relatedness and relationship . Mol Ecol Notes . 2006 ; 6 : 576 ± 579 .
69. Malenfant RM , Davis CS , Cullingham CI , Coltman DW . Circumpolar genetic structure and recent gene flow of polar bears: a reanalysis . PLOS ONE . 2016 ; 11 ( 3 ): e0148967. https://doi.org/10.1371/journal. pone. 0148967 PMID: 26974333
70. Overland JE , Wang M , Walsh JE , Stroeve JC. 2013 . Future Arctic climate changes: Adaptation and mitigation time scales . Earth's Future . 2013 : 2: https://doi.org/10.1002/2013EF000162
71. Lyrholm T , Gyllensten U . Global matrilineal population structure in sperm whales as indicated by mitochondrial DNA sequences . Proc. R. Soc. Lond. B , 1988 ; 265 : 1679 ± 1684 .
72. de March BGE , Maiers LD , Friesen MK . An overview of genetic relationships of Canadian and adjacent populations of belugas (Delphinapterus leucas) with emphasis on Baffin Bay and Canadian eastern Arctic populations . NAMCO Sci. Publ . 2002 ; 4 : 17 ± 38 .
73. Pearce JM , Talbot SL , Petersen MR , Rearick JR . Limited genetic differentiation among breeding, molting, and wintering groups of the threatened Steller's eider: the role of historic and contemporary factors . Conserv Gen . 2005 ; 6 : 743 ± 757 .
74. O 'Corry-Crowe G , Taylor BL , Gelatt T , Loughlin TR , Bickham J , Basterretche M , et al. Demographic independence along ecosystem boundaries in Steller sea lions revealed by mtDNA analysis: implications for management of an endangered species . Can. J. Zool . 2006 ; 84 : 1796 ± 1809 .
75. Turgeon J. , Duchesne P , Colbeck GJ , Postma LD , Hammill MO . Spatiotemporal segregation among summer stocks of beluga (Delphinapterus leucas) despite nuclear gene flow: implication for the endangered belugas in eastern Hudson Bay . Conserv Gen . 2012 ; 13 : 419 ± 433 .
76. Dutton PH , Roden SE , Stewart KR , LaCasella E , Tiwari M , Formia A et al. Population stock structure of leatherback turtles (Dermochelys coriacea) in the Atlantic revealed using mtDNA and microsatellite markers . Conserv Gen . 2013 ; 14 : 625 ± 636 .
77. Baker CS , Steel D , Calambokidas J , Falcone E , GoncaÂlez-Peral U , Barlow J et al. Strong maternal fidelity and natal philopatry shape genetic structure in North Pacific humpback whales . Mar Ecol Prog Ser . 2013 ; 494 : 291 ± 306 .
78. Forschler MI , del Val E , Bairlein F. Extraordinary high natal philopatry in a migratory passerine . J Ornithol . 2010 ; 151 : 745 ± 748 .
79. Bonanomi S , Therkildsen NO , Retzel A , Hedeholm RB , Pedersen MW , Meldrup D et al. Historical DNA documents long-distance natal homing in marine fish . Mol Ecol . 2016 ; 25 : 2727 ± 2734 . https://doi.org/ 10.1111/mec.13580 PMID: 26859133
80. Hoffman JI . Forcada J. Extreme natal philopatry in female Antarctic fur seals (Arctocephalus gazella) . Z Saugetierkd . 2013 ; 77 : 71 ± 73 .
81. Colbeck GJ , Duchesne P , Postma LD , Lesage V , Hammill MO , Turgeon J . Groups of related belugas (Delphinapterus leucas) travel together during their seasonal migrations in and around Hudson Bay . Proc R Soc B . 2013 ; 280 : 20122552 . https://doi.org/10.1098/rspb. 2012 .2552 PMID: 23222451
82. Matthews CJD , Ferguson SH . Weaning age variation in beluga whales (Delphinapterus leucas) . J. Mamm . 2015 ; 96 : 425 ± 437 .
83. O 'Corry-Crowe G , Mahoney A , Suydam R , Quakenbush L , Whiting A , Lowry L , et al. Genetic profiling links changing sea ice to shifting beluga whale migration patterns . Biol. Lett . 2016 ; 12 : 20160404 . https://doi.org/10.1098/rsbl. 2016 .0404
84. Pusey AE . Sex-biased dispersal and inbreeding avoidance in birds and mammals . Trends Ecol Evol . 1987 ; 2 : 295 ± 299 . https://doi.org/10.1016/ 0169 - 5347 ( 87 ) 90081 - 4 PMID: 21227869
85. Seaman G , Barger E , Lee R Jr. Buckland beluga whale traditional ecological knowledge project . Final report to the Native Village of Buckland . 2015 ; 109pp.
Wade PR , Angliss RP . Guidelines for assessing marine mammal stocks: report of the GAMMS workshop , April 3 ± 5 , 1996 , Seattle, Washington. U.S. Dept. Commer. NOAA Tech. Memo. NMFS-OPR-12.
1997 ; 93p .
87. Kleinenberg SE , Yablokov AV , Bel'kovich BM , Tarasevich MN . Beluga (Delphinpaterus leucas): Investigations of the Species . Academy of Sciences of the USSR , Moscow. 1964 . Translated by the Israel Program for Scientific Translations , 1969 . 375pp.
88. Grebmeier JM , Bluhm BA , Cooper LW , Danielson SL , Arrigo KR , Blanchard AL et al. Ecosystem characteristics and processes facilitation persistent microbenthic biomass hotspots and associated benthivory in the Pacific Arctic . Prog Oceanogr . 2015 ; 136 : 92 ± 114 .
89. Harwood LA , Innes S , Norton P , Kingsley MCS . Distribution and abundance of beluga whales in the Mackenzie estuary, southeast Beaufort Sea, and west Amundsen Gulf during late July 1992 . Can J Fish Aquatic Sci . 1996 ; 53 : 2262 ± 2273 .
90. Lowry LF , Kingsley MCS , Hauser DDW , Clarke J , Suydam R . Aerial survey estimates of abundance of the Eastern Chukchi Sea stock of beluga whales (Delphinapterus leucas ) in 2012 . Arctic. 2017 ; 70 : 273 ± 286 .
91. Grebmeier JM , Overland JE , Moore SE , Farley EV , Carmack EC , Frey KE et al. A major ecosystem shift in the northern Bering Sea . Science . 2006 ; 311 : 1461 ± 1464 . https://doi.org/10.1126/science. 1121365 PMID: 16527980