Scholarly article on topic 'Holarctic genetic structure and range dynamics in the woolly mammoth'

Holarctic genetic structure and range dynamics in the woolly mammoth Academic research paper on "Biological sciences"

Share paper

Academic research paper on topic "Holarctic genetic structure and range dynamics in the woolly mammoth"




^^^ CrossMark

Cite this article: Palkopoulou E, Dalen L, Lister AM, Vartanyan S, Sablin M, Sher A, Edmark VN, Brandstrom MD, Germonpré M, Barnes I, Thomas JA. 2013 Holarctic genetic structure and range dynamics in the woolly mammoth. Proc R Soc B 280: 20131910.

Received: 23 July 2013 Accepted: 16 August 2013

Subject Areas:

evolution, ecology, molecular biology Keywords:

ancient DNA, Mammuthus primigenius, extinction, refugia, climate

Holarctic genetic structure and range dynamics in the woolly mammoth

Eleftheria Palkopoulou1'2, Love Dalén19, Adrian M. Lister3, Sergey Vartanyan4, Mikhail Sablin5, Andrei Sher6+, Veronica Nystrom Edmark1, Mikael D. Brandstrom7, Mietje Germonpré8, Ian Barnes9 and Jessica A. Thomas9,*

department of Bioinformatics and Genetics, Swedish Museum of Natural History, 10405 Stockholm, Sweden

2Department of Zoology, Stockholm University, Stockholm 10691, Sweden

3Department of Earth Sciences, Natural History Museum, London SW7 5BD, UK

4Northeast Interdisciplinary Research Institute, Far East Branch, Russian Academy of Sciences, Magadan 685000, Russia

5Zoological Institute of Russian Academy of Sciences, Saint-Petersburg 199034, Russia

6Institute of Ecology and Evolution, Russian Academy of Sciences, Moscow 119071, Russia

7Department of Forest Mycology and Pathology, Swedish University of Agricultural Sciences, 10691 Uppsala, Sweden

8Operational Direction 'Earth and History of Life', Royal Belgian Institute of Natural Sciences, Vautierstraat 29, 1000 Brussels, Belgium

9School of Biological Sciences, Royal Holloway University of London, Egham, Surrey TW20 0EX, UK

Ancient DNA analyses have provided enhanced resolution of population histories in many Pleistocene taxa. However, most studies are spatially restricted, making inference of species-level biogeographic histories difficult. Here, we analyse mitochondrial DNA (mtDNA) variation in the woolly mammoth from across its Holarctic range to reconstruct its history over the last 200 thousand years (kyr). We identify a previously undocumented major mtDNA lineage in Europe, which was replaced by another major mtDNA lineage 32-34 kyr before present (BP). Coalescent simulations provide support for demographic expansions at approximately 121 kyr BP, suggesting that the previous interglacial was an important driver for demography and intra-specific genetic divergence. Furthermore, our results suggest an expansion into Eurasia from America around 66 kyr BP, coinciding with the first exposure of the Bering Land Bridge during the Late Pleistocene. Bayesian inference indicates Late Pleistocene demographic stability until 20-15 kyr BP, when a severe population size decline occurred.

Author for correspondence:

Eleftheria Palkopoulou


deceased 11 August 2008. ^Present address: Department of Biology, University of York, York Y010 5DD, UK.

Electronic supplementary material is available at or via

1. Introduction

One of the greatest strengths of an ancient DNA (aDNA) approach is that it enables the study of genetic change through time. Analyses of samples across particular geographical regions through time have revealed unexpected patterns of local population extinction and recolonization [1-6]. However, while such studies are invaluable for investigating the interaction between population dynamics and local changes in the environment, it is not always clear how different lineages evolved and where recolonizing populations originated. A more comprehensive approach, encompassing the full geographical extent of a species' distribution, is needed to fully understand its biogeographic history [7,8]. Furthermore, a large sample size is important to provide good spatio-temporal coverage and sufficient detail for the reconstruction of evolutionary events.

The woolly mammoth (Mammuthus primigenius) is one of the best-studied taxa in the field of aDNA. Adapted to the cold and arid steppe-tundra, mammoths were widespread during the Late Pleistocene (ca 116-12 kyr BP), with a range that extended from Western Europe to the northern part of North America [9]. However, the availability of well-preserved permafrost samples

Royal Society Publishing

Informing the science of the future

© 2013 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License, which permits unrestricted use, provided the original author and source are credited.

in Siberia and Alaska has meant that genetic studies have predominantly focused on the region of Beringia [6,10-13]. By contrast, little is known about the population structure at the western end of mammoth distribution in Europe, although aDNA extracted from a single European specimen indicated a high level of sequence divergence from other mammoth populations [6]. Furthermore, while morphological data revealed that woolly mammoths were present in Europe from around 200 kyr BP until the end of the Pleistocene [14,15], little is known about their population dynamics within Europe, nor about the extent of gene flow between European and Asian populations.

Previous genetic studies on Beringian mammoths identified two deeply divergent monophyletic mitochondrial DNA (mtDNA) lineages, clades I and II, hypothesized to have evolved in isolation on either side of the Bering Strait [6,10]. Clade I (haplogroups C, D and E in Debruyne et al. [11]) had a widespread distribution during the later stages of the Late Pleistocene, but appears to have originated in North America and dispersed into Eurasia during the Middle or early Late Pleistocene [6,11]. Clade II (haplogroup A in Debruyne et al. [11]) had a much more limited geographical distribution in eastern Siberia and is thought to have a Siberian origin [6,10]. It appears that these two genetic lineages coexisted in northeast Siberia for thousands of years before clade II disappeared at approximately 40 kyr BP [6,10]. Clade I, however, survived well into the Holocene. The last mainland populations of clade I mammoth persisted in areas of northern Siberia until ca 11 kyr BP [16] but, in contrast to most megafaunal species that went extinct around the Pleistocene/Holocene transition, small populations of woolly mammoth survived to the mid-Holocene, until ca 6 kyr on St Paul Island [17] and 4 kyr on Wrangel Island [18].

