Ancestral Polymorphisms and Sex-Biased Migration Shaped the Demographic History of Brown Bears and Polar Bears

Recent studies have reported discordant gene trees in the evolution of brown bears and polar bears. Genealogical histories are different among independent nuclear loci and between biparentally inherited autosomal DNA (aDNA) and matrilineal mitochondrial DNA (mtDNA). Based on multi-locus genomic sequences from aDNA and mtDNA, we inferred the population demography of brown and polar bears and found that brown bears have 6 times (aDNA) or more than 14 times (mtDNA) larger population sizes than polar bears and that polar bear lineage is derived from within brown bear diversity. In brown bears, the effective population size ratio of mtDNA to aDNA was at least 0.62, which deviated from the expected value of 0.25, suggesting matriarchal population due to female philopatry and male-biased migration. These results emphasize that ancestral polymorphisms and sex-biased migration may have contributed to conflicting branching patterns in brown and polar bears across aDNA genes and mtDNA.


Introduction
Genealogical discordance is highlighted in phylogenetic relationships between brown bears (Ursus arctos) and polar bears (U. maritimus). Phylogenetic analyses based on mitochondrial DNA (mtDNA) indicate that the root of polar bear lineage is within brown bear diversity, representing a paraphyletic relationship [1,2]. In contrast, two recent studies based on autosomal DNA (aDNA) found early divergence of polar bears from the common ancestor (i.e. monophyly) of brown bears by exploring one species tree from multi-locus data [3,4]. Hailer et al. [3] sequenced 14 independent aDNA loci in geographically diverse samples from brown bears and polar bears and estimated that polar bears split off 603 thousand years ago (ka), followed by brown bear divergence at 125 ka. The authors concluded that discrepancies in phylogenies between aDNA and mtDNA were the result of mtDNA introgression from brown bears to polar bears due to hybridization. Based on genome scale sequence data, Miller et al. [4] suggested the ancient divergence of polar bears from brown bears with an estimated split time of 4 to 5 million years ago (MYA). The authors concluded that extensive gene flow had maintained the genetic relationship between brown bears and polar bears until 200 to 100 ka and that polar bear mtDNA was replaced with brown bear mtDNA.
In contrast to previous reports on discordant gene trees between aDNA and mtDNA, Nakagome et al. [5] pointed out that genealogical histories from aDNA genes were compatible with mtDNA history and that the majority of discrepancies could be traced to stochastic processes in which ancestral polymorphisms and subsequent incomplete lineage sorting shaped a unique gene genealogy. Based on sequence data from Hailer et al. [3], Nakagome et al. [5] simulated the stochastic process under the coalescent model [6] and tested whether genealogy at each locus significantly deviated from the paraphyletic model of mtDNA phylogeny. The results were completely different from Hailer et al. [3], and none of the aDNA genes could reject the paraphyletic model. Indeed, we observed conflicting gene genealogies among the aDNA loci (see Figure 1 and Figure S1), which likely result from different patterns of ancestral polymorphisms and incomplete lineage sorting. Gene trees from eight and four out of 13 loci analyzed in the previous study [5] showed paraphyletic and multifurcating topology, respectively. Although one locus, AZIN1, depicted monophyletic topology, only one mutation defined the brown bear clade. The origin of brown bears and polar bears remains unresolved, and further studies are required to reconcile these dramatically different results.
There are two types of discrepancies in the evolution of brown bears and polar bears, one among unlinked aDNA genes and the other between aDNA and mtDNA. The first discrepancy can be attributed to ancestral polymorphisms and incomplete lineage sorting [7]. Given enough times (,5N e generations, where N e is the effective number of chromosomes), DNA sequences in a population are likely to coalesce within each population and any loci can be expected to depict a consistent gene tree with the species tree. Because N e determines the amount of ancestral polymorphisms and the probability of coalescence in a population, gene tree discordance among loci depends on N e . In contrast, the second discrepancy may correlate with behavioral traits of female philopatry and male-biased dispersal in brown bears [8,9]. A recent study pointed out that sex-biased migration is a plausible cause of inconsistent bear phylogenies among aDNA, mtDNA, and sex-linked DNA [10]. To address these issues, we inferred a demographic history of brown and polar bears using aDNA and mtDNA.

