Skip to main content
  • Research article
  • Open access
  • Published:

Late Quaternary loss of genetic diversity in muskox (Ovibos)



The modern wildherd of the tundra muskox (Ovibos moschatus) is native only to the New World (northern North America and Greenland), and its genetic diversity is notably low. However, like several other megafaunal mammals, muskoxen enjoyed a holarctic distribution during the late Pleistocene. To investigate whether collapse in range and loss of diversity might be correlated, we collected mitochondrial sequence data (hypervariable region and cytochrome b) from muskox fossil material recovered from localities in northeastern Asia and the Arctic Archipelago of northern North America, dating from late Pleistocene to late Holocene, and compared our results to existing databases for modern muskoxen.


Two classes of haplotypes were detected in the fossil material. "Surviving haplotypes" (SHs), closely similar or identical to haplotypes found in modern muskoxen and ranging in age from ~22,000 to ~160 yrbp, were found in all New World samples as well as some samples from northeastern Asia. "Extinct haplotypes" (EHs), dating between ~44,000 and ~18,000 yrbp, were found only in material from the Taimyr Peninsula and New Siberian Islands in northeastern Asia. EHs were not found in the Holocene muskoxen specimens available for this study, nor have they been found in other studies of extant muskox populations.


We provisionally interpret this evidence as showing that genetic variability was reduced in muskoxen after the Last Glacial Maximum but before the mid-Holocene, or roughly within the interval 18,000-4,000 yrbp. Narrowing this gap further will require the recovery of more fossils and additional genetic information from this interval.


The large-mammal (megafaunal) diversity of North America and Eurasia was considerably reduced by extinctions that occurred around the time of the Pleistocene/Holocene transition (i.e., ca. 10,000 years ago) [1]. Interestingly, some megafaunal species (e.g., woolly mammoth, Mammuthus primigenius) became extinct at this time in North American and Eurasian continental areas but survived on nearby islands for thousands of years thereafter [2, 3]. Some other species whose ranges likewise straddled the Old and New Worlds became extinct on one continent, but not the other, at the end of the Pleistocene (e.g., horse, Equus caballus; saiga, Saiga tatarica). A third pattern, involving a much lengthier decline, characterizes Ovibos moschatus (tundra muskox), the subject of this paper (fig. 1, 2).

Figure 1
figure 1

Sketch map of Holarctic region in polar view, noting places mentioned in text and specimens sampled for mtDNA analysis. Dark gray area, muskox range ca. AD 1870 [6]; C, cape; Cr, creek; I, island; L, lake; R, river. During the late Pleistocene, Ovibos ranged from western Europe to North America via Beringia. In the Old World, presence of muskox fossils on islands north of mainland Asia presumably indicates that Ovibos was able to extend its range during full glacial times to the subaerial parts of the continental shelf (dashed line). In the New World, in addition to central Alaska, unglaciated areas in the western Arctic Archipelago may have continuously supported muskoxen also, although the finite 14C fossil record for these islands does not begin until about 6800 yrbp [40].

Figure 2
figure 2

Finite 14C date distributions for fossil specimens of Ovibos. (A) northern North America and Greenland (N = 72; X = 3222 yr, SD = 4313 yr) and (B) Arctic periphery of Asia east of Kara Sea (N = 38; X = 22,522 yr, SD = 11,145 yr), collated from the literature [32,40,41]. In addition to a notable difference in dating intensity, with many more published 14C dates being available for New World Ovibos, there is also a considerable difference in sample statistics. Skewed date distribution for New World Ovibos reflects long-term effect of former icecaps on amount of habitable area [42]. In northeastern Asia, where only small, local (cordilleran) icecaps formed, potential range was less affected and Ovibos was presumably continuously present during the last 40–50,000 years. Star in B denotes 7,000 yr gap between "last" Pleistocene and "first" Holocene dated occurrences of Ovibos in Asia. Slight offset in "period of heightened extinction" in New World vs. Old World, amounting to approximately 500–1000 yr, reflects a moderate difference in "last" occurrence dates for mainland mammoths in the two areas [32].

According to the available fossil record, the lineage terminally represented by O. moschatus arose in Eurasia but entered North America at least as early as the early Pleistocene [46]. At roughly the same time several other closely-related Asian ovibovine lineages (e.g., Soergelia, Praeovibos) also entered North America, where they evidently prospered before disappearing later in the Pleistocene. By the start of the Holocene, muskox diversity had been reduced to surviving O. moschatus; the species is now native only to the islands of the Arctic Archipelago and Greenland (fig. 1), if late 20th-century efforts to reintroduce muskoxen into various parts of their ancient range are ignored [6].

This lengthy history of muskox diversification in Quaternary Holarctica stands in sharp contrast to the extremely limited genetic diversity characterizing modern O. moschatus. Groves' [7, 8] analysis of sequence data for cytochrome b (cyt b) and the mitochondrial hypervariable region in 37 individuals of O. moschatus from various localities in Alaska, Canada, and Greenland revealed only marginal genetic differentiation within extant populations. Groves [7] identified a total of 8 hypervariable region haplotypes, but these differed individually by only 1 or 2 nucleotide substitutions in approximately two-thirds of the cases, for a maximum of 9 differences – a total of 1.3% variation across 697 bp (base pairs), including insertions and deletions. Extreme genetic homogeneity is also indicated by microsatellite analysis [[9], but see also [10]] and the nuclear Mhc DRB 3 gene [11]. As several different populations were sampled for these studies, in the absence of any other plausible explanation it seems reasonable to conclude that living muskoxen simply lack genetic diversity. The nuclear and mitochondrial data therefore complement each other and imply that muskoxen either passed through a long period of very small effective population sizes, or experienced one or more genetic bottlenecks in very recent times. The investigation reported here was prompted by the possibility that ancient DNA (aDNA) analyses might shed some light on the complex history of extirpation, replacement, and differential survival that have evidently affected muskoxen populations during the late Quaternary. Ancient DNA analysis has been successfully applied to late Pleistocene and Holocene (including the early 20th century) remains by various groups [e.g., [1222]]. Because of existing uncertainty concerning species boundaries among some late Quaternary members of tribe Ovibovini, the terms "ovibovine" and "muskox" will be used to make general reference to such taxa, without endorsing any particular view as to their taxonomic distinctiveness or geographical distribution.


