Polar bears exhibit genome-wide signatures of bioenergetic adaptation to life in the arctic environment.

Polar bears (Ursus maritimus) face extremely cold temperatures and periods of fasting, which might result in more severe energetic challenges than those experienced by their sister species, the brown bear (U. arctos). We have examined the mitochondrial and nuclear genomes of polar and brown bears to investigate whether polar bears demonstrate lineage-specific signals of molecular adaptation in genes associated with cellular respiration/energy production. We observed increased evolutionary rates in the mitochondrial cytochrome c oxidase I gene in polar but not brown bears. An amino acid substitution occurred near the interaction site with a nuclear-encoded subunit of the cytochrome c oxidase complex and was predicted to lead to a functional change, although the significance of this remains unclear. The nuclear genomes of brown and polar bears demonstrate different adaptations related to cellular respiration. Analyses of the genomes of brown bears exhibited substitutions that may alter the function of proteins that regulate glucose uptake, which could be beneficial when feeding on carbohydrate-dominated diets during hyperphagia, followed by fasting during hibernation. In polar bears, genes demonstrating signatures of functional divergence and those potentially under positive selection were enriched in functions related to production of nitric oxide (NO), which can regulate energy production in several different ways. This suggests that polar bears may be able to fine-tune intracellular levels of NO as an adaptive response to control trade-offs between energy production in the form of adenosine triphosphate versus generation of heat (thermogenesis).

Introduction

The polar bear (Ursus maritimus) is often thought of as a prime example of adaptive evolution in response to the extremes of life in the high Arctic. Polar bears face many challenges that result in high energetic demands. In the winter, males and nonpregnant females, which do not hibernate, must maintain a constant body temperature in an environment where external temperatures may regularly be as low as −50 °C. This is further compounded by winds, which may lead to convective losses of greater than 75% of the metabolic heat produced (Best 1982). Polar bear fur provides relatively poor insulation during extreme cold conditions (Øritsland 1970), and it has been suggested that the adipose tissue of polar bears is an adaptation for increased energy storage (Pond et al. 1992), potentially making thermal regulation difficult. In contrast, brown bears (U. arctos) and American black bears (U. americanus) at northern latitudes enter dens in the fall, providing thermal protection during the winter. The ground temperature in a black bear den near Fairbanks, AK, was found to be −9 °C when the external temperature was −46 °C (Folk et al. 1972). Similarly, brown bears were found to den for between 37 and 45 days longer in northern than in middle/southern locations (Manchi and Swenson 2005). Therefore denning behavior may help mitigate the challenges of thermoregulation in brown and black bears.

Another important energetic challenge that polar bears face relates to fasting during periods of food scarcity. In some areas of the Arctic, summer sea ice melt forces bears to move to shore where access to their primary high-fat marine mammal diet is limited. The length of these periods and the extent of fasting differ across their range, with relatively little known about fasting in polar bears that stay on multiyear ice in the northern portion of their distribution. In western Hudson Bay, Canada, near the southern extent of their range, polar bears may fast for periods of 4 months or longer (Stirling et al. 1999). These challenges are intensified for pregnant females, which may fast up to 8 months, beginning in the summer and extending through the winter, when they give birth and begin nursing young, all before returning to the sea ice in the spring to hunt again (Atkinson and Ramsay 1995; Thiemann et al. 2006; Robbins et al. 2012). Brown and black bears also go through seasonal cycles of decreased food availability, though these cycles are more predictable. Across most of their distribution, brown bears consume a diet composed primarily of carbohydrates supplemented by terrestrial sources of meat, but on the coast, salmon make a major contribution to their diet (Rode and Robbins 2000; Mowat and Heard 2006). Food availability is limited during winter, and bears enter dens where they fast during hibernation (Hilderbrand et al. 2000). Recent work suggests that hibernating brown and black bears have similar energy expenditures to hibernating pregnant female polar bears; however, fasting polar bears that remain active have much higher energy demands, with metabolic rates similar to the basal metabolic rates of mammals in general (Robbins et al. 2012). The degree to which polar bears can derive sufficient energy during prolonged periods of fasting affects individual survival and whether females are able to provide proper nutrition for their cubs or if they may abort reproductive efforts altogether (Derocher and Stirling 1996; Robbins et al. 2012).

At the cellular level, respiration (particularly oxidative phosphorylation [OxPhos]) provides the primary source of energy, directly influencing metabolic performance. Therefore, energetic challenges, such as thermal regulation and fasting, may lead to strong selection on the function of the mitochondrial and nuclear genes involved in this pathway. Polar bears, which inhabit colder environments and may undergo long periods of metabolically inefficient fasting, may have higher energetic demands when compared with their sister species, the brown bear. We hypothesize that the mitochondrial and nuclear genomes of polar bears may exhibit evidence of molecular adaptation in genes involved in cellular respiration to tune energy and heat production in response to demands of life in the Arctic.

Evidence for molecular adaptation has been found in many different organisms that demonstrate unique energy demands. Positive selection on mitochondrial genes involved in cellular respiration has been found in carnivorous plants (Jobson et al. 2004), reptiles (Castoe et al. 2008), and fish (Garvin et al. 2011). In mammals, evidence of adaptation has been found for several species, such as bats, which require relatively large amounts of energy for flight (Shen et al. 2010), killer whales inhabiting Antarctic pack ice (Foote et al. 2011), and primates, which have energetically costly brain activity (Grossman et al. 2004). Adaptations for energy metabolism, including cytochrome c oxidase activity, have been identified in the nuclear genome of the yak (Bos grunniens), which lives in cold, high-altitude environments (Qiu et al. 2012). Additional mechanisms for molecular adaptation have been proposed in rodents and hibernating ground squirrels. Proliferation of mitochondria in brown adipose tissue (BAT) and increased mitochondrial activity may be an adaptation for maintaining and quickly increasing body temperature above lower limits during hibernation (Boyer and Barnes 1999; Lowell and Spiegelman 2000).

Molecular adaptations in many of the species noted earlier have been identified in the protein complexes associated with the OxPhos pathway, which produces the majority of cellular energy in the form of adenosine triphosphate (ATP). OxPhos takes place in the mitochondria where electrons are transferred successively from one protein complex to the next until they reach cytochrome c oxidase, where O2 is reduced to H2O. As the electrons proceed down the chain, protons are pumped into the intermembrane space, resulting in an electrochemical gradient. Finally, ATP synthase uses the energy of this gradient to synthesize ATP (Ludwig et al. 2001). The mammalian mitochondrial genome encodes 13 proteins that make up critical subunits of the OxPhos complexes. The remaining approximately 1,500 genes necessary for the regulation and production of energy are encoded in the nuclear genome (Wallace 2005).

OxPhos is a complex process, and its efficiency may be regulated at many different levels. Mutations in the DNA sequences encoding the proteins of the OxPhos complexes may directly impact their function, such as in bats and primates. The mitochondrially encoded subunits form critical catalytic reaction centers and therefore are particularly sensitive to mutations (Fernández-Vizarra et al. 2009). Such mutations may be compensated for in nuclear encoded genes to allow for correct function (Rand et al. 2004; Osada and Akashi 2012). Regulation may also occur through the action of nuclear-encoded transcription factors and other proteins involved in the transcription and translation of the protein subunits. Additionally, nuclear-encoded subunits are thought to play a role in assembly of the complexes (Ludwig et al. 2001; Fernández-Vizarra et al. 2009; Pierron et al. 2012). In addition to the presence and proper assembly of the OxPhos complexes, necessary substrates, such as oxygen and adenosine diphosphate, and the correct mitochondrial membrane potential are required for OxPhos to occur (Hüttemann et al. 2008). Additionally, the cytochrome c oxidase complex can be inhibited by allosteric binding of ATP (Arnold and Kadenbach 1997). It can also be inhibited by high levels of nitric oxide (NO), which compete with O2 for the binding site. This inhibition decreases O2 consumption, alters the electrochemical gradient across the inner-membrane, and decreases ATP production (Ghafourifar and Cadenas 2005). In contrast to this, moderate levels of NO have been shown to increase mitochondrial biogenesis and activity, as well as influence expression levels of uncoupling protein 1 (UCP1) (Nisoli et al. 2003), which uncouples the production of energy during OxPhos to produce heat instead of ATP, a process referred to as nonshivering thermogenesis (Lowell and Spiegelman 2000). Up-regulation of UCP1 has been demonstrated to play a role in maintaining stable body temperature of small mammals occupying colder climates in China (Li et al. 2001), as well as hibernating arctic ground squirrels (Spermophilus parryii) from the northern portion of their range in Alaska (Yan et al. 2006). Given the complex process of cellular respiration and OxPhos, molecular adaptations in polar and brown bears may be manifested in several different ways.

