Translator Disclaimer
14 October 2011 Phylogeography of Steller sea lions: relationships among climate change, effective population size, and genetic diversity
Author Affiliations +

The biology of the Steller sea lion (Eumetopias jubatus) has been the subject of intense scientific investigation. This is primarily due to the rapid decline of population size in the western part of the species' range since the 1970s and the subsequent Threatened and Endangered species listings that had direct impact on the management of one of the world's largest fisheries. The Steller sea lion has emerged as an indicator species representing the environmental health of the North Pacific Ocean and Bering Sea. In this study, to better understand the historical processes that have culminated in the extant populations of E. jubatus, a large genetic data set consisting of 3 mitochondrial regions for >1,000 individuals was analyzed from multiple phylogeographic and demographic perspectives. The results describe the role of climate change in shaping the population structure of E. jubatus. Climatically associated historical processes apparently involved differential demographic responses to ice ages (and putative glacial vicariance) dependent on population size. Ice ages during times of small effective population size promoted restricted gene flow and fragmentation, and ice ages occurring during times of large population size promoted gene flow and dispersal. These results illustrate that effective population size has a profound effect on how species respond to climate change, an observation with obvious implications for large mammals and endangered species under the present conditions of imminent anthropogenically caused climate change. In addition, the results confirm previous observations of strongly biased historic and contemporary gene flow involving dispersal from west to east. Furthermore, phylogenetic patterns in combination with available fossil data suggest the potential of an Asian origin of E. jubatus. The results of this study provide a detailed scenario for the history that has shaped contemporary populations of E. jubatus.

An intensive research effort into the biology of the endangered Steller sea lion (Eumetopias jubatus) has resulted in significant advances in our understanding of the evolutionary history, population structure, demographic trends, physiology, ecology, and evolution of this species (Baker et al. 2005; Bickham et al. 1996, 1998; Harlin-Cognato et al. 2006; Hoffman et al. 2006, 2009; Kenyon and Rice 1961; Loughlin et al. 1987, 1992; O'Corry-Crowe et al. 2006; Phillips et al. 2009a, 2009b; Trites and Donnelly 2003). For example, population genetic and phylogenetic investigations have demonstrated clear population structure at rookery, region, and stock levels, with phylogenetic delineation increasing at broader geographic groupings (Baker et al. 2005; Bickham et al. 1996, 1998; O'Corry-Crowe et al. 2006). In addition, a population decline of >80% resulted in listing the western stock as endangered under the United States. However, the eastern stock has maintained stable numbers over recorded time but still retains a threatened listing. The determination of major coincident patterns in genetics, morphology, and demography recently has led to the formal description of these 2 forms as the subspecies E. j. jubatus (consisting of the Asian and western stocks) and E. j. monteriensis (eastern stock—Phillips et al. 2009a; Fig. 1). Refining our understanding of the processes that ultimately led to the current population structure can enhance our existing knowledge of mammalian evolution and guide the future management of Steller sea lions in the North Pacific and Bering Seas. As a result of extensive research conducted to support Endangered Species Act listing, we have a relatively broad understanding of many aspects of the biology and life-history characteristics of E. jubatus and a preliminary understanding of how the evolutionary history of the species has influenced observed contemporary genetic patterns. Given the distribution of E. jubatus around the North Pacific Rim (within a region where glacial movement has been particularly dynamic throughout the Pleistocene—Mann and Peteet 1994; Fig. 1) and the dependence of E. jubatus on linearly distributed rookeries for breeding and birthing (many of which are within areas known to have been covered by ice during previous glacial maxima—Grosswald and Hughes 2002; Mann and Peteet 1994; Westgate et al. 2001), the hypothesis of Harlin-Cognato et al. (2006) that genetic patterns are in part the result of phylogeographic occurrences driven by climatic change (changes in the extent of ice cover, glacial vicariance, and sea-level changes) is reasonable. Furthermore, the presence of a detailed and accurate phylogeographic history encompassing the entire range of E. jubatus is necessary to provide evolutionary context for the life-history characteristics documented in other studies.


Map of the distribution of Eumetopias jubatus and recognized subspecies. Abbreviations for the regions are as follows: OKH  =  Sea of Okhotsk, KUR  =  Kuril Islands, KAM  =  Kamchatka Peninsula, COM  =  Commander Islands, WAL  =  western Aleutian Islands, CAL  =  central Aleutian Islands, EAL  =  eastern Aleutian Islands, WGA  =  western Gulf of Alaska, BER  =  Bering Sea, CGA  =  central Gulf of Alaska, PWS  =  Prince William Sound, SEA  =  southeastern Alaska, BRC  =  British Columbia, ORE  =  Oregon, NCA  =  northern California. Rookeries included in this study are indicated by arrows and numbered following Appendix I.


The well-defined population model of E. jubatus is characterized by rookery dependence with high female philopatry. This information, in combination with a thorough genetic sampling of individuals through years of fieldwork, is useful for interpreting genetic patterns from a population genetics or phylogeographic perspective. In addition, a detailed working hypothesis of global climatic history based on ice core data (Loulergue et al. 2008) provides logical a priori phylogeographic expectations and a potential to build accurate phylogeographic inferences about E. jubatus.

We used a large genetic database of E. jubatus of >1,000 individuals (approximately 1.63% of the female population from 79% of known rookeries). Sequence data are taken from 2 mitochondrial genes, cytochrome b (Cytb) and nicotinamide adenine dinucleotide dehydrogenase 1 (ND1), and the hypervariable region 1 (HVRI) of the control region. Evolutionary signal is maximized from this data set by recovering cryptically recurrent substitutions from HVRI data (Phillips et al. 2009b). The end result is a robust and highly resolved genetic data set. Multiple analytical approaches have been used to test the major hypothesis that climate change has shaped contemporary genetic patterns of diversity and divergence in E. jubatus. These analyses help develop a long-term perspective on gene flow and elucidate the geographic center of origin of the species.

Materials and Methods

Sampling and molecular methodologies