Our results permit empirical recognition of two haplotype groupings within late Quaternary muskoxen. The first consists of haplotypes recovered from Pleistocene samples collected at various localities in northeastern Asia, including islands north of the mainland (fig. 1, 3). These haplotypes, which differ from all New World material so far examined as well as from the few known Asian specimens that are of Holocene age, appear to have completely disappeared and may therefore be grouped as "EHs" (extinct haplotypes). All other known haplotypes, which include those found in all living muskoxen, are very similar inter se and thus comprise a second grouping, "SHs" (surviving haplotypes).

Figure 3
figure 3

Timeline showing distribution of specimens yielding sequence information. Symbols indicate whether a given specimen exhibited an "extinct" (EH) or a "surviving" (SH) haplotype. Specimens yielding substantial sequence data are in bold. SHs (i.e., haplotypes represented in extant O. moschatus) were found in all New World samples tested, as well as some Old World samples. EHs have been found so far only in Asian samples.

Sequence verification

A total of 9 fossil samples yielded reproducible sequence covering cyt b and the hypervariable region of the mtDNA genome; with the exception of OMTai14, each fragment from these samples was successfully amplified at least twice (additional fig. 1, 2, 3; table 1). One other sample (OMBolL15564) yielded consistent results for cyt b but the HV sequence could not be accurately scored (table 1, additional fig. 4). Sequences were then determined via multiple clones produced from each PCR product, a standard aDNA protocol for insuring correct sequence determination [12, 22] (additional fig. 1, 2, 3). In addition, several samples were independently extracted by different individuals, amplified, and cloned. Clones were sequenced in one or another of the three separate laboratories involved in this investigation (AMNH, New York; Medigenomix GmbH, Munich; SR&D GmbH, Martinsried).

Table 1 Information on bone samples utilized in aDNA studies, with 14C age estimates

For example, to determine the cyt b sequence for sample OMTai38, two PCRs were performed for the first fragment and second overlapping fragment and five clones sequenced for each PCR product (additional fig. 2). The sequences summarized in tables 3, 4, 5 are thus the consensus of independent replications and multiple-clone sequences for each individual PCR fragment (with a few exceptions), and are unlikely to be laboratory artifacts (additional fig. 1, 2, 3). For this particular sample, all clones from the 3' fragments were identical with one exception: in clone 1.2, there is a C-T transition which most likely represents a polymerase error. In the 5' fragment, all five clones of the first PCR were identical. In the second fragment three classes of sequence were obtained, represented by clones 2.3 and 2.4; these were identical to the clones from the first PCR and matched the 3' fragment in the overlap. The second and third class of sequences (2.1, 2.2, and 2.5) were different from all other clones obtained in the study, unique to the second PCR amplification, and did not match the 3' fragment in the overlap. Blast searching the sequence revealed perfect or near perfect homology for both divergent sequences with other bovid species such as the domestic sheep, Ovis aries. These likely represent nuclear mitochondrial DNA sequences (Numts) and thus the sequence represented in the two independent fragments was scored as the bona fide mtDNA sequence, the consensus of which was used for further study.

Table 3 DNA sequence variation in the cyt b gene
Table 4 Sequence comparisons between Holocene and Pleistocene samples
Table 5 Summary of variable sites for 1,180 bp of cyt b sequence, OMTai39 compared to 5 modern muskox

This example represents a relatively extreme case, as most PCRs yielded identical or nearly identical sequences that differed by minor among-clone variations that can be explained by DNA polymerase/sequencing errors or DNA damage in the templates (for example, note the clone sequences for OMTai39 in additional fig. 2). Where a position or several positions were not determinable after two independent amplifications, additional PCRs were performed, clones sequenced, and the majority base scored as the bona fide mtDNA sequence for each position of each sequence. In other words, each base of each sequence was covered by multiple PCRs and multiple clone sequences to determine the majority base for a given position as in Krings et al. [12].

In addition, six other samples (OMTai95, OMWra2, OMTai31670, OMYuk02, OMArA09, and OMYak12) yielded some sequence information, but the amount of DNA recovered was insufficient to permit two independent replications for each fragment and/or retrieve all four fragments of interest in this study (see table 1).

Overlapping PCR reactions were used to identify possible Numts and to prevent mischaracterizing Numts as organellar mtDNA [14, 23]. In some species DNA extracted from hair may yield Numts at a higher frequency than do other tissues [24], making it a useful proxy for gauging the possibility of Numt contamination. A hair sample (OMModCl, table 1) from a living muskox was amplified, products cloned, and multiple sequences determined using several of the primer combinations used in this study. Sequences were consistent in the overlapping regions and matched DNA from other modern muskoxen, indicating that the primers used in this study do not favor Numts over organellar mtDNA (additional fig. 1, 2).

Cytochrome b