Investigation of lineage-specific molecular adaptations in polar and brown bears is somewhat complicated by the fact that mitochondrial and nuclear genomes demonstrate different evolutionary histories in these species. Therefore, careful interpretation is required. In phylogenetic trees constructed from mitochondrial DNA, brown bears are paraphyletic, with brown bears from the Admiralty, Baranof, and Chichagof Islands of the Alexander Archipelago in southeastern Alaska (hereafter referred to as ABC brown bears) being more closely related to polar bears than to other brown bears (Talbot and Shields 1996), contrary to expectations based on their phenotype (fig. 1). Phylogenetic analyses incorporating the mitochondrial genome sequence from an ancient polar bear jawbone, dated to approximately 110,000–130,000 years before present, suggests that ABC brown bears and polar bears diverged approximately 150,000 years ago (Lindqvist et al. 2010). However, analysis of nuclear intron sequences (Hailer et al. 2012) and whole-genome sequences (Miller et al. 2012) of brown and polar bears indicate monophyly of both species, as expected. Genome sequences suggest a divergence time possibly as long as millions of years ago, followed by periods of hybridization and introgression (Miller et al. 2012). This hybridization is likely an underlying cause of the discordance between the mitochondrial and nuclear phylogenetic trees, with the time since divergence of polar bears and ABC brown bears being much shorter for mitochondrial DNA than nuclear DNA.

F

Discordance among mitochondrial and nuclear phylogenies for brown and polar bears. The dashed line is a diagrammatic representation of the mitochondrial gene tree, in which brown bears are paraphyletic. Thick black lines represent the nuclear gene tree, in which brown and polar bears form reciprocally monophyletic groups.

Polar bears are emblematic of adaptation and human-induced climate change, and here, we describe analysis of the mitochondrial and nuclear genomes of both polar and brown bears to test the hypothesis that polar bears have evolved lineage specific molecular adaptations in genes related to cellular respiration in response to environmental conditions in the high Arctic. To date, most studies of adaptive evolution in the cellular respiration pathway, particularly those in nonmodel, wildlife species, have focused solely on the mitochondrial genes involved in the process. However, 99% of genes involved are encoded in the nuclear genome (Wallace 2005), necessitating a genomic approach. Although sequencing costs have steadily declined in recent years (Wetterstrand 2013), genomic resources of the required magnitude are still unavailable for the majority of wildlife species. Further, most species are difficult to study in the wild and may be impossible to study through breeding, knock-out, or common garden experiments. Therefore, genome-wide analyses of positive selection and functional divergence provide a tractable and robust method for examining adaptations. We use extensive, recently developed genomic resources of brown and polar bears (Miller et al. 2012) to investigate genes that demonstrate signatures of functional divergence, including positive selection.

Materials and Methods

Mitochondrial Genome

Data Sets

We investigated all 13 genes of the mitochondrial genome, which are critical for energy production. We examined three data sets composed of varying numbers of taxa, representing different distributions of branch lengths. For Data Set 1, we examined sequences from 13 bear taxa and 13 carnivore outgroups (supplementary table S1, Supplementary Material online). Bear taxa included all extant species, additional members of distinct lineages within brown bears, plus sequences from an extinct polar bear lineage and both the extinct cave (U. spelaeus) and American giant short-faced (Arctodus simus) bears. For Data Set 2 we examined the subset of only the 13 bear taxa, and for Data Set 3 we investigated a further subset of bears belonging to the subfamily Ursinae, which includes all bears of the genus Ursus, including the sun bear (U. malayanus) and sloth bear (U. ursinus).

Phylogenetic Trees

Mitochondrial genome sequences were downloaded from GenBank (supplementary table S1, Supplementary Material online) and the control region removed. The remaining sequence (15,625 bp), consisting of the 13 protein-coding genes, 22 tRNAs, and two rRNAs, was aligned using the program MUSCLE (Edgar 2004) with default parameter settings (http://www.ebi.ac.uk/Tools/msa/muscle/, last accessed February 15, 2014). Maximum likelihood phylogenetic analyses were conducted using the RAxML BlackBox web-server (http://phylobench.vital-it.ch/raxml-bb/, last accessed February 15, 2014) (Stamatakis et al. 2008) assuming the general time reversible substitution model with rate heterogeneity modeled by a proportion of invariant sites and a proportion of sites with rates approximated by a gamma distribution having four discrete rate categories (GRT + I + G). The statistical support for each clade was assessed through 1,000 bootstrap replicates. The topology of the resulting tree (supplementary fig. S1, Supplementary Material online) was consistent with previously published results concerning bear phylogenetic relationships (Krause et al. 2008; Lindqvist et al. 2010), except we found a lack of support for the sister relationship between the American black (U. americanus) and the Asian black bear (U. thibetanus). All other branches received high support, and the phylogeny was well resolved.

Selection Analyses

To detect changes in evolutionary rate and signatures of positive selection, we analyzed the alignments of codon sequences and their corresponding tree under a maximum likelihood framework using the program CODEML in the PAML software package (Yang 2007). CODEML estimates the nonsynonymous/synonymous substitution ratio (dN/dS, ω), with ω < 1 indicating purifying selection, ω = 1 signifying neutral evolution, and ω > 1 suggestive of positive selection. In general, purifying selection acts across genes to preserve their function, which can mask evidence for positive selection at one or a few sites. To avoid this issue, we performed site-specific analyses, which take into account heterogeneity in selective pressures by defining different codon site classes and allowing different ω ratios for each site class (Nielsen and Yang 1998; Yang et al. 2000). We compared 1) model M3 (three classes of sites) with M0 (a single class of sites), 2) model M2a (positive selection) with M1a (nearly neutral), and 3) model M8 (rates at sites follow a beta distribution with a category for positive selection) with M7 (beta distribution only). Each model attaches a log-likelihood (ln L) value to each examined alignment and tree topology. To determine whether there was statistically significant support for each alternative model being tested over the null model, their ln L were compared using a log-likelihood ratio test (LRT, twice the difference between their ln L values). The comparison of M3 versus M0 is a test for variable rates of evolution, and not an explicit test for positive selection like the comparisons of M2a versus M1a and M8 versus M7, but may provide additional evidence for positive selection if the third category of sites has ω > 1.

The site-specific models described earlier are conservative because they examine evidence for accelerated evolutionary rates and positive selection across all the lineages in the tree. Therefore, to detect differences in evolutionary rates and evidence of positive selection on specific sites along a prespecified lineage, we also performed branch-site tests (Yang and Nielsen 2002; Bielawski and Yang 2004; Zhang et al. 2005). These models split the branches of the tree into a foreground branch of primary interest, with the rest of the branches designated as background branches. For these analyses, we compared 1) clade model D to sites model M3, which allows a class of sites to be under divergent selective pressure between foreground branches and the rest of the tree, and 2) modified Model A, which includes an extra class of sites under positive selection with ω > 1 in the foreground branches, to the null model A, in which the last site class has ω = 1. Again the LRT was used to evaluate the models. For all explicit tests of positive selection, we estimated the Bayes Empirical Bayes (BEB) posterior probability for sites to be under positive selection (Yang et al. 2005).

Tests for Functional Divergence

Because the codon-based analyses in PAML do not take into account the magnitude of change in the physiochemical properties of the amino acid resulting from a nonsynonymous substitution, we further analyzed genes with significant evidence of increased evolutionary rate and positive selection using the program TreeSAAP (McClellan and McCracken 2001; Woolley et al. 2003). This program compares the distribution of observed changes inferred from the phylogenetic tree to the distribution expected at random under neutral conditions, with significance determined through a goodness-of-fit test under a χ2 distribution. We included the default 31 physiochemical properties in the analysis with a sliding window size of 20 codons. Substitutions were considered significant if they were assigned to the highest categories representing extreme changes in physiochemical properties (categories 6–8) and if they were significant at P < 0.01. Structural modeling was performed in Visual Molecular Dynamics (VMD; Humphrey et al. 1996) using the crystal structures available in the Protein Data Bank.