While several hypotheses have been proposed to explain the genetic changes that took place during the Late Pleistocene, a full picture is yet to emerge regarding the origin of different genetic lineages and the timing of changes in demography and genetic variation. This study has three aims. First, to establish the woolly mammoth's Late Pleistocene genetic structure across the whole of its Holarctic distribution. Second, to examine different hypotheses regarding the mammoth's evolutionary history, including levels of genetic diversity and the timing of local population turnover events, range expansions and contractions. Third, to improve the resolution of phylogenetic and demographic analyses, since previous studies on mammoth genetics have been shown to potentially lack significant molecular signal [19]. In this study, we extract and analyse mtDNA from woolly mammoth specimens from across the Holarctic, expanding the genetic sampling both spatially and temporally. We include specimens from Europe as well as from Siberia identified as Middle Pleistocene in age. Combined with previously published mtDNA sequences, the dataset comprises more than 300 mammoth specimens, thus enabling a thorough reconstruction of the species' population history from the Late Middle Pleistocene (LMP) up until its extinction.

mammoth (see electronic supplementary material, table S1). MtDNA amplification was performed as in Barnes et al. [6], targeting a 741 bp region, including the 3' end of the cytochrome b gene (CytB), two tRNA genes (tRNA-Thr and tRNA-Pro), and the first hypervariable part of the control region (CR1). Pre-PCR laboratory work was performed in dedicated aDNA laboratories at Royal Holloway, University of London and the Swedish Museum of Natural History in Stockholm, following standard protocols and procedures (for details, see the electronic supplementary material). Radiocarbon dating was performed at the Oxford Radiocarbon Accelerator Unit using accelerator mass spectrometry. Radiocarbon dates were calibrated in OxCal v. 4.1 [20] with the IntCal09 calibration curve [21] (see the electronic supplementary material, tables S1 and S2).

(b) Phylogenetic analyses

Our ancient mtDNA sequences were aligned with homologous woolly mammoth sequences available on GenBank [6,10-13, 22-25] (see the electronic supplementary material, table S2 for accession numbers) in GENEIOUS v. 5.0.1 [26]. We used Partition Finder [27] to select the best-fit partitioning scheme and DNA substitution model (see electronic supplementary material, table S3). Bayesian phylogenies were generated using MrBayes v. 3.2.1 [28], with two African elephant (Loxodonta cyclotis, Loxodonta africana) and one Asian elephant (Elephas maximus) sequences as outgroups (accession nos: AY359274, NC000934 and EF588275). Four chains were run for 20 million generations and sampled every 2000. BEAST v. 1.7 [29] was used to construct dated phylo-genies and to estimate the dates of three specimens that were stratigraphically dated to the LMP, using the tip-dating method. For the dating analyses, the temporal signal of the data was first assessed using the date randomization test [19]. BEAST analyses were run under three different population models, for 200 million generations and sampled every 20000. Convergence for both phylogenetic analyses was assessed in Tracer [30]. A median-joining network was constructed using Phylonet5 v. 1.0.0 (A. Helgason 2009, unpublished). Changes in phylogeographic patterns over time were visualized in GenGIS [31].

(c) Coalescent simulations and approximate Bayesian computation

Serial coalescent simulations were run with Bayesian Serial SimCoal [32] and analysed in an approximate Bayesian computation (ABC) framework [33] using the ABC package in R [34]. Summary statistics (see electronic supplementary material, table S4) were calculated in Arlequin v. 3.5 [35]. Using the age of first reproduction as a proxy for generation time [36], we assumed a generation time of 15 years [37]. The mean mutation rate (9.56% per site/106 years), transition bias (0.98) and shape parameter of gamma distribution (0.107) estimated from BEAST were used in all simulations. One million iterations were run for each scenario. Two regression methods were employed, local linear regression [38] and the neural networks algorithm [39] with 1% acceptance ratio. The latter method is recommended when highly dimensional summary statistics are used in order to transform the number of possibly correlated variables into a smaller number of variables [40]. Additional simulations with a generation time of 20 years were run to assess the effect of generation time on the outcome of the analysis.

2. Material and methods

(a) DNA analysis and radiocarbon dating

We recovered DNA from specimens of bone, tooth and tusk (n = 88) collected from most of the Holarctic range of the woolly

(i) Models 1A and 1B

Two scenarios were simulated to examine whether a history of population isolation in multiple refugia during the previous interglacial could explain the observation of different mtDNA clades: (i) three populations (representing the three mtDNA clades) that split from each other with constant effective population size

(Nef) and a uniform prior for the split time versus (ii) three populations that split from each other (as above) and experienced a simultaneous bottleneck, followed by exponential growth. Nef, split time and start time of expansions were sampled from uniform priors (for details, see the electronic supplementary material, figure S1).

(ii) Models 2A and 2B

We estimated the split time between the North American and Eurasian populations carrying mtDNA clade I, which could represent the time when North American individuals crossed over the Bering Land Bridge, leading to the introduction and expansion of clade I woolly mammoths in Eurasia [6,11]. In addition to the hypothesized dispersal time, subsequent gene flow in the opposite direction was assessed. Two alternative models were thus evaluated with two populations (representing North America and Eurasia) that diverged from each other: (i) without gene flow after the divergence event or (ii) with subsequent gene flow from Eurasia to North America. Uniform priors were used for Nef, split time and migration rate (for details, see the electronic supplementary material, figure S2).

