Demographic History of the Brown Bear (Ursus arctos) on Hokkaido Island, Japan, Based on Whole-Genomic Sequence Analysis.

Previous studies of the brown bear (Ursus arctos) on Hokkaido Island, Japan, have detected three geographically distinct subpopulations representing different mitochondrial lineages and shown that gene flow between subpopulations has occurred due to male-biased dispersal. In this study, we determined whole-genomic sequences for six Hokkaido brown bears and analyzed these data along with previously published genomic sequences of 17 brown bears from other parts of the world. We found that the Hokkaido population is genetically distinct from the other populations, keeping genetic diversity higher than the endangered populations in western Europe but lower than most populations on the continents. A reconstruction of historical demography showed no increase in population size for the Hokkaido population during the Eemian interglacial period (130,000-114,000 years ago). In a phylogenetic analysis of the autosomal data, the Hokkaido population formed a clade distinct from North American and European populations, showing that it has maintained genetic diversity independently from continental populations following geographical isolation on the island. This autosomal genetic similarity contrasts with the geographically separate mitochondrial lineages on Hokkaido and indicates the occurrence of male-driven gene flow between subpopulations.

Significance

The demographic history of brown bears, which are distributed widely in the northern hemisphere is complex, so it is still unclear although the previous studies using various DNA sequences have been done. This study examined the bear population on the Hokkaido Island, Japan, with the whole-genome analysis, and compared it with the continental ones. As a result, it was clarified that this population which, has been geographically isolated by the small island in the Far East, has experienced a different history of migration and demography from continental ones. Our results can contribute to further understanding of the population demography of brown bears.

Introduction

The brown bear (Ursus arctos) is widely distributed in the northern hemisphere (McLellan et al. 2017) and has become a model system (Davison et al. 2011) for studying mammalian evolution at the population level. Brown bear populations show complex demographic histories that reflect geographical isolation (Randi et al. 1994), bottlenecks caused by human activity (Waits et al. 2000), hybridization between related species (Pongracz et al. 2017), and phenotypic adaptation to a variety of environments (Sato et al. 2011; Colangelo et al. 2012). Various previous studies have shed light on these processes, though not always reaching the same conclusions. Studies based on mitochondrial cytochrome b and control region sequences (Davison et al. 2011) and a complete mitogenome analysis (Hirata et al. 2013) have shown that several brown bear lineages across the northern hemisphere. Analyses of Y-chromosome sequences, however, showed no phylogeographic correlation with the mitochondrial DNA (mtDNA) lineages (Bidon et al. 2014; Hirata et al. 2017). Based on whole-genome sequences, Cahill et al. (2013, 2015, 2018) and Liu et al. (2014) reported gene flow from polar bears (Ursus maritimus) included within the same mtDNA lineages as brown bears from the Admiralty, Baranof, and Chichagof (ABC) islands to brown bears, whereas Miller et al. (2012) proposed gene flow from brown bears to polar bears, during speciation. Similarly, Barlow et al. (2018) showed that gene flow occurred from the now-extinct cave bear (Ursus spelaeus), which is placed to the outgroup of the brown and polar bears, to brown bears. In an analysis of several nuclear genes, Hailer et al. (2012) found that brown bears are divided into two main clades, the European and North American lineages. These various data suggest that four major factors have contributed to brown bear evolution: 1) expansion several times across the entire northern hemisphere, 2) male dispersal contributing to the genetic structure of current populations, 3) hybridization with related species (polar bear and cave bear) after dispersal, and 4) relative isolation of local populations leading to genetic homogeneity within these populations.

Three different mtDNA lineages, clades 4, 3a2, and 3b appear to have migrated from the Eurasian Continent to Hokkaido Island, Japan, at different times (Hirata et al. 2013) and become allopatrically distributed in southern, central, and eastern Hokkaido, respectively (Matsuhashi et al. 1999; Hirata et al. 2013). Some gene flow, however, occurs between subpopulations within Hokkaido due to males dispersing more widely than females (Hirata et al. 2017). On Hokkaido, allopatric distribution and admixture between subpopulations have thus both been involved in structuring the brown bear population.