All tissue samples used in this study were hind-flipper punches taken from pups biopsied at natal rookeries a few weeks after birth and were obtained in accordance with guidelines of the American Society of Mammalogists (Gannon et al. 2007) under Marine Mammal Protection Act permit 782-1532-02. Sampling included specimens from 48 of 61 total known breeding rookeries distributed across the North Pacific Rim (Fig. 1). Because of strong natal philopatry observed in E. jubatus, this sampling strategy leads to a high probability that pup haplotypes observed at rookeries are the result of breeding at that particular rookery.

Sequence data for the complete Cytb (1,140 base pairs [bp]) and HVRI (238 bp) for 1,021 individuals were available from previous studies (Baker et al. 2005; Harlin-Cognato et al. 2006; Phillips et al. 2009b). In addition to these Cytb and HVRI data, at least 1 individual representing each observed HVRI/Cytb haplotype combination also was sequenced for the entire ND1 gene (957 bp). Because a large proportion of total observed HVRI/Cytb composite haplotypes were observed at a low frequency (69% ≤ 5 observations), the sample subsets of individuals for ND1 likely provided an accurate overall depiction of the distribution of ND1 haplotypes and provided a reasonable alternative to sequencing all 1,021 samples for ND1. Polymerase chain reactions for ND1 were performed as described by Baker et al. (2005) using primers LGL 287 (CCTACGTGATCTGAGTTCAGACC) and LGL 563 (GGTATGAGCCCGATAGCTTA) with thermal profiles consisting of an initial denaturing cycle of 95°C for 5 min followed by 32 cycles of 95°C for 30 s, 50°C for 30 s, and 72°C for 2.5 min. For sequencing, internal sequencing primers Ejub-ND1-434seq(F) (CCATTATTCTCCTGTCAGTAC) and Ejub-ND1-636(R) (GGCCTGCTGCATATTCTACG) were used with a 50°C annealing temperature. Reactions were cleaned by centrifugation through Sephadex G-50 columns and then sequenced on an ABI 3730 platform (Applied Biosystems, Warrington, United Kingdom).

The ND1 sequences were aligned in Sequencher 4.8 (GeneCodes Corporation, Ann Arbor, Michigan) and examined for the presence of insertions–deletions and stop codons in the coding region to check for any signs that nuclear inserts had been sequenced incidentally. A neighbor-joining tree was constructed in PAUP* version 4.0b10 (Swofford 2003) to identify all unique ND1 haplotypes.

Data matrices construction and selection

To maximize phylogeographic signal from the mitochondrial DNA (mtDNA) data set an approach involving the concatenation of HVRI and coding region sequences, followed by the preferential weighting of coding region data and subsequent phylogeny construction, was used—similar to that of Bandelt et al. (2002)—to characterize HVRI recurrent substitution. Characterized recurrent substitutions were incorporated into the data matrix as pseudo positions at the end of the file, and the original position was removed. Phillips et al. (2009b) demonstrated little indication of recurrent substitution in mtDNA coding regions of E. jubatus, justifying the application of this strategy. Using this method basal tree structure was determined largely by coding region data, areas unlikely to be subject to homoplasy at the taxonomic level and lineage age considered here, and resolution at more terminal branching was improved by the disclosure of recurrent substitution within HVRI. To confirm the appropriateness of this method an Approximately Unbiased (AU) test was performed. In this test phylogenies constructed from genetic data matrices that were either extended to describe recurrent substitution at HVRI, were not extended to express this variation, or did or did not include the ND1 data (included as a variable here due to the nonexhaustive sampling strategy for this gene) were produced and then compared for explanatory power given each data matrix, following the methods outlined by Shimodaira (2002). The input phylogenies for this assessment were generated in MrBayes version 3.0 following all program author implementation specifications (Ronquist and Huelsenbeck 2003) and using the best fit model of evolution for each gene determined by the Akaike information criterion in ModelTest 3.7 (Posada and Crandall 1998; Table 1). For the AU test site-wise log-likelihoods for each nucleotide position in each data matrix were calculated independently given each phylogeny. Log files containing site-wise log-likelihood values served as input for the program CONSEL (Shimodaira and Hasegawa 2001) where the AU test was performed and the confidence set (α  =  0.05) was obtained.

Table 1.

Results of Akaike information criterion selected best-fit models of nucleotide substitution for each gene. α  =  gamma distribution shape parameter. HVRI(Cytb) refers to the hypervariable region 1 (HVRI) sequence extended to express recurrent substitution identified by cytochrome b (Cytb) weighting during neighbor-joining tree construction. Similarly, HVRI(Cytb/ND1) refers to the HVRI sequence extended to express recurrent substitution identified by Cytb and nicotinamide adenine dinucleotide dehydrogenase 1 (ND1) weighting during neighbor-joining tree construction.


Phylogeographic methods

The complete single-locus nested clade phylogeographical analysis (NCPA) procedure conducted in this study followed the methodologies outlined in several other studies as follows. A haplotype network was constructed in TCS version 1.3 (Clement et al. 2000) following the rules of parsimony using the data matrix previously selected through AU testing. Next, the network was nested into clades following the nesting algorithm described by Templeton et al. (1987), Templeton and Sing (1993), and Crandall (1996). Testing for significant associations (χ2) of clades with geography was performed in the program GeoDis by conducting 10,000 random permutations of clades (genetic variation) among sampling locations (rookeries—Posada et al. 2000; Templeton and Sing 1993). Geographic distances used were calculated as dispersal distances along the continental shelf, rather than great circle distances, to depict more accurately probable dispersal patterns in E. jubatus. Along-shelf dispersal distances (measured as distance between rookeries passing through any intermediate rookeries) are likely most appropriate because dispersal among rookeries by E. jubatus is thought to largely be confined to productive continental shelf waters, and great circle distances (point-to-point measurements around a sphere) would, in some instances, imply long dispersal routes through deep oceanic waters. Significant results from the GeoDis analysis were interpreted using the most recent inference key (made available 15 December 2008).