3. Results

(a) Mitochondrial DNA diversity

The complete 741 bp mtDNA sequence was successfully amplified for 56 out of 88 woolly mammoth specimens. Owing to low DNA preservation, the remaining specimens yielded partial sequences. Of these, only 16 could be sequenced for a short 79 bp fragment that contains polymorphic sites informative for clade identification (see electronic supplementary material, table S1). Twenty-nine novel haplotypes were identified (figure 1a and electronic supplementary material, figure S3). These sequences together with previously published homologous mtDNA sequences (see electronic supplementary material, table S2) comprised a total dataset of 320 sequences.

(b) Genetic structure and demographic change

Bayesian phylogenetic analyses using MRBAYES and BEAST produced very similar topologies (figure 1b and electronic supplementary material, figures S4 and S5). Three major mtDNA clades were identified, each supported by moderate to high posterior probability. Of the three clades, clades I and II have been described earlier [6,10], while we here show that clade III is a lineage containing specimens from Europe. Three specimens from northern Yakutia (northeastern Siberia) that were stratigraphically dated to the LMP were found to be basal to the three clades (figure 1b and electronic supplementary material, figures S4 and S5). The date randomization test indicated that the dataset contained sufficient temporal signal for meaningful estimation of molecular rates and dating [19]. Tip-dating in BEAST produced molecular dates for these individuals between ca 175 and 200 kyr BP (see electronic supplementary material, table S5).

Coalescent simulations indicated that the model including a bottleneck followed by exponential growth for the three populations fits the observed data better than the model without growth (BF = 21-70 and electronic supplementary material table S6). Analysis of posterior distributions with the neural networks regression algorithm suggested that the three populations split at ca 196 kyr BP (95% credibility interval (CI): 156-261 kyr BP) and went

through a contemporaneous demographic expansion at ca 121 kyr BP (95% CI: 80-148 kyr BP; figure 2 and electronic supplementary material, table S7). Although the 95% CI of the expansion time is relatively wide, its median is close to the end of the last interglacial period (Marine Isotope Stage (MIS) 5e), suggesting population isolation in interglacial refugia as a likely explanation for the observed divergent mitochondrial lineages. Analysis of additional simulations with a generation time of 20 years (instead of 15 years) produced very similar posterior distributions for split and expansion time ( for details, see the electronic supplementary material, table S8).

(c) Range dynamics of clade I woolly mammoths

Most North American sequences clustered within a distinct sub-clade of clade I (haplogroup C in figure 1 and electronic supplementary material, figures S3-S5). Earlier studies had suggested a North American origin of this clade [6,11]. According to this hypothesis, the star-like patterns in the main subclade of clade I (haplogroups D and E in figure 1a) could suggest that founder effects took place in Eurasia after the dispersal of North American clade I mammoths across Beringia. We estimated the split time between the North American and Eurasian populations carrying mtDNA clade I to ca 66 kyr BP (95% CI: 4996 kyr BP; figure 2c; electronic supplementary material, table S7), which could correspond to the time when North American mammoths dispersed into Eurasia. It is notable that the two oldest dated clade I specimens found in eastern Siberia, with radiocarbon estimates of greater than 60 kyr BP (hap 2 and hap 5 in electronic supplementary material, figure S3 and see tables S1 and S2 for details), carry haplogroups identical to, and two mutational steps from, the modal haplotype in hap-logroup D, which is what would be expected if this modal haplotype represents the founding lineage of the expansion that took place after the dispersal event.

The presence of North American sequences outside haplo-group C, in the main subclade of clade I (figure 1) could be indicative of post-colonization gene flow in an eastwards direction [11]. The simulations, however, provided higher support (BF > 200; electronic supplementary material, table S6) for the model without gene flow from the Eurasian population to the North American population after their split (model 2A in electronic supplementary material, figure S2). Thus, an alternative explanation could be that these lineages represented polymorphisms that were already present in the ancestral population in North America.

(d) Genetic turnover events

Following the expansion from North America to Eurasia, clade I appears to have been sympatric with clade II in Central and East Siberia until the demise of the latter (figure 3; see also [6,10]). Based on finite radiocarbon dates from clade II specimens, it seems that this clade disappeared around 45 kyr BP (although it should be noted that one clade II specimen has yielded an infinite date of more than 33 000 14C years, roughly corresponding to more than 37.5 kyr BP). Furthermore, clade I woolly mammoths continued expanding to the west, all the way to Europe, where they replaced the endemic population carrying mtDNA clade III ( figure 3). Clade III disappeared from the fossil record at ca 34 kyr BP, whereas clade I seems to have made its first appearance in Europe at ca 32 kyr BP (see electronic supplementary material, table S1).

clade I

haplogroup E

□ Siberia

□ Europe

□ North America

□ Chukotka/Kamtchatka

□ Wrangel Island after isolation

clade I

clade III






Figure 1. Median-joining haplotype network and Bayesian phylogeny of woolly mammoth mtDNA sequences. (a) Haplotype colours indicate their geographical location. Shaded areas correspond to haplogroups as in Debruyne et al. [11]. Black dots represent missing haplotypes. Haplotype size is proportional to its frequency within the dataset except for the three most frequent haplotypes within haplogroups D and E that have frequencies above 10. (b) The labels at the tips of the phylogeny are coloured according to the geographical origin of the sequences. The timescale on the x-axis is in calendar years before present. Bayesian posterior probabilities of major internal nodes above 0.8 are shown.

g 0.00015