The goal of this study is to further clarify the demographic history of brown bears, especially those in Hokkaido by analyses of whole-genomic sequences obtained by next-generation sequencing technology. Here we report on the genetic diversity, genetic structure and genealogy, demographic dynamics of the Hokkaido brown bear based on whole-genomic sequences; compare our results with those from previous studies on other populations; and discuss the demographic history between the Hokkaido bears and continental populations.

Results

Whole-Genome Sequencing and Genetic Diversity

Our analysis included the whole-genomic sequences we obtained from six brown bears from Hokkaido (accession numbers: DRR276774–DRR276779) and previously reported sequences of 17 individuals from other areas (fig. 1 and supplementary table S1, Supplementary Material online). The proportion of the genome covered was 96.6–99.4% (supplementary table S1, Supplementary Material online). The number of single nucleotide polymorphisms (SNPs) found per individual ranged from 3,726,121 to 4,748,290, with a depth range of 8.61 to 54.87 (supplementary table S1, Supplementary Material online). Supplementary figures S1 and S2, Supplementary Material online, show the heterozygosity and nucleotide diversity, respectively, in 50-kb sliding windows. Both values were lower for the Hokkaido individuals than for most individuals from other populations but higher than in endangered populations in western Europe.

Fig. 1.

Sampling locations for brown bears included in this study. The size of circles shows the number of samples from various regions. For Hokkaido, sample numbers from previous studies (Matsuhashi et al. 1999; Hirata et al. 2017) are shown. The colors of the circles for Hokkaido individuals indicate the mtDNA lineage: white, southern Hokkaido (clade 4); dark gray, central Hokkaido (clade 3a2); and light gray, eastern Hokkaido (clade 3b).

Phylogenetic Tree Based on the Mitogenome

A phylogenetic tree (fig. 2) reconstructed using complete mtDNA sequences shows two main clades, one comprising bears from western Europe, a polar bear, and bears from southeastern Alaska (ABC islands), and the other comprising bears from Alaska (Kenai and Denali), Hokkaido, and eastern Europe. Within the latter clade, the Hokkaido individuals comprise three distinct lineages, clade 4, 3b, and 3a2, from southern, eastern, and central Hokkaido, respectively, and an individual from Georgia was included in neither clades 3a1 nor 3a2.

Fig. 2.

Phylogenetic tree based on the mitogenome (A). Numbers near nodes are bootstrap values in percent. Clade designations to the right refer to Hirata et al. (2013). The Alaskan brown bear was removed from the result because its sequences data (SRR7758718) do not include mitochondrial DNA data. Neighbor-joining tree based on genetic distances for autosomal genomic data (B).

Genetic Structure and Genealogy

A principal component analysis (PCA) for all brown bears examined (fig. 3) showed separation among Hokkaido, European, and North American individuals. In particular, the first principal component (horizontal axis) indicated strong separation between the Hokkaido population and other populations. In a PCA for only the Hokkaido individuals (supplementary fig. S3, Supplementary Material online), the first principal component (horizontal axis) showed separation of the three known subpopulations. The second principal component (vertical axis) separated southern Hokkaido 1 from southern Hokkaido 2, whereas in each case the two individuals from central Hokkaido and those from eastern Hokkaido were closely related to one another.

Fig. 3.

Result of the PCA for all brown bears. The horizontal axis shows the first principal component, and the vertical axis the second principal component.

Figure 4 shows the genetic structure of the bears analyzed, based on ADMIXTURE. At K = 3, which was the most likely number of clusters (fig. 4), three clusters, Hokkaido, Europe, and North America are evident. At K = 2 (fig. 4), which was the second most likely number of clusters, the Hokkaido individuals are separate from the others. Figure 4 shows the probability that the individual from Russia belongs to both the European population and other populations.

Fig. 4.

Result of ADMIXTURE analyses at K = 3 (A) and K = 2 (B). Populations are labeled at the bottom.

The Fst value between Europe (Slovakia, Georgia, southern Sweden 1, southern Sweden 2, Slovenia, northern Italy, central Italy, and Spain) and North America (Alaska, Denali, Kenai, Admiralty 1, Admiralty 2, Baranof, Chichagof 1, and Chichagof 2) was lower than that between Hokkaido and the others (table 1), supporting the results shown in figures 3 and 4. Fst values among the Hokkaido subpopulations (supplementary table S2, Supplementary Material online) were consistent with a closer relationship between southern and central Hokkaido, as indicated by the results of the first principal component of PCA (supplementary fig. S3, Supplementary Material online).