Nuclear Genome

We implemented three complementary approaches to target genes or regions containing single-nucleotide polymorphisms (SNPs) in polar and brown bears. Then, independently for each lineage, we tested for significant enrichment or depletion of gene ontology (GO) terms and KEGG pathways associated with cellular respiration, and compared and further analyzed the results to examine lineage-specific differences in these species. For our analyses, we started with approximately 12 million nuclear SNPs that were previously identified by Miller et al. (2012) from whole-genome alignments of 23 polar bears, 3 brown bears (including 2 ABC and 1 non-ABC brown bears), and 1 black bear. These data are publicly available in the Genome Diversity shared data library on the Galaxy platform (https://usegalaxy.org/, last accessed February 20, 2014). The dog genome (Canis familiaris assembly version 2.0) was used as a reference for aligning bear sequences. To avoid redundancy, coding SNPs were mapped to the transcript with the longest coding sequence of each gene (i.e., canonical transcripts) (Hubbard et al. 2009; Flicek et al. 2012). GO terms and KEGG pathways were transferred to bear genes from the annotation of the dog assembly in the Ensembl database (Flicek et al. 2012). Data were analyzed using the Genome Diversity tools available in the Galaxy platform (Bedoya-Reina et al. 2013).

Using a maximum parsimony approach, we classified each SNP as being derived in either brown or polar bears. In this approach, the nucleotide present in two of the three species (brown, polar, and black bears) was considered to be the ancestral state. By using the tools available in Galaxy, we implemented three different methods to select three subsets of SNPs (fig. 2). For the first method (Method 1), we selected all SNPs mapped to coding regions where at least one allele was derived in either the brown or polar bear lineages. For example, SNPs derived in the polar bear lineage were required to have at least one copy of an allele that differed from the alleles in both brown and black bears. For the second and third methods, we only included derived, fixed SNPs mapped to coding regions and then examined genes with SNPs predicted to lead to a functional change (Method 2) or genes potentially under selection (Method 3). For a SNP to be considered fixed, all individuals of the species were required to have two copies of the same allele. By using this methodology, we avoided the inclusion of possible introgressed genes or loci with incomplete lineage sorting and selected the markers possibly involved in the adaptation of each clade. To investigate the functional effect of each nonsynonymous mutation in polar and brown bear lineages (Method 2), functional changes were predicted with Polyphen2 (Adzhubei et al. 2010). In this program, amino acid changes are classified as benign, possibly damaging (i.e., possibly causing a functional change), or probably damaging. We analyzed the change between the amino acid predicted in the recent common ancestor of polar and brown bears, and the amino acid fixed in each lineage. We also used the same data set to investigate genes under positive selection (Method 3). We estimated the ratio of nonsynonymous to synonymous mutations, ω, for each gene using the Jukes–Cantor model (Jukes and Cantor 1969). This calculation was done on a genome-wide scale using a custom script that implemented a divide-and-conquer algorithm, and which is available upon request. First, we investigated genes with ω > 1 in the brown or polar bear lineages (regardless of the value of ω in the other lineage, Method 3A), and then we investigated genes specifically that were under positive selection (ω > 1) in one lineage but under purifying selection (ω < 1) in the other lineage (Method 3B).

F

Overview of methods used to investigate the nuclear genomes of brown and polar bears for adaptations related to cellular respiration. The allele present in two of the three bear species was assumed to be ancestral. For example, a derived allele in polar bear would differ from the alleles found in brown and black bears. Dark gray indicates analysis of function, and light gray indicates analysis for positive selection. Bb is brown bear and Pb is polar bear.

Enrichment Tests

Genes were assigned functions based on GO terms and KEGG pathways associated with their annotation from the dog genome. Genes may have more than a single function, and therefore, the number of GO terms is larger than the number of genes. We analyzed these genes using tools available in Galaxy to test whether they were enriched or depleted in GO terms and KEGG pathways related to cellular respiration (i.e., those containing “oxi,” “ATP,” “energ,” or “mito”) when compared with all other genes in the genome in all other GO categories or KEGG pathways. Statistical significance was assessed with two-tailed Fisher exact tests, as suggested by Rivals et al. (2007). To test for statistically significant differences in ω for genes grouped in GO terms and KEGG pathways related to cellular respiration, we conducted Mann–Whitney U tests. Workflows of analyses are available at https://usegalaxy.org/u/talesdemileto/h/polar-bear (last accessed February 22, 2014). Also, the galaxy compressed files are available in the web page http://www.bx.psu.edu/∼oscar/ (last accessed February 15, 2014).

Additional Selection Analyses

Additional branch-site PAML analyses were conducted for genes in significant GO terms and KEGG pathways associated with cellular respiration identified in Methods 2 and 3. Because we only used SNPs that were fixed in each species, as well as fixed between species, only a single sequence each was used to represent polar, brown, and black bears. Orthologous genes for panda (Ailuropoda melanoleuca) and dog were obtained from the Ensembl database (Flicek et al. 2012). We used the program MUSCLE (Edgar 2004) to obtain codon sequence alignments and made corrections manually.

Results

Selection and Functional Divergence Analyses in Bears

In site-specific analyses of each of the three mitochondrial data sets, there was strong evidence of heterogeneous evolutionary rates among codon sites (P < 0.00001). As expected, a large proportion of sites were highly conserved, and the majority of remaining sites either demonstrated a relaxation of purifying selection (0 < ω < 1) or evolved in neutral manner (ω ∼ 1) (data not shown). In branch-site analyses of Data Set 1 (the full data set), comparison of the evolutionary rates between the bear clade and the clade composed of all the other carnivores suggested that 8 of 13 genes have significantly accelerated rates of evolution in bears (P ≤ 0.01) (supplementary fig. S2, Supplementary Material online). However, tests for positive selection were not significant for those genes except for ND5, which was found to be under positive selection in both bears and other carnivores (P = 0.032), suggesting that it may be important for carnivores in general. Within the data set composed of just bears (Data Set 2), site-specific analyses, which averaged across all bear lineages, suggested evidence for positive selection on two genes: COXIII (P = 0.025) and ND4 (P = 0.00001). For COXIII, ω = 1.384 at 1.3% of sites: the BEB analysis suggested evidence for positive selection at two positions, one of which (position 155 in the carnivore alignment, Asp ←→ Ser) was also predicted by TreeSAAP to lead to functional changes in the protein and occurred near interaction sites with two nuclear-encoded subunits (Tsukihara et al. 1996) (supplementary table S2 and fig. S3, Supplementary Material online). For ND4, ω = 1.474 at 2.5% of sites: the BEB analysis suggested the potential for positive selection at eight codon positions, one of which (codon position 251, Asp → Asn) had a high posterior probability of 0.997, was also detected as functionally divergent in TreeSAAP, and occurred in the region of the structure thought to be involved in proton pumping (Efremov et al. 2010; Efremov and Sazanov 2011) (supplementary table S2 and fig. S4, Supplementary Material online). COXIII or ND4 were not found to be under positive selection in branch-site analyses of the brown, ABC brown, or polar bear lineages.

Selection and Functional Divergence Analyses in Polar and Brown Bears

In branch-site analyses of the full data set (including outgroups, Data Set 1), significant differences (P = 0.006) in evolutionary rates of the COXI gene were detected between polar bears and all other lineages, with approximately 0.2% of the sites having ω ≈ 117 for polar bears, whereas ω = 0.52 for the other taxa (table 1). These results seem to be primarily driven by sequence divergence of the modern polar bear, as the test for heterogeneous rates of evolution for the branch to the ancient polar bear was not significant (P = 1.000). Additionally, the test comparing the clade containing the polar bears plus the ABC brown bear to all other taxa was marginally significant (P = 0.0421), whereas tests including non-ABC brown bears and other bears were not significant (P > 0.254). Explicit tests for positive selection were also not significant. Branch-site analyses on the data set of bears only (excluding outgroups, Data Set 2) recovered similar results as seen with the full data set (table 1), with two exceptions. In contrast to the analysis that included the outgroups, with the reduced data set the analysis comparing the clade containing the two polar bears plus the ABC brown bear to the clade containing all other bears did not detect any significant differences in evolutionary rates (P = 0.165). Also, in these analyses, significant differences in evolutionary rate of COXI were found between the group that included black, cave, brown, ABC brown, and polar bears and the group that included all other bears (P = 0.012), although evolutionary rates suggest strong purifying selection in the other bears and a relaxation of purifying selection in the black/brown/polar bear group (ω = 0.0669 vs. ω = 0.505), as opposed to the high ω values (e.g., ω ≈ 117) found along the branch leading to the polar bear (table 1). A similar pattern of increased evolutionary rates in polar bears was found in analyses that included an even more restricted taxon sampling of only bears in the subfamily Ursinae (Data Set 3). These restricted analyses provide support that polar bears have an accelerated rate of evolution in the COXI gene, even when compared with closely related bear species.

Table 1

Comparison of Evolutionary Rates for the Mitochondrial COXI Gene

Bear Taxa in Ingroup	Bears + Outgroups				Bears Only
	Proportion of Sites	ω (Ingroup)	ω (Others)	P	Proportion of Sites	ω (Ingroup)	ω (Others)	P
Ancient polar	0.2%	0.979	0.669	1.000	6.3%	1.227	0.172	0.763
Modern polar	0.2%	117.800	0.521	0.006	0.8%	28.692	0.438	0.026
All polar	0.2%	117.546	0.521	0.006	0.8%	28.639	0.437	0.026
ABC/polar	0.2%	14.023	0.525	0.042	0.8%	4.087	0.430	0.165
Brown/ABC/polar	0.2%	2.669	0.556	0.254	1.5%	1.530	0.285	0.114
Cave/brown/ABC/polar	0.2%	1.152	0.612	0.624	0.5%	1.349	0.518	0.457
Black/cave/brown/ABC/polar	6.5%	0.203	0.061	1.000	2.5%	0.505	0.067	0.012

Note.—Data sets were composed of 13 bear and 13 other carnivore taxa (left) and 13 bear taxa only (right). Bolded P values indicate comparisons for which the rates of the ingroup bear taxa differ significantly from the rest of the taxa.

BEB analysis of the COXI gene, which codes for a critical subunit of cytochrome c oxidase, suggests two potential sites under positive selection: codon site 57 (Val → Ile) with a posterior probability of 0.754 and codon site 512 (Met → Thr) with a posterior probability of 0.598. Functional divergence analysis of this gene suggests that substitutions at both of these sites may lead to changes in amino acid properties (P < 0.001, supplementary table S2 and fig. S5, Supplementary Material online). At position 57, this change was found to have occurred along four branches in the tree: the branch to the extinct short-faced bear (A. simus), to the sun bear (U. malayanus), to the American black bear, and to the non-ABC brown bear from Kodiak Island, Alaska, although these lineages did not demonstrate significantly increased evolutionary rates in branch-site tests. Structural analysis demonstrates that the amino acid at position 57 is located in the second alpha helix region and occurs close to several predicted heme-binding sites (positions 54, 58, and 61–62; fig. 3A). At position 512, the potential change in function was found to have occurred along the branch leading to the modern polar bear but not to the ancient polar bear. The amino acid in the modern polar bear sequence was found in 25 other modern polar bears from across their range, including the Chukchi and southern Beaufort seas near Alaska, and from Svalbard, Norway, but not in the ancient polar bear, which was discovered in Svalbard. Structural analysis of the bovine COXI protein suggests that amino acid 512 occurs in a portion of the protein relatively near the subunit I/subunit Vb interface (position 495) (fig. 3B).

F

Structure of the bovine cytochrome c oxidase complex. (A) Amino acid position 57 of COXI is highlighted in orange. This position occurs near predicted heme-binding sites at positions 54, 58, and 61–62. (B) Amino acid position 512 of COXI protein is highlighted in green and surrounded by beta-pleated sheets of the nuclear-encoded subunit Vb protein.

Derived Substitutions in Brown and Polar Bears

For the first analysis (Method 1), we selected all SNPs mapped to coding regions where at least one allele was derived in either the brown or polar bear lineages. About 7 million SNPs in brown and/or polar bears had a fixed variant in black bears. A total of 1,973,677 SNPs were found to have at least one allele originating in the polar bear lineage, and of these, 7,429 were nonsynonymous substitutions. These nonsynonymous SNPs occurred in 4,729 genes that were associated with a total of 5,638 GO terms. We also found that 3,573,736 SNPs had at least one allele originating in brown bears. Of these, 7,162 represented nonsynonymous substitutions, which occurred in 4,405 genes that were associated with 5,666 GO terms.

Of the 4,729 genes with nonsynonymous SNPs originating in the polar bear lineage, a total of 2,057 genes were significantly enriched or depleted for 257 GO terms (supplementary table S3, Supplementary Material online). Of the 4,405 genes with nonsynonymous SNPs originating in brown bears, a total of 1,723 genes were significantly enriched or depleted for 271 GO terms. Over both analyses, a total of eight of the significantly enriched/depleted terms were related to cellular respiration (table 2). Four of these were enriched in both brown bear and polar bears, two terms related to ATP binding and ATPase activity were enriched solely in brown bears, and the term “ATPase activity coupled to transmembrane movement of substances” was significantly enriched for derived nonsynonymous substitutions exclusively in polar bears. Among the genes found annotated with this GO term was CFTR encoding for the cystic fibrosis transmembrane conductance regulator, which mediates the pH concentration of the mitochondria and other organelles (Chandy et al. 2001). One GO term, “ATP synthesis couple proton transport” was significantly depleted in polar bears. Fifty-six KEGG pathways were significantly enriched or depleted in polar bears, whereas 53 were significantly differentially distributed in brown bears. Among these, the OxPhos pathway, which is strongly associated with cellular respiration, was significantly depleted for nonsynonymous, derived SNPs in both brown and polar bear lineages (table 2). No pathway associated with cellular respiration was significantly enriched.

Table 2

Significant GO Terms and KEGG Pathways Associated with Cellular Respiration for Genes that Have At Least One Derived Allele (Method 1)

GO Term/KEGG Pathway	Enriched/Depleted	Total Genes	Polar Bears	Brown Bears
ATPase activity	+	68	32**	26*
ATPase activity, coupled to transmembrane movement of substances	+	16	10**	6
ATP binding	−	575	183	180**
ATP catabolic process	+	66	32***	27*
ATP synthesis coupled proton transport	−	16	0**	2
Mitochondrial nucleoid	+	24	14**	12*
OxPhosa	−	114	14**	11***
Phospholipid-translocating ATPase activity	+	6	3	4*
Transcription from mitochondrial promoter	+	3	3*	3*

aKEGG pathway.

*P < 0.05.

**P < 0.01.

***P < 0.001.

Functional Effects of Fixed, Derived Nonsynonymous Substitutions

For the second analysis (Method 2), we only considered derived, fixed SNPs mapping to coding regions that were predicted to lead to a functional change. About 1 million SNPs with an allele fixed in black bears were found to have different alleles fixed in polar and brown bears; 7,635 of these were mapped to the coding regions of 4,156 genes, and 2,224 SNPs were found to cause nonsynonymous amino acid changes. In the polar bear lineage, a total of 409 nonsynonymous, fixed, derived substitutions were predicted to have a functional effect and these occurred in 381 genes. A total of 200 GO terms and 9 KEGG pathways were significantly enriched or depleted for these genes, of which 5 GO terms were associated with cellular respiration (table 3, supplementary table S4, Supplementary Material online). In comparison, 170 substitutions in 159 genes were fixed and predicted to have a functional effect in the brown bear lineage. A total of 136 GO terms and 6 KEGG pathways were significantly enriched or depleted, and of these, 7 GO terms were associated with cellular respiration (table 3). No KEGG pathways associated with cellular respiration were enriched or depleted for damaging substitutions in either of the two lineages.

Table 3

Significant GO Terms Associated with Cellular Respiration for Genes with Fixed, Derived, Nonsynonymous Substitutions Predicted to Lead to Functional Change (Method 2)

GO Term	Enriched/ Depleted	Total Genes	Polar Bears	Brown Bears
Extrinsic to mitochondrial inner membrane	+	2	0	1*
Heme-transporting ATPase activity	+	1	1*	0
Mitochondrial translation	+	3	1	1*
NO metabolic process	+	1	1*	1*
NO production involved in inflammatory response	+	1	1*	0
Nitric-oxide synthase activity	+	2	1*	0
Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen	+	39	3*	1
Positive regulation of glucose import in response to insulin stimulus	+	4	0	1*
Positive regulation of glucose metabolic process	+	5	0	1*
Positive regulation of fatty acid beta-oxidation	+	4	0	1*
Positive regulation of nitric-oxide synthase biosynthetic process	+	4	1	1*

*P < 0.05.

Brown bears demonstrated significant enrichment for the GO terms “positive regulation of glucose import in response to insulin stimulus” and “positive regulation of glucose metabolic process,” which are both related to the uptake of glucose into cells. Genes with substitutions causing functional impacts originating in the polar bear lineage were found to be overrepresented in GO terms involved in the regulation of NO. TLR4 (encoding the toll-like receptor 4 precursor), annotated with the term “nitric oxide production involved in inflammatory response,” was found to have a significantly higher rate of evolution in polar bears than brown bears (P = 0.04, polar bear = 999.0, brown bear = 3.304; supplementary table S7, Supplementary Material online). This gene was not found to be under positive selection in either lineage. The gene NOS3 was identified in the term “nitric oxide synthase activity,” which encodes for the endothelial NO synthase (eNOS). Additionally, the gene CPS1, associated with the enriched GO term “nitric oxide metabolic process,” was found to have damaging mutations in both the brown and polar bear lineages. Although still under debate, this gene may encode the mitochondrial NO synthase (mtNOS) in liver cells, or, through its control of the urea cycle, it may indirectly influence cellular NO levels (Scaglia et al. 2004; Ghafourifar and Cadenas 2005). Neither of these genes exhibited evidence for significantly increased evolutionary rate or positive selection when subjected to additional selection analyses using PAML.

Selective Pressures on Fixed, Derived Nonsynonymous Substitutions

For the final analysis (Method 3), we again only considered derived, fixed SNPs mapped to coding regions and then examined genes potentially under selection as assessed through ω. We first investigated if there was significant enrichment of any GO terms and KEGG pathways for genes potentially under positive selection (ω > 1) in each lineage, regardless of the value of ω in the other lineage (Method 3A). A total of 721 genes in polar bears and 427 genes in brown were found to have a ω greater than 1. For polar bears, a total of 240 GO terms were significantly enriched in these genes, 10 of which were associated with cellular respiration (table 4, supplementary table S5, Supplementary Material online). For brown bears, 274 significant GO terms were identified, 13 of which were associated with cellular respiration. Polar bears were significantly enriched for genes in two terms: “Mitochondrial alpha-ketoglutarate dehydrogenase complex,” which is related to the citric acid cycle, and “NADH oxidation,” which functions in glycolysis and gluconeogenesis and includes the gene ALDOB. ALDOB was found to have a higher evolutionary rate in polar bears, but this was not significant (supplementary table S7, Supplementary Material online). Eighteen KEGG pathways were enriched for genes in polar bears, as were 15 in brown bears, but none was associated with cellular respiration. One KEGG pathway, the OxPhos pathway, was significantly depleted for these genes in polar bears, but no evidence of significant depletion was detected in brown bears (table 4).

Table 4

Significant GO Terms and KEGG Pathways Associated with Cellular Respiration in Genes under Positive Selection (ω > 1, Method 3A)

GO Term/KEGG Pathway	Enriched/ Depleted	Total Genes	Polar Bears	Brown Bears
Extrinsic to mitochondrial inner membrane	+	2	0	1*
Heme-transporting ATPase activity	+	1	1*	0
Intrinsic to mitochondrial inner membrane	+	1	1*	0
Mitochondrial alpha-ketoglutarate dehydrogenase complex	+	1	1*	0
Mitochondrial degradosome	+	2	0	1*
Mitochondrial mRNA catabolic process	+	2	0	1*
Mitochondrial mRNA polyadenylation	+	1	0	1*
Mitochondrial outer membrane	+	28	4*	0
Mitochondrial RNA 3′-end processing	+	2	0	1*
Mitochondrial RNA 5′-end processing	+	1	0	1*
Mitochondrial RNA catabolic process	+	1	0	1*
NADH oxidation	+	1	1*	0
NO metabolic process	+	1	1*	1*
NO production involved in inflammatory response	+	1	1*	0
OxPhosa	−	114	0*	2
Oxidoreductase activity, acting on single donors with incorporation of molecular oxygen, incorporation of two atoms of oxygen	+	7	2*	0
Positive regulation of mitochondrial RNA catabolic process	+	2	0	1*
Regulation of cellular respiration	+	4	0	2**
RNA import into mitochondrion	+	1	0	1*
rRNA import into mitochondrion	+	2	0	1*
tRNA aminoacylation for mitochondrial protein translation	+	1	0	1*
Vacuolar proton-transporting V-type ATPase complex assembly	+	1	1*	0

aKEGG pathway.

*P < 0.05.

**P < 0.01.

Next, we investigated genes potentially under positive selection (ω > 1) in one lineage but under purifying selection (ω < 1) in the other lineage (Method 3B). We identified 114 GO terms with significant differences in ω between these species. For 44 of these terms, the average ω ratio was greater than 1 in polar bears and lower than 1 in brown bears (supplementary table S6, Supplementary Material online). Three terms associated with cellular respiration were enriched (table 5), including “negative regulation of nitric-oxide synthase activity.” Additionally, the term “gluconeogenesis” was also enriched. A total of five KEGG pathways showed a higher mean ω in polar than in brown bears, but none of them was associated with cellular respiration. A total of 70 GO terms and 5 KEGG pathways showed a mean ω higher than 1 in brown bears and lower than 1 in polar bears, but none of them were associated with cellular respiration.

Table 5

GO Terms Associated with Cellular Respiration in Genes with a Significantly Higher ω in Polar than in Brown Bears (Method 3B)

GO Term	Number of Genes with ω Polar > 1 > Brown
Gluconeogenesis	5*
Negative regulation of NO synthase activity	3*
Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen	11*
Oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, reduced flavin or flavoprotein as one donor, and incorporation of one atom of oxygen	4*

Note.—These genes may be under positive selection in polar bears and under purifying selection in brown bears. No genes with higher ω in brown bears were found to be significantly enriched or depleted in functions related to cellular respiration.

*P < 0.05.

We observed that for genes potentially under positive selection in the polar bear lineage, there was an enrichment of GO terms associated with NO regulation. For instance, we found the GO terms “nitric oxide production involved in inflammatory response” (Method 3A, table 4) and “negative regulation of nitric-oxide synthase activity” (Method 3B, table 5) to be enriched. Genes with GO terms associated with NO regulation that have a higher ω ratio in polar than in brown bears include CAV3 (encoding caveolin-3), ENG (encoding endoglin), and TLR4.

Comparison of the Genes Identified in the Polar and Brown Bear Lineages

Overall, the genes found to be associated to the overrepresented GO terms and KEGG pathways, as identified by the three methods implemented here, were largely unique to each lineage, although there was some overlap between species (supplementary table S8, Supplementary Material online). A total of 179 genes were uniquely identified in brown bears, 44 uniquely in polar bears, and 32 were identified in both lineages. Shared genes represented 15% of the total for brown bears and 42% for polar bears. Although there was overlap for the broader functions identified by each method (e.g., NO regulation), the genes identified within those functions differed. Method 1 (shared, derived alleles) uniquely identified 176 genes in brown bear, 30 in polar bear, and 32 in both. Method 2 (lineage specific SNPs predicted to result in functional change) uniquely identified three genes in brown bears, 6 in polar bears, and 1 in both. Finally, Method 3 (lineage specific SNPs demonstrating evidence for positive selection) identified 3 genes in brown bears, 10 in polar bears, and 1 in both. This also demonstrates that Method 1 is the broadest and least restrictive method to target SNPs associated with cellular respiration in the two bear species.

Some genes for each lineage were identified through multiple approaches. Three genes were found in brown bears by two different methods: EARS2 (an aminoacyl-tRNA synthetase) by Methods 1 and 3, KDR (which forms a receptor for vascular endothelial growth factor) by Methods 1 and 2, and NOA1 (NO-associated protein 1) by Methods 2 and 3. The TLR4 and ABCB6 (mitochondrial ABC transporter) genes were found in polar bears by Methods 2 and 3, whereas CPS1 was found in both bear lineages by all the three methods used.

Further, PAML tests on the genes belonging to significantly enriched GO terms and KEGG pathways failed to demonstrate significant evidence for positive selection along the polar or brown bear lineages. However, seven genes displayed significantly increased rates of evolution in the polar bear lineage (some discussed earlier): TLR4, ALOX5, CYP7A1, OGFOD1, PCK1, PPOX, and MSTO1 (supplementary table S7, Supplementary Material online).

Discussion

Effective Population Size and Molecular Evolution in Bears

We found evidence in the mitochondrial and/or nuclear genomes for accelerated evolutionary rates, positive selection, and functional divergence in bears, with particularly strong signals in the polar bear lineage. Interpretation of accelerated evolutionary rates can be challenging because several factors can explain the results, including positive selection, mutational biases, and genetic drift (Yang and Bielawski 2000). Carnivores, such as bears, are often top predators and may have low effective population sizes (Ne). When effective population sizes are low, genetic drift dominates and selection can be relatively ineffective, allowing slightly deleterious mutations to accumulate (Allendorf and Luikart 2007), potentially leading to increased evolutionary rates. This may be particularly relevant for mitochondrial genomes, which are haploid and demonstrate maternal inheritance, and therefore have a lower effective population size than nuclear genes. Here, we identified several genes that exhibit increased evolutionary rates in bears in general, and in brown and polar bears specifically. The higher evolutionary rate in bears could potentially be associated with a lower Ne in these species when compared with other carnivores.

Comparisons of Ne among species and between studies is inherently difficult because estimates of effective population size are often made using a wide variety of techniques that may be relevant at different time scales and that depend on different assumptions (Luikart et al. 2010). Worldwide Ne estimates of extant bears from the literature are approximately 4,000 for polar bears (Cronin et al. 2009), 20,000 for brown bears (Paetkau et al. 1998), and 171,000 for American black bears (Hellgren and Vaughan 1989), based on a census to effective population size ratio of about 0.18, which was determined independently for all three species. Recent whole-genome sequencing of giant pandas (Ai. melanoleuca) estimated a present effective population size of approximately 5,000 (Zhao et al. 2013). Effective population size estimates in other carnivores range broadly. An estimate for tigers (Panthera tigris) is approximately 1,300 individuals (Smith and McDougal 1991), whereas for leopards (P. pardus) in Tanzania, the estimate may be closer to 40,000 (Spong et al. 2000). An estimate for the European mink (Mustela lutreola) is about 2,200 (Lodé and Peltier 2005), but estimates for Antarctic seals are orders of magnitude higher at 151,000 for Weddell seal (Leptonychotes weddellii), 255,000 for Ross seal (Ommatophoca rossii), and 880,000 for crabeater seal (Lobodon carcinophaga) (Curtis et al. 2011). It should also be noted that it has been shown that effective population sizes may fluctuate considerably over the course of millions of years (Miller et al. 2012; Zhao et al. 2013). Therefore, it is difficult to say whether or not the overall Ne of bears is significantly lower than for other carnivores and whether increased evolutionary rates in bears are primarily due to genetic drift. Polar bears have relatively small effective size, and population sizes may have fluctuated in the past during periods of climatic change (Miller et al. 2012), suggesting a potential role for drift, but they also have likely experienced particularly strong selective pressures associated with the energetic demands of living in the extreme Arctic environment.

Positive Selection and Functional Divergence in Bears

Branch-site analyses comparing the evolutionary rates of mitochondrial protein-coding genes in bears and other carnivores suggest that bears have higher evolutionary rates in 8 of 13 genes. Beyond the influence of genetic drift, differences in the branch lengths between bears and outgroups may also influence analyses of evolutionary rates, as saturation of substitutions may have occurred on the longer branches to the outgroups. Simulations suggest that as divergence increases, saturation does occur and power of the branch-site test decreases (Gharib and Robinson-Rechavi 2013). However, no corresponding increase in false positives has been found (Gharib and Robinson-Rechavi 2013). Selection pressures are also important to consider. In general, bears may face different energetic challenges than most other carnivores, because the diets of several species (particularly in the largest genus, Ursus) rely on carbohydrates from vegetation, which may lack protein and other vital nutrients (Hilderbrand et al. 1999; Hilderbrand et al. 2000; Rode and Robbins 2000). Therefore, it is possible that natural selection and adaptation may have played a role in their responses to meet these challenges.

Results here suggest that the mitochondrial COXIII and ND4 genes exhibit some evidence for both positive selection and functional divergence in bears (but not other carnivores), which may be suggestive of molecular adaptation. COXIII is part of the cytochrome c oxidase complex, which is an important regulator of OxPhos and subsequent production of ATP. The amino acid change in position 155 in bears could potentially influence the ability of this subunit to interact with two nuclear-encoded subunits, which in general are thought to regulate conformation and activity of the complex. Although subunits I and II (coded for by COXI and COXII, respectively) are thought to be the primary units involved in electron transfer, subunit III may play a role in proton pumping (Fernández-Vizarra et al. 2009). ND4 is part of the first respiratory complex, NADH-quinone reductase, which is also involved in electron transport, and it has been suggested to form part of the proton-pumping mechanism (Efremov et al. 2010; Efremov and Sazanov 2011). Changes in this gene may lead to changes in the ability of this complex to pump protons across the mitochondrial membrane, although further work is necessary to gain a better understanding of the structure of this complex.

Evolution of Cytochrome c Oxidase I in Brown and Polar Bears

Despite a lack of explicit evidence for positive selection, we identified several potentially interesting evolutionary signatures in the COXI gene of brown and polar bears. As mentioned earlier, COXI is one of two core subunits of the cytochrome c oxidase complex and is critical to ATP production. Codon position 57 of COXI had a low BEB posterior probability of 0.75, which should be interpreted cautiously because the explicit tests for differential rates of molecular evolution and positive selection were not significant. However, positive selection can be difficult to detect, even for relatively more powerful branch-site tests and BEB analyses (Zhang et al. 2005). Structural analyses suggest that this site is near important binding regions of the cytochrome c oxidase complex (fig. 3A), and functional analyses in TreeSAAP suggest that the substitution observed here may lead to functional changes. However, in brown bears, the substitution only occurred in an Alaskan brown bear lineage. Additional analyses, particularly of a larger set of brown bears, are necessary to determine whether this substitution may be under selection in a portion of their range, which could be indicative of local adaptation in cellular respiration.

A potentially interesting evolutionary signal was also found at codon position 512 of the COXI gene in the polar bear lineage. This site also had a low BEB posterior probability (of ∼0.60), so results should be interpreted cautiously. Computational predictions did suggest that the substitution at this site could lead to a significant change in the physiochemical properties of the amino acid (supplementary table S2, Supplementary Material online) and potentially in the function of the enzyme. This substitution was found to be fixed in a larger set of 25 polar bears from across their range, including bears captured in Alaska and Svalbard, Norway. This may not be surprising as the low Ne of polar bears could lead to relatively quick fixation of an allele when compared with bear species with higher Ne. Further, the low Ne of polar bears would suggest that drift would dominate evolutionary changes, as opposed to selection, unless the changes were strongly beneficial or deleterious. The consequence of the amino acid change in position 512 is not immediately clear because it is not directly involved in polypeptide or ligand interactions (fig. 3B), and because this position is relatively variable among carnivores. It is relatively near the site of interaction with a nuclear-encoded subunit, and therefore, it potentially could affect regulation of the complex. A similar substitution pattern has been noted in other mammal species as well (da Fonseca et al. 2008).

Examination of the molecular evolution of mitochondrial genes in polar bears is also somewhat complicated by their evolutionary history. Nuclear and mitochondrial gene trees demonstrate incongruent topologies (fig. 1), the exact cause of which remains a contentious topic. In particular, the origin of the mitochondrial genomes of polar and ABC brown bears is still unclear. Analyses have suggested that the mitochondrial genomes of ABC brown and polar bears diverged approximately 150,000 years ago (Lindqvist et al. 2010), and if brown bear mtDNA was transferred into polar bears (Edwards et al. 2011) around this time, then there would have been relatively little time for independent evolution to occur in polar bears. However, if ABC brown bears harbor a mitochondrial genome that originated in the polar bear lineage (Cahill et al. 2013), we expect that ABC brown bear mitochondrial DNA would demonstrate similar patterns of molecular evolution to polar bears. There is little or no evidence for higher evolutionary rates of the COXI gene on the branch leading to either the ABC brown bear or the extinct ancient polar bear, which existed 110,000–130,000 years ago (Alexanderson et al. 2013), shortly after the mitochondrial divergence date, but there is some evidence for increased evolutionary rate in modern polar bears suggesting that an increase in rate may be relatively recent. The short length of the branch to the modern polar bear (supplementary fig. S1, Supplementary Material online) may hinder accurate detection of positive selection because short branches contain little information, subsequently leading to low power (Zhang et al. 2005; Nozawa et al. 2009; Gharib and Robinson-Rechavi 2013). Investigation of three data sets composed of taxa with varying branch lengths were consistent, though, providing some measure of the robustness of the results. Overall, given the evolutionary history of the mitochondrial genome and the effective population size of polar bears, it may be difficult to gain a full understanding of the role of positive selection in shaping the mitochondrial genome in this lineage.

Previous work has suggested that on a global scale, rates of evolution in the genome of bears may be similar to those of other mammals (Zhao et al. 2010). However, evidence for potentially adaptive molecular evolution has been uncovered as well. In the giant panda, a large group of genes that are potentially under positive selection demonstrate functions related to the sensory system, including taste and olfaction. These genes may have a beneficial effect on diet (e.g., through increased bamboo consumption) and mate location (Zhao et al. 2013). Expressed sequence tag approaches have identified genes that appear to evolve more quickly in the American black bear. Evidence for positive selection was found for one of these, the CSRP3 gene, which also had a substitution that may impact protein function (Zhao et al. 2010). This gene helps regulate the force of heart contraction, and it was suggested that this could be an adaptation for lower heart rate, and therefore energy preservation, during hibernation (Zhao et al. 2010). Here, we examine whether polar bears may demonstrate lineage-specific adaptations in genes related to cellular respiration due to their unique energy requirements in the high Arctic environment.

We have extensively relied on GO term and KEGG pathway analysis to identify the functional implications of our results, which may lead to overstatement of the biological significance of our findings. Nevertheless, the approach we implemented is a very efficient way to extract major biological and evolutionary implications from large gene data sets, and it has been proven to be useful in multiple studies (Fedorov et al. 2011; Qiu et al. 2012; Zhao et al. 2013). Additionally, the robustness of our results is supported by consistently finding related categories and pathways across the results obtained by different methods.

Nuclear Genes Related to Cellular Respiration in Brown and Polar Bears

The overwhelming majority of proteins required for regulation and production of energy during the process of cellular respiration are encoded in the nuclear genome. Our investigation of the nuclear genomes of brown and polar bears revealed that genes with substitutions that were derived in either the brown bear or polar bear lineages were significantly enriched for functions related to mitochondrial activity in general, and cellular respiration in particular.

Some functional categories were shared between brown and polar bears. In the analyses utilizing the presence of at least one derived allele (not necessarily fixed in either brown or polar bears, Method 1), we identified the shared terms “Mitochondrial nucleoid” and “Transcription from a mitochondrial promoter” (table 2). Genes in these terms include POLG, the mitochondrial DNA polymerase (Spelbrink et al. 2000), DNA2, which is a helicase/nuclease that acts to stabilize mitochondrial and nuclear DNA (Duxin et al. 2009), C10orf2, which is a helicase important for mtDNA replication (Longley et al. 2010), and TFAM, which codes for an important mitochondrial transcription factor and is also involved in both mtDNA replication and repair (Shi et al. 2012). The gene PPARGC1B, an important transcription factor and regulator of mitochondrial activation and OxPhos (Handschin and Spiegelman 2006), was also annotated with these terms. Because substitutions in these genes have not become fixed yet between the brown and polar bear lineages, they may have been either recently derived or perhaps are important for dealing with energetic challenges in both lineages.

From the analysis of derived alleles (Method 1), we also found that for both brown and polar bears, the KEGG pathway “Oxidative Phosphorylation” had fewer genes with derived substitutions than would be expected by chance (table 2). This would suggest that, in general, genes involved directly in OxPhos are highly conserved in brown and polar bears, consistent with the hypothesis that energy production is particularly important in these lineages. In analyses of genes that may be under positive selection (ω > 1, Method 3A), we found, again, that the KEGG pathway for OxPhos was depleted, but this time only in polar bears (table 4). This suggests that the proteins involved in OxPhos are strongly conserved in polar bears. In American black bear, it was noted that genes potentially involved in physiological processes tightly linked to hibernation are often more highly conserved than other genes (Zhao et al. 2010). This suggests a pattern of particularly strong constraints on genes with important physiological functions related to energy production and conservation.

In addition to these, polar bears also demonstrated significant enrichment of genes under positive selection in two other terms related to cellular respiration: “mitochondrial alpha-ketoglutarate dehydrogenase complex” and “NADH oxidation.” The gene BCKDHB was identified under the term “mitochondrial alpha-ketoglutarate dehydrogenase complex.” This complex is known to be an important control point in the citric acid cycle. The gene ALDOB was identified under the term “NADH oxidation” and is involved in both glycolysis and gluconeogenesis. Although not involved directly in OxPhos, substitutions in these genes may still influence the overall process of cellular respiration in polar bears.

Glucose Uptake and Metabolism in Brown Bears

In brown bears, two GO terms related to glucose were found to be significantly enriched for genes with substitutions predicted to cause functional change (Method 2): “Positive regulation of glucose import in response to insulin stimulus” and “Positive regulation of glucose metabolic process” (table 3). The same gene, IRS1 (Insulin receptor substrate 1), was associated with both GO terms. Insulin signals cells to take up glucose from the blood, and if necessary convert it to glycogen or triglycerides for energy storage. The diet of brown bears varies considerably across their distribution, but in most areas, it is largely dominated by carbohydrates and supplemented with meat when available, although salmon make a major contribution in coastal areas (Rode and Robbins 2000; Mowat and Heard 2006). A functional change in the IRS1 protein could aid brown bears in regulating the uptake and use of glucose from their diet (Tamemoto et al. 1994), which would be especially relevant during times of hyperphagia before hibernation. An increase in insulin has been noted in black bears during hyperphagia (Palumbo et al. 1983), and a functional change in IRS1 may allow brown bears to store more energy for later use during hibernation.

Alternatively, in recent studies a decrease in glucose catabolism has been suggested during hibernation in ground squirrels (Yan et al. 2008) and black bears (Fedorov et al. 2011), as evidenced by a reduction in transcription of genes involved in carbohydrate catabolism in hibernating versus active animals. Instead, transcriptional levels of genes involved in beta oxidation and lipid catabolism are increased during hibernation (Yan et al. 2008; Fedorov et al. 2011), as fat reserves act as the primary energy source during this time. Although changes in expression levels may allow genomic plasticity in immediate responses to the environment, changes at the nucleotide sequence level of the IRS1 gene, as documented here, could represent a longer term evolutionary adaptation. A functional change in IRS1 could help prevent glucose uptake (i.e., insulin resistance) or catabolism during brown bear hibernation and aid in the metabolic switch to lipid catabolism. The substitution in IRS1 identified here was fixed in the three brown bears for which nuclear genome sequences are available. Sequencing this gene in brown bears from across their range will yield additional insights into its importance in brown bear feeding and physiology.

Enrichment in the Polar Bear Lineage for Functions Related to NO

Investigations of genes with fixed, derived substitutions exhibiting functional divergence (Method 2) or evidence of positive selection (Method 3) demonstrated enrichment in more GO terms related to NO production and regulation in polar bears than in brown bears. In mammals, NO is known to have pleiotropic effects, and the role of NO is complex because it depends on the cell type, concentration, and the length of action. However, there is increasing evidence that NO is important in regulating energy metabolism in mammals (Dai et al. 2013). NO is synthesized from l-arginine by three isoforms of NOS that were initially named in relation to their tissue localization (Förstermann and Kleinert 1995): neuronal NOS (nNOS or NOS I), inducible NOS (iNOS or NOS II), and endothelial NOS (eNOS or NOS III). There is some evidence that an NOS isoform may also be present in mitochondria (mtNOS) (Ghafourifar and Cadenas 2005).

In the polar bear lineage, the genes associated with significant GO terms related to NO include NOS3, CPS1, TLR4, CAV3, and ENG. The NOS3 gene, which encodes eNOS (Janssens et al. 1992), demonstrated a fixed substitution, which computational predictions suggested has a significant impact on function (Method 2). Therefore, in polar bears, the levels of NO production by one of the three NO synthesis proteins may be directly altered.

The gene CPS1 was identified in both functional analyses and analyses of selection (Methods 2 and 3). There is some evidence that CPS1 may encode the mitochondrial NOS and therefore directly regulate NO concentrations, although this is contentious (reviewed in Ghafourifar and Cadenas 2005). On the other hand, because CPS1 acts during the first step of the urea cycle, which produces l-arginine, it indirectly regulates the availability of the substrate for NOS, and hence NO levels (Pearson et al. 2001; Scaglia et al. 2004). A previous study has shown down-regulation of the expression of CPS1 in hibernating Arctic ground squirrels, which may be an adaptation for using amino acids as a source for gluconeogenesis to fuel tissues with high energy demands, such as the brain and BAT (where adaptive thermogenesis takes place; see below) (Yan et al. 2008). Therefore, the identification of CPS1 may indicate its importance in either directly or indirectly controlling levels of intracellular NO or in energy production for nonshivering thermogenesis.

The genes TLR4, CAV3, and ENG are known to influence the activity of NOS enzymes. Stimulation of the TLR4 gene, which was identified during functional analyses (Method 2), has been shown to increase expression of iNOS in response to lipopolysaccharides/toxins and subsequently increases the concentrations of NO (Heo et al. 2008). CAV3 and ENG, both detected during analyses of selection (Method 3B), regulate the activity of eNOS: CAV3 modulates activity of eNOS and is expressed predominantly in striated muscle cells (Williams and Lisanti 2004), and ENG stabilizes eNOS and increases activity (Toporsian et al. 2005). Overall, these findings suggest lineage-specific evolution in polar bears in genes involved in the regulation and production of NO.

NO and Adaptation in Polar Bears

NO has several functions related to cellular respiration and the control of energy production. At high levels, NO interacts with superoxide producing peroxynitrite, which irreversibly inhibits cytochrome c oxidase and other respiratory chain complexes, inducing apoptosis (Ghafourifar and Cadenas 2005). At low concentrations, NO competes directly with oxygen at cytochrome c oxidase, and reversibly negatively regulates this complex (Brown 2001). Thus, NO can directly regulate cellular respiration, oxygen consumption, and energy production. This could potentially be beneficial to polar bears if they selectively decrease ATP production during times of rest or fasting, which would help them maintain energy stores. Also, by decreasing energy production while resting, polar bears may minimize oxygen consumption and therefore decrease heat lost through respiration in subfreezing temperatures.

Adaptive, or nonshivering, thermogenesis is a regulated production of heat in response to environmental temperature or feeding and can protect against cold exposure or regulate energy balance after changes in diet (Lowell and Spiegelman 2000). It has recently been shown that NO plays a role in this process. Environmental signals trigger the release of the hormone noradrenaline, stimulating β-adrenergic receptors, which activate eNOS, and subsequently lead to the production of NO (Bossy-Wetzel and Lipton 2003; Nisoli et al. 2003). As indicated earlier, high concentration of NO inhibits mitochondrial respiration. However, through activation of cyclic guanosine monophosphate, moderate levels of NO induce a master regulator of mitochondrial biosynthesis, leading to an increased expression of genes encoding components of the respiratory chain complexes and ultimately increasing mitochondrial biogenesis (Bossy-Wetzel and Lipton 2003; Nisoli et al. 2003). A primary molecule involved in cold-induced thermogenesis is the UCP1, the activation of which allows protons to flow back through the mitochondrial membrane, uncoupling ATP production, and resulting in the generation of heat (Lowell and Spiegelman 2000).

BAT is the primary site for adaptive thermogenesis. BAT is present in small mammals, such as rodents, and it has been demonstrated that an increase in UCP1 synthesis in BAT plays a key role in thermogenesis for hibernating small mammals occupying colder climates (Milner et al. 1989; Li et al. 2001; Yan et al. 2006). Increased mitochondrial activity may be an adaptation for quickly increasing and maintaining body temperature above lower limits during hibernation (Boyer and Barnes 1999; Lowell and Spiegelman 2000). In hibernating Arctic ground squirrels, expression of genes related to fatty acid catabolism and gluconeogenesis were increased in BAT, and it was suggested that these processes provide the energy to sustain adaptive thermogenesis or other high energy organs, like the brain (Yan et al. 2008). Glycerol, which arises from the breakdown of triglycerides, provides the substrate for gluconeogenesis, and thus can supply additional energy for increased thermogenesis (Yan et al. 2008). Increased expression of genes in the gluconeogenesis pathway, including PCK1, a main control point for the pathway, has also been found in hibernating black bears (Fedorov et al. 2011).

In larger mammals, BAT is present at birth but becomes sparse during development, and it has been suggested that BAT may be lacking in adult bears (Davis et al. 1990; Jones et al. 1999). However, thermogenesis has also been found to occur in white adipose tissue (Ye et al. 2013) and may potentially occur in muscle through more complex mechanisms (Wu et al. 1999; Lowell and Spiegelman 2000; Meyer et al. 2010). Interestingly, our analyses identified five genes related to gluconeogenesis that are potentially under positive selection in polar bears, including PCK1. If adaptive thermogenesis does occur in polar bears, and is fueled by gluconeogenesis similar to Arctic ground squirrels, it would not be surprising for genes in the gluconeogenesis pathway to be under selection, because the diet of polar bears is composed primarily of blubber and is essentially devoid of any substantial source of carbohydrates (Ramsay and Hobson 1991; Grahl-Nielsen et al. 2003). Recent work suggests that gluconeogenesis does increase in mice fed a high-fat diet (Obici et al. 2012). If polar bears utilize adaptive thermogenesis, it is possible that regulation of NO concentrations provides a mechanism to control energy balance at different times of the year and in response to individual needs. When polar bears have relatively more abundant food supplies, they could utilize heat production through thermogenesis to maintain their body temperature. During the summer, when temperatures are warmer and prey are less abundant, adaptive thermogenesis would not be necessary, thus allowing polar bears to minimize their energy expenditures and retain more of their fat reserves.

Considering the signals of positive selection and functional divergence observed in a large group of NOS genes and their effectors (i.e., NOS3 or eNOS, CPS1, and TLR4), it is plausible that substitutions in these genes may allow fine tuning of NO production and provide an adaptive response for polar bears in controlling the trade-offs between oxygen consumption, heat production, and energy production in the form of ATP (supplementary fig. S6, Supplementary Material online). Such adaptations would be beneficial during times of low energy demand, to make energy production more efficient for dealing with periods of low food availability, to prevent heat loss from respiration in cold environments, or to uncouple ATP production in favor of releasing heat for warmth. It is unknown whether this may represent a more generalized mechanism in high-Arctic species, but it may be an adaptive response in polar bear compared with its lower latitude relatives, the brown bear, and American black bear. This study represents one of the first genome-wide analyses for nonmodel, nonprimate organisms and may provide insights into the genetic basis of adaptation in other Arctic species, which face an uncertain future in a rapidly changing climate (Stroeve et al. 2007; Overland and Wang 2013).

Supplementary Material

Supplementary figures S1–S6 and tables S1–S8 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).