DNA Sequence Data
We used two different datasets of population-based sequences from autosomal DNA (aDNA) genes and complete mitochondrial DNA (mtDNA) ( Table S1). All sequence data were downloaded from the European Molecular Biology Laboratory (EMBL) Nucleotide Sequence Database (http://www.ebi.ac.uk/embl/) and GenBank (http://www.ncbi.nlm.nih.gov/genbank/). aDNA sequences included brown bears (2n = 36) and polar bears (2n = 36) from 14 loci [3]. Brown bear sequences were derived from broad geographic regions in Norway, Romania, Sweden, Bulgaria, U.S.A., and Canada (see Table S1 in Hailer et al. [3]). mtDNA sequences were from brown bears (n = 9) and polar bears (n = 26). To infer an ancestral allele at each segregating site among brown and polar bears, we used orthologous sequences from the giant panda for 14 aDNA loci [11] and the complete mtDNA sequences of one giant panda (the accession number: EF212882) and two American black bears (JX196366; NC003426).
In the mtDNA genome (except for the D-loop region), some of the sequences contained multiple substitutions, which did not follow the infinite site model in population genetic analyses. Under the infinite site model, one mutation event can occur at a previously un-mutated site and can be disrupted by back mutation or recurrent mutation events. We detected incompatibilities between mutated sites using the four-gamete test implemented in DnaSP version 5.1 [12,13]. If mutations are incompatible, we can infer recurrent mutation at the sites because mtDNA is the haploid, non-recombination sequence. To minimize the effects of multiple substitutions on population genetic analyses, we generated two mtDNA datasets by excluding (1) sequences (Set-I) or (2) multiple substitution sites (Set-II). (1) We excluded sequences with multiple substitutions from brown bears or polar bears and aligned the mtDNA sequences. However, the aligned sequences still contained multiple substitutions spanning from the 8,269 nucleotide position (np) to 15,439 np, and approximately half of the mtDNA sequences were trimmed. As a result, the 8,268 bp mtDNA dataset (the sequence from 1 np to 8,268 np) was comprised of data from four brown bears and ten polar bears (Table S1). (2) We excluded multiple substitution sites from alignment of 35 mtDNA sequences (nine from brown bears and 26 from polar bears), with a total sequence length of 15,403 bp.

Reconstruction of Gene Trees for aDNA Loci
We estimated gene genealogy for 13 aDNA loci from brown bears (2n = 36) and polar bears (2n = 36) using GENETREE [14]. Because the LRGUK locus included a high frequency of recombinant sequences and did not follow the infinite site model, this locus was excluded from GENETREE analysis. Genealogy was deduced from the path of mutations to the most recent common ancestor (MRCA) under the infinitely many-site model [15]. Under this model, a gene is thought to be an infinite sequence of completely linked sites where mutations occur at sites that have never experienced mutations before. Given the gene tree, we computed maximum likelihood estimates for the  population mutation rate and generated the empirical distribution of time to MRCA (T MRCA ), as described in a previous study [5].

Bayesian Inference for Demographic History
To estimate demographic parameters, we applied kernel approximate Bayesian computation (kernel-ABC) [16,17] to the sequence data from aDNA or mtDNA. ABC is a statistical framework to obtain approximation of a posterior estimate without likelihood function. Sequence data are summarized into a set of summary statistics. The advantage of kernel-ABC is that the cost of computing time is independent of the dimension of data, and thus high-dimensional summary statistics can be incorporated into the analyses. In this study, we used site frequency spectrum (SFS) and haplotype frequency spectrum (HFS) as the set of summary statistics. SFS consists of the site frequency, j i , which is the number of segregating sites in which the mutant nucleotide is present on i sequences in the sample. SFS reflects a pattern of mutations among segregating sites and refines a partition of data determined by conventional summary statistics, such as the number of segregating sites, nucleotide diversity [18], and Tajima's D [19]. A previous study showed that SFS can improve approximation of the posterior estimate given the full  Table S2). doi:10.1371/journal.pone.0078813.g002 data compared with conventional statistics [16]. However, HFS consists of the haplotype frequency, d j , which is the number of haplotypes observed in j sequences in the samples. HFS can account for recombination patterns at a locus. Because sequence data from aDNA include many recombinants, the combination of SFS and HFS provides more detailed information than SFS alone. We summarized the sequence data into SFS for brown bears (SFS uar ) or polar bears (SFS uma ) and two-dimensional HFS (2D-HFS) in which each haplotype was shared between brown bears and polar bears or specific to the population. For aDNA data, we calculated SFS uar , SFS uma , and 2D-HFS for each locus ( Figures  S2a-S2n) and merged the data into a set of summary statistics across 14 loci. The allelic state (ancestral/derived allele) at each segregating site was determined by alignment with an orthologous giant panda sequence. As mtDNA is a haploid genome, we generated 2D-SFS in which each derived allele was shared between or specific to either brown bears or polar bears ( Figures  S2o and S2p). Because it is generally difficult to determine the allelic state using the giant panda sequence due to a higher mutation rate in mtDNA than aDNA, we used two mtDNA sequences from American black bears with phylogenetic positions that are closer to brown and polar bears compared with the giant panda [20]. If both American black bear sequences had the same allele for the segregating site as brown and polar bears, it was defined as the ancestral allele. Otherwise, we used the allele that was consistent with the giant panda sequence. Although some of the sites may contain back or recurrent mutations, most of the sites likely follow the infinite site model. Based on the observed summary statistics for aDNA and mtDNA, we estimated posterior means of parameters by kernel-ABC.
The demographic model represents population divergence between brown bears and polar bears, and the parameters in the model are the effective population size in brown bears (N uar ), polar bears (N uma ), and divergence time between brown bear and polar bear populations (T). Prior distributions for each parameter were represented as a log-normal distribution (LN) with a variance as a square of the mean (m). We used different sets of prior distributions between aDNA and mtDNA as follows:  , where n is the number of simulations) were selected based on previously described 10-fold cross validation [16,17]. This algorithm was repeated 100 times and the mean and standard deviation (S.D.) of the posterior mean estimate for each parameter were calculated.
All simulations were performed using the program package ms, which generates samples from the coalescent model [21]. In coalescent simulation of aDNA genes, we assumed that a recombination rate was equal to a mutation rate. The mutation rate was calculated from the average number of substitutions between the giant panda and brown and/or polar bears assuming that divergence of the giant panda and brown/polar bears is 12 MYA, which represents the oldest remains from a member of the giant panda lineage and is compatible with the molecular clock estimate [22,23]. For aDNA genes, we averaged mutation rates across 14 loci to derive an estimate of 1.314610 28 bp/site/ generation. The number of substitutions in mtDNA was corrected by the Tajima-Nei distance model [24], and the mutation rate was estimated to be 7.036610 28 bp/site/generation.