A 114-bp segment of cyt b was successfully retrieved from most of the working fossil samples (table 3); in three cases (OMTai14, 39, and 46) it was possible to amplify a 376 bp fragment. All specimens identified as SHs that yielded cyt b sequence were identical. Only partial sequence data could be recovered for OMWra2 and OMYak12 (table 3). However, for comparable runs of sequence, these specimens were identical to all other SHs.

Owing to institutional sampling restrictions and the low quality of retrieved DNA, analysis of OMTai95 had to be based on one repetition of each of the two overlapping fragments. Although results for this specimen are therefore tentative, its sequence does not differ from that of OMTai38 and OMTai 46 (not shown). Thus, despite some individual variation, in all cases samples could be classed unambiguously as either EH or SH.

Specimens in the EH group differed uniformly from those classed as SH by three fixed differences (table 3, positions 324, 346, and 366). There were also some internal differences within the EH group. OMTai14 and OMTai39 possessed an additional base substitution not seen in other members (position 436), and OMBolL15564 expressed a unique A to G substitution at position 425. Given the difficulties encountered in determining the mitochondrial HV region sequence for this last sample, this result was excluded from analysis.

The Holocene samples (N = 43) exhibited no variations. Only one haplotype was retrieved (table 4), which is consistent with Groves and Shields' [8] evidence for limited among-individual variation at the cyt b locus in living muskoxen. By contrast, the Pleistocene samples yielded 2 haplotypes among 4 sequences and thus greater among-individual haplotype and nucleotide diversity in a much smaller sample.

As a proxy for resolving the magnitude of differences among cyt b haplotypes, 1,180 bp of sequence in overlapping fragments was determined for one of the best-performing specimens, OMTai39 (table 5) and aligned by eye against 5 modern muskox cyt b sequences deposited in GenBank (accession numbers U17862 and U90300-U90303, here designated as OMMod01-05). Among modern muskoxen, within-group differences in cyt b sequences range from 1 to 4 base substitutions (table 5). By comparison, OMTai39 differed from all modern samples by 5–7 substitutions. Domestic sheep (Ovis aries) [GenBank: NC001941] and goat (Capra hircus) [GenBank: NC005044] were used as outgroups for phylogenetic analysis. The takin Budorcas taxicolor was not included in the analysis as Groves and Shields [8] have shown that the relationship between this species and muskox is not demonstrably closer than that of other caprines. When submitted to phylogenetic testing using maximum parsimony, neighbor joining, and maximum likelihood methods, OMTai39 placed in all analyses outside modern muskoxen cyt b sequences (bootstrap support, 80–98%) supporting the distinctiveness of Pleistocene vs. Holocene haplotypes (fig. 4).

Figure 4
figure 4

Cyt b (1,180 bp) from a representative EH specimen (OMTai39) compared to modern Ovibos moschatus (OMMod01-05). Maximum parsimony tree using goat and sheep as outgroups. Calculated bootstrap support is shown for parsimony/neighbor joining/maximum likelihood at each node. Although analysis returned an evaluation of "ns" (no statistical support) for a given grouping in only one instance, significance of results should not be overemphasized. Nevertheless, it is parsimonious to conclude that (1) all muskoxen form a cluster as against sheep and goat, and (2) moderns form a structureless subcluster as against the extinct haplotype represented by OMTai39.

Hypervariable region

Samples from Yukon and islands in the Canadian Arctic Archipelago (OMYuk01, OMArA03, OMArA04, and OMArA05) displayed minor haplotype variation affecting 2–4 positions in the studied portion of the hypervariable region (additional fig. 1). No fixed differences were found between these samples and sequence data on modern muskoxen in GenBank [7](cf. OMMod01-08, OMModCl, table 5).

One position in the GenBank samples was not observed in any of the fossil sequences, nor was it found in a modern muskox sequence (OMModC1) that was acquired specifically to test results from this study (table 6, position 168). Because all of the sequences presented here are consensus results, generated from multiple PCRs and multiple clones sequenced on both strands, the difference in the database for modern samples is best interpreted as a sequencing artifact in the published sequences.

Table 6 Summary of variable positions in the hypervariable region of modern and ancient muskoxen

OMTai23658, one of two Holocene samples from Taimyr, yielded two unique haplotypes for the hypervariable region, identified here as A and B (table 6, position 144). Multiple PCRs, cloning, and sequencing of clones by two independent laboratories established that these haplotypes differ by 1 base, the difference likely being due to heteroplasmy (additional fig. 1; see also [26]). Haplotype A differed by 0–2 positions and haplotype B by 1–3 positions from sequences recovered for other fossil muskoxen from the Canadian Arctic Archipelago. The difference between both haplotypes of OMTai23658 and modern muskoxen was 1–3 substitutions.

All of the other successfully sequenced samples from Taimyr (OMTai14, OMTai38, OMTai39, and OMTai46) are of Pleistocene age (> 10,000 years ago). They were found to differ by 0–3 substitutions from one another but by 2–6 positions from modern muskoxen and the Canadian fossil samples. In addition, all the Pleistocene Taimyr samples presented 1 fixed difference not observed in any of the other samples studied (table 6, position 104), while OMTai38 and OMTai46 shared 2 base changes not found in any other specimen group (table 6, positions 148 and 166). OMTai14 and OMTai39 shared one position not found in any other specimen group (table 6, position 2).

The DNA yield from one other sample, OMBolL15564 from the New Siberian Islands, was low. However, it yielded reproducible EH sequence for cyt b. This is consistent with the sample's origin (near Taimyr) and age (late Pleistocene). Unfortunately, the HV region yielded several different sequences per amplification, including sequences that were similar to SH sequences, EH sequences, and ones that diverge from all others encountered in the course of this study (additional fig. 4). Since results for this specimen were highly inconsistent, we suspect some form of sample contamination but cannot identify the likely source. This required that we remove it form further analysis.