Table 1

Autosomal Fst Values Among Groups of Bears Representing the Hokkaido, European, and North American Populations

	Hokkaido	Europe
Hokkaido	—	—
Europe	0.228 (0.0012)	—
North America	0.262 (0.0014)	0.156 (0.0008)

Note.—Europe includes eight individuals (Slovakia, Georgia, southern Sweden 1, southern Sweden 2, Slovenia, northern Italy, central Italy, and Spain). North America includes eight individuals (Alaska, Denali, Kenai, Admiralty 1, Admiralty 2, Baranof, Chichagof 1, and Chichagof 2). Hokkaido includes six individuals, two each from central, eastern, and southern Hokkaido. Standard deviation (SD) values are in parentheses.

In the neighbor-joining tree (fig. 2) constructed by using identity-by-state distances between autosomal genomes, the Hokkaido and North American individuals comprise distinct clades embedded within a paraphyletic European group, with the individual from Russia most closely related to the Hokkaido clade which was also suggested in figure 4.

A plot of ƒ3 statistics (supplementary fig. S5, Supplementary Material online) against geographical locality showed that the Russia, Kenai, and Denali individuals are more closely related to the Hokkaido population than are the other populations examined, which was also supported by the result of ƒ4 (Asiatic black bear, individuals from Denali, Kenai, and Russia; others from Europe or North America, individuals from Hokkaido), although the relationship between bears from Denali and Kenai, and other North American individuals was closer (supplementary table S3, Supplementary Material online).

In trees constructed by using TreeMix (Supplementary fig. S4Supplementary Material online), the three groups appear as distinct clades. Both trees in supplementary figure S4, Supplementary Material online, assuming no genetic flow and supplementary figure S4, Supplementary Material online, assuming that one gene flow occurred showed that the Russian individual lies outside the clade including the European individuals. In addition, it indicates gene flow between the Asiatic black bear and the brown bears of Georgia.

We applied ƒ4 statistics to clarify the genetic relationship between the continental population and the Hokkaido population. Values of ƒ4(Asiatic black bear, European populations; North American populations, individuals from Hokkaido) were calculated to investigate the genetic relationship between populations in Europe (Slovakia, Georgia, Sweden, Slovenia, northern Italy, central Italy, and Spain) and the Hokkaido population (supplementary table S4, Supplementary Material online), and values of ƒ4(Asiatic black bear, North American populations; the Russian individual, individuals from Hokkaido) were calculated to detect the genetic relationship between populations in North America (Alaska, Denali, Kenai, Admiralty, Baranof, and Chichagof) and those in Hokkaido (supplementary table S5, Supplementary Material online). Although the values of ƒ4(Asiatic black bear, European populations; North American populations except Denali and Kenai, Hokkaido) indicated slightly close relatedness between Hokkaido individuals and European populations (ƒ4 = 0.0003500–0.001510, |Z score| > 3 except in two cases), the relatedness between European populations and individuals from Denali and Kenai (ƒ4 = −0.001644 to −0.000436, |Z score| > 3) was closer. Bears from Hokkaido and that from North America were not closely related to each other (ƒ4 = −0.0019050 to −0.001230, |Z score| > 3).

Values for ƒ4(Asiatic black bear, X; Hokkaido 1, Hokkaido 2), where X indicates brown bears other than those from Hokkaido, were calculated to estimate allele-sharing bias within the Hokkaido populations. The analysis suggested that southern Hokkaido 1 shares more alleles with European brown bears than with other Hokkaido individuals (fig. 5 and supplementary table S6, Supplementary Material online) although only ƒ4(Asiatic black bear, Spain; South Hokkaido 1, eastern Hokkaido 1) and ƒ4(Asiatic black bear, Sweden; South Hokkaido 1, eastern Hokkaido 1) were statistically significant (|Z score| > 3) and both Z scores were slightly above 3, indicating that the statistical significance was marginal.

Fig. 5.