Test of Natural Selection in mtDNA
To detect signals of natural selection in mtDNA sequences (Set I and Set II) in brown bears, we calculated Tajima's D [19], Fu and Li's D and F [25], and Fay and Wu's H [26] using DnaSP version 5.1 [12]. American black bear sequences were used as outgroup to determine ancestral states of alleles at segregating sites. The significance of these statistics was evaluated by the null distribution of 10,000 coalescent simulations under the constant size model using DnaSP [12].

Estimation of T MRCA
We calculated T MRCA in the genealogy between each of the geographic lineages in brown bears (Kenai, Kodiak, France, or ABC Islands) and polar bear lineages by kernel-ABC. mtDNA sequence data for Set-I and Set-II were summarized into 2D-SFS ( Figures S3 and S4). The demographic parameter is T MRCA = T n +T n-1 +Á Á Á+ T 2 , where T i follows an exponential distribution with the parameter of i(i -1)/2 under the coalescent model. Based on the algorithm described in the previous section,

Results and Discussion
The demographic model assumed in this study is population divergence between brown bears and polar bears, including three parameters of effective population size in brown bears (N uar ), polar bears (N uma ), and divergence time between brown bear and polar bear populations (T). To encompass the genetic diversity of widely distributed brown bears, we considered 14 aDNA loci (18 brown bears and 18 polar bears) from Hailer et al. [3] and complete mtDNA sequences (9 brown bears and 26 polar bears) from Miller et al. [4] and GenBank (see Table S1). In this study, we used kernel-ABC to incorporate high dimensional summary statistics and improve approximation of posterior estimate given data [16,17]. mtDNA genomes included multiple substitutions and did not follow the infinite site model in which each mutation occurs at a previously un-mutated site. We generated two mtDNA datasets by excluding sequences (Set-I) or multiple substitution sites (Set-II) and analyzed both datasets to evaluate the effects of sequenceediting methods on demographic inferences.