As may be seen in table 4, the Pleistocene sequences had more transversions (2 of 5 mutations) and slightly higher haplotype and nucleotide diversity per site measured than did Holocene sequences. However, given the almost 9-fold higher representation of Holocene compared to Pleistocene sequences, the statistics should be treated with caution. Higher among-sequence diversity (i.e., 2 haplotypes from 4 sequences as opposed to 4 from 44) for Holocene samples supports the concept that genetic diversity in muskoxen has been restricted. This is also consistent with the results obtained for cyt b (table 4).

To further clarify among-sequence relationships, we depict the sequence network suggested by our data in figure 5: nodes reflect missing intermediates and correlate with the number of mutational steps between sequences. As expected, all haplotypes of the SH group are very closely related. Most of the extinct haplotypes recovered differ somewhat among themselves, and all Taimyr EHs differ more substantially from SHs than do SHs among themselves. The distance between the extinct haplotypes and the closest SH sample is two to three mutational steps. The network supports the distinctiveness of most of the HV haplotypes of Pleistocene age.

Figure 5
figure 5

Minimum spanning network for hypervariable region sequences. The network is derived from the alignment shown in table 6; though indicated, the indel was omitted as it is a likely sequencing artifact in the database. Identical sequences are grouped together; EH samples are indicated by shading. The mutational change and sample mean radiocarbon age ranges (in yrbp) are shown next to each branch.


In view of the small number of specimens analyzed, conclusions have to be tentative. The conclusion best supported by the data is that, compared to modern muskoxen, late Quaternary Ovibos populations were less homogeneous for the loci sampled in this investigation. In particular, one group of haplotypes (EHs) occurring in fossil specimens from northeastern Asia differs from haplotypes in existing databases (SHs). The actual number of differences is not great (see Results), but their existence is thrown into sharp relief by the virtual lack of variation in the modern wildherd. In view of the multiple safeguards utilized to insure accuracy and repeatability (extraction of several samples in independent laboratories, sequencing of multiple independent clones from two or more amplifications, using overlapping PCR fragments to determine sequences, and rejection of suspect or ambiguous results), we feel that we can exclude laboratory artifacts or DNA damage as the source of these differences.

When the extinct haplotypes might have disappeared is not known. In our sample, the latest accepted record (of 5) is 18,310 ± 70 yrbp (radiocarbon years before present), which tells us that populations in which these haplotypes occurred survived at least into the Late Glacial Maximum, but not whether they lasted until the period of intensified megafaunal extinction at the end of the Pleistocene approximately 10,000 years ago (although in view of the extensive end-Pleistocene radiocarbon record for Ovibos in Asia [32], we regard this as more likely than not).

Another conclusion that is supported by the limited evidence currently available is that the SH group was present in mid-Holocene Asian muskoxen. This inference is based on the fact that OMTai23658 from Taimyr yielded sequence data indistinguishable from modern. (We had less success with the only other Eurasian Holocene sample available to us, OMTai31670. For the sequences covered it was initially found to be identical to modern muskox, but repeated efforts to replicate results were unsuccessful.) Given available data, however, it is premature to conclude that EHs were definitely not present as late as the middle Holocene. If in the future other specimens of Holocene age come to light in Asia, it would be well worthwhile to sample them for aDNA studies.

A last point of interest is that Ovibos, one of the few high-latitude megafaunal mammals to have survived into recent times, has clearly done so with reduced genetic variability. As already noted, at what point before the present this variability was lost cannot be satisfactorily established with existing data. For species for which aDNA sequence information is substantially better, such as bears and bison, there is molecular evidence for repeated loss and turnover of haplotype diversity during the late Quaternary (cf. [19, 20, 27]). Although no aDNA investigation of the giant Irish deer, Megaloceros, has yet been published, Stuart et al. [28] have recently shown that the radiocarbon record for this taxon reveals dramatic changes in population distributions during the late Pleistocene, followed ultimately by total collapse in the mid-Holocene. Thus, the data presented here, in combination with studies performed on other megafaunal species, is consistent with a scenario in which genetic diversity in muskoxen was reduced prior to the Pleistocene/Holocene transition.

It would also be of great interest to know whether late Quaternary small mammals (few of which became extinct) were affected by population crashes at the same time as larger species were. Extensive work on the mitochondrial genetics and phylogeography of extant lemmings by Fedorov and coworkers [29, 30] suggests that, in these rodents at least, bottlenecking and catastrophic loss of genetic diversity did not occur in late glacial/postglacial time. The robustness of this inference could be tested with aDNA methods similar to the ones utilized here, as there are large samples of reasonably well-dated fossil lemming material in the paleontological collections of various museums (cf. [31]).


In summary, the combined genetic and chronometric evidence presented in this paper establishes that Ovibos moschatus was genetically more diverse in the late Pleistocene than it is today. Precisely when haplotypes not found in extant muskoxen were lost, and under what conditions, are matters that cannot be well constrained at present. However, the data decisively show that now-extinct haplotypes persisted in Asia to a point only a few thousand years earlier than the accepted onset of the end-Pleistocene extinctions. It is reasonable to predict that, with better sampling, it may eventually be possible to show that populations of Ovibos moschatus, far from passing through this period of heightened extinction unscathed, were in fact brought down to small numbers over much of their original range, perhaps almost to the vanishing point [32]. Why muskoxen should have survived nevertheless in these circumstances, while many other high-latitude species did not, is yet another puzzle within the overarching mystery of late Quaternary extinctions.