Range of values for ƒ4 (Asiatic black bear, X; Hokkaido 1, Hokkaido 2). For each boxplot, the boxes range from the 25th to 75th percentiles and the lines extend to 1.5 times the distance from the 25th to 75th percentile. The lines inside the boxes indicate the median value. The dots indicate ƒ4 values, with different shading indicating the population from which X derives: black, European population (Slovakia, Georgia, Sweden, Slovenia, northern Italy, central Italy, and Spain); dark gray, North American population (Alaska, Denali, Kenai, Admiralty, Baranof, and Chichagof); light gray, Russian population. Hokkaido 1 and Hokkaido 2 indicate individuals from Hokkaido. Lower values mean that more alleles are shared between Hokkaido 1 and X. Asterisks indicate mean the range of values of ƒ4(Asiatic black bear, the European population; southern Hokkaido 1, Hokkaido 2).

Demographic Dynamics

Historical demography was reconstructed based on the pairwise sequentially Markovian coalescent (PSMC). This analysis indicated that the effective population size (Ne) of most European and North American brown bears increased about 120 ka (thousands of years ago) (fig. 6 and supplementary fig. S6, Supplementary Material online). In contrast, none of the six Hokkaido individuals indicated an increase in population size at that time (fig. 6).

Fig. 6.

PSMC estimates of brown bear effective population size (Ne) through time. Comparison of results including individual central Hokkaido 1 and two individuals from Europe and North America (Slovakia and Baranof) (A). Results for all Hokkaido individuals (B). The vertical gray shading indicates the Eemian interglacial period (130–114 ka).

Discussion

Genetic Diversity

Heterozygosity and nucleotide diversity values both indicated that the Hokkaido brown bears are higher in genetic diversity than the endangered brown bears in Europe, although the Hokkaido population has long been isolated from populations in continental Eurasia. This result is incongruent with previous studies based on major histocompatibility complex class-II DQA genes (Goda et al. 2009), DNA fingerprinting (Tsuruga et al. 1994), and protein polymorphism (Tsuruga et al. 1996), in which the Hokkaido brown bears showed lower genetic diversity than other mammals. We conclude that the diversity of the whole genome has been maintained at a higher level than that of only coding regions as shown in the previous studies, as has also been reported for macaques (Osada et al. 2015). Previous studies have shown that the genetic diversity of the whole genome declines in cases of 1) long-term population decline and inbreeding (Xue et al. 2015; Benazzo et al. 2017), and 2) strong recent inbreeding (Prufer et al. 2014; Benazzo et al. 2017). Although the Hokkaido bears have experienced a population decline since the Last Glacial Period (70,000–10,000 years ago; fig. 6), our results indicated no events leading to a decline of genetic diversity in the whole genome.

Demography of Brown Bears

Matsuhashi et al. (1999) and Hirata et al. (2013) showed that three mtDNA lineages are distributed allopatrically in southern, central, and eastern Hokkaido. Our study also showed the same three lineages and the corresponding allopatric distributions. Most previous mtDNA analyses found that clade 3a included only two clades, the East European and Alaskan lineage (clade 3a1) and the central Hokkaido lineage (clade 3a2; Hirata et al. 2013; Tumendemberel et al. 2019). However, an APLP analysis by Hirata et al. (2014) detected another haplotype from western Iran and the Caucasus (eastern 2) in clade 3a that was distinct from clades 3a1 and 3a2. Consistent with this result, the haplotype from Georgia was included in clade 3a but not in clade 3a1.

Our PCA and the genetic cluster analyses revealed apparent genetic differentiation among Hokkaido, European, and North American bears, with the Hokkaido individuals separate from both the European and North American groups. The first principal component of PCA, the clustering of ADMIXTURE at K = 2, and Fst values all suggest that the nuclear genome of brown bears on Hokkaido differs from that in the continental populations. Analyses by Hirata et al. (2013) based on the whole mitogenome indicated that clade 3a2 in central Hokkaido was the last of the three lineages (clades 4, southern; 3b, eastern; 3a2, central) to migrate to Hokkaido, and that clade 3a1 (the sister group to clade 3a2) spread widely from eastern Europe to Alaska via the Bering region after that. The fossil record is congruent with this result because no fossil record of clade 3a brown bears has been confirmed prior to 10 ka in North America (Leonard et al. 2000; Barnes et al. 2002), which could be later than the isolation of Hokkaido Island from the Eurasian Continent (Ohshima 1990).