To cross-validate NCPA statistical significances an analysis of variance (ANOVA) was performed for each clade that returned a significant NCPA test statistic. For this analysis each haplotype/nested clade was treated as a group and the nesting clade as the population to which the groups belonged. Covariance components were calculated from the observed distribution of haplotypes/nested clades among each rookery. Significance was determined by randomly permuting the geographic distribution of groups among each respective population for 1,000 iterations and recalculating the test statistics for each iteration to obtain the null distributions.

Additionally, data on methane concentration recovered from ice cores with dense coverage extending 800,000 years ago were incorporated for comparison with dated significant single-locus NCPA inferences. This data set consisted of 2,103 data points with an average time resolution of approximately 380 years (Loulergue et al. 2008). Because methane is a globally mixed and long-lived greenhouse gas, it is considered a valuable indicator of climatic oscillation (Houghton et al. 2001). In the North Pacific Ocean climatic oscillations clearly have promoted the advance and retreat of glacial bodies that likely have modified biotic distributions (Grosswald and Hughes 2002; Mann and Peteet 1994; Westgate et al. 2001). If demographic events inferred by single-locus NCPA within E. jubatus are largely the result of glacial vicariance, we would expect global methane values (i.e., climate) at the time of these events to correlate with single-locus NCPA derived inferences in a meaningful way. Because preliminary analysis verified the normality of the methane data (Kolmogorov–Smirnov test; D  =  0.068, P > 0.99), and because sample size was large, significance of this relationship was assessed through a 2-tailed Z-test by treating methane concentrations at estimated dates of demographic events as random samples from a normally distributed population of values.

Because preliminary analysis of single-locus NCPA results indicated sequestering of the same demographic inference type into distinct temporal clusters, the probability of such a pattern arising by chance was assessed as a binomial probability. Eighteen unique inference outcomes exist in the most recently revised inference key. Therefore, the probability of success (i.e., obtaining any 1 specific inference in a single trial) was 1/18  =  0.0556, and the probability of failure to obtain any 1 specific inference in a single trial was 17/18  =  0.9444. The exact binomial formula is P(k out of n)  =  n!/[k!(nk!)](pk)[q(nk)], where n is the number of trials, k is the number of successes, p is the probability of success on a given trial, and q is the probability of failure on a given trial.

To date historical events inferred by single-locus NCPA for this assessment divergence dating using BEAST version 1.4.8 (Bayesian evolutionary analysis sampling trees—Drummond and Rambaut 2007) was carried out as described below for the Bayesian phylogeny, except with monophyly constraints added to ensure tree topology matched that observed in the haplotype network. Monophyly constraints were expressed in the starting newick tree used to initiate the analysis.

Divergence dating and Bayesian phylogeny estimation within the lineage of E. jubatus were performed simultaneously in the program BEAST. The basal position of Callorhinus within the family Otariidae, the sister relationship of the Arctocephaline clade to the Otaria, Zalophus, Eumetopias clade, and the ZalophusEumetopias split are well supported (Arnason et al. 2006; Higdon et al. 2007) and are the 3 node date priors included in this study. Multiple studies have assigned divergence dates within the family Otariidae, with dates estimated from fossil evidence representing a hard lower bound (Marshall 1990) and multiple molecular date estimates producing a range of values. For example, although the oldest available Eumetopias specimen (potentially an extinct form) recovered suggests a minimum ZalophusEumetopias divergence of 2 million years ago (mya—Repenning 1976), molecular dates for this split are 4.5, 6, and 8 mya, depending on the estimation algorithm and the loci included (Arnason et al. 2006; Higdon et al. 2007). As a conservative measure—conservative in that older molecular-based divergence estimates imply a greatly reduced rate of sequence evolution in otariids, which is not likely given the observed sequence diversity within Eumetopias—the dates estimated by Higdon et al. (2007), based on a 50-gene supertree that provided the most recent molecular dates, were used as date priors in this study (Callorhinus divergence, 8.2 mya; Arctocephaline–Otaria, Zalophus, Eumetopias clade split, 5.2 mya; and ZalophusEumetopias split, 4.5 mya).

Initial divergence dating within the Eumetopias lineage was performed using the 2 coding regions (Cytb and ND1), with each partition receiving its own model of DNA evolution as previously determined through model testing. In this analysis all codon positions were included. However, dates estimated from the inclusion of only 1st and 2nd codon positions gave very similar node dates, indicating that homoplasy at coding regions within the family Otariidae was not heavily influencing divergence estimation (data not shown). All program operations for this analysis strictly followed the program author guidelines. Dates estimated for the time to most recent common ancestor and the next most basal set of divergences within the Eumetopias lineage were recorded and retained as node date priors in a 2nd BEAST analysis to estimate dates of the more terminal nodes using only HVRI data. The objective of dating basal nodes with coding region data and terminal nodes with HVRI data was to partition substitutions from the 3 genes into areas of the tree where they retain the most phylogenetic information. From the final tree file the maximum clade credibility tree was obtained using TreeAnnotator (part of the BEAST package). In addition, because previous studies have demonstrated that the Eumetopias phylogeny contains substantial information about population history, a Bayesian skyline plot allowing for 5 discrete changes in population size was constructed using standard Markov chain Monte Carlo sampling procedures to estimate posterior distributions of theta (θ  =  Nefτ, where Nef is female effective population size and τ is generation length) through time using a flexible demographic model and directly from the sample of gene sequences (Drummond et al. 2005). Five demographic changes were allowed to identify historic demographic changes while at the same time not overparameterizing the analysis. A generation time for E. jubatus of 10 years was assumed for this estimation (Calkins and Pitcher 1982; York 1994).


Gene variability, recurrent substitution, and data matrix selection

Approximately 1.6 million sequenced base pairs (bp) were included in this study. We identified 18 Cytb and 82 HVRI haplotypes from 1,021 individual E. jubatus (Appendix I). The Cytb network (not shown) contained no reticulations and consisted of 2 halves that largely describe the 2 previously identified subspecies. The halves of the cladogram are separated by an inferred haplotype and 2 amino acid changes. Conversely, the network produced from the 82 observed HVRI haplotypes (not shown) contained numerous reticulations. We sequenced 202 individuals, including at least 1 representative of all HVRI/Cytb composite haplotypes, for the 957 nucleotides of ND1, revealing 11 haplotypes (see Appendix II for accession numbers). A parsimony network constructed from these haplotypes included 1 inferred haplotype and no reticulations (not shown). In summary, relationships among haplotypes within the Cytb and ND1 data sets were resolved although they exhibited only broadscale geographic resolution provided by the moderate number of observed variable sites (19 in Cytb and 12 in ND1). Conversely, the HVRI data set contained 41 variable sites within 238 bp; however, the frequency of reticulations throughout the network prevented an accurate depiction of the evolutionary relationships among haplotypes. The frequencies of all composite 3-gene mitochondrial haplotypes are described in detail in Appendix II.

Character changes mapped onto a neighbor-joining tree constructed from preferential weighting of Cytb and ND1 recovered 86 substitutions at 19 sites within HVRI, a 30% increase in the number of recurrent substitutions that otherwise would be detected by character mapping using only the HVRI data. Furthermore, results of the AU test based on site-wise log-likelihood scores selected the phylogeny constructed from Cytb, ND1, and HVRI extended to express all of these 86 substitutions as the optimal for describing 3 of 4 data matrices (Table 2). Notable is the observation that the phylogeny produced from the most simplistic data set (that including Cytb and HVRI) was best explained by the most complex data set (that including Cytb, ND1, and HVRI extended to express identified recurrent substitutions, termed HVRI(x)/Cytb/ND1 in Table 2 ). As a result, this data matrix was used for downstream phylogeographic analyses.

Table 2.

Approximately Unbiased (AU) test (Shimodaira 2002) of topological congruence based on site-wise log-likelihoods of 4 molecular data matrices for reciprocal gene trees. The tree receiving 1st-position ranking for each comparison is demarcated by bold and italicized numbering. See text for definition of acronyms.


Phylogeographic and phylogenetic patterns

The haplotype network constructed from the final data matrix generally maintains a network structure consisting of basal haplotypes that are frequent and widely distributed and terminal haplotypes that are sequestered geographically and occur in lower frequency relative to interior haplotypes (Fig. 2). From this network the divergence between subspecies is clearly demarked, although E. j. monteriensis is distributed across 2 parts of the network that are separated from each other by 2 widely distributed and 3 inferred haplotypes. An additional characteristic of this network is that, except for 1 instance, range-wide haplotypes (those observed in all 3 stocks) are sequestered to regions of the network otherwise consisting exclusively of haplotypes of E. j. jubatus. Of these range-wide haplotypes, their occurrence in E. j. monteriensis (eastern stock) is rare (6 of 7 of these haplotypes occur fewer than 4 times in E. j. monteriensis) and is restricted to southeastern Alaska (the westernmost region of E. j. monteriensis).


Nested design of the statistical parsimony network constructed from the HVRI(x)/Cytb/ND1 data matrix. Circles represent observed haplotype linkages, and the distributions and frequencies are indicated by coloring and size, respectively. Lines connecting haplotypes represent substitutions defining their relationships. Small black circles indicate inferred haplotypes. Reference Table 3 for information about dating and single-locus nested clade phylogeographical analysis (NCPA) results. See text for definitions of gene acronyms.


Table 3.

Summary of significant results for the nested clade phylogeographical analysis (NCPA). Values within the Dc and Dn columns are distance (km), with an indication of whether each value was significantly large (>) or small (<), and its associated P-value. n.s.  =  not significant at P  =  0.05. I − T  =  interior minus tips distance. Inferences were drawn from the key made available 15 December 2008. CRE  =  contiguous range expansion; IC  =  inconclusive outcome; RGF w/ IBD  =  restricted gene flow with isolation by distance; RGF/LDD  =  restricted gene flow or dispersal but with some long-distance dispersal; PF/LDC  =  past fragmentation or long-distance colonization (or both). Estimated dates are expressed as millions of years ago (mya).


Table 3.

Summary of significant results for the nested clade phylogeographical analysis (NCPA). Values within the Dc and Dn columns are distance (km), with an indication of whether each value was significantly large (>) or small (<), and its associated P-value. n.s.  =  not significant at P  =  0.05. I − T  =  interior minus tips distance. Inferences were drawn from the key made available 15 December 2008. CRE  =  contiguous range expansion; IC  =  inconclusive outcome; RGF w/ IBD  =  restricted gene flow with isolation by distance; RGF/LDD  =  restricted gene flow or dispersal but with some long-distance dispersal; PF/LDC  =  past fragmentation or long-distance colonization (or both). Estimated dates are expressed as millions of years ago (mya).


Eighteen demographic events were detected by single-locus NCPA, consisting of 5 contiguous range expansions, 5 instances of restricted gene flow with isolation by distance, 5 inferences of restricted gene flow or dispersal but with some long-distance dispersal, 2 past fragmentation or long-distance colonization events (or both), and 1 inconclusive outcome (Table 3). The inference of restricted gene flow with isolation by distance obtained at the total-cladogram level describes the overall demographic phenomenon of population structure of E. jubatus. The inference of past fragmentation or long-distance colonization (or both) for clade 4-2 captures, in part, the ancient divergence between E. j. jubatus and E. j. monteriensis. The relationship between subspecies within clade 4-2 and at the total cladogram level corresponds to the general pattern observed in the Bayesian phylogeny describing the subspecies divergence (Fig. 3). Tests performed for clades 4-1 and 4-3 both returned an inference of restricted gene flow or dispersal but with some long-distance colonization. These 2 clades are the largest subspecies-specific groupings in the analysis, and their associated inferences correspond with previously identified dispersal trends in E. jubatus described by dispersal generally occurring among adjacent rookeries with potential for longer distance migration due to their high vagility. The remaining inferences were recovered throughout portions of the network that describe further variation within subspecies.


Maximum clade credibility Bayesian phylogeny. Haplotypes are color coded by their observed distributions, and the black outgroup haplotype is Zalophus californianus according to Wilson and Reeder (2005). Nodes receiving a posterior probability support value of 0.5 or higher are indicated by small red stars. Significant nodes containing clades, where the sequestering of haplotypes into distinct geographic regions is observed, are numbered, and the estimated divergence date and 95% highest probability density intervals (HPDs) for these nodes are listed. Although posterior support for node 6 was below 0.5, this node is numbered for reference in the text.


As would be expected under the assumption of a low error rate of NCPA, of the 18 clades indicated as having statistically significant geographic associations through NCPA, 10 also yielded significant ANOVA test statistics (Table 4). Of these 10 clades, all but 1 received inferences involving some type of long-range movement or dispersal of haplotypes (contiguous range expansion, past fragmentation or long-distance colonization [or both], or restricted gene flow with some long-distance dispersal; events promoting significant among-group variation). Conversely, of the 8 clades returning nonsignificant ANOVA test statistics, all but 2 of these (inconclusive outcome and contiguous range expansion) involved restricted gene flow often with isolation by distance (events that would hinder accumulation of intergroup variation).

Table 4.

Results of ANOVAs for each clade receiving significant test statistics through the nested clade phylogeographical analysis (NCPA) statistical testing procedure. Descriptions of NCPA inference abbreviations are found in Table 3. P-values are for global FST calculated from the analysis, and significant statistics are bold and italicized. Significance was determined by 1,000 permutations of rookery occurrences of haplotypes among clades.


The mean methane concentration at dates estimated for the time to most recent common ancestor of all significant single-locus NCPA inferences was 468.44 parts per billion by volume (ppbv). Although this value is irrespective of the 95% highest probability density intervals (HPDs) surrounding the mean divergence dates, this value was significantly different from a mean methane concentration over the last 800,000 years of 519.97 ppbv (Z  =  2.41, P  =  0.015). By contrast the mean methane concentration at estimated times for all other clades was not significant (  =  499.16; Z  =  1.64, P  =  0.10). Furthermore, dates for significant demographic inferences show a pattern of temporal clustering that generally correlates with periods of low methane concentrations and time periods previously characterized as periods of glacial maxima. The inference returned at the total-cladogram level dates to within the 5th glacial maximum and all other inferences form 2 distinct groupings most closely associated with the 2nd and 3rd glacial maxima (Rohling et al. 1998; Fig. 4). The HPDs for the date estimates are generally large and span the time frame of methane-inferred climatic oscillations. However, 1 SD from the mean dates generally resides within major climatic periods in which the mean is positioned. This important nuance justifies the overall confidence in the datings.


Histogram of methane concentration (in parts per billion by volume) over the last 500,000 years. The nonlinearity of the x axis is a result of heterogeneity of methane sampling over time (Loulergue et al. 2008) and should be considered when interpreting this figure. Dates of time to most recent common ancestor for significant single-locus nested clade phylogeographical analysis (NCPA) clades are plotted onto this histogram, and the type of inference obtained for each clade is color coded. Small black hash marks above the histogram define point estimates for nonsignificant clades. Gray rectangles along the bottom of the histogram indicate time spans of 5 previous glacial periods, as defined by Rohling et al. (1998).


Comparing the patterns of demographic events with climatological cycling relative to estimated historic population sizes provided information on how the relationship between climate (glaciation) and population size potentially interacted to influence demography (Figs. 4 and 5). Notable is the observation that all inferences (6) clustering around the time of the 3rd glaciation describe processes involving restriction of gene flow. These 6 inferences refer to deeper clade nestings (3- and 4-step clades) and date to periods of low ancestral population sizes. The most common inference within this cluster was restricted gene flow or dispersal but with some long-distance dispersal, occurring 4 times. The binomial probability that 4 of 6 inferences within this time period would be the same by chance is 0.00013. In contrast, inferences clustering around the 2nd glaciation or earlier all refer to shallower nestings (1- and 2-step clades) corresponding to a time of population size increase in E. jubatus. In this cluster several inferences describe some sort of restricted gene flow; however, the most common type of inference within this cluster was contiguous range expansion, occurring 5 times. The probability of such a chance occurrence is 0.00017. No significant demographic inferences were recovered in more recent times, around the last glacial maximum, a period in which population sizes were estimated to be highest.


Bayesian skyline plot with mean population size and 95% highest probability density intervals (HPDs) plotted against time. As in Fig. 4, significant NCPA inferences are plotted and color coded.


The time to most recent common ancestor for E. jubatus (Fig. 3, node 1), based on coding region data and outgroup node priors from Higdon et al. (2007), was estimated at 0.360 mya (95% HPD  =  0.145–0.876). The initial divergence from node 1 resulted in the sequestering of all exclusively E. j. monteriensis haplotypes into a single, yet polyphyletic, grouping (node 2). Terminal to node 2 were all 40 haplotypes exclusive to E. j. monteriensis and 22 exclusive to E. j. jubatus, with 8 and 12 of these exclusive to the western and Asian stocks, respectively. Within this grouping 2 major clades were identified, 1 of which (node 4) was comprised exclusively of haplotypes of E. j. monteriensis, except for 4 haplotypes observed only in the Asian stock. Node 6, being sister to node 3 (Fig. 3), retained the remaining haplotypes observed exclusively in E. j. monteriensis in addition to haplotypes occurring across the species' distribution; however, the posterior probability of this node was <0.5. The basal clades terminal to nodes 4 and 6 consisted of haplotypes of E. j. jubatus.

On the other half of the phylogeny terminal to node 3 all but 5 range-wide haplotypes and 1 haplotype observed in the western and eastern stocks—3 of 48 occurrences of this haplotype were in the eastern stock—are haplotypes exclusive to E. j. jubatus. Although within this clade 11 haplotypes were observed throughout E. j. jubatus, 17 and 13 haplotypes were restricted to the western and Asian stocks, respectively. The divergence observed at clade 5 leads to 2 lineages, 1 of which contains haplotypes found either exclusively in the Asian stock or shared between the Asian and western stock, but none exclusive to the western stock.

Although posterior support for some nodes in this phylogeny was <0.5, the major phylogenetic pattern depicted supports previously hypothesized relationships among the 3 stocks and the 2 subspecies. However, the current results both describe the early divergence in the species history between the 2 extant and currently recognized subspecies and also indicate a complex history leading to the current distribution of haplotypes. Moving from the base to the tips of the phylogeny of E. jubatus, clear patterns of lineage sorting at different stages of completion are observable. Although terminal to the initial divergence all haplotypes of E. j. monteriensis (eastern stock) are sequestered to node 2 (Fig. 3), it is not until about 200,000 years later at node 4 that a large group of haplotypes of E j. monteriesis approach monophyly—1 Asian stock haplotype is within the eastern stock clade that is itself terminal to an exclusively Asian stock clade. The relationships among haplotypes terminal to node 6 present a pattern of less complete lineage sorting, with some tip clades exclusively E. j. monteriensis, with others still retaining many haplotypes of E. j. jubatus. Although node 6 was unsupported, a comparison of the completion of sorting between nodes 4 and 6 illustrates the stochastic nature of the lineage sorting process; although nodes 6 and 4 likely have similar coalescence times, sorting is closer to completion terminal to node 4.


In this study single-locus NCPA was used to draw inferences on population history. Multiple independent statistical approaches were used to validate the inferences returned by single-locus NCPA. ANOVA returned patterns of group (clade) significance that correspond to what would be expected given the types of demographic events returned for each clade; demographic events that promote significant among-group variance allocation (those events that distribute genetic diversity among rookeries) were usually significant, and those demographic events that promote within-group variance retention were found most often to be nonsignificant. Furthermore, results from the single-locus NCPA, in combination with reconstructions of historical Nef and ice-core inferred climatological cycling, suggest that the patterns of demographic events are related to the climate and population size at the time. The chance probability of such a pattern was assessed in this study by using Z-tests and binomial probabilities and found to be low. These combined statistical observations describe a pattern of phylogeography of E. jubatus that includes an interaction between population size, climate, and their dependence on linearly distributed rookeries.

Our results support the hypothesis that climate change has a major role in shaping the observed distribution of genetic diversity in E. jubatus. One of the major inferences pertaining to this hypothesis is that glacial vicariance has largely shaped the evolutionary history of E. jubatus. The observation that historic events related to restricted gene flow generally occurred during a time of low Nef suggests that vicariance, in conjunction with genetic drift acting on low ancestral numbers of haplotypes, promoted geographic sequestering throughout portions of the distribution. In addition, inferences of occasional long-distance dispersal during this time potentially describe the species' ability to establish new territory via dispersal during a time of low Nef and reduced competition for breeding territory. Conversely, all inferences of contiguous range expansion correlate with a population expansion as described by the Bayesian skyline plot (Fig. 5). This indicates that as Nef increased the response to glacial vicariance changed. During this time the population of E. jubatus as a whole was expanding both in numbers and in range. Individuals would have been forced to disperse greater distances to find breeding territory as Nef increased and glacial cover reduced the number of ice-free rookeries. Finally, the lack of significant demographic inferences recovered around the time of the last glacial maximum is likely the result of large, sustained Nef, previous saturation of all suitable habitats, and perhaps a lack of phylogeographic resolution at this timescale.

In general, these results can be interpreted as descriptive of how abiotic factors can influence the evolutionary fate of a species. Specifically for E. jubatus, one of the results of this study is that relatively warm periods have promoted population expansion and dispersal by increasing available rookery habitat. Perhaps the most important aspect of the current results is the implication that future climate change has the potential to affect distributions and demography of contemporary populations of species in ways similar to that documented in this study. It is clear that climate change causes stress, including the necessity of finding new breeding sites, adapting to new food sources, and so on. It can be inferred from the present study that population size influences the capacity for dealing with such stressors, because the genetic signatures resulting from similar climatic patterns are distinctly different for large and small populations.

Although the Bayesian phylogeny presented in Fig. 3 describes relationships among haplotypes similar to those presented in the haplotype network, this perspective provided additional insight into both historical and contemporary processes. This tree describes the ancient divergence between subspecies, with node 1 representing the base of the lineage of E. jubatus describing the initial documentable divergence 0.360 mya. Following from this estimation, during the 5th glacial period when population sizes were low, the distribution of E. jubatus was split by glacial vicariance into at least 2 disjunct populations initiating subspecific differentiation. The observation that 3 Asian stock haplotypes form the base of the large E. j. monteriensis clade terminal to node 4 describes an ancient “out of Asia” relationship between the extreme ends of the distribution of E. jubatus and indicates that Asian portions of the species' range served as the source population for the eastward expansion. Although node 6 is unsupported, the relationship observed between sister clades here is similar to that observed at node 4, with an E. j. jubatus clade forming the base of a clade consisting largely of E. j. monteriensis. That no haplotypes were observed in the Asian stock and E. j. monteriensis, but not the western stock, indicates the lack of contemporary dispersal between the Asian stock and E. j. monteriensis. However, the observation that 6 of 7 range-wide haplotypes belong to E. j. jubatus clades does indicate recent directional dispersal from west to east. The oldest Eumetopias fossil was recovered from the Onma formation in Japan dated to the Pliocene about 2 mya (thought to represent an extinct member of the genus EumetopiasKaseno 1951; Mitchell 1968; Repenning 1976). This information, paired with the findings from the current study, indicates a potential geographic center of origin for Eumetopias.

In this study a novel method was used to account for generally cryptic recurrent substitutions. This improved resolution by decreasing the phylogenetic obscurity produced by homoplasy while improving confidence in branching order. Potentially, this approach is applicable to any phylogenetic analysis in which coding genes, or otherwise conservative loci, are used in combination with a more rapidly evolving locus in which recurrent substitutions produce homoplasy. In the event of evidence for coding region homoplasy, these regions should not be weighted to resolve hypervariable region recurrent substitutions. To determine the appropriateness of using this method several aspects of the data set should be considered. First, the estimated time to most recent common ancestor is an important parameter, because the greater the age of the lineage, the more time for homoplasy to accumulate in all regions. By using this information in combination with knowledge about the mutation rate of coding regions, network construction and evaluation (whether coding region–based networks contain reticulations), and consistency and homoplasy index calculations, a logical decision can be made regarding the use of this approach.

In summary, by the implementation of multiple statistical approaches to maximize genetic signal and draw evolutionary inferences, a strong, multifaceted perspective of evolutionary history of E. jubatus has been developed. This history reflects a major influence of climate change and glacial vicariance, with demographic response contingent upon population size at the time. In addition, phylogenetic patterns indicate the direction of dispersal over evolutionary timescales and, paired with limited fossil data, present a hypothesis for the geographic center of origin of E. jubatus located in Asia. The results of this study can serve as a working hypothesis for future studies to develop a better understanding of the evolution of E. jubatus.


We thank the many people who helped in obtaining tissue samples—especially V. Burkanov (coordinator for the collecting trips in Russia), T. Loughlin, and D. Calkins—and also the National Oceanic and Atmospheric Administration cruises in the Aleutians and the Gulf of Alaska and the Alaska Department of Fish and Game. K. M. Nichols, O. E. Rhodes, and J. A. DeWoody provided valuable insights and suggestions during the course of the doctoral studies of CDP. All sea lion tissue samples were collected under authorization of United States Marine Mammal Permits 358–1888 (Alaska Department of Fish and Game) and 782–1532 (National Marine Mammal Lab). Comments from 2 anonymous reviewers were very constructive and greatly improved this manuscript. This study was conducted as part of the doctoral dissertation of CDP at Purdue University. Funding was provided by the National Marine Fisheries Service, the Lily Endowment, and internal funds provided by the Department of Forestry and Natural Resources, Purdue University.

Literature Cited


U Arnason et al. 2006. Pinniped phylogeny and a new hypothesis for their origin and dispersal. Molecular Phylogenetics and Evolution 41:345–354 Google Scholar


A. R Baker et al. 2005. Variation of mitochondrial control region sequences of Steller sea lions, Eumetopias jubatus: the three-stock hypothesis. Journal of Mammalogy 86:1075–1084 Google Scholar


H. J Bandelt L Quintana-Murci A Salasand V Macaulay 2002. The fingerprint of phantom mutations in mitochondrial DNA data. American Journal of Human Genetics 71:1150–1160 Google Scholar


J. W Bickham T. R Loughlin J. K Wickliffeand V. N Burkanov 1998. Geographic variation in the mitochondrial DNA of Steller sea lions: haplotype diversity and endemism in the Kuril Islands. Biosphere Conservation 1:107–117 Google Scholar


J. W Bickham J. C Pattonand T. R Loughlin 1996. High variability for control-region sequences in a marine mammal: implications for conservation and biogeography of Steller sea lions (Eumetopias jubatus). Journal of Mammalogy 77:95–108 Google Scholar


D. G Calkinsand K. W Pitcher 1982. Population assessment, ecology and trophic relationships of Steller sea lions in the Gulf of Alaska. Pp. 445–546 in Environmental assessment of the Alaskan continental shelf: final reports of principal investigators. United States Department of Commerce, National Oceanic and Atmospheric Administration, Juneau, Alaska 19:1–565 Google Scholar


M. D Clement D Posadaand K. A Crandall 2000. TCS: a computer program to estimate gene genealogies. Molecular Phylogenetics and Evolution 3:102–113 Google Scholar


K. A Crandall 1996. Multiple interspecies transmissions of human and simian T-cell leukemia/lymphoma virus type I sequences. Molecular Biology and Evolution 13:115–131 Google Scholar


A. J Drummondand A Rambaut 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7:214 Google Scholar


A. J Drummond A Rambaut B Shapiroand O. G Pybus 2005. Bayesian coalescent inference of past population dynamics from molecular sequences. Molecular Biology and Evolution 22:1185–1192 Google Scholar


W. L Gannon R. S Sikes and the Animal Care and Use Committee of the American Society of Mammalogists. 2007. Guidelines of the American Society of Mammalogists for the use of wild mammals in research. Journal of Mammalogy 88:809–823 Google Scholar


M. G Grosswaldand T. J Hughes 2002. The Russian component of an arctic ice sheet during the last glacial maximum. Quaternary Science Reviews 21:121–146 Google Scholar


A Harlin-Cognato J. W Bickham T. R Loughlinand R. L Honeycutt 2006. Glacial refugia and the phylogeography of Steller's sea lion (Eumatopias jubatus) in the North Pacific. Journal of Evolutionary Biology 19:955–969 Google Scholar


J. W Higdon O. R. P Bininda-Emonds R. M. D Beckand S. H Ferguson 2007. Phylogeny and divergence of the pinnipeds (Carnivora: Mammalia) assessed using a multigene dataset. BMC Evolutionary Biology 7:216 Google Scholar


J. I Hoffman K. K Dasmahapatra W Amos C. D Phillips T. S Gelattand J W Bickham 2009. Contrasting patterns of genetic diversity at three different genetic markers in a marine mammal metapopulation. Molecular Ecology 18:2961–2978 Google Scholar


J. I Hoffman C. W Matson W Amos T. R Loughlinand J. W Bickham 2006. Deep genetic subdivision within a continuously distributed and highly vagile marine mammal, the Steller's sea lion (Eumetopias jubatus). Molecular Ecology 15:2821–2832 Google Scholar


J. T Houghton et al. 2001. Climate change 2001: the scientific basis. Contributions of working group I to the third assessment report of the intergovernmental panel on climate change. Cambridge University Press, Cambridge, United Kingdom Google Scholar


Y Kaseno 1951. Pliocene pinniped remains from Kanazawa, Ishikawa Prefecture, Japan. Translations of the Proceedings of the Palaeontological Society of Japan 2:57–64 Google Scholar


K. W Kenyonand D. W Rice 1961. Abundance and distribution of the Steller sea lion. Journal of Mammalogy 42:223–234 Google Scholar


T. R Loughlin M. A Perezand R. L Merrick 1987. Eumetopias jubatus. Mammalian Species 283:1–7 Google Scholar


T. R Loughlin A. S Perlovand V. A Vladimirov 1992. Range-wide survey and estimation of total number of Steller sea lions in 1989. Marine Mammal Science 8:220–239 Google Scholar


L Loulergue et al. 2008. Orbital and millennial-scale features of atmospheric CH4 over the past 800,000 years. Nature 453:383–386 Google Scholar


D. H Mannand D. M Peteet 1994. Extent and timing of the last glacial maximum in southwestern Alaska. Quaternary Research 42:136–142 Google Scholar


C. R Marshall 1990. The fossil record and estimating divergence times between lineages: maximum divergence times and the importance of reliable phylogenies. Journal of Molecular Evolution 30:400–408 Google Scholar


E. D Mitchell 1968. The Enaliarctinae, a new group of extinct aquatic Carnivora and a consideration of the origin of the Otariidae. Bulletin of the American Museum of Natural History 151:203–284 Google Scholar


Crowe, G O'Corry- et al. 2006. Demographic independence along ecosystem boundaries in Steller sea lions revealed by mtDNA analysis: implications for management of an endangered species. Canadian Journal of Zoology 84:1796–1809 Google Scholar


C. D Phillips J. W Bickham J. C Pattonand T. S Gelatt 2009a. Systematics of Steller sea lions (Eumetopias jubatus): subspecies recognition based on concordance of genetics and morphometrics. Occasional Papers, Museum of Texas Tech University 283:1–15 Google Scholar


C. D Phillips et al. 2009b. Assessing substitution patterns, rates, and homoplasy at HVRI of Steller sea lions, Eumetopias jubatus. Molecular Ecology 18:3379–3393 Google Scholar


D Posadaand K. A Crandall 1998. Modeltest: testing the model of DNA substitution. Bioinformatics 14:817–818 Google Scholar


D Posada K. A Crandalland A. R Templeton 2000. GeoDis: a program for the cladistic nested analysis of the geographical distribution of genetic haplotypes. Molecular Ecology 9:487–488 Google Scholar


C. A Repenning 1976. Adaptive evolution of sea lions and walruses. Systematic Zoology 25:375–390 Google Scholar


E. J Rohling M Fenton F. J Jorissen P Betrand G Ganssenand J. P Caulet 1998. Magnitudes of sea-level lowstands of the past 500,000 years. Nature 394:162–165 Google Scholar


F Ronquistand J. P Huelsenbeck 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19:1572–1574 Google Scholar


H Shimodaira 2002. An approximately unbiased test of phylogenetic tree selection. Systematic Biology 51:492–508 Google Scholar


H Shimodairaand M Hasegawa 2001. CONSEL: for assessing the confidence of phylogenetic tree selection. Bioinformatics 17:1246–1247 Google Scholar


D. L Swofford 2003. PAUP*: phylogenetic analysis using parsimony (*and other methods). Version 4.0b 10. Sinauer Associates, Inc., Publishers, Sunderland, Massachusetts Google Scholar


A. R Templeton E Boerwinkleand C. F Sing 1987. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. I. Basic theory and an analysis of alcohol dehydrogenase activity in Drosophila. Genetics 117:343–351 Google Scholar


A. R Templetonand C. F Sing 1993. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping. IV. Nested analyses with cladogram uncertainty and recombination. Genetics 134:659–669 Google Scholar


A. W Tritesand C. P Donnelly 2003. The decline of Steller sea lions Eumetopias jubatus in Alaska: a review of the nutritional stress hypothesis. Mammal Review 33:3–28 Google Scholar


J. A Westgate S. J Preece D. G Froese R. C Walter A. S Sandhuand C. E Schweger 2001. Dating early and middle (Reid) Pleistocene glaciations in central Yukon by tephrochronology. Quaternary Research 56:335–348 Google Scholar


D. E Wilsonand D. M Reeder (eds.). 2005. Mammal species of the world: a taxonomic and geographic reference. 3rd ed. Johns Hopkins University Press, Baltimore, Maryland Google Scholar


A. E York 1994. The population dynamics of northern sea lions, 1975–1985. Marine Mammal Science 10:38–51 Google Scholar


Appendix I

A Sample sizes (n) for Eumetopias jubatus categorized at 4 geographic scales: subspecies, genetically identifiable stock, region, and deme (rookery). Rookery numbers correspond to Fig. 1.


Appendix II

Frequency of HVRI/Cytb/ND1 haplotype linkages across the 3 genetically identifiable stocks. The 1st series of letters in the name refers to the hypervariable region 1 (HVRI) haplotype name, the following number refers to the cytochrome b (Cytb) haplotype name, and the number after the N in each name refers to the nicotinamide adenine dinucleotide dehydrogenase 1 (ND1) haplotype name. GenBank accession numbers for HVRI haplotypes: AY340876–AY340937, FJ948491–FJ948546. GenBank accession numbers for Cytb haplotypes: DQ144995–DQ145021, FJ948486–FJ948490. GenBank accession numbers for ND1 haplotypes: GQ477068–GQ477078.


Appendix II.


American Society of Mammalogists
C. D. Phillips, T. S. Gelatt, J. C. Patton, and J. W. Bickham "Phylogeography of Steller sea lions: relationships among climate change, effective population size, and genetic diversity," Journal of Mammalogy 92(5), 1091-1104, (14 October 2011).
Received: 1 September 2010; Accepted: 1 April 2011; Published: 14 October 2011

Back to Top