Despite their wide holarctic distribution during the late Quaternary, ovibovines are comparatively rare as fossils except in those few areas (e.g., Taimyr Peninsula, parts of Yukon; fig. 1) where conditions were evidently optimal for them during much of the late Pleistocene. We obtained our specimens from two broad areas, designated "New World" and "Old World" (respectively, northern North America including the Arctic Archipelago, and the Arctic periphery of Asia east of the Kara Sea). A total of 31 muskoxen samples were subjected to DNA extraction and PCR amplification using various primer combinations. Of these, 16 samples worked with varying degrees of success; table 1 and fig. 3 provide information on these specimens regarding locality, age, and other relevant data (see also additional figs. 1, 2, 3, 4). The 15 remaining samples, from Alaska, Russia, and the Canadian Arctic Archipelago, did not produce amplifiable DNA on extraction and are not shown. It is particularly unfortunate that none of the Alaskan material yielded DNA. There may be an age effect: Mateus et al. [33] showed that muskoxen were common in Alaska earlier in the Pleistocene, but became very rare after 40,000 yrbp (which is the upper practical limit of radiocarbon dating). Designations combine taxon acronym/geographical origin/sample number: thus OMArA05 stands for Ovibos moschatus/Arctic Archipelago/sample 5.

Radiocarbon dating record

All specimens listed in table 1 had either been subsampled for 14C dating as part of this study or had been dated previously by other workers. Age estimates are cited in radiocarbon years before present without δ13C correction. All Beta Analytic estimates are accelerator dates on the collagen fraction of submitted specimens and typically have low associated standard errors. Other dates quoted in the table are based on older liquid scintillation technology which required large samples. However, accuracy and precision are not necessarily poorer in the latter case: OMTai23658 has been dated by both methods, and mean age estimates are statistically indistinguishable at 2 sigmas (table 1).

All samples yielded 13C/12C ratios within normal ranges for the materials dated (usually bone), and no difficulties with any samples were reported by the radiocarbon laboratory. Radiocarbon age estimates in excess of 40,000 yrbp are generally considered to be of questionable accuracy, and specimens so dated could be considerably older [37]. Although EH samples are generally older than SH samples, there is overlap around the time of the Last Glacial Maximum 21,000 calendar years ago [38], indicating that haplotype groups were coeval. Greater temporal overlap would probably have been found had sample sizes been larger.

Fig. 2 collates published finite 14C records for the study areas in the Old and New Worlds. Note that the distribution of dates in the New World record is quite asymmetrical, with few early dates and a heavy concentration in the late Holocene. A similar asymmetry is seen in fig. 3, which depicts the date distribution of specimens used in this study. The explanation concerns substantial differences in the amount of habitable range available to muskoxen during the period covered by radiocarbon dating. Because almost all of Canada and Greenland were covered by icecaps in the last (late Wisconsinan) cold phase, few early muskox dates would be expected from this area. Following icecap retreat, the muskoxen record increases markedly, especially after the mid-Holocene with return of a significant biomass of vascular plants at high latitudes [39]. By contrast, the Arctic periphery of Asia was not significantly glaciated, and muskox were evidently present throughout this zone during the latest Pleistocene, but rare or absent thereafter until ~4,000 yrbp when they made a brief reapperance.

DNA extraction, PCR, and sequencing

All fossil samples were collected by the first author using an electric drill with stainless steel "keyhole" coring drill bits that were cleaned and air dried after each use. Samples were placed in individual plastic containers, labeled, and stored while in the field in cool conditions and on return in a -20°C freezer in an isolated lab. Protocols for DNA extractions, PCR conditions, PCR reamplification, cloning strategy, colony PCR and sequencing, have been described in detail in Krings et al. [12] and Greenwood [34], with the variant that (for most samples) DNA extractions were conducted with NucleoSpin® DNA-Trace Kits in combination with NucleoSpin® funnel columns (Macherey & Nagel), Gene Clean Kit (Qbiogene), or GOLD Forensic DNA Kit (peqlab) according to manufacturer protocols and eluted in 55 to 60μl 10 mM Tris-HCl, pH 8.5. Reamplifications of PCR products when necessary were done either by direct amplification from of an aliquot of the first PCR or subsequent to gel extraction of the primary amplification band using a Qiagen Gel Extraction Kit. A buffer control was taken through each extraction procedure with no sample added and used as a PCR negative control in each PCR reaction to control for contamination. PCR primers, combinations used, and sequences are shown in table 2. Modern hair and blood DNA were prepared in a separate laboratory and never brought to laboratories conducting aDNA work so as to avoid potential contamination. Both hair and blood were prepared using Qiagen QIAamp DNA Mini Kits and PCR performed using Roche Expand Taq or Promega Taq polymerase according to manufacturer instructions. PCR conditions used were 30 cycles of 94°C 30 seconds, 50°C 30 seconds and 72°C 30 seconds. Products were cloned and clones sequenced as described in Greenwood (2002).

Table 2 Primers used in reconstructing cyt b and HV sequences

Negative water controls and in many cases mock extraction controls were included in experiments to control for contamination. Extractions, amplification, and sequencing for samples OMTai23658, OMTai39, OMTai46, OMYuk01, OMArA03, OMArA04, and OMArA05 were performed independently at SR&D and Medigenomix; in each case the same sequences were obtained. Each amplified region was covered at least twice if possible; multiple clones were sequenced to determine the correct sequence and to avoid mis-scoring errors due to DNA damage [35]. In the few cases in which available material was insufficient for double coverage, the sequence result is regarded as tentative. The number of replications for each PCR for each sample is indicated in table 1. Sequences were also derived from overlapping PCR fragments as previously described [12] (see additional fig. 1, 2, 3). [All consensus sequences have been deposited with the following accession numbers in GenBank: AY839538-AY839551 (cyt b); AY839552-AY839563 (hypervariable region); and AY83956 (the only substantially complete cyt b sequence retrieved in this study for OMTai39).]