In contrast, the neighbor-joining tree and ƒ3 statistics showed that the Hokkaido brown bears are most closely related to the individuals from Russia, Denali, and Kenai, reflecting some genetic affinity between the Hokkaido population and individuals occurring relatively close by Hokkaido. Moreover, ƒ4 statistics suggested the Hokkaido bears are more closely related to European bears than North American bears except for Denali and Kenai. The results infer clade 3a migration route from the Eurasia continent onto Hokkaido, Denali and Kenai via the land bridges (supplementary fig. S7, Supplementary Material online). Matsuhashi et al. (1999) and Hirata et al. (2013) concluded based on mtDNA data that ancestral bears migrated to Hokkaido from continental Eurasia via Sakhalin Island. The sea depth between Hokkaido and Sakhalin is relatively shallow, so Hokkaido had been a part of the Eurasia continent rather than an island before 12 ka (Ono 1990 ), which is consistent with the result in this study. Mizumachi et al. (2021) analyzed mtDNA sequences from ancient and modern samples of brown bears on Sakhalin and found that only clade 3a1, which is distributed widely on the continent, occurs on Sakhalin. Although the mtDNA haplogroups detected on Hokkaido are clearly different from those in continental Eurasia, admixture by male-biased dispersal between continental Eurasia and Hokkaido could have occurred prior to the opening of Soya Strait between Sakhalin and Hokkaido ∼12 ka (Ohshima 1990). In addition, ancestors of the Denali and Kenai bears are estimated to have migrated to North America via Bering land bridge after migration onto Hokkaido because clade 3a1, which is not found in Hokkaido, distributes around Denali and Kenai (Leonard et al. 2000; Barnes et al. 2002). The closer relationship between Denali–Kenai and the European bears from Eurasia than that between Hokkaido and European bears is reflected to the ƒ4 statistics.

Our neighbor-joining tree based on autosomal genomic data showed brown bears divided into three groups (Hokkaido, Europe, and North America), all of which include mtDNA lineage 3a. This suggests that populations having mitochondrial 3a haplotypes have dispersed across a wide range and contributed to current populations by hybridization among local populations. Although Matsuhashi et al. (2001) and Davison et al. (2011) found a clear allopatric distribution of mtDNA lineages, our results from the whole autosomal genome indicated that individuals from the same population may be closely related even if their mtDNA haplogroups are different, which is consistent with the results of Hailer et al. (2012). This pattern could have resulted from wide male-biased dispersal (McLellan and Hovey 2001). Bidon et al. (2014) analyzing Y-chromosome sequences indicated no allopatric distribution as shown in mtDNA analysis (Davison et al. 2011) and male-biased gene flow caused by sex-biased dispersal. It suggests that autosomal homogenization has occurred by male-biased gene flow within the groups after their dispersal.

More locally, values for ƒ4(Asiatic black bear, X; Hokkaido 1, Hokkaido 2) showed that southern Hokkaido 1 individual shared more alleles with European individuals than with other Hokkaido individuals. This suggests that the southern Hokkaido population has been isolated by a stronger geographic barrier than the other Hokkaido populations; the sampling location (Kikonai) for individual southern Hokkaido is located at the southernmost end of Hokkaido, where admixture with other Hokkaido populations may have been restricted. As a result, some genetic factors from the continent may be maintained. More samples are needed, however, to elucidate the genetic features of the brown bears in southern Hokkaido.

The results of PSMC detected no increase in Ne at around 120 ka for brown bear individuals from Hokkaido. The recovery in Ne observed for continental brown bears is generally correlated with increased temperatures in the Eemian interglacial period 130–114 ka (Miller et al. 2012; Cahill et al. 2013; Benazzo et al. 2017). The Hokkaido population appears to have undergone two demographic events. In the first event, ancestors of the Hokkaido population migrated to Hokkaido before the interglacial period and underwent a different demographic history from the continental population. Hirata et al. (2013) estimated that the southern Hokkaido (clade 4) and eastern Hokkaido (clade 3b) lineages derived from mtDNA continental lineages about 190 and 160 ka, respectively. Their population demography could have differed from that of continental populations after those times. The second demographic event was gene flow involving the autosomal genome into populations having different mtDNA haplogroups, as might be expected with male-biased dispersal; in this context, Hirata et al. (2017) demonstrated that male dispersal does occur between populations with different mtDNA haplogroups. The central Hokkaido population (clade 3a2) likewise did not show the Ne recovery, although Hirata et al. (2013) estimated that it migrated to Hokkaido 50 ka. By dispersal, alleles of the southern and eastern Hokkaido lineages have been shared with the central Hokkaido population.

