Circumpolar Genetic Structure and Recent Gene Flow of Polar Bears: A Reanalysis

<div><p>Recently, an extensive study of 2,748 polar bears (<i>Ursus maritimus</i>) from across their circumpolar range was published in PLOS ONE, which used microsatellites and mitochondrial haplotypes to apparently show altered population structure and a dramatic change in directional gene flow towards the Canadian Archipelago—an area believed to be a future refugium for polar bears as their southernmost habitats decline under climate change. Although this study represents a major international collaborative effort and promised to be a baseline for future genetics work, methodological shortcomings and errors of interpretation undermine some of the study’s main conclusions. Here, we present a reanalysis of this data in which we address some of these issues, including: (1) highly unbalanced sample sizes and large amounts of systematically missing data; (2) incorrect calculation of <i>F</i>_<i>ST</i> and of significance levels; (3) misleading estimates of recent gene flow resulting from non-convergence of the program BayesAss. In contrast to the original findings, in our reanalysis we find six genetic clusters of polar bears worldwide: the Hudson Bay Complex, the Western and Eastern Canadian Arctic Archipelago, the Western and Eastern Polar Basin, and—importantly—we reconfirm the presence of a unique and possibly endangered cluster of bears in Norwegian Bay near Canada’s expected last sea-ice refugium. Although polar bears’ abundance, distribution, and population structure will certainly be negatively affected by ongoing—and increasingly rapid—loss of Arctic sea ice, these genetic data provide no evidence of strong directional gene flow in response to recent climate change.</p></div

Hierarchical analysis of molecular variance (AMOVA) for nuclear microsatellites among management units within the six genetic clusters identified in this paper and shown in Table 2.

<p>For this analysis, we followed Peacock <i>et al</i>., 2015 [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.ref020" target="_blank">20</a>] by including northern Davis Strait in the Eastern Archipelago cluster and southern Davis Strait in the Hudson cluster. Western Laptev was included in the Western Basin cluster and Eastern Laptev was included in the Eastern Basin cluster. However, results did not differ significantly when the Laptev Sea and Davis Strait MUs were excluded entirely.</p

Genetic diversity statistics for 18 management units of polar bears.

<p>For microsatellite data, a maximum of 30 individuals have been retained from each management unit from the original dataset of 2,748 individuals, and only the 14 loci indicated in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.g001" target="_blank">Fig 1</a> have been used. Molecular diversity indices for mitochondrial DNA were calculated in Arlequin using pairwise differences with no gamma correction.</p

CLUMPAK-averaged Structure outputs for 20 independent runs of <i>K</i> = 2–7, which were clustered and averaged using CLUMPAK.

<p>Numbers under each K-value indicate the proportion of runs that converged to the solution presented. Minority modes supported by at least two runs are provided in Fig A in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.s001" target="_blank">S1 File</a>. Management unit abbreviations are as in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.t001" target="_blank">Table 1</a>. </p

Circos plots of gene flow using 14 or 21 loci among: three clusters corresponding to our Structure results for <i>K</i> = 4 (excl. Norwegian Bay samples), four clusters identified in Peacock <i>et al</i>., 2015 [20] (excl. Norwegian Bay samples), and five clusters corresponding to our Structure results for <i>K</i> = 6 (excl. Norwegian Bay samples).

<p>Segment colours are as in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.g004" target="_blank">Fig 4</a> and are sized proportionally to the population size estimates in Table D of <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.s001" target="_blank">S1 File</a>, though polar bear population sizes are estimated with very broad confidence intervals, particularly in the Polar Basin, where reliable estimates are not available for most MUs. The width of each ribbon where it meets a segment on the circumference indicates the proportion of migrants into (but not out of) each region. Black ribbons are significantly directional based on non-overlapping 95% CIs of immigration rates; grey ribbons are not significantly directional.</p

Sampling locations for 476 of the 495 polar bears used in this analysis; the remainder did not have lat–long coordinates.

<p>Individuals are colour-coded by genetic cluster similarly to the colour scheme for <i>K</i> = 6 in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.g002" target="_blank">Fig 2</a>. Black samples are unassigned (i.e., <i>Q</i>_<i>max</i> < 0.5). Uncoloured individuals are those that were used in the original study but were not included in our random subset of 30 individuals per management unit; their predicted cluster memberships based on BAPS trained clustering are shown in Fig C of <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.s001" target="_blank">S1 File</a>. Management unit abbreviations are as in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.t001" target="_blank">Table 1</a>. Approximate sea ice extent during the breeding season is shown using measurements for April 15, 2008 [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.ref078" target="_blank">78</a>], though there is great spatial heterogeneity in sea ice thickness and concentration, as well as great intra-seasonal and inter-annual variability. Note that this map (and the data) does not reflect a 2014 boundary change between NB and SB made by the territorial governments and the co-management boards with management authority for these two subpopulations, because it has not yet been evaluated by the IUCN Polar Bear Specialist Group.</p

Circumpolar Genetic Structure and Recent Gene Flow of Polar Bears: A Reanalysis - Fig 3

<p>(a) Structure output of mean likelihood ± SD calculated from 20 independent runs for each value of <i>K</i> from 1 to 20. (b) <i>ΔK</i> calculated using the Evanno method in CLUMPAK. (c) Probability of K calculated using the Pritchard method in CLUMPAK.</p

Comparisons of pairwise <i>F</i>_<i>ST</i> values for nuclear microsatellites with <i>Φ</i>_<i>ST</i> values from mitochondrial DNA; the line indicates the expectation of under isolation as given in [56].

<p>(a) S4 Fig from Peacock <i>et al</i>., 2015 [<a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.ref020" target="_blank">20</a>], reproduced here under the terms of its Creative Commons CC0 license. (b) A recreated version of this graph, generated using our recalculated <i>F</i>_<i>ST</i> and <i>Φ</i>_<i>ST</i> values. In (b), M’Clintock Channel, Norwegian Bay, and Viscount Melville Sound were excluded because of inadequate sample sizes for mitochondrial DNA (<i>N</i> ≤ 3) and the Laptev Sea was excluded because this MU was significantly out of Hardy–Weinberg equilibrium. Points for brown bears were not recalculated and are not shown. Coloured points indicate intra-cluster MU pairs (coloured as in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.g004" target="_blank">Fig 4</a>); grey points indicate inter-cluster MU pairs.</p

Missing data in Peacock <i>et al</i>., 2015 [20].

<p>The size of the square at each management unit–locus intersection is proportional to the amount of missing data at that locus in that management unit. Management unit abbreviations are as in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.t001" target="_blank">Table 1</a>. Asterisks denote loci that were retained for the reanalysis presented in this paper.</p

Hierarchical analysis of molecular variance (AMOVA) for mitochondrial DNA among management units within the six genetic clusters identified in this paper and shown in Table 2.

<p>Note that many management units (incl. the entire Norwegian Bay cluster) were excluded entirely from this AMOVA because of inadequate sampling. Because we lacked sample location information for downloaded haplotypes, we were unable to split Davis Strait or the Laptev Sea into northern/southern or eastern/western samples; therefore, these MUs were removed for this calculation in addition to the MUs that were removed for low sample sizes in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0148967#pone.0148967.t002" target="_blank">Table 2</a>.</p