Phylogenetic analysis

Maximum parsimony, neighbor joining, and maximum likelihood analyses were executed using PAUP 4.0b10 [36] with Ovis (sheep) and Capra (goat) as outgroups for the nearly-complete cyt b sequence of OMTai39 (fig. 4). Although 100 bootstrap replicates were performed for each tree generated, the limited variability present in the data requires that no great weight be placed on the results of such experiments (i.e., chances are high that the few variable sites will be included in the bootstrap and will therefore nearly always replicate the same tree). Nevertheless, a strictly parsimonious interpretation of the distribution of certain fixed differences (see text) indicates that sequence data can be used to diagnostically distinguish extinct from surviving haplotypes. Networks for the hypervariable region sequences derived from the alignment shown in table 6 were drawn by hand (fig. 5).


  1. Martin PS, Steadman DW: Prehistoric extinctions on islands and continents. Extinctions in Near Time: Causes, Contexts, and Consequences. Edited by: MacPhee RDE. 1999, New York, Kluwer, 17-56.

    Chapter  Google Scholar 

  2. Vartanyan SL, Arslanov KA, Tertychnaya TV, Chernov SB: Radiocarbon dating evidence for mammoths on Wrangel Island, Arctic Ocean, until 2000 BC. Radiocarbon. 1995, 37: 1-6.

    Google Scholar 

  3. Guthrie RD: Radiocarbon evidence of mid-Holocene mammoths stranded on an Alaskan Bering Sea island. Nature. 2004, 429: 746-749. 10.1038/nature02612.

    Article  PubMed  Google Scholar 

  4. McDonald JN, Ray CE: The autochthonous North American musk oxen Bootherium, Symbos, and Gidleya (Mammalia: Artiodactyla: Bovidae). Smithsonian Contributions to Paleobiology. 1989, 66: 1-77.

    Article  Google Scholar 

  5. McDonald JN, Ray CE, Harington CR: Taxonomy and zoogeography of the musk ox genus Praeovibos Staudinger, 1908. Illinois State Museum Scientific Papers. 1991, 23: 285-314.

    Google Scholar 

  6. Lent PC: Muskoxen and their Hunters. 1999, Norman, University of Oklahoma Press

    Google Scholar 

  7. Groves P: Intraspecific variation of mitochondrial DNA of muskoxen based on control region sequences. Can J Zool. 1997, 75: 568-575.

    Article  CAS  Google Scholar 

  8. Groves P, Shields GF: Cytochrome B sequences suggest convergent evolution of the Asian takin and Arctic muskox. Mol Phylogenet Evol. 1997, 8: 363-374. 10.1006/mpev.1997.0423.

    Article  CAS  PubMed  Google Scholar 

  9. Holm LE, Forchhammer MC, Boomsma JJ: Low genetic variation in muskoxen (Ovibos moschatus) from western Greenland using microsatellites. Mol Ecol. 1999, 8: 675-679. 10.1046/j.1365-294x.1999.00615.x.

    Article  CAS  PubMed  Google Scholar 

  10. Van Coeverden PJ, Boag P: Optimization of novel polymorphic microsatellites in muskox (Ovibos moschatus) leads to an increased estimate of microsatellite diversity. Mol Ecol Notes. 2004, 4: 713-715. 10.1111/j.1471-8286.2004.00768.x.

    Article  Google Scholar 

  11. Mikko S, Roed K, Schmutz S, Anderson L: Monomorphism and polymorphism at the Mhc DRB loci in domestic and wild ruminants. Immunol Rev. 1999, 167: 169-178.

    Article  CAS  PubMed  Google Scholar 

  12. Krings M, Stone A, Schmitz RW, Krainitzki H, Stoneking M, Pääbo S: Neandertal DNA sequences and the origin of modern humans. Cell. 1997, 90: 19-30. 10.1016/S0092-8674(00)80310-4.

    Article  CAS  PubMed  Google Scholar 

  13. Hadly EA, Kohn MH, Leonard JA, Wayne RK: A genetic record of population isolation in pocket gophers during Holocene climatic change. Proc Natl Acad Sci U S A. 1998, 95: 6893-6896. 10.1073/pnas.95.12.6893.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  14. Greenwood AD, Castresana J, Feldmaier-Fuchs G, Pääbo S: A molecular phylogeny of two extinct sloths. Mol Phylogenet Evol. 2001, 18: 94-103. 10.1006/mpev.2000.0860.

    Article  CAS  PubMed  Google Scholar 

  15. Greenwood AD, Lee F, Capelli C, DeSalle R, Tikhonov AN, Marx PA, MacPhee RDE: Evolution of endogenous retrovirus-like elements of the woolly mammoth (Mammuthus primigenius) and its relatives. Mol Biol Evol. 2001, 18: 840-847.

    Article  CAS  PubMed  Google Scholar 

  16. Lambert DM, Ritchie PA, Millar CD, Holland B, Drummond AJ, Baroni C: Rates of evolution in ancient DNA from Adélie penguins. Science. 2002, 295: 2270-2273. 10.1126/science.1068105.

    Article  CAS  PubMed  Google Scholar 

  17. Orlando L, Leonard JA, Thenot A, Laudet V, Guerin C, Hänni C: Ancient DNA analysis reveals woolly rhino evolutionary relationships. Mol Phylogenet Evol. 2003, 28: 485-499. 10.1016/S1055-7903(03)00023-X.

    Article  CAS  PubMed  Google Scholar 

  18. Matisoo-Smith E, Robins JH: Origins and dispersals of Pacific peoples: evidence from mtDNA phylogenies of the Pacific rat. Proc Natl Acad Sci U S A. 2004, 101: 9167-9172. 10.1073/pnas.0403120101.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  19. Shapiro B, Drummond A, Rambaut A, Wilson MC, Matheus PE, Sher AV, Pybus OG, Gilbert MT, Barnes I, Binladen J, Willerslev E, Hansen AJ, Baryshnikov GF, Burns JA, Davydov S, Driver JC, Froese DG, Harington CR, Keddie G, Kosintsev P, Kunz ML, Martin LD, Stephenson RO, Storer J, Tedford R, Zimov S, Cooper A: Rise and fall of the Beringian steppe bison. Science. 2004, 306: 1561-1565. 10.1126/science.1101074.

    Article  CAS  PubMed  Google Scholar 

  20. Leonard JA, Wayne RK, Cooper A: Population genetics of ice age brown bears. Proc Natl Acad Sci U S A. 2000, 97: 1651-1654. 10.1073/pnas.040453097.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  21. Leonard JA, Vilà C, Wayne RK: Legacy lost: genetic variability and population size of extirpated US grey wolves (Canis lupus). Mol Ecol. 2004, 14: 9-17. 10.1111/j.1365-294X.2004.02389.x.

    Article  Google Scholar 

  22. Hofreiter M, Serre D, Rohland N, Rabeder G, Nagel D, Conard N, Munzel S, Pääbo S: Lack of phylogeography in European mammals before the last glaciation. Proc Natl Acad Sci U S A. 2004, 101: 12963-12968. 10.1073/pnas.0403618101.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  23. Greenwood AD, Capelli C, Possnert G, Pääbo S: Nuclear DNA sequences from late Pleistocene megafauna. Mol Biol Evol. 1999, 16: 1466-1473.

    Article  CAS  PubMed  Google Scholar 

  24. Greenwood AD, Pääbo S: Nuclear insertion sequences of mitochondrial DNA predominate in hair but not in blood of elephants. Mol Ecol. 1999, 8: 133-137. 10.1046/j.1365-294X.1999.00507.x.

    Article  CAS  PubMed  Google Scholar 

  25. Thalmann O, Hebler J, Poinar HN, Pääbo S, Vigilant L: Unreliable mtDNA data due to nuclear insertions: a cautionary tale from analysis of humans and other great apes. Mol Ecol. 2004, 13: 321-335. 10.1046/j.1365-294X.2003.02070.x.

    Article  CAS  PubMed  Google Scholar 

  26. Ivanov PL, Wadhams MJ, Roby RK, Holland MM, Weedn VW, Parsons TJ: Mitochondrial DNA sequence heteroplasmy in the Grand Duke of Russia Georgij Romanov establishes the authenticity of the remains of Tsar Nicholas II. Nat Genet. 1996, 12: 417-420. 10.1038/ng0496-417.

    Article  CAS  PubMed  Google Scholar 

  27. Orlando L, Bonjean D, Bocherens H, Thenot A, Argant A, Otte M, Hänni C: Ancient DNA and the population genetics of cave bears (Ursus spelaeus) through space and time. Mol Biol Evol. 2002, 19: 1920-1933.

    Article  CAS  PubMed  Google Scholar 

  28. Stuart AJ, Kosintsev PA, Higham TFG, Lister AM: Pleistocene to Holocene extinction dynamics in giant deer and woolly mammoth. Nature. 2004, 431: 684-689. 10.1038/nature02890.

    Article  CAS  PubMed  Google Scholar 

  29. Fedorov VB, Goropashnaya AV, Jarrell GH, Fredga K: Phylogeographic structure and mitochondrial DNA variation in true lemmings (Lemmus) from the Eurasian Arctic. Biological Journal of the Linnaean Society. 1999, 66: 357-371. 10.1006/bijl.1998.0271.

    Article  Google Scholar 

  30. Fedorov VB, Goropashnaya AV, Jaarola MA, Cook JA: Phylogeography of lemmings (Lemmus): no evidence for postglacial colonization of Arctic from the Beringian refugium. Mol Ecol. 2003, 12: 725-731. 10.1046/j.1365-294X.2003.01776.x.

    Article  CAS  PubMed  Google Scholar 

  31. Smirnov NG, Fedorov VB: Holarctic collared lemmings: traces of their spread as related to the history of the arctic biota. Russian Journal of Ecology. 2003, 34: 370-376. 10.1023/A:1027304332042.

    Article  Google Scholar 

  32. MacPhee RDE, Tikhonov AN, Mol D, de Marliave C, van der Plicht H, Greenwood AD, Flemming C, Agenbroad L: Radiocarbon chronologies and extinction dynamics of the late Quaternary megafauna of the Taimyr Peninsula, Russian Federation. Journal of Archaeological Science. 2002, 29: 1017-1042. 10.1006/jasc.2001.0802.

    Article  Google Scholar 

  33. Matheus P, Kunz M, Guthrie RD: Using frequency distributions of radiocarbon dates to detect relative changes in Pleistocene mammal populations – a test case from northern Alaska. Proceedings of the 3rd International Mammoth Conference, May 24–29. 2003, [] , Dawson City, Yukon, Abstracts and Program

    Google Scholar 

  34. Greenwood AD: Late Pleistocene DNA extraction and analysis. Techniques in Molecular Systematics and Evolution. Edited by: DeSalle R, Giribet G, Wheeler W. 2002, Basel, Birkhauser-Verlag, 288-301.

    Chapter  Google Scholar 

  35. Hofreiter M, Jaenicke V, Serre D, Haeseler A, Pääbo S: DNA sequences from multiple amplifications reveal artifacts induced by cytosine deamination in ancient DNA. Nucleic Acids Res. 2001, 29: 4793-4799. 10.1093/nar/29.23.4793.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Swofford DL: PAUP*: Phylogenetic Analysis using Parsimony (* and Other Methods); Version 4. 2002, Sunderland, MA, Sinauer

    Google Scholar 

  37. Geyh M, Schleicher H: Absolute Age Determination. 1991, New York, Springer-Verlag

    Google Scholar 

  38. Siegert M: Ice Sheets and Late Quaternary Environmental Change. 2001, New York, Wiley

    Google Scholar 

  39. Pielou EC: After the Ice Age: The Return of Life to Glaciated North America. 1991, Chicago, University of Chicago Press

    Book  Google Scholar 

  40. Harington CR: Annotated Bibliography of Quaternary Vertebrates of Northern North America, with Radiocarbon Dates. 2003, Toronto, University of Toronto Press

    Google Scholar 

  41. Sher AV, Kuzmina SA, Kuznetsova TV, Sulerzhitsky LD: New insights into the Weichselian environment and climate of the East Siberian Arctic, derived from fossil insects, plants, and mammals. Quaternary Science Reviews. 2005, 24: 533-569. 10.1016/j.quascirev.2004.09.007.

    Article  Google Scholar 

  42. Dyke AS: An outline of North American deglaciation with emphasis on central and northern Canada. Quaternary Glaciations: Extent and Chronology, Part 2: North America. Edited by: Ehlers J, Gibbard PL. 2004, Amsterdam, Elsevier, 373-424.

    Chapter  Google Scholar 

  43. Vereshchagin NK, Baryshnikov GF: Quaternary mammalian extinctions in northern Eurasia. Quaternary Extinctions: a Prehistoric Revolution. Edited by: Martin PS, Klein RG. 1984, Tucson, University of Arizona Press, 483-516.

    Google Scholar 

  44. Harington CR: A postglacial muskox (Ovibos moschatus) from Grandview, Manitoba, and comments on the zoogeography of Ovibos. Canadian National Museum of Natural Sciences Publications in Paleontology. 1970, 2: 1-13.

    Google Scholar 