0.0003 -

75 150

time (kyr BP)

(b) 0.00020


40 60 80 100 time (kyr BP)

Figure 2. Posterior distributions for the time parameters of models IB, 2A and relative sea-level (RSL) variation in the Bering Strait. (a) Expansion time (T0) and (b) split time (T-|) of the three populations in model 1B. (c) Split time (T^ between the North American and Eurasian population carrying mtDNA clade I in model 2A. Time is given in calendar years before present. The thick red and blue lines show the posterior density curve, and the dotted lines show the prior distribution. (d) Figure redrawn from Hu et al. [41]. The thick black line represents the RSL and the horizontal grey line shows the current Bering Strait depth. RSL below the grey line indicates exposure of the Bering Strait. The blue vertical dashed line indicates the median of the posterior distribution of the split time (T1) between the North American and Eurasian population in model 2A.

(e) Bayesian analysis of changes in population size through time

The Bayesian skyline plot (BSP) showed a severe and sudden decline in Nef during MIS 2, starting at ca 20-15 kyr BP and ending early in MIS 1, at ca 10 kyr BP (figure 4). The observed reduction in Nef was at least 10-fold from ca 20 000 individuals (95% highest posterior density (HPD): 40 000-10 000) to ca 1000 individuals (95% HPD: 3000-400) assuming an average generation time of 15 years. Even though the HPD intervals are wide for the periods preceding and following the bottleneck, Nef seems to have remained relatively constant until this decline.

4. Discussion

The last interglacial, known in Europe as the Eemian and equated with MIS 5e (130-116 kyr BP), was a period characterized by climate at least as warm as today [42]. Tundra habitats in mainland Europe and Siberia were partly replaced by mixed or more open forests and sea levels were as high as, or even higher than that at present [42,43]. Species adapted to the cold, dry and open conditions of the steppe-tundra, such as the woolly mammoth, were likely affected by these environmental changes through range contractions and demographic declines [44]. Previous research employing back-casted climate models has proposed that habitat suitable for the woolly mammoth in Eurasia was severely reduced during the Eemian interglacial, restricting the mammoth's geographical range to a few northern areas [45]. Our

44 + kyr BP • • êÊÈÊv'' ^d y

34 - 44 kyr BP • V

• •

24 - 34 kyr BP • • J * • v H

14 - 24 kyr BP • * Ë • « • • • •V '."I,

4 - 14 kyr BP • * ^ •

Figure 3. Spatial distribution of radiocarbon-dated and genetically analysed mammoth specimens. Dates are given in calendar years before present. Colours indicate clade membership of the specimens: clade I; purple, clade II; pink, clade III; green.

genetic results appear to be consistent with this hypothesis, as they indicate that demographic expansions, presumably preceded by bottlenecks, took place around the end of the last interglacial (ca 121 kyr BP; figure 2a). This suggests that woolly mammoth populations were constrained in refugia

1 x 106

1 x 103

0 25000 50000 75000 100000 125000 calendar years before present

Figure 4. BSP of all dated woolly mammoth specimens. Effective population size (Nef) on the y-axis is biased by a factor of generation time assumed to 15 years. The thick solid line shows the average product of effective population size and generation time and the shaded area represents the 95% HPD. (Online version in colour.)

during the Eemian interglacial, and subsequently expanded in population size as well as geographically across Eurasia and North America as the climate became cooler and drier. This implies that severe bottlenecks at the onset of interglacials may have been a recurrent feature in the woolly mammoth's evolutionary history, and consequently lends support for the hypothesis that the end-Pleistocene retraction of the woolly mammoth from most of its range was driven by climate change.

From a genetic perspective, isolation in interglacial refugia and the subsequent demographic expansions indicated in our analyses could also provide an explanation for the existence of three distinct mtDNA clades during the Late Pleistocene. Divergent mtDNA lineages have been detected in many other taxa, and it has been argued that these most likely evolved owing to long-term allopatric isolation in glacial refugia [46,47]. As a cold-adapted species, the opposite pattern is expected for the woolly mammoth [48]. The presence of three distinct and highly diverged woolly mammoth mtDNA lineages (figure 1b) suggests that at least three separate interglacial refugia may have existed during the Eemian. The European distribution of clade III mammoths points to the existence of an interglacial refugium in western Eurasia, whereas the restricted distribution of clade II mammoths (figure 3) and American origin of clade I [6,11] imply two additional refugia, in northern Siberia and North America, respectively. However, the observed phylogenetic pattern should be interpreted with caution because it is based on mtDNA, which is a single non-recombining locus [49].

Our estimate of split time for the three populations, representing the three mtDNA clades, (ca 196 kyr BP) appears to be much younger than the tip calibrated divergence time estimated by Gilbert et al. [10] (ca 1 Ma). However, coalescence times of gene lineages (e.g. estimated in BEAST) typically predate population split times [50]. Moreover, it should be noted that tip calibration of the molecular clock requires a sufficient temporal signal in the data [19,51]. In contrast to our dataset, it has been shown that the full mtDNA genomes presented in Gilbert et al. [10] lacked such signal [19], likely owing to small sample size and limited temporal coverage. Thus, more complete mtDNA genomes need to be sequenced in order to further resolve the divergence times among the mtDNA clades.

Three specimens in our dataset from Yakutia (eastern Siberia) were stratigraphically dated to the LMP. These were

phylogenetically basal to, and thus fall outside the diversity of, the three mtDNA clades (figure 1b and electronic supplementary material, figures S4 and S5). The ages of these specimens, as estimated in BEAST, were approximately 175-200 kyr BP, and thus fall within the LMP (see electronic supplementary material, table S5). These specimens therefore predate the time to the most recent common ancestor of each mtDNA clade (figure 1b) as well as the expansion time of the three populations belonging to each clade, inferred from the simulations (approx. 121 kyr BP; figure 2a), and could therefore represent ancestral variation that existed before the Eemian population bottlenecks.

The decrease in temperatures and retreat of forests that occurred during the early stages of the last glaciation likely allowed for fragmented mammoth populations to expand their range [45]. Following the refugial scenario discussed above, woolly mammoths carrying mtDNA clade I would have expanded across northern North America, whereas clade II mammoths spread into Central and East Siberia and clade III mammoths colonized most of Europe. These regions were thus initially inhabited by discrete woolly mammoth populations, at least with respect to maternal gene flow, for thousands of years.

The onset of the last glaciation led to a gradual accumulation of ice in the Scandinavian and North American ice sheets, which in turn induced a decrease in sea levels [52]. The Bering Land Bridge that had been inundated owing to elevated sea levels during the Eemian became exposed, providing a corridor that connected northeast Siberia and Alaska. Previous studies have suggested that woolly mammoths belonging to clade I originated in North America and colonized Eurasia by crossing over the Bering Land Bridge but the timing of this expansion has been debated [6,11]. Using coalescent simulations, we estimated that the population split between North American and Eurasian mammoths carrying mtDNA clade I occurred at ca 66 kyr BP (figure 2c), a timing that could represent the westward dispersal of North American woolly mammoths into Eurasia. This coincides with the first exposure of the Bering Land Bridge since the penultimate glacial period [41] (figure 2d). Therefore, climate-induced sea-level changes could have promoted the expansion of clade I mammoths into Eurasia. It should be noted that our simulations did not explicitly test the hypothesis that woolly mammoths carrying mtDNA clade I from North America colonized Eurasia [6,11] but under this assumption inferred the timing of the presumable dispersal event.

Gene flow in the opposite direction, from Eurasia to North America after the populations split, was not supported by the simulations (see electronic supplementary material, table S6). This implies that the existence of sequences belonging to haplogroups D and E in North America [11] is owing to incomplete lineage sorting in the North American population. Assuming that founder effects took place in Siberia during the initial colonization, higher genetic drift would be expected in the Eurasian population compared with the North American population. This could be the reason why we observed an asymmetry in lineage sorting among the two populations.

If North American woolly mammoths expanded into Siberia, they would have encountered and likely interbred with genetically distinct individuals (carrying mtDNA clade II) that resided in Central and East Siberia. The two

mtDNA clades subsequently coexisted in sympatry for more than 20 kyr until clade II disappeared at around 40 kyr BP ( figure 3), either owing to selection or, more likely, genetic drift [6,10]. Furthermore, our data show that the geographical distribution of clade I continued to expand westwards, reaching European Russia at ca 32 kyr BP (figure 3). This first appearance of mammoths carrying mtDNA clade I in Europe appears to coincide with the disappearance of the endemic population in Europe (carrying mtDNA clade III) at ca 34 kyr BP (figure 3). However, in contrast to the genetic turnover in northeast Siberia, there is no evidence for a temporal overlap between clades I and III in Europe, which could indicate that the European population became extinct, and that Europe subsequently was recolonized by Siberian mammoths carrying clade I. There is size variation among mammoths in the last glaciation of Europe, but it does not map clearly to clade membership and may be owing to environmental effects, for example, vegetational productivity [53-55]. In Western Europe, large body size was maintained up to the last populations in the late-glacial [56].

It is interesting to note that similar population turnovers have been identified in cave bears from Central Europe at ca 31 kyr BP [2] and collared lemmings from northwestern Europe between ca 43-32 kyr BP [5]. The possible synchrony of these events could indicate that the turnovers were driven by regional changes in the environment. These population turnovers also appear to coincide approximately with the earliest findings of Gravettian culture in northwest Europe, around 34 kyr BP [57,58]. In the case of the mammoth, the extinction of clade III might even be associated with the emergence of the Gravettian culture itself, although this seems less likely as this would not explain how Siberian mammoths could have immediately recolonized Europe, and then survived there for at least another 15 kyr.

Spatio-temporal changes in fossil abundance suggest that northern populations of woolly mammoth generally declined between about 27 and 18 kyr BP (broadly the LGM), whereas populations in Central and South Siberia initially increased [16,59], and subsequently declined during the B0lling-Aller0d interstadial. The Younger Dryas (ca 12.9-11.7 kyr BP) was associated with extirpation in North America and southern Siberia, but temporary expansion in northern Siberia and into northeast Europe where mainland populations survived into the early Holocene (until ca 11 kyr BP) [9,16]. Such complex range contractions and asynchronous local extinctions resulted in the final demise of the mainland populations [9,60]. Previous studies have attempted to recover a signal of population decrease or the final extinction using genetic data but failed to do so ([6,11,59] but see [13]). Lorenzen et al. [36] simulated and compared different demographic models and found support for a population increase before the LGM (ca 26 kyr BP) in Eurasian mammoths. However, no significant changes in population size were observed

when looking at temporal dynamics of global Nef, or could be inferred from modelling the potential inhabitable range size through time. In contrast to these previous studies, we identified a marked end-Pleistocene population reduction (figure 4). Larger sample size with thorough temporal coverage, together with the inclusion of mammoth sequences from the post-bottleneck population on Wrangel Island in our dataset, offered the resolution required to capture significant changes in Nef [59]. We suggest that the revealed drastic drop in female effective population size most probably reflects the cumulative effects of the Late Pleistocene decline in most of the geographical range of woolly mammoths as well as the population bottleneck on Wrangel Island [13].

The results from this study reveal that the Late Pleistocene history of the woolly mammoth was characterized by a complex series of demographic changes, range expansions and clade replacements. Thus, while the high prevalence of mammoths in the fossil record might imply a stable and abundant species, populations of the woolly mammoth appear to have been highly dynamic. Both genetic data and the radiocarbon record indicate a dramatic final demographic decline at the end of the last glaciation. However, our results suggest that this decline was mirrored by a similar decline during the previous interglacial, a pattern that has also previously been observed in other cold-adapted taxa, such as reindeer [61], arctic fox [62] and polar bear [63]. It thus seems likely that environmental changes played a significant role in shaping the woolly mammoth's demographic history, with warm periods restricting the amount of available habitat and cold periods leading to population expansions, both owing to increases in the amount of steppe-tundra and through sea-level-driven exposure of the Bering Land Bridge. Resolving why the woolly mammoth survived in refugia during earlier interglacials, but not during the Holocene, may thus provide the key to understand the mechanism behind its final extinction.

Acknowledgements. The authors thank David Nogues-Bravo and three anonymous referees for comments on a previous version of this paper. The authors acknowledge the Canadian Museum of Nature (CMN), University of Alaska Fairbanks (UAF), Pirkko Ukkonen, Dale Guthrie, Martin Street, E.N. Mashchenko, Neil Clark, Peter Kurhy, Haowen Tong and Dustin White for providing samples. The authors are grateful to Beth Shapiro for samples and laboratory analyses, and Tony Stuart for samples and discussion. The authors also thank Ludovic Orlando, Christian Anderson, Tom Gilbert, Edson Sandoval-Castellanos and Pontus Skoglund for discussion and advice.

Funding statement. Funding for this study was provided by the EU FP6 ERA-NET project CLIMIGRATE, the Swedish Research Council (VR) and Marie Curie Actions FP6 grant no. 041545. Radiocarbon dating was funded through grant no. NF/2008/1/17 from the NERC radiocarbon facility. E.P. acknowledges funding from Stiftelsen Lars Hiertas Minne and the project 'IKY scholarships' financed by the operational program 'Education and Lifelong Learning' of the European Social Fund (ESF) and the NSRF 2007-2013. A.L. was funded by NERC grant nos. NE/D003105 and NE/G005982.


1. Leonard JA, Vila C, Fox-Dobbs K, Koch PL, Wayne RK, Van Valkenburgh B. 2007 Megafaunal extinctions and the disappearance of a specialized wolf ecomorph. Curr. Biol. 17,1146-1150. (doi:10. 1016/j.cub.2007.05.072)

2. Hofreiter M, Munzel S, Conard NJ, Pollack J, Slatkin M, Weiss G, Pääbo S. 2007 Sudden replacement of cave bear mitochondrial DNA in the late Pleistocene. Curr. Biol. 17, R122-R123. (doi:10.1016/j.cub.2007. 01.026)

3. Dalen L, Nystrom V, Valdiosera C, Germonpre M, Sablin M, Turner E, Angerbjorn A, Arsuaga JL, Gotherstrom A. 2007 Ancient DNA reveals lack of postglacial habitat tracking in the arctic fox. Proc. Natl Acad. Sci. USA 104, 6726-6729. (doi:10.1073/pnas.0701341104)

4. Prost S, Smirnov N, Fedorov VB, Sommer RS, Stiller M, Nagel D, Knapp M, Hofreiter M. 2010 Influence of climate warming on arctic mammals? New insights from ancient DNA studies of the collared lemming Dicrostonyx torquatus. PLoS ONE 5, e10447. (doi:10.1371/journal.pone.0010447)

5. Brace S et al. 2012 Serial population extinctions in a small mammal indicate Late Pleistocene ecosystem instability. Proc. Natl Acad. Sci. USA 109, 20 53220 536. (doi:10.1073/pnas.1213322109)

6. Barnes I, Shapiro B, Lister A, Kuznetsova T, Sher A, Guthrie D, Thomas MG. 2007 Genetic structure and extinction of the woolly mammoth, Mammuthus primigenius. Curr. Biol. 17, 1072-1075. (doi:10. 1016/j.cub.2007.05.035)

7. Campos PF et al. 2010 Ancient DNA analyses exclude humans as the driving force behind late Pleistocene musk ox (Ovibos moschatus) population dynamics. Proc. Natl Acad. Sci. USA 107, 5675-5680. (doi:10.1073/pnas.0907189107)

8. Edwards Ceiridwen J et al. 2011 Ancient hybridization and an Irish origin for the modern polar bear matriline. Curr. Biol. 21, 1251-1258. (doi:10.1016/j.cub.2011.05.058)

9. Stuart AJ, Sulerzhitsky LD, Orlova LA, Kuzmin YV, Lister AM. 2002 The latest woolly mammoths (Mammuthus primigenius Blumenbach) in Europe and Asia: a review of the current evidence. Q. Sci. Rev. 21, 1559-1569. (doi:10.1016/S0277-3791(02)00026-4)

10. Gilbert MTP et al. 2008 Intraspecific phylogenetic analysis of Siberian woolly mammoths using complete mitochondrial genomes. Proc. Natl Acad. Sci. USA 105, 8327-8332. (doi:10.1073/pnas. 0802315105)

11. Debruyne R et al. 2008 Out of America: ancient DNA evidence for a New World origin of Late Quaternary woolly mammoths. Curr. Biol. 18, 1320-1326. (doi:10.1016/j.cub.2008.07.061)

12. Nystrom V, Dalen L, Vartanyan S, Liden K, Ryman N, Angerbjorn A. 2010 Temporal genetic change in the last remaining population of woolly mammoth. Proc. R. Soc. B 277, 2331-2337. (doi:10.1098/rspb. 2010.0301)

13. Nystrom V et al. 2012 Microsatellite genotyping reveals end-Pleistocene decline in mammoth autosomal genetic variation. Mol. Ecol. 21, 3391-3402. (doi:10.1111/j.1365-294X.2012. 05525.x)

14. Lister AM, Sher AV. 2001 The origin and evolution of the woolly mammoth. Science 294, 1094-1097. (doi:10.1126/science.1056370)

15. Lister AM, Sher AV, van Essen H, Wei GB. 2005 The pattern and process of mammoth evolution in Eurasia. Q. Int. 126 - 28, 49 - 64. (doi:10.1016/j. quaint.2004.04.014)

16. Nikolskiy PA, Sulerzhitsky LD, Pitulko VV. 2011 Last straw versus Blitzkrieg overkill: climate-driven changes in the Arctic Siberian mammoth population and the Late Pleistocene extinction problem. Q. Sci. Rev. 30, 2309-2328. (doi:10.1016/j.quascirev.2010. 10.017)

17. Veltre DW, Yesner DR, Crossen KJ, Graham RW, Coltrain JB. 2008 Patterns of faunal extinction and paleoclimatic change from mid-Holocene mammoth and polar bear remains, Pribilof Islands, Alaska.

Q. Res. 70, 40 - 50. (doi:10.1016/j.yqres.2008. 03.006)

18. Vartanyan SL, Garutt VE, Sher AV. 1993 Holocene dwarf mammoths from Wrangel Island in the Siberian Arctic. Nature 362, 337-340. (doi:10.1038/ 362337a0)

19. Ho SYW, Lanfear R, Phillips MJ, Barnes I, Thomas JA, Kolokotronis S-O, Shapiro B. 2011 Bayesian estimation of substitution rates from ancient DNA sequences with low information content. Syst. Biol. 60, 366-375. (doi:10.1093/sysbio/syq099)

20. Ramsey B. 2009 Bayesian analysis of radiocarbon dates. Radiocarbon 51, 337-360.

21. Reimer P et al. 2009 IntCal09 and Marine09 radiocarbon age calibration curves, 0-50,000 years cal BP. Radiocarbon 51, 1111-1150.

22. Krause J et al. 2006 Multiplex amplification of the mammoth mitochondrial genome and the evolution of Elephantidae. Nature 439, 724-727. (doi:10. 1038/nature04432)

23. Rogaev EI, Moliaka YK, Malyarchuk BA, Kondrashov FA, Derenko MV, Chumakov I, Grigorenko AP. 2006 Complete mitochondrial genome and phylogeny of Pleistocene mammoth Mammuthus primigenius. PLoS Biol. 4, 403 - 410. (doi:10.1371/journal.pbio. 0040073)

24. Gilbert MTP et al. 2007 Whole-genome shotgun sequencing of mitochondria from ancient hair shafts. Science 317, 1927-1930. (doi:10.1126/ science.1146971)

25. Enk JM, Yesner DR, Crossen KJ, Veltre DW, O'Rourke DH. 2009 Phylogeographic analysis of the mid-Holocene Mammoth from QagnaxA Cave, St Paul Island, Alaska. Palaeogeogr. Palaeoclimatol. Palaeoecol. 273, 184-190. (doi:10.1016/j.palaeo. 2008.12.019)

26. Drummond AJ et al. 2012 Geneious v5.6. See http:// wwwgeneiouscom.

27. Lanfear R, Calcott B, Ho SYW, Guindon S. 2012 PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29, 1695-1701. (doi:10. 1093/molbev/mss020)

28. Ronquist F et al. 2012 MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 61, 539-542. (doi:10.1093/sysbio/sys029)

29. Drummond A, Rambaut A. 2007 BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7, 214. (doi:10.1186/1471-2148-7-214)

30. Rambaut A, Drummond AJ. 2007 Tracer v1.4: MCMC trace analyses tool. See software/tracer/.

31. Parks DH, Porter M, Churcher S, Wang S, Blouin C, Whalley J, Brooks S, Beiko RG. 2009 GenGIS: a geospatial information system for genomic data. Genome Res. 19, 1896-1904. (doi:10.1101/gr. 095612.109)

32. Anderson CNK, Ramakrishnan U, Chan YL, Hadly EA. 2005 Serial SimCoal: a population genetics model for data from multiple populations and points in time. Bioinformatics 21, 1733-1734. (doi:10.1093/ bioinformatics/bti154)

33. Beaumont MA. 2010 Approximate Bayesian computation in evolution and ecology. Annu. Rev. Ecol. Evol. Syst. 41, 379 - 406. (doi:10.1146/ annurev-ecolsys-102209-144621)

34. Csillery K, Francois O, Blum MGB. 2012 abc: an R package for approximate Bayesian computation (ABC). Methods Ecol. Evol. 3, 475-479. (doi:10. 1111/j.2041-1210X.2011.00179.x)

35. Excoffier L, Lischer HEL. 2010 Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564-567. (doi:10.1111/j.1755-0998.2010.02847.x)

36. Lorenzen ED et al. 2011 Species-specific responses of Late Quaternary megafauna to climate and humans. Nature 479, 359-364. (doi:10.1038/ nature10574)

37. Sukumar R. 1989 The Asian elephant: ecology and management. Cambridge, UK: Cambridge University Press.

38. Beaumont MA, Zhang WY, Balding DJ. 2002 Approximate Bayesian computation in population genetics. Genetics 162, 2025-2035.

39. Blum MGB, Francois O. 2010 Non-linear regression models for approximate Bayesian computation. Stat. Comput. 20, 63-73. (doi:10.1007/s11222-009-9116-0)

40. Csillery K, Blum MGB, Gaggiotti O, Francois O. 2010 Approximate Bayesian computation (ABC) in practice. Trends Ecol. Evol. 25, 410-418. (doi:10. 1016/j.tree.2010.04.001)

41. Hu A, Meehl GA, Otto-Bliesner BL, Waelbroeck C, Han W, Loutre M-F, Lambeck K, Mitrovica JX, Rosenbloom N. 2010 Influence of Bering Strait flow and North Atlantic circulation on glacial sea-level changes. Nat. Geosci. 3, 118-121. (doi:10.1038/ ngeo729)

42. Kukla GJ et al. 2002 Last interglacial climates. Q. Res. 58, 2-13. (doi:10.1006/qres.2001.2316)

43. Sher AV. 1991 Problems of the last interglacial in Arctic Siberia. Q. Int. 10-12, 215-222. (doi:10. 1016/1040-6182(91)90053-q)

44. Lister AM. 2004 The impact of Quaternary ice ages on mammalian evolution. Proc. R. Soc. Lond. B 359, 221 -241. (doi:10.1098/rstb.2003.1436)

45. Nogues-Bravo D, Rodriguez J, Hortal J, Batra P, Araujo MB. 2008 Climate change, humans, and the extinction of the woolly mammoth. PLoS Biol. 6, e79. (doi:10.1371/journal.pbio.0060079)

46. Hewitt GM. 1996 Some genetic consequences of ice ages, and their role in divergence and speciation. Biol. J. Linnean Soc. 58, 247 - 276.

47. Taberlet P, Cheddadi R. 2002 Quaternary refugia and persistence of biodiversity. Science 297, 2009 - 2010. (doi:10.1126/science.297.5589.2009)

48. Stewart JR, Lister AM, Barnes I, Daleun L. 2010 Refugia revisited: individualistic responses of species

in space and time. Proc. R. Soc. B 277, 661-671. (doi:10.1098/rspb.2009.1272)

49. Ballard JWO, Whitlock MC. 2004 The incomplete natural history of mitochondria. Mol. Ecol. 13, 729-744. (doi:10.1046/j.1365-294X.2003.02063.x)

50. Edwards SV, Beerli P. 2000 Perspective: gene divergence, population divergence, and the variance in coalescence time in phylogeographic studies. Evolution 54, 1839-1854. (doi:10.2307/2640530)

51. Drummond AJ, Pybus OG, Rambaut A, Forsberg R, Rodrigo AG. 2003 Measurably evolving populations. Trends Ecol. Evol. 18, 481-488. (doi:10.1016/ S0169-5347(03)00216-7)

52. Lambeck K, Chappell J. 2001 Sea level change through the last glacial cycle. Science 292, 679-686. (doi:10.1126/science.1059549)

53. Germonpre M. 1993 Osteometric data on Late Pleistocene mammals from the Flemish Valley. Belgium: Documents de Travail de l'Institut royal des Sciences Naturelles de Belgique.

54. Germonpre M, Sablin M, Khlopachev GA, Grigorieva GV. 2008 Possible evidence of mammoth hunting

during the Epigravettian at Yudinovo, Russian Plain. J. Anthropol. Archaeol. 27, 475 - 492. (doi:10.1016/j. jaa.2008.07.003)

55. Maschenko EN. 2002 Individual development, biology and evolution of the woolly mammoth Mammuthus primigenius (Blumenbach, 1799). Cranium. 19, 4-120.

56. Lister AM. 2009 Late-glacial mammoth skeletons (Mammuthus primigenius) from Condover (Shropshire, UK): anatomy, pathology, taphonomy and chronological significance. Geol. J. 44,

447 - 479. (doi:10.1002/gj.1162)

57. Jacobi RM, Higham TFG. 2008 The 'Red Lady' ages gracefully: new ultrafiltration AMS determinations from Paviland. J. Hum. Evol. 55, 898-907. (doi:10. 1016/j.jhevol.2008.08.007)

58. Pesesse D, Flas D. 2012 The Maisierian, at the Edge of the Gravettian. Proc. Prehistoric Soc. 78, 95-109. (doi:10.1017/S0079497X00027122)

59. MacDonald GM, Beilman DW, Kuzmin YV, Orlova LA, Kremenetski KV, Shapiro B, Wayne RK, Van Valkenburgh B. 2012 Pattern of extinction of the

woolly mammoth in Beringia. Nat. Commun. 3, 893. (doi:10.1038/ncomms1881)

60. Lister AM, Stuart AJ. 2008 The impact of climate change on large mammal distribution and extinction: evidence from the last glacial/interglacial transition. C. R. Geosci. 340, 615-620. (doi:10. 1016/j.crte.2008.04.001)

61. Flagstad O, Roed KH. 2003 Refugial origins of reindeer (Rangifer tarandus L.) inferred from mitochondrial DNA sequences. Evolution 57, 658-670.

62. Dalen L, Fuglei E, Hersteinsson P, Kapel CMO, Roth JD, Samelius G, Tannerfeldt M, Angerbjorn A. 2005 Population history and genetic structure of a circumpolar species: the arctic fox. Biol. J. Linnean Soc. 84, 79-89. (doi:10.1111/j.1095-8312.2005. 00415.x)

63. Miller W et al. 2012 Polar and brown bear genomes reveal ancient admixture and demographic footprints of past climate change. Proc. Natl Acad. Sci. USA 109, E2382-E2390. (doi:10.1073/pnas. 1210506109)