This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

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.

Results

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.

Conclusion

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.

Background

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.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.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
[4-6]. 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 DRB3 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., [12-22]]. 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.

Results

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.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).

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.

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.

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.

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.

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.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.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.

Discussion

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]).

Conclusion

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.

Methods

Samples

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).

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).

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.

Acknowledgements

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.

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