Download references


The authors gratefully acknowledge financial support by grants from the National Science Foundation OPP 0117400 (RDEM, ADG), the Niarchos Fund (RDEM), and the Evelyn Stefansson Nef Foundation (RDEM). The authors wish to thank Lars Giesen (Medigenomix GmbH, Germany), Simone Watzinger and Jörg Hauf (Scientific Research and Development GmbH, Germany) for excellent technical assistance. For access to specimens and other favors, we thank Richard Harington, Alison Murray, Kieran Shepherd (all of the Canadian Museum of Nature, Ottawa), Meng Jin (American Museum of Natural History, New York), and Bernard Buigues (Cerpolex/Mammuthus Project, Khatanga). We are grateful to Mike Bunce, Eske Willerslev, Hendrik Poinar, and especially Claudia Englbrecht for helpful discussions. Patricia Wynne executed the base map used for figure 1 and Laura Carter helped with the editing of the final manuscript.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Ross DE MacPhee.

Additional information

Authors' contributions

RDEM participated in the design of the study, coordinated the drafting of the paper with the other authors, acquired and interpreted radiocarbon dates, and helped with phylogenetic analysis; ANT participated in the drafting of the paper and provided data and information from Russian sources; DM managed acquisition of specimens for sequence analysis in Russia and Canada, and participated in the drafting of the paper; and ADG participated in the design of the study, oversaw all lab work, sequence alignment, and phylogenetic analysis, and participated in drafting the paper. All authors read and approved the final manuscript.

