Greg O’Corry-Crowe, Robert Suydam, Lori Quakenbush, Brooke Potgieter, Lois Harwood, Dennis Litovka, Tatiana Ferrer, John Citta, Vladimir Burkanov, Kathy Frost, Barbara Mahoney

PLoS ONE JOURNAL

Abstract

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, R_st = 0.09–0.13), and likely demographic independence of (F_st-mtDNA = 0.02–0.66), and in many cases limited gene flow (F_st-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 (P_BAYES = 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 (r_winter-F_st-mtDNA = 0.9, r_summer-F_st-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.

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. These, in turn, have implications for gene-culture coevolution and create novel challenges for management and conservation, including the identification of management units, the assignment of migrating animals to population of origin, and the assessment of risk at the population level.

The primary factors influencing population subdivision in migratory species include; 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. 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. species with extended periods of parental care natal philopatry and thus population divergence may involve parentally directed learning of migration routes termed migratory culture. Structuring in migratory species also depends on which sex disperses with sex-biased dispersal potentially leading to the demographic isolation of population units even in the face of gene flow and seasonal overlap. Patterns of dispersal and heterogeneity in migratory species can be shaped by geography and other environmental factors (e.g., water temperature, weather patterns, sea-ice cover) that directly dictate the path and timing of migration, for example, and thus the level and duration of population overlap. And finally, glacial history, especially at higher latitudes, where the historical location of refugial populations and the sequence in which new habitat emerged following deglaciation can have a lasting effect on contemporary patterns of distribution and migration, and thus dispersal and population structure.Disentangling the contributions of such factors to migration, dispersal and structuring in migratory species is challenging. Genetic investigations offer huge potential in resolving the interplay between behaviors, the environment, and demographic history in shaping seasonal movements, dispersal and population structure in these species.

The beluga whale, Delphinapterus leucas, is a mid-sized cetacean of Arctic and sub-Arctic waters where all these factors likely interplay to shape contemporary migration, dispersal, breeding, and population structure patterns. In this study we use spatial and temporal patterns of genetic variation and individual relatedness to test a series of hypotheses about these aspects of belugas in the Pacific.

Across much of their range belugas migrate between wintering areas at or near the southern margin of the sea-ice and summering grounds in seasonally more open water farther north where they feed, molt and raise their young. Many Arctic populations are migratory with some completing annual circuits in excess of 6,000km while subarctic populations are less so, often exhibiting substantial overlap between winter and summer ranges. Highly gregarious, beluga whales congregate by the thousands at several geographically discrete nearshore locations following ice-breakup in summer. Breeding is believed to occur primarily in winter and early spring, and adult males may select different habitats than females and younger animals at certain times of year. These attributes, combined with the potential for intersecting migration routes and for separate summering groups sharing the same wintering ground during the period of peak mating, contribute to complex patterns of dispersal and population structure in beluga whales, raise the possibility of stable migratory cultures, and create unique challenges for species management, including the identification and monitoring of stocks.

In the Pacific, beluga whales inhabit three geographically distinct regions; the Gulf of Alaska, the Sea of Okhotsk, and the Bering, Chukchi and Beaufort Seas, termed the BCB region. Within the BCB region, there are six major summer aggregation areas near shore that are labeled in as follows: Bristol Bay and Norton Sound in the eastern Bering Sea, Kotzebue Sound and off Kasegaluk Lagoon in the eastern Chukchi Sea, Mackenzie Delta-Amundsen Gulf region of the eastern Beaufort Sea, and Anadyr Bay in the western Bering Sea (Figs and). There are also three discrete summering locations in the Okhotsk Sea, and a small resident population in Cook Inlet and isolated group in Yakutat Bay in the Gulf of Alaska.

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. 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. 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. Considering the peak breeding season for beluga whales is during winter and early spring, 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 leading to their identification as demographically distinct management stocks. A few studies have documented lower levels of nDNA (microsatellite) heterogeneity among geographic strata compared to mtDNA which has been taken as evidence of extensive male-mediated gene flow among summering groups, possibly on shared wintering grounds.

These earlier studies, however, have a number of limitations with respect to resolving patterns of migration, dispersal and gene flow of Pacific beluga whales basin-wide. All were of too short a duration to directly test for natal philopatry and migratory culture in this long-lived mammal, and to assess the stability of migration and dispersal patterns over ecological time frames. All had a regional focus that risked incomplete assessments of the molecular and behavioral ecology of this highly vagile species in an environment with few geographic barriers. The failure to include all populations, for example, can lead to inaccurate cluster analyses and erroneous assignment tests. Low marker number and sample size in many cases limited interpretations of observed patterns of differentiation. Furthermore, differing laboratory conditions prevented the comparison of nDNA (microsatellite) data sets among studies and populations, while few whales were sampled on migration or outside traditional coastal summering areas.

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 for a total of 1,647 whales from across their Pacific range. We test hypotheses on how the spatial proximity of summering and wintering areas may influence levels of dispersal and interbreeding among discrete summering aggregations, some of the findings of which have recently been referred to in the tracking paper by Citta et al. We include previously unsampled and undersampled locations, greatly expand the analysis of migrating whales, and extend the sampling time frame in many areas to encompass more than three decades from 1978–2010. With this extended time frame we assess patterns of kinship within groups and among years, determine the stability of migration patterns and structure over ecological time frames, and further investigate the population origins of beluga whales in one area (Kotzebue Sound, Alaska) that has witnessed dramatic changes in the pattern of annual return of beluga whales in summer. Finally, we explore how model-based cluster analysis of population genetic structure and Bayesian inference of recent rates of dispersal may be influenced by sample size, uneven sample numbers, locus number and differing degrees of subdivision.

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. A more detailed summary of sample numbers across the four distinct sampling periods of this study are provided in. Tissues were preserved in 20% dimethyl sulfoxide (DMSO) saturated with NaCl. Total DNA was extracted from each sample using standard cell lysis-high-salt extraction protocols, and each sample was screened for variation within 410bp of the mtDNA control region according to previously published methods. Total DNA was also extracted from a number of dried tissue and teeth samples using standard silica-based ancient DNA methods. The gender of each sample was determined by PCR-based methods and all samples were analyzed for variation within eight microsatellite loci typed on beluga whales or other cetacean species. Earlier studies including an analysis of 288 belugas in the current study (including 9 known mother-calf pairs), demonstrated that all eight loci were inherited in a Mendelian fashion and were not sex-linked. Comparisons of observed to expected genotypic frequencies using the Micro-Checker program (version 2.2.3) for the four best sampled areas found only one of 32 tests consistent with null alleles (i.e., homozygote excess), and no evidence of large allelic dropout (0/32 tests) or scoring errors due to stuttering (0/32 tests).

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 software and by estimating both haplotypic (H) and nucleotide (π) diversity using Arlequin 3.5 software. Estimation of genetic diversity (H_e, H_o and number of alleles) and probabilities of identity (I) for nuclear loci was performed using Cervus 3.0. Tests for Hardy-Weinberg (H-W) and linkage equilibria were performed in Genepop 4.1 with P-values estimated from 500,000 iterations of the data using the Mark chain method of Guo and Thompson. An analysis of both mtDNA and nDNA data revealed that the populations were not in mutation-drift equilibrium.

We used both frequency-based (F_st) and distance-based statistics (Φ_st, R_st) to calculate the degree of genetic differentiation among spatial and temporal strata. Fixation indices were estimated using an analysis of variance framework, and homogeneity tests based on 50,000 permutations of the data were performed in Arlequin. The model used for estimating the evolutionary distances between pairs of mtDNA sequences was simple pairwise differences. In the case of nuclear loci, genetic distances were based on stepwise mutation models. Confidence intervals for F_st 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.

The Bayesian model-based clustering method, Structure 2.3.4, was used to infer population structure and the likely number of populations (K), and to assign individuals to population of origin. We conducted 28 distinct analyses with unique parameter settings: analyses were run both with and without admixture, with and without sampling location included as prior information, and they contained sample sets of varying size that were randomly sampled from the complete sample set using R (n = 15, 30, 50, 100, and all) that included individuals successfully genotyped at a range of loci (n ≥ 6 loci, n ≥ 7 loci, n = all 8 loci, see. For each parameter set we used a burn-in period of 50,000 iterations followed by 2x10⁵ iterations to collect data. Multiple runs (n = 5) were conducted for each value of K (n = 8) to ensure convergence for a total of 1,120 separate runs. The models of Hubisz et al. use sampling location as prior information to reveal further underlying structure without risking detecting structure that is not present. To identify which individuals were likely immigrants (or descendants of recent immigrants) to their populations, we defined prior probabilities that each individual was an immigrant (v = 0.01–0.10) and incorporated information on the geographic origin of each sample before running the clustering analysis. We used Clumpak to integrate results from multiple runs for different values of K and to implement methods for choosing K.

To further explore the origins of individual whales, we assigned individuals to populations based on estimated likelihoods of their genotype and/or mtDNA haplotype arising in each of the sampled populations under assumptions of random assortment of alleles and independence of loci using the programs Doh and whichrun 4.1. A log₁₀ ratio (LOD) score of the most likely allocation to the second most likely of ≥1 (i.e., ratio of ≥10) denoted confidence in the assignment to a single population.

Differences in allelic diversity among populations can result in genotypes with higher calculated likelihoods in populations other than their source (i.e., nominal) population, even though the genotype may still be relatively common in the source population. In such cases an exclusion test is also required to ascertain likely population of origin. We thus examined the estimated likelihoods of genotypes relative to those of other genotypes within each baseline population. Using the Assignment Calculator in the Doh program and assuming Hardy-Weinberg Expectations (HWE), 1,000 new individuals were generated from the gene pools of each population, and the natural log of the estimated probability of an individual’s genotype was compared to the likelihood distribution of the generated genotypes.

Dispersal patterns over ecological time scales were also investigated using the Bayesian approach of Wilson and Rannala  in the program BayesAss 1.3. 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 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 (F_st) at mtDNA and levels of differentiation and average relatedness at nuclear markers (F_st, F_is, r, mAIc, vAIc) for males and females. The tests rely on sampling individuals post-dispersal in order to exclude the homogenizing effects of past immigration and interbreeding by either sex through the sampling of their offspring. We assumed dispersal occurs primarily at the juvenile and sub-adult stage and so only adult cohorts were used. The randomization approach of Goudet et al. was used to test whether there was significant contemporary sex-biased dispersal. Using Fstat 2.9.3 we compared observed differences to a null distribution based on 10,000 randomizations of the data.

We applied the Bayesian stock-mixture method of Pella and Masuda to assess the population composition of migrating herds. This method incorporates uncertainty in the estimation of stock proportions of the ‘mixture’, or in this case migrant samples, and unlike conditional maximum likelihood methods improves stock determination by using information in the mixture sample to improve estimates of baseline relative frequencies. Using the program Bayes, we ran multiple chains with different starting parameters (including different prior stock proportions) and used the univariate statistic, the shrink factor, to monitor convergence of chains to the posterior probability. Analyses were conducted separately on mtDNA and nDNA data, and on mtDNA-and nDNA data combined.

Finally, because philopatry should promote high levels of relatedness among individuals sampled across years and generations within an area we tested this by comparing estimates of average relatedness and proportions of first and second-order relationships (based on nDNA) within and between years using the programs coancestry and ml-relate.

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

A total of 1,116 individual whales were successfully screened for polymorphism at six or more microsatellite loci. 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. Average expected heterozygosity (H_e) for the 8 loci ranged from 0.511 to 0.799 across the geographic strata. 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 Delta-Amundsen Gulf, Canada, and the probabilities of identity (P(ID)) for each area were on the order of 10⁻⁵ to 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 (Φ_st = 0.105–0.321) and frequency-based (F_st = 0.079–0.285) statistics. 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 (F_st = 0.106–0.507), the only exception being the low level of mtDNA heterogeneity observed between Kotzebue Sound and Mackenzie-Amundsen (F_st = 0.022). The three primary summer concentration areas within the Sea of Okhotsk also displayed substantial mtDNA heterogeneity (F_st = 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 F_st 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 and among most summer concentration areas and resident populations. As with mtDNA, Cook Inlet was the most distinct of all geographic strata for microsatellite variation ( = 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. 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.

There was widespread agreement across the various parameter settings used in the model-based cluster analysis. 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 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 ≈1.0). These population clusters corresponded to Cook 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. Likewise, Anadyr Bay did not form a discrete population cluster in these analyses with individuals assigned to 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.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. 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. Stronger correlations with genetic heterogeneity were found among wintering areas (r_mtDNA-F_st = 0.9, p = 0.007) compared to summering areas (r_mtDNA-F_st = 0.1, p = 0.351) and for mtDNA (r = 0.1–0.9) compared to nDNA (r = 0.12–0.27).

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. Of the likely migration destinations, Whichrun assigned the majority of Point Hope whales (95.4%) and whales from the other three areas (85.7%) to the Mackenzie-Amundsen concentration area (i.e. eastern Beaufort Sea) for mtDNA. 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). Breaking assignments out by marker type revealed that mtDNA had a greater influence than any other locus on cumulative likelihoods.

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. 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). 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. 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. 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.

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. 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 of dispersal within BCB over such a short time frame (a few generations) using this method.

In the assessment of sex-biased gene flow through comparison of relative ratios of nDNA to mtDNA differentiation, estimates of microsatellite heterogeneity (F_st) for those pair-wise comparisons involving Cook Inlet were, on average, 7.1 times lower than corresponding estimates for mtDNA. 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 (F_st 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 (F_st ratio = 55.6) due to the very low levels of microsatellite differentiation observed between these two areas, a finding also reflected in the clustering analysis.

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 (F_{st males} = 0.218 v. F_{st females} = 0.235). Similarly, using fstat allelic frequency differences at nuclear markers between the eastern Beaufort and eastern Chukchi were significantly lower for adult males compared to adult females (F_{st males} = 0.010 v. F_{st females} = 0.022, p = 0.043, while estimated average relatedness within each area was significantly lower for males compared to females (r_males = 0.0204 v. r_females = 0.0419, p = 0.042). These sex differences were more pronounced in large (i.e., older) adults.

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. Of all the individual whales identified as possible migrants or of having recent migrant ancestry, the majority (74%) were adult males. While the number of males sampled in each population tends to exceed that of females, comparison of the sex ratio for the most thoroughly sampled population, the eastern Chukchi Sea (M:F– 1.56:1), to that for likely migrants or migrant descendants identified here (M:F– 5.67:1) revealed a significant male bias (χ² = 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 ( = 0.017 ± 0.1, p >0.05), and was significantly lower than corresponding mean values among areas ( = 0.308 ± 0.15; t-test: p = 0.0003). Results were similar for nDNA, while Structure assigned whales from each summering concentration area sampled in different time periods into the same population clusters. As with the original combined dataset, Norton Sound exhibited low differentiation (F_st) from, and clustered with, both Bristol Bay and the eastern Beaufort Sea. The one exception was Kotzebue Sound, where substantial mtDNA differentiation (F_st = 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. 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. 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, , within years did not differ significantly from 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. 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).

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 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, and documented stability in migration and dispersal behavior over ecological (i.e., decadal) time frames with notable exceptions. A much expanded nDNA dataset revealed a pattern of limited interbreeding among many distinct subpopulations within the BCB region where the spatial proximity of winter ranges and spring migration routes appears to have a greater influence on the degree of both individual and genetic exchange than summer ranges. This study supports recent satellite telemetry evidence of a series of mostly distinct, if overlapping, wintering areas within the BCB region, and contrasts with earlier genetic studies 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 and gene flow of this long-lived, highly vagile species than previously possible. The model-based cluster analysis and assignment tests were greatly improved by larger sample sizes and locus number, informed run assumptions, and the inclusion of all potential population clusters. The limited information content in the microsatellite data within regions severely limited the utility of the program BayesAss to estimate recent rates of dispersal, something recognized by a number of other recent studies. Below we discuss the investigation’s findings and their implications for beluga whale behavioral ecology and management in more detail.

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. These basin-wide findings expand on earlier regional-level studies and are supported by recent satellite telemetry studies that have mapped winter as well as summer movements. Together, these investigations indicate that the Alaska and Kamchatka Peninsulas have been effective barriers to dispersal and gene flow over time. Projections for continued climate warming across the Arctic in concert with inherent natal homing tendencies of beluga whales (see below) will likely act to maintain distinct regional populations in the 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 Lagoon emerged as discrete population clusters early on in the various cluster analyses, that is, at low values of K 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 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, and see, it has rarely been explicitly tested. It relies on the assumptions of a simple drift-dispersal model to explain predominantly frequency-based (F_st) differentiation. However, multiple factors can complicate interpretations of mtDNA subdivision in terms of contemporary dispersal patterns: local populations may not be in equilibrium, sampling may be inadequate, and movement patterns in migratory species may change across years and generations (see below). Using large sample sizes collected across many years we identified very low numbers of females and males as likely migrants or of migrant ancestry and recorded consistent patterns of mtDNA heterogeneity among most areas over decades of sampling. With age of first reproduction estimated at 8.3yr for Nearctic belugas this time series likely spanned multiple generations and supports philopatry to natal population or summering group as a dominant trait in beluga whales. The discovery that substantial numbers of closely related whales, including parent-offspring pairings, were sampled across years, and even decades, within summering areas is direct proof of natal homing.

A recent study summarizing findings from multiple satellite telemetry investigations revealed that beluga whales from separate summering grounds across the BCB region return to somewhat geographically discrete wintering areas in the northern Bering Sea. 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 or level of interbreeding (see below).

Natal philopatry has been documented or inferred in a wide range of migratory species and long-distance travelers, including birds, fish, reptiles, and mammals. 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. Whitehead recently declared that the evidence for culture in non-humans is rare and that documenting cultural stability (in non-humans) is difficult. Furthermore, while genetic differentiation may reflect the faithful transmission of migratory culture, it is not direct evidence of it.

Colbeck et al. recently found that groups of closely related beluga whales often migrate together in the Hudson Bay region of Canada and posited that this could facilitate the learning of migration routes. Matthews and Ferguson confirmed that beluga whale calves are dependent on their mothers for two or more years, and thus multiple migratory cycles. Such findings when viewed in conjunction with our findings of limited dispersal among areas, and substantial numbers of close kin not only within herds and years at discrete summering areas 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, and thus their degree of overlap, has likely also been shaped by physical factors that directly shape migration pathways (e.g., sea ice, geography) and by historical factors such as the routes of postglacial colonization. Any shift in environmental parameters could presumably disrupt migration routes and seasonal habitat use patterns (see below), and/or the level of philopatry to those demographically distinct circuits. A recent study revealed that beluga whale movements are resilient to inter-annual variation in sea ice conditions which may be related to cultural inheritance. However, the study did observe that anomalies in beluga whale migration patterns were correlated with anomalous ice years.

Dispersal and gene flow

The increased statistical power of the expanded data set that we used allowed a more quantitative assessment of sex-biased dispersal than earlier studies. While both sexes tend to be philopatric, we found that when dispersal does occur it is predominantly by adult males. Male-biased dispersal is typical in mammals and has been interpreted to be a strategy to reduce competition for mates and avoid inbreeding. 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. In light of the present study’s findings it appears 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. We also, however, rejected panmixia on a common wintering ground. This differs from an earlier microsatellite study in this region by Brown Gladden et al. based on a much smaller dataset (i.e., less locations, samples, years and loci). It also differs from recent findings on beluga whales in Canada and the Russian Far East. Turgeon et al. found no convincing evidence of nDNA differentiation among demographically discrete summering assemblages in the Hudson Bay region concluding this likely reflected near random mating on a common breeding area. Similarly, Meschersky et al. found no evidence of nDNA differentiation among two of the three demographically distinct summering aggregations in the Sea of Okhotsk (also see mtDNA findings in 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 and the time of year when these discrete subpopulations are in closest proximity. For example, while patterns of mtDNA variation in Norton Sound were consistent with a basin-wide 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. While the summering grounds of these two subpopulations are far from Norton Sound (1,000km and 2,000km, respectively) a recent telemetry study revealed that the wintering range of the Norton Sound subpopulation may overlap to some degree with that of the eastern Beaufort Sea and Bristol Bay whales. Similarly, low nDNA differentiation and clustering might suggest the same for the Anadyr Bay subpopulation and beluga whales in both the eastern Beaufort and the Sea of Okhotsk. However, low sample size from the latter precludes a well-supported inference at this time.

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, 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. This agrees with earlier genetic investigations and recent telemetry studies where the spring migration of eastern Beaufort whales occurs earlier and through denser sea ice than eastern Chukchi Sea belugas. The discovery that a few individual whales at some of these spring locations had higher likelihoods of eastern Chukchi Sea ancestry or of mixed-ancestry, however, indicates that the Bering Strait region is also an area of mixing in spring. This is also evident in recent movement data where Citta et al. observed that for tagged eastern Beaufort Sea whales to migrate north in spring through the Bering Strait earlier than the eastern Chukchi belugas they had to pass through the latter’s primary wintering area.

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. This was followed by decades of inconsistent returns. Earlier genetic studies found evidence that the pre-decline whales were likely a demographically distinct subpopulation and that more recent occurrences of belugas in the Sound in the 1990s and 2000s most likely involved whales from the Beaufort Sea. These studies, however, did not include all potential source subpopulations in the BCB region and so the origins of Kotzebue whales remained provisional. The current study included all major summer aggregations in the BCB region and confirmed that the pre-decline summer assemblage in Kotzebue Sound was likely a demographically discrete subpopulation and that beluga whales entering the Sound post-decline in the 1990s and 2000s (at least for those years where we had large sample sizes) were different from those we originally sampled and were most likely comprised of whales primarily from the eastern Beaufort Sea sub population.

Evolution, ecology and management

The macro-geographic patterns of genetic differentiation in beluga whales across the North Pacific suggest divergence of the evolutionary trajectories of regional populations on millennial timescales. Regional adaptation of beluga populations is likely over such time frames while rescue of small and/or declining populations like that in Cook Inlet by whales from other regions is probably remote. The stability of heterogeneity among spatial strata on decadal time scales within regions and the evidence for generational-scale kinship within strata may indicate the kind of cultural stability in beluga whale populations that is necessary for gene-culture coevolution where distinct cultures drive the evolution of neutral and functional genes over much faster time scales.

The study revealed that all spatially discrete summering aggregations within regions are demographically independent subpopulations and as such should be defined as separate stocks under management objectives to maintain healthy, functioning populations across the species range. The greatly increased dataset over previous studies now means that five of the six summering aggregations within the BCB region can be more definitively defined as the following stocks: Bristol Bay, Eastern Bering Sea (Norton Sound), eastern Chukchi Sea (Kasegaluk Lagoon), eastern Beaufort Sea (Mackenzie-Amundsen), and Gulf of Anadyr (Anadyr Bay). 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. We define stocks, for example, by their primary coastal summering location.

This study’s findings, in concert with new satellite telemetry data, force us to alter this perception. The demographically distinct summering aggregations within the BCB region, for example, in most cases: do not appear to interbreed extensively, return to discrete wintering areas, and 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 (N_e) and how to model risk, recovery and population viability. Furthermore, when the migratory circuits of these subpopulations are considered in light of the heterogeneity of Pacific Arctic and subarctic marine ecosystems the individual subpopulations likely also have quite distinct ecologies and perhaps other unique aspects to their behavior and population biology. It follows then that different migratory populations are likely exposed to different threats.

Finally, the study confirms the critical importance of the Northern Bering Sea and Bering Strait region to beluga whales. This relatively small area is home to tens of thousands of beluga whales from multiple subpopulations, some among the largest in the world, for example, for much of the year. This region is also important for other marine mammal species (Citta et al. in review.) and is currently undergoing dramatic environmental and ecosystem change. 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.