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 (χ²) 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!(n − k!)](p^k)[q^(n−k)], 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 Zalophus–Eumetopias 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 Zalophus–Eumetopias 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 Zalophus–Eumetopias 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 (θ  =  N_efτ, where N_ef 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).

Fig.2.

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 D_c and D_n 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 D_c and D_n 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.

Fig.3.

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 F_ST 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 (X¯  =  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.

Fig.4.

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.

Fig.5.

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.