Demographic History based on Autosomal DNA
To obtain posterior estimates of the demographic parameters given the aDNA data, we summarized the sequence data into SFS and HFS. Summary statistics were calculated for each locus and merged into a set of observed summary statistics across 14 loci ( Figure S2). Based on observed data and prior conditions described above, we estimated posterior means of the parameters by kernel-ABC (Table 1). N uar was estimated at 46,434 (standard deviation (S.D.), 520), which is approximately six times larger than N uma = 7,867 (S.D. 252). Divergence time was estimated to be 136,597 (S.D. 2,264) generations ago, which was translated into 1.37 MYA based on an assumed generation time of 10 years. Nakagome et al. [5] showed that the time to most recent common ancestor in brown bears (T MRCA-uar ) was more than T = 1.37 MYA in eight out of 13 loci, while T MRCA in polar bears (T MRCAuma ) was less than the divergence time in most of these loci. Our population genetic analysis for aDNA genes supports the paraphyletic model of brown and polar bear lineages (Fig. 1).
We examined the robustness of posterior estimates relative to prior distributions. We used a set of 10 different prior distributions for each demographic parameter and compared posterior means of the parameter among prior conditions ( Fig. 2 and Table S2).

Demographic History based on mtDNA
To resolve gene tree discordance between aDNA and mtDNA, we estimated demographic parameters using mtDNA sequences. We used two mtDNA datasets (Set-I and Set-II) as the observed data in kernel-ABC. Set-I consists of 8,268 bp sequences from four brown bears and 10 polar bears, while Set-II includes 15,403 bp from nine brown bears and 26 polar bears (Table S1). Each of the datasets was summarized into a matrix of SFS (2D-SFS) ( Figure  S2). The estimated value of N uar /2 or N uma /2 was defined as the effective population size of mtDNA (N mtDNA-uar or N mtDNA-uma ). The posterior estimate of N mtDNA-uar (Set-I: mean 6 S.D. = 67,13261,405; Set-II: 57,7776726) was approximately 29 or 14 times larger than N mtDNA-uma (Set-I: 2,3516113; Set-II: 4,232641) ( Table 1). The divergence time between brown bear and polar bear populations (Set-I: 30,9926607 generations ago; Set-II: 32,8476337 generations ago) was estimated to be about 0.31 MYA and 0.33 MYA, respectively. These results suggest that posterior estimates from Set-I are compatible with those from Set-II and that our estimations do not likely depend on sequenceediting methods. Under the equilibrium assumptions of random mating, an equal number of males and females, and Poisson distribution of offspring, the effective population size based on mtDNA (N mtDNA ) is expected to be a quarter of that based on aDNA (N aDNA ). The ratio of N mtDNA to N aDNA (N mtDNA /N aDNA ) in polar bears was 0.149 and 0.269 in Set-I and Set-II, respectively, which was relatively close to the expected value of 0.250 (Table 1). In contrast, the N mtDNA /N aDNA ratio in brown bears was 0.723 and 0.622 in Set-I and Set-II, respectively, which was much larger than the expected value.

Effects of Sex-biased Migration
Deviation from the expected ratio could have resulted from the following three forces: (i) differences in mutation rates between aDNA and mtDNA; (ii) natural selection on mtDNA; and (iii) sexbiased demographic events leading to different effective population sizes of males and females. For (i), we calculated mutation rates for aDNA and mtDNA based on the average number of substitutions between the giant panda and brown/polar bears. The number of substitutions in mtDNA was corrected using the Tajima-Nei distance model [24]. As expected, the mutation rate was higher in mtDNA (Set-I: 7.036610 28 bp/site/generation; Set-II: 7.838610 28 bp/site/generation) than in aDNA (1.314610 28 ), although this difference was already incorporated into estimation of N mtDNA and N aDNA . For (ii), we tested natural selection of mtDNA in brown bears using neutrality statistics (Table S3). However, none of the statistics showed signals indicating natural selection. Evolutionary forces of (i) and (ii) are unlikely to result in a high N mtDNA-uar /N aDNA-uar ratio. For (iii), male-biased dispersal and female philopatry are known behavioral traits in brown bears [8,9]. A recent study proposed the presence of sex-biased migration in ancestral populations of bears, which led to high or low genetic divergence of mtDNA or Y-linked DNA in bear speciation [10]. Female philopatry has been shown to increase N mtDNA relative to N aDNA [27], and male-biased migration has contributed to higher genetic differentiation in mtDNA than aDNA [28]. Hence, the observed N mtDNA-uar /N aDNA-uar value can be explained by sex-biased migration in brown bears.
Female philopatry may lead to matriarchal population stratification in geographically different regions (Kenai, Kodiak, France, and ABC Islands). Based on Set-I and Set-II datasets, we considered mtDNA genealogy between each of the geographic lineages in brown bears and polar bears. First, we counted the number of mutations specific to brown bear, ancestral polar bear, or polar bear lineages ( Figures S3 and S4). More mutations in the brown bear lineage have accumulated in Kenai or Kodiak, followed by France and the ABC Islands. Next, we estimated T MRCA for each genealogy by kernel-ABC based on 2D-SFS. The  D. 1,576), which was lower than estimates from brown/polar bear comparisons. All of the T MRCA estimates based on Set-II were generally consistent with those based on Set-I ( Figure S4). Therefore, the T MRCA Kenai-uma , T MRCA France-uma , and T MRCA ABC-uma estimates likely reflect split times in the mtDNA genealogy between each of geographic lineages in brown bears and the ancestral lineage of polar bears. These results represent the highly differentiated mtDNA lineages of brown bears in each geographic region.

Evolutionary History of Brown Bears and Polar Bears
Biparentally inherited aDNA and matrilineal mtDNA show that the polar bear lineage derived from diverse brown bear lineages (Fig. 3). However, the divergence times between brown bear and polar bear populations are different between aDNA (1.37 MYA) and mtDNA (Set-I: 0.31 MYA; Set-II: 0.33 MYA) demographic histories. This discrepancy likely reflects different effects of sexbiased migration on aDNA and mtDNA genealogies. The aDNA genealogies are products of female philopatry and male-biased migration, while the mtDNA genealogy has only been shaped by female philopatry, and evolutionary traces of brown bears and polar bears have been recorded in mtDNA lineages. Geographic differentiations of brown bear lineages were shown in the mtDNA genealogy, although these patterns were likely homogenized in aDNA genealogies by male-biased migration between geographic regions, as suggested in a previous study [10]. Therefore, the mtDNA genealogy represents patterns of population divergence in the speciation history of polar bears from brown bears. The divergence time based on aDNA is close to T MRCA Kenai-uma (Set-I: 1.12 MYA; Set-II: 1.20 MYA), showing that the ancestral lineage of polar bears in mtDNA genealogy split from the most distant brown bear lineage at this time. In contrast, the divergence time from mtDNA is almost equal to T MRCA ABC-uma (Set-I: 0.33 MYA; Set-II: 0.42 MYA), suggesting that the ancestral lineage of polar bears split from the closest brown bear lineage at this time. A plausible explanation for the evolution of brown bears and polar bears is as follows: the initial divergence of brown bear lineages occurred at 1.12 MYA; geographic lineages of brown bears gradually diverged from the ancestral lineage of other brown bears and polar bears and expanded into broad geographic regions; after that, the polar bear lineage appeared from the ancestral lineage at 0.33 MYA and differentiated into the arctic region. This recent divergence is compatible with the oldest fossil record of polar bears, which is estimated to be between 110 ka and 130 ka [29]. Hence, the genetic diversity of brown bears and polar bears can be shaped by ancestral polymorphisms and sex-biased migration.
In conclusion, Bayesian inference based on the coalescent model highlights the fact that population demography underlies patterns of genetic variation across independent loci even though stochastic effects shape a unique genealogy at each locus. Furthermore, sex-biased migration has a significant effect on population demography [30,31], and we observed that female philopatry and male-biased dispersal generated discordant genealogy in brown and polar bears between aDNA and mtDNA. aDNA, mtDNA, and sex-linked DNA have different modes of inheritance and thus reflect their own history of evolutionary demographics and natural selection. As population-based genomic sequencing becomes available to analyze a variety of wild animals, population genetic approaches can be applied to phylogenetic problems to model stochastic processes and inheritance characteristics involved in speciation history.  . The number of bins in Set I or Set II is defined by the full spectra in brown bears (n = 4) and polar bears (n = 10) or by the full spectrum in brown bears (n = 9) and the given intervals in polar bears (0, 1-3, 4-6, 7-9, 10-12, 13-15, 16-18, 19-21, 22-24, 25-26). The bins with a frequency of ''0'' (Set I or Set II) and a fixed frequency (Set I) were excluded (black).