Ross DE MacPhee and Alex D Greenwood contributed equally to this work.

Electronic supplementary material


Additional File 1: Alignment of hypervariable region clones from overlapping PCR fragments used to determine the sequences for each specimen analyzed. The consensus sequence of OMModCl was used as a reference. Dots indicate identity to the consensus sequence. Differences are indicated. Base pair sequences that could not be determined for a given clone are marked "N". Gaps are indicated by "-". The first clone provides the name of each sample according to the labels provided in table 1. The numbers indicate PCR replication + clone number: thus 1.1 indicates PCR 1, first clone sequence. (PDF 163 KB)


Additional File 2: Alignment of cyt b clones from overlapping PCR fragments used to determine the sequences for each specimen analyzed. Symbols as in additional fig. 1. (PDF 160 KB)


Additional File 3: Alignment of clones from overlapping PCR fragments used to determine 1,180 bp of OMTai39 cyt b. Symbols as in additional fig. 1. The reference is a modern muskox sequence taken from GenBank. Four additional modern muskox sequences are included in the alignment. (PDF 164 KB)


Additional File 4: Clone sequences derived from multiple amplifications of sample OMBolL15564. All symbols are the same as in additional fig. 1, 2, 3. The modern DNA consensus sequence from additional fig. 1 is used as a reference sequence. (PDF 158 KB)

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

MacPhee, R.D., Tikhonov, A.N., Mol, D. et al. Late Quaternary loss of genetic diversity in muskox (Ovibos). BMC Evol Biol 5, 49 (2005).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: