Migratory culture, population structure and stock identity in North Pacific beluga whales (Delphinapterus leucas)

PLOS ONE, Mar 2018

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 Bering-Chukchi-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 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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

http://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0194201&type=printable

Migratory culture, population structure and stock identity in North Pacific beluga whales (Delphinapterus leucas)

March 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. wixsite.com/pbbe/. Funding: A research grant was provided by the Alaska Beluga Whale Committee and administered by the North Slope Borough Department of Wildlife Management (http://www.northslope.org/ 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 Fisheries Service. 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. Introduction 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 [4] 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, 3 ]. (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 [ 8, 9 ] with sex-biased dispersal potentially leading to the demographic isolation of population units even in the face of gene flow and seasonal overlap [ 10 ]. (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 [ 11 ], 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 [12].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 [ 13, 14 ]. 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 [ 20, 21 ]. Breeding is believed to occur primarily in winter and early spring [ 22 ], 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 [ 26 ] and for separate summering groups sharing the same wintering ground during the period of peak mating [ 27 ], 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 overlap [ 31 ]. Considering the peak breeding season for beluga whales is during winter and early spring [ 22 ], 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 [ 37, 38 ]. 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 [ 27, 39 ]. 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 areas. 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 literature [ 39 ] 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. [ 31 ]. 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 [ 32, 33 ]. Total DNA was also extracted from a number of dried tissue and teeth samples using standard silica-based 5 / 32 1 2 3 4 5 6 7 8 9 10 11 A, B C D E, F Cook Inlet Cook Inlet Bristol Bay E. Bering Sea E. Chukchi Sea E. Beaufort Sea ancient DNA methods [ 40, 41 ]. The gender of each sample was determined by PCR-based methods [42] and all samples were analyzed for variation within eight microsatellite loci (S1 Table) typed on beluga whales [ 43 ] or other cetacean species [ 44, 45 ]. 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 [ 27, 34, 43 ]. Comparisons of observed to expected genotypic frequencies using the Micro-Checker program (version 2.2.3) [46] 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). Data analysis 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 6.0 [ 47 ] software and by estimating both haplotypic (H) [ 48 ] and nucleotide (π) [ 49 ] diversity using ARLEQUIN 3.5 software [ 50 ]. Estimation of genetic diversity (He, Ho and number of alleles) and probabilities of identity (I) for nuclear loci was performed using CERVUS 3.0 [ 51 ]. Tests for Hardy-Weinberg (H-W) and linkage equilibria were performed in GENEPOP 4.1 [ 52 ] with P-values estimated from 500,000 iterations of the data using the Mark chain method of Guo and Thompson [ 53 ]. An analysis of both mtDNA and nDNA data revealed that the populations were not in mutation-drift equilibrium (see [ 34 ]). 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 [ 54, 55 ], 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 [ 56, 57 ], 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. [58] 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 [ 59 ] 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 [ 61 ]. 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 [62] 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 [63] 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 ([ 34 ], 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 [ 64 ]. 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. [ 64 ] 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 [ 65 ] 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 [ 66 ], 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 [67] and ML-RELATE [ 68 ] (See S1 Supporting information). Results 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 Supporting information). 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 n = n = Gulf of Alaska 141 9 / 32 Migratorycultureandphilopatryinbelugawhales t r C ( .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 o n g e ia th P re d y dn tta e b a p h t ed ie l z ra w t o n en o l e 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 l ,itc ,sn i tea ittsa itao taum red s t u tsr o re on eh i n f e d t B th se d r a n d b ,a n fo , a s d 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 w a o 0 g h n 2 ± u r f l o o 88 e f e ize 19 b s c if u s i l c a le om a v p r P - f p m d in g sa i o ) in a r (B d e h p c on it e i p w d lo s te re taa cea i r 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 : d e so vo ed isr r A D e g 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 c a 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 5 t A A 0 i . e ( N 0 n n a m i r MA e k u S lu n B sa aL C K a 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 i r e B 11/32 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) (S3 Table). Seasonal migration 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 stratum. Dispersal 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 east Chukotka east Chukotka east Chukotka east Chukotka east Chukotka east Chukotka east Chukotka east Chukotka east Chukotka Little Diomede Little Diomede Little Diomede Little Diomede Little Diomede Little Diomede Little Diomede Little Diomede Little Diomede Little Diomede Little Diomede Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope Point Hope LOD² ²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 [ 22 ], comparison of the sex ratio for the most thoroughly sampled population, the eastern Chukchi Sea (M:F± 1.56:1) [ 22 ], 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). Discussion 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 [ 31 ], and contrasts with earlier genetic studies [ 27 ] that 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 25 GLGs. 25 GLGs. 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 [ 34, 69 ]. Below we discuss the investigation's findings and their implications for beluga whale behavioral ecology and management in more detail. Demographic history 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 movements [ 17, 31 ]. 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 [70] in concert with inherent natal homing tendencies of beluga whales (see below) will likely act to maintain distinct regional populations in the future. 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 [4], 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 [ 22 ] 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 [ 31 ] (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 [ 2, 78 ], fish [79], reptiles [ 76 ], and mammals [ 80 ]. In 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 [4] 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 [4], it is not direct evidence of it. Colbeck et al. [ 81 ] 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 [ 82 ] 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 [ 83 ] 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. [ 27 ] 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. [ 75 ] 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. [ 39 ] 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 [ 22 ] and the time of year when these discrete subpopulations are in closest proximity [ 31 ] (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 [ 31 ]. 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. Seasonal migration 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 [32] and recent telemetry studies [ 31 ] 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. [ 31 ] 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. Temporal anomalies 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 [ 20 ]. This was followed by decades of inconsistent returns [ 85 ]. 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 [ 33, 83 ]. 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 [4] 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 [86]. 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 [ 18, 31 ], 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 [88] 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 threats. 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 [ 89, 90 ], 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 change [ 88, 91 ]. 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. Supporting information S1 Supporting information. This is a file that contains supporting information text. (DOCX) 25 / 32 S1 Fig. Median Joining network of beluga whale mtDNA haplotypes across their North Pacific Ocean. (TIF) 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). (TIF) 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. (TIF) 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. (TIF) 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. (XLSX) 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. (XLSX) 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. (XLSX) 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. (XLSX) 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. (XLSX) 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. (XLSX) 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. (XLSX) 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. (2013). (DOCX) Acknowledgments 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. Author Contributions 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. 2002. 63. 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


This is a preview of a remote PDF: http://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0194201&type=printable

Greg O’Corry-Crowe, Robert Suydam, Lori Quakenbush, Brooke Potgieter, Lois Harwood, Dennis Litovka, Tatiana Ferrer, John Citta, Vladimir Burkanov, Kathy Frost, Barbara Mahoney. Migratory culture, population structure and stock identity in North Pacific beluga whales (Delphinapterus leucas), PLOS ONE, 2018, DOI: 10.1371/journal.pone.0194201