Demographic History Scenario of Brown Bears in Hokkaido

Our study aimed to clarify the demographic history of brown bears in Hokkaido from four points of view: genetic diversity, genetic structure and genealogy, demographic dynamics, and gene flow. From our results, we concluded the following scenario and showed it in supplementary fig. S7, Supplementary Material online.

Brown bears migrated to Hokkaido before the Eemian interglacial period (supplementary fig. S7, Supplementary Material online) because the PSMC estimates indicated no historical increase in population size on Hokkaido. The divergence time of the southern and eastern Hokkaido lineages (clades 4 and 3b) estimated by Hirata et al. (2013) supports this conclusion.

Hybridization between different lineages has occurred following migration because the central Hokkaido lineage (clade 3a2), which is estimated to have migrated to Hokkaido after 50 ka (Hirata et al. 2013), also showed no Ne recovery in the interglacial period. The Y-chromosomal DNA analysis by Hirata et al. (2017) found no correlation between the genetic affinity of paternal brown bear haplotypes and the geographical distribution on Hokkaido. This suggests that populations having different mtDNA lineages underwent admixture due to male dispersal, contributing to the current structure of the Hokkaido population.

Isolation by distance (IBD) and geographic barriers have also affected the Hokkaido population, leading to hybridization bias between individuals. Southern Hokkaido 1 shared more alleles with the individuals from Europe than with other Hokkaido individuals, which reflects some maintenance of genetic factors from the continent. In Scandinavia, IBD has led to a strong correlation between genetic and geographic distances (Kopatz et al. 2014; Schregel et al. 2018), and the Hokkaido population could likewise be affected. Separation of the Hokkaido population into three subpopulations as indicated by the first principal component in the PCA analysis, and higher Fst values between eastern and central-southern Hokkaido individuals, indicate not only IBD but also some effect by a geographic barrier across the eastern–southern/central boundary, although what is the geographic barrier is still unknown. This difference in Fst values was also observed by Hirata et al. (2017), although the cause remains unknown.

Materials and Methods

Samples and Genome Sequencing

Our study utilized muscle tissue samples from six Hokkaido brown bears (one individual per sex per mtDNA lineage; fig. 1 and and supplementary table S1, Supplementary Material online) that had been collected in previous studies (Matsuhashi et al. 1999; Hirata et al. 2017). Total genomic DNA was extracted with the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) and stored at 4°C or −20°C until use.

Paired-end libraries were prepared by using the TruSeq DNA Nano Sample Prep Kit (Illumina, San Diego, USA), and whole-genome shotgun sequences were generated by Macrogen Japan (Tokyo, Japan), using an Illumina HiSeq X sequencer. Samples were sequenced to an average depth of 30×. FastQC ver. 0.11.8 (Babraham Bioinformatics 2011) and trimomatic ver.0.39 (Bolger et al. 2014) were used to check the quality of the sequences and remove Illumina adapters.

Mapping and Variant Calling

Whole-genomic sequence data from 17 brown bears (Miller et al. 2012; Cahill et al. 2013, 2015; Benazzo et al. 2017; Barlow et al. 2018; Taylor et al. 2018) were downloaded for analysis. Figure 1 and supplementary table S1, Supplementary Material online, indicate the sampling locations of these bears. One individual from Alaska (sample name: Alaska) is not shown in figure 1 because no information on the sampling location was available. The reads from 23 brown bears and 1 Asiatic black bear (Kumar et al. 2017; supplementary table S1, Supplementary Material online) were mapped to the polar bear reference genome (Liu et al. 2014; UrsMar_1.0, accession number: GCF_000687225.1) by using the BWA-MEM algorithm implemented in the BWA package v.0.7.17 (Li and Durbin 2009) with the “-M” command option. The number of reads aligned to the reference genome and the proportion of the genome covered was summarized by using CollectAlignmentSummaryMetrics in GATK ver. 4. 1. 2. 0 (McKenna et al. 2010). Single nucleotide variants were called by using HaplotypeCaller in GATK, and then further hard-filtered under the following parameters: significant fisher strand test > 60.0; variant confidence/quality by depth < 2.0; an RMS mapping quality (MQ) < 40.0; strand odds ratio > 9.0; MQRankSum < −12.5, and significant read position bias (ReadPosRankSum) < −8.0.

The reference polar bear genome includes thousands of short fragments (<200 bp), which are not suitable for many population genetic analyses, and so only SNPs on 356 scaffolds longer than 100 kb and no sex-linked scaffolds (Cahill et al. 2013; Bidon et al. 2014), and one mtDNA scaffold, were analyzed. Genmap ver. 1.3.0 (Pockrandt et al. 2020) was used to calculate genome mappability and exclude SNPs with (30, 2)-mappability < 1 (i.e., the 30-mer starting from the site is unique in the genome even allowing for 2 mismatches). To assess the reliability of sequence data, the mean depth of SNPs per individual was calculated by using vcftools ver. 0. 1. 17 (Danecek et al. 2011).

To assess genetic diversity, heterozygosity and nucleotide diversity were calculated in 50 kb sliding windows by using plink ver. 1.9 (Purcell et al. 2007) under the options “–het” and vcftools, respectively.

A phylogenetic tree based on the whole mitogenome was reconstructed with the maximum likelihood method implemented in MEGA X (Kumar et al. 2018). The optimal substitution model determined under the Bayesian information criterion was HKY + G (gamma distributed rate variation). The Asiatic black bear (supplementary table S1, Supplementary Material online) was included as an outgroup taxon because previous studies (Talbot and Shields 1996; Hirata et al. 2013) have found that polar bears group within clade 2, which includes brown bears from the ABC islands.

PCA and Bayesian clustering analyses were performed by using EIGENSOFT version 7.2.1 (Patterson et al. 2006) and ADMIXTURE ver. 1.3.0 (Alexander et al. 2009), respectively. SNPs under linkage disequilibrium (LD) were excluded before analyses by using plink options “–indep-pairwise 50 5 0.5.”

To construct a neighbor-joining tree, a distance matrix was calculated based on the identity-by-state plink option “–distance.” The tree was plotted by using the R package phangorn ver. 2. 5. 5 (Schliep 2011), with the Asiatic black bear (supplementary table S1, Supplementary Material online) as an outgroup. TreeMix ver. 1.13 (Pickrell and Pritchard 2012) was run to estimate patterns of population splits and mixtures in brown bear populations, with the polar bear reference genome (Liu et al. 2014) as the outgroup.

To infer which individuals were closely related to the Hokkaido individuals, ƒ3 statistics were calculated with the AdmixTools software package (Patterson et al. 2012). Fst values (Hudson et al. 1992; Keinan et al. 2007) were calculated by using EIGENSOFT.

To estimate the relationship among populations, ƒ4 statistics were calculated with the AdmixTools software package. The smaller the value of ƒ4(A, B; C, D) calculated using the allele frequencies of A, B, C, and D, the more alleles are shared between B and C (Patterson et al. 2012). Because polar bears admixed with brown bears in the past (Miller et al. 2012; Cahill et al. 2013; Liu et al. 2014), we used the sequence of the Asiatic black bear as the outgroup.

We tested three situations, ƒ4 (Asiatic black bear, Denali, Kenai, and Russia; other Europe or North America, Hokkaido), ƒ4(Asiatic black bear, Europe; North America, Hokkaido), ƒ4(Asiatic black bear, North America; Russia, Hokkaido), and ƒ4(Asiatic black bear, X; Hokkaido 1, Hokkaido 2), where X indicates brown bears other than those from Hokkaido, to detect the genetic relationship between bears from Europe and North America, and allele sharing bias between other brown bears on the continent and Hokkaido.

PSMC methods (Li and Durbin 2011) were applied to estimate past demography. The diploid consensus sequence for each sample was created from a bam file by using SAMtools ver. 1.9 (Li et al. 2009), with the “mpileup” command set to the “-C50” option. Each consensus sequence was transformed into the required format by using the fq2psmcfa tool in the PSMC package, and was analyzed with PSMC under the option “-N25 -t15 -r5 -p ‘4 + 25*2 + 4 + 6’.” The mutation rate and generation time were set at 1.82*10−8 mutations per site per generation and 11 years, respectively, as indicated by previous studies (Liu et al. 2014; Benazzo et al. 2017).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.