# Ancient DNA suggests modern wolves trace their origin to a Late Pleistocene expansion from Beringia

Liisa Loog, Olaf Thalmann, Mikkel-Holger Sinding, Verena Schuenemann, Angela Perri, Mietje Germonpré, Herve Bocherens, Kelsey Witt, Jose Samaniego Castruita, Marcela Velasco, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Liisa Loog, Olaf Thalmann, Mikkel-Holger Sinding, Verena Schuenemann, Angela Perri, et al. *An- cient DNA suggests modern wolves trace their origin to a Late Pleistocene expansion from Beringia*. Molecular Ecology, 2020, 29 (9), pp.1596-1610. 10.1111/mec.15329. hal-03353285

HAL Id: hal-03353285

https://hal.science/hal-03353285v1

Submitted on 22 Aug 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Ancient DNA suggests modern wolves trace their origin to a Late Pleistocene expansion from Beringia

Liisa Loog<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> ID

Olaf Thalmann<fnref n="5" />

Mikkel-Holger S. Sinding<fnref n="6" />,<fnref n="7" />,<fnref n="8" /> ID

Verena J. Schuenemann<fnref n="9" />,<fnref n="10" />,<fnref n="11" />

Angela Perri<fnref n="12" />

Mietje Germonpré<fnref n="13" />

Herve Bocherens<fnref n="10" />,<fnref n="14" />

Kelsey E. Witt<fnref n="15" />

Jose A. Samaniego Castruita<fnref n="6" />

Marcela S. Velasco<fnref n="6" />

Inge K. C. Lundstrøm<fnref n="6" />

Nathan Wales<fnref n="6" />,<fnref n="16" />

Gontran Sonet<fnref n="17" /> ID

Laurent Frantz<fnref n="2" />

Hannes Schroeder<fnref n="6" />

Jane Budd<fnref n="18" />

Elodie-Laure Jimenez<fnref n="13" />

Sergey Fedorov<fnref n="19" />

Boris Gasparyan<fnref n="20" />

Andrew W. Kandel<fnref n="21" />

Martina Lázničková-Galetová<fnref n="22" />,<fnref n="23" />,²⁴

Hannes Napierala²⁵

Hans-Peter Uerpmann<fnref n="9" />

Pavel A. Nikolskiy²⁶,²⁷

Elena Y. Pavlova²⁷,²⁸

Vladimir V. Pitulko²⁷

Karl-Heinz Herzig<fnref n="5" />,²⁹

Ripan S. Malhi³⁰

Eske Willerslev<fnref n="2" />,³¹,³²

Anders J. Hansen<fnref n="8" />,³¹

Keith Dobney³³,³⁴,³⁵

M. Thomas P. Gilbert<fnref n="6" />,³⁶

Johannes Krause<fnref n="9" />,³⁷

Greger Larson<fnref n="1" />

Anders Eriksson<fnref n="2" />,³⁸,³⁹ ID

Andrea Manica<fnref n="2" />

<footnote n="1">¹Research Laboratory for Archaeology and History of Art, University of Oxford, Oxford, UK</footnote>

<footnote n="2">²Department of Zoology, University of Cambridge, Cambridge, UK</footnote>

<footnote n="3">³Manchester Institute of Biotechnology, School of Earth and Environmental Sciences, University of Manchester, Manchester, UK</footnote>

<footnote n="4">⁴Department of Genetics, University of Cambridge, Cambridge, UK</footnote>

<footnote n="5">⁵Department of Pediatric Gastroenterology and Metabolic Diseases, Poznan University of Medical Sciences, Poznan, Poland</footnote>

<footnote n="6">⁶EvoGenomics, GLOBE Institute, University of Copenhagen, Copenhagen, Denmark</footnote>

<footnote n="7">⁷Natural History Museum, University of Oslo, Oslo, Norway</footnote>

<footnote n="8">⁸The Qimmeq Project, University of Greenland, Nuussuaq, Greenland</footnote>

<footnote n="9">⁹Institute for Archaeological Sciences, University of Tübingen, Tübingen, Germany</footnote>

<footnote n="10">¹⁰Senckenberg Centre for Human Evolution and Palaeoenvironment, University of Tübingen, Tübingen, Germany</footnote>

<footnote n="11">¹¹Institute of Evolutionary Medicine, University of Zurich, Zurich, Switzerland</footnote>

<footnote n="12">¹²Department of Human Evolution, Max Planck Institute for Evolutionary Anthropology, Leipzig, Germany</footnote>

<footnote n="13">¹³OD Earth and History of Life, Royal Belgian Institute of Natural Sciences, Brussels, Belgium</footnote>

<footnote n="14">¹⁴Department of Geosciences, Palaeobiology, University of Tübingen, Tübingen, Germany</footnote>

<footnote n="15">¹⁵School of Integrative Biology, University of Illinois at Urbana-Champaign, Urbana, IL, USA</footnote>

<footnote n="16">¹⁶BioArch, Department of Archaeology, University of York, York, UK, USA</footnote>

<footnote n="17">¹⁷OD Taxonomy and Phylogeny, Royal Belgian Institute of Natural Sciences, Brussels, Belgium</footnote>

<footnote n="18">¹⁸Breeding Centre for Endangered Arabian Wildlife, Sharjah, United Arab Emirates</footnote>

<footnote n="19">¹⁹Mammoth Museum, Institute of Applied Ecology of the North of the North-Eastern Federal University, Yakutsk, Russia</footnote>

<footnote n="20">²⁰Institute of Archaeology and Ethnography, National Academy of Sciences of the Republic of Armenia, Yerevan, Republic of Armenia</footnote>

<footnote n="21">²¹Heidelberg Academy of Sciences and Humanities: The Role of Culture in Early Expansions of Humans, Tübingen, Germany</footnote>

<footnote n="22">²²Department of Anthropology, University of West Bohemia, Pilzen, Czech Republic</footnote>

<footnote n="23">²³Moravian museum, Brno, Czech Republic</footnote>

<footnote>Thalmann, Sinding, and Schuenemann contributed equally to this work.</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>

<footnote>© 2019 The Authors. *Molecular Ecology* published by John Wiley & Sons Ltd</footnote>
<sup>24</sup>Hrdlička Museum of Man, Faculty of Science, Charles University, Praha, Czech Republic

<sup>25</sup>Institute of Palaeoanatomy, Domestication Research and History of Veterinary Medicine, Ludwig-Maximilians-University Munich, Munich, Germany

<sup>26</sup>Geological Institute, Russian Academy of Sciences, Moscow, Russia

<sup>27</sup>Institute for Material Culture History, Russian Academy of Sciences, St Petersburg, Russia

<sup>28</sup>Arctic and Antarctic Research Institute, St Petersburg, Russia

<sup>29</sup>Institute of Biomedicine and Biocenter of Oulu, Medical Research Center and University Hospital, University of Oulu, Oulu, Finland

<sup>30</sup>Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign, Urbana, IL, USA

<sup>31</sup>Centre for GeoGenetics Globe Institute, University of Copenhagen, Copenhagen, Denmark

<sup>32</sup>Wellcome Trust Sanger Institute, Cambridge, UK

<sup>33</sup>Department of Archaeology, Classics and Egyptology, University of Liverpool, Liverpool, UK

<sup>34</sup>Department of Archaeology, University of Aberdeen, Aberdeen, UK

<sup>35</sup>Department of Archaeology, Simon Fraser University, Burnaby, BC, Canada

<sup>36</sup>Norwegian University of Science and Technology, University Museum, Trondheim, Norway

<sup>37</sup>Max Planck Institute for the Science of Human History, Jena, Germany

<sup>38</sup>Department of Medical & Molecular Genetics, King's College London, Guys Hospital, London, UK

<sup>39</sup>cGEM, Institute of Genomics, University of Tartu, Tartu, Estonia

**Correspondence**

Liisa Loog, Department of Genetics,
University of Cambridge, Cambridge, UK.
Email: liisaloog@gmail.com

Greger Larson, Research Laboratory for
Archaeology and History of Art, University
of Oxford, Oxford, UK.
Email: greger.larson@arch.ox.ac.uk

Anders Eriksson, cGEM, Institute of
Genomics, University of Tartu, Tartu,
Estonia.
Email: aeriksson75@gmail.com

Andrea Manica, Department of Zoology,
University of Cambridge, Cambridge, UK.
Email: am315@cam.ac.uk

**Funding information**

Russian Science Foundation, Grant/Award
Number: N16-18-10265 RNF; Grantová
Agentura České Republiky, Grant/Award
Number: GAČR15-064465; Polish National
Science Center, Grant/Award Number:
2015/19/P/NZ7/03971; Lundbeckfonden,
Grant/Award Number: R52-5062;
Natural Environment Research Council,
Grant/Award Number: NE/K003259/1
and NE/K005243/1; H2020 European
Research Council, Grant/Award Number:
339941-ADAPT, 647787-LocalAdaptation,
681396-Extinction Genomics and ERC-
2013-StG 337574-UNDEAD

**Abstract**

Grey wolves (*Canis lupus*) are one of the few large terrestrial carnivores that have maintained a wide geographical distribution across the Northern Hemisphere throughout the Pleistocene and Holocene. Recent genetic studies have suggested that, despite this continuous presence, major demographic changes occurred in wolf populations between the Late Pleistocene and early Holocene, and that extant wolves trace their ancestry to a single Late Pleistocene population. Both the geographical origin of this ancestral population and how it became widespread remain unknown. Here, we used a spatially and temporally explicit modelling framework to analyse a data set of 90 modern and 45 ancient mitochondrial wolf genomes from across the Northern Hemisphere, spanning the last 50,000 years. Our results suggest that contemporary wolf populations trace their ancestry to an expansion from Beringia at the end of the Last Glacial Maximum, and that this process was most likely driven by Late Pleistocene ecological fluctuations that occurred across the Northern Hemisphere. This study provides direct ancient genetic evidence that long-range migration has played an important role in the population history of a large carnivore, and provides insight into how wolves survived the wave of megafaunal extinctions at the end of the last glaciation. Moreover, because Late Pleistocene grey wolves were the likely source from which all modern dogs trace their origins, the demographic history described in this study has fundamental implications for understanding the geographical origin of the dog.

**KEYWORDS**

Approximate Bayesian Computation, ancient DNA, coalescent modelling, megafauna,
Pleistocene, population structure, population turnover, wolves

# 1 | INTRODUCTION

The Pleistocene epoch harboured a large diversity of top predators,
although most became extinct during or soon after the Last Glacial

Maximum (LGM), ~21,000 years ago (Barnosky, Koch, Feranec, Wing,
& Shabel, 2004; Clark et al., 2012). The grey wolf (*Canis lupus*) was one
of the few large carnivores that survived and maintained a wide geo-
graphical range throughout the period (Puzachenko & Markova, 2016),
and both the palaeontological and archaeological records attest to the
continuous presence of grey wolves across the Northern Hemisphere
for at least the last 300,000 years (Sotnikova & Rook, 2010) (reviewed
in Appendix S1). This geographical and temporal continuity across the
Northern Hemisphere contrasts with analyses of complete modern
genomes, which have suggested that all contemporary wolves and
dogs descend from a common ancestral population that existed as
recently as 20,000 years ago (Fan et al., 2016; Freedman et al., 2014;
Skoglund, Ersmark, Palkopoulou, & Dalén, 2015).

These analyses point to a bottleneck followed by a rapid radia-
tion from an ancestral population around or just after the LGM. The
geographical origin and dynamics of this radiation remain unknown.
Resolving these demographic changes is necessary for understand-
ing the ecological circumstances that allowed wolves to survive the
Late Pleistocene megafaunal extinctions. Furthermore, because
dogs were domesticated from Late Pleistocene grey wolves (Larson
et al., 2012), a detailed insight into wolf demography during this time
period would provide an essential context for reconstructing the his-
tory of dog domestication.

Reconstructing past demographic events solely from modern ge-
nomes is challenging because multiple demographic histories can lead
to similar genetic patterns in present-day samples (Groucutt et al.,
2015). Analyses that incorporate ancient DNA sequences can elim-
inate some of these alternative histories by quantifying changes in
population genetic differences through time. While nuclear markers
provide greater power relative to mitochondrial DNA (mtDNA), the
latter is more easily retrievable and better preserved in ancient sam-
ples due to its higher copy number compared to nuclear DNA, thus
allowing for the generation of data sets with greater geographical and
temporal coverage. In particular, analysing samples dated to before,
during and after the demographic events of interest greatly increases
the power to infer past demographic histories. Furthermore, the nu-
clear mutation rate in canids is poorly understood, leading to wide
date ranges for past demographic events reconstructed from panels
of modern whole genomes (e.g., Fan et al., 2016; Freedman et al.,
2014). Having directly dated samples from a broad time period al-
lows us to estimate mutation rates with higher accuracy and precision
compared to alternative methods (Drummond, Nicholls, Rodrigo, &
Solomon, 2002; Rambaut, 2000; Rieux et al., 2014).

Demographic processes, such as range expansions and contrac-
tions, that involved space as well as time are particularly challeng-
ing to reconstruct as they often lead to patterns that are difficult
to interpret intuitively (Groucutt et al., 2015). Hypotheses involving
spatial processes can be formally tested using population genetic
models that explicitly represent the various demographic processes
and their effect on genetic variation through time and across space
(Eriksson et al., 2012; Eriksson & Manica, 2012; Posth et al., 2016;
Raghavan et al., 2015; Warmuth et al., 2012). The formal integration
of time and space into population genetics frameworks allows for
the analysis of sparse data sets, a common challenge when dealing
with ancient DNA (Loog et al., 2017).

Here, we use a spatially explicit population genetic framework to
model a range of different demographic histories of wolves across

the Northern Hemisphere that involve combinations of population
bottlenecks, turnover and long-range migrations as well as local gene
flow. To estimate model parameter and formally test hypotheses of
the origin and population dynamics of the expansion of grey wolves
during the LGM, we assembled a substantial data set (Figure 1; Table
S1), spanning the last 50,000 years and the geographical breadth
of the Northern Hemisphere. This data set consists of 90 modern
and 45 ancient wolf whole mitochondrial genomes (38 of which are
newly sequenced). In the following, we first present a phylogenetic
analysis of our sequences and a calibration of the wolf mitochondrial
mutation rates. We then perform formal hypothesis testing using
Approximate Bayesian Computation (ABC) with our spatiotempo-
rally explicit models. We conclude with a discussion of how our find-
ings relate to earlier studies and implications for future research.

2 | MATERIALS AND METHODS

2.1 | Data preparation

We sequenced whole mitochondrial genomes of 40 ancient wolf
samples. Sample information, including geographical locations,
estimated ages and archaeological context information for the an-
cient samples, is provided in Table S1 and Appendix S1. Of the 40
ancient samples, 24 were directly radiocarbon dated for this study
and calibrated using the IntCal13 calibration curve (see Table S1 for
radiocarbon dates, calibrated age ranges and accelerator mass spec-
trometry [AMS] laboratory reference numbers). DNA extraction,
sequencing and quality filtering, and mapping protocols used are
described in Appendix S2.

We included 16 previously published ancient mitochondrial wolf
genomes (Table S1 and Appendix S2). To achieve a uniform data set,
we reprocessed the raw reads from previously published samples
using the same bioinformatics pipeline as for the newly generated
sequences.

We subjected the aligned ancient sequences to strict quality cri-
teria in terms of damage patterns and missing data (Figures S3–S5).
First, we excluded all whole mitochondrial sequences that had more
than one-third of the whole mitochondrial genome missing (excluding
the mitochondrial control region—see below) at minimum three-fold
coverage. Second, we excluded all ancient whole mitochondrial se-
quences that contained more than 0.1% of singletons showing signs
of deamination damage typical for ancient DNA (C to T or A to G sin-
gletons). After quality filtering, we were left with 32 newly sequenced
and 13 published ancient whole mitochondrial sequences (Table S1).

We also excluded sequences from archaeological specimens that post-date the end of the Pleistocene and that have been identified as dogs (Table S1), because any significant population structure resulting from a lack of gene flow between dogs and wolves could violate the assumption of a single, randomly mating canid population. Some of the Pleistocene specimens used in the demographic analyses (TH5, TH12, TH14) have been argued to show features commonly found in modern dogs and have therefore been suggested to
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** Geographical distribution of modern (<500 years old, circles) and ancient (>500 years old, triangles) samples (a) and temporal distribution of ancient samples (b) used in the analyses. The geographical locations of the samples have been slightly adjusted for clarity (see Table S1 for exact sample locations). *Samples dated by molecular dating [Colour figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

represent Palaeolithic dogs (e.g., Druzhkova et al., 2013; Germonpré, Lázničková-Galetová, Losey, Räikkönen, & Sablin, 2015; Germonpré, Lázničková-Galetová, & Sablin, 2012; Germonpré et al., 2009; Sablin & Khlopachev, 2002). Here, we disregard such status calls because of the controversy that surrounds them (Crockford & Kuzmin, 2012; Drake, Coquerelle, & Colombeau, 2015; Morey, 2014; Perri, 2016), and because early dogs would have been genetically similar to the local wolf populations from which they derived. This reasoning is supported by the close proximity of these samples to other wolf specimens confidently described as wolves in the phylogenetic tree (see Figure S10).

Finally, we sequenced six samples from modern wolves and added 66 modern published wolf sequences from NCBI, two sequences from Freedman et al. (2014), 13 sequences from Sinding et al. (2018) and three sequences from Gopalakrishnan et al. (2018) (Table S1). Data from Sinding et al. (2018) and Gopalakrishnan et al. (2018) were newly assembled following the same bioinformatics protocols as were used for newly sequenced modern wolf samples (see Appendix S2). This resulted in a final data set of 135 complete wolf mitochondrial genome sequences, of which 45 were ancient and 90 were modern. We used the CLUSTALW alignment tool (version 2.1) (Larkin et al., 2007) to generate a joint alignment of all genomes.
To avoid the potentially confounding effect of recurrent mutations in the mitochondrial control region (Excoffier & Yang, 1999) in pairwise difference calculations, we removed this region from all subsequent analyses. This resulted in an alignment of sequences 15,466 bp in length, of which 1,301 sites (8.4%) were variable. The aligned data set is given in Appendix S1.

## 2.2 | Phylogenetic analysis

We calculated the number of pairwise differences between all samples (Figure S6) and generated a neighbour-joining tree based on pairwise differences (Figure S7). This tree shows a clade consisting of samples exclusively from the Tibetan region and the Indian subcontinent that are deeply diverged from all ancient and other modern wolf samples (see also Aggarwal, Kivisild, Ramadevi, & Singh, 2007; Sharma, Maldonado, Jhala, & Fleischer, 2004). A recent study of whole genome data showed a complex history of South Eurasian wolves (Fan et al., 2016) that is beyond the scope of our study. While their neighbouring phylogeny grouped South Eurasian wolves with East and North East Asian wolves (Fan et al., 2016: Figure 3), they cluster outside of all other grey wolves in a principal component analysis (Fan et al., 2016: Figure 4), and also show a separate demographic history within a Pairwise Sequentially Markovian Coalescent analysis (PSMC) (Fan et al., 2016: Figure 5). Because our study did not possess sufficient samples from the Himalayas and the Indian subcontinent to unravel their complex demography, we excluded samples from these regions and focused on the history of North Eurasian and North American wolves, for which we have good coverage through time and space.

We used PARTITIONFINDER (Lanfear, Calcott, Ho, & Guindon, 2012) and BEAST (version 1.8.0) (Drummond, Suchard, Xie, & Rambaut, 2012) to build a tip calibrated wolf mitochondrial tree (with a strict global clock, see Appendix S1 for full details) from modern and directly dated ancient samples, and to estimate mutation rates for four different partitions of the wolf mitochondrial genome (see Tables S3 and S4 for results).

We used BEAST to molecularly date seven sequences from samples that were not directly radiocarbon dated (TH4, TH6, TH14, TU15) or that had been dated to a period beyond the limit of reliable radiocarbon dating (>48,000 years ago) (CGG12, CGG29, CGG32). We estimated the ages of the samples by performing a BEAST run where the mutation rate was fixed to the mean estimates from the previous BEAST analysis and all other parameter settings were set as described in Appendix S1. We cross-validated this approach through a leave-one-out analysis where we sequentially removed a directly dated sample and estimated its date as described above. We find a close fit ($R^2$ = 0.86) between radiocarbon and molecular dates (Figure S9). We combined the seven undated samples with the 110 ancient and modern samples from the previous run and used a uniform prior ranging from 0 to 100,000 years to estimate the ages of the seven undated samples (see Table S5 for results).

Finally, in order to estimate the mitochondrial divergence time between the South Eurasian (Tibetan and Indian) and the rest of

our wolf samples, we performed an additional BEAST run in which we included all modern and ancient grey wolves (N = 129) as well as five Tibetan and one Indian wolf, and used parameters identical to those described above. The age of the ancient samples was set as the mean of the calibrated radiocarbon date distribution (for radio-carbon-dated samples) or as the mean of the age distribution from the BEAST analyses (for molecularly dated samples).

## 2.3 | Isolation by distance analysis

We performed isolation by distance (IBD) analyses to see the extent to which wolf mitochondrial genetic variation shows population structure. To this end, we regressed the pairwise geographical distances between 84 modern wolf samples (Table S1) against their pairwise genetic (mitochondrial) distances. The geographical distance between all sample pairs was calculated in kilometres as the great circle distance from geographical coordinates, using the Haversine formula (Sinnott, 1984) to account for the curvature of the Earth as follows:

$$G_{ij} = 2r \arcsin(\sqrt{\sin^2((\varphi_j - \varphi_i)/2)^2 + \cos(\varphi_i) \cos(\varphi_j) \sin^2((\lambda_i - \lambda_j)/2)^2})$$

where G is the distance in kilometres between individuals i and j; $\varphi_i$ and $\varphi_j$ are the latitude coordinates of individuals i and j, respectively; $\lambda_i$ and $\lambda_j$ are the longitude coordinates of individuals i and j, respectively; and r is the radius of the earth in kilometres. The pairwise genetic distances were calculated as the proportion of sites that differ between each pair of sequences (excluding the missing bases), using dist.dna function in the R package APE (Paradis, Claude, & Strimmer, 2004).

## 2.4 | Geographical deme definitions

We represented the wolf geographical range as seven demes, defined by major geographical barriers through time.

1. The European deme is bordered by open water from the north and the west (the Arctic and the Atlantic oceans, respectively); the Ural Mountains from the east; and the Mediterranean, the Black and the Caspian Sea and the Caucasus mountains from the south.

2. The Middle-Eastern deme consists of the Arabian Peninsula, Anatolia and Mesopotamia and is bordered by the Black Sea, the Caspian Sea and the Aral Sea in the north; the Indian Ocean in the south; the Tien Shen mountain range, the Tibetan Plateau and the Himalayas from the east; and the Mediterranean Sea in the west.

3. The Central North Eurasian deme consist of the Siberian Plateau and is bordered by the Arctic Ocean from the north; the Ural Mountains from the west; the Lena River and mountain ranges of northeastern Siberia (Chersky and Verkhoyansk ranges) from the east; and the Tien Shen mountain range, the Tibetan Plateau and the Gobi Desert from southeast.
4. The *East Eurasian* deme is bordered by the Tien Shen mountain range, the Tibetan Plateau and Gobi desert from the west; the Pacific Ocean from the east; and the Lena river and the mountain ranges of northeastern Siberia (Chersky and Verkhoyansk ranges) from the north.

5. The *Beringia* deme spans the Bering Strait, which was a land bridge during large parts of the Late Pleistocene and the Early Holocene. It is bordered to the west by the Lena River and mountain ranges of northeastern Siberia (Chersky and Verkhoyansk ranges), and to the south and east by the extent of the Cordillerian and Laurentide ice sheets during the LGM.

6. The *Arctic North America* deme consists of an area of the North American continent east of the Rocky Mountains and west of Greenland, that was covered by ice during the last glaciation and is at present known as the Canadian Arctic Archipelago.

7. The *North America* deme consists of an area in the Northern American subcontinent up to and including the area that was covered by the Cordillerian and Laurentide ice sheets during the last glaciation (Raghavan et al., 2015).

## 2.5 | AMOVAS

To quantify the extent that our geographical demes capture genetic variation in the data we performed analyses of molecular variance (AMOVAS) (Excoffier, Smouse, & Quattro, 1992). We calculated the pairwise genetic distance between all modern wolf (n = 84, Table S1) sample pairs as described above (Isolation by distance analysis) and partitioned the samples, based on their geographical locations, into seven populations corresponding the geographical demes, as described in Section 2.4, Geographical deme definitions. We used these demes as the level of analyses and performed 1 million permutations using the *amova* function in the R package pegas (version 0.10). We found strong support for our geographical demes (p < 10⁻⁶) with 24.4% of the variance within the data set explained by the chosen demes.

## 2.6 | Demographic scenarios

We tested a total of 16 demographic scenario combinations, from four different kinds of demographic scenarios (illustrated in Figure 4a):

1. Static model (the null hypothesis)—neighbouring demes exchange migrants, no demographic changes.

2. Bottleneck scenarios—demes exchange migrants as in the static model but populations have different size in different time periods. We consider three time periods: 0–15,000 years ago, 15,000–40,000 years ago, and >40,000 years ago.

3. Expansion scenarios—demes exchange migrants as in the static model but a single deme (which itself has a continuous population through time) experiences an expansion starting between 5,000 and 40,000 years ago (at a minimum rate of 1,000 years per deme, so the whole world could be colonized within 3,000 years

or faster). The deme of origin has a continuous population through time while native populations in all other demes experience replacement—allowing us to formally test both the continuity and replacement hypotheses in each of the demes.

4. Combinations of scenarios 2 and 3.

## 2.7 | Population genetic coalescent framework

We implemented coalescent population genetic models for the different demographic scenarios to sample gene genealogies.

In the static scenario, we simulated local coalescent processes (Kingman, 1982) within each deme (scaled to rate 1/K per pair of lineages, where K is the mean time to the most recent common ancestor (TMRCA) in a deme and is thus proportional to the effective population size). In addition, we moved lineages between demes according to a Poisson process with rate m per lineage. To match the geographical and temporal distribution of the data, we represented each sample with a lineage from the corresponding deme and date.

The bottleneck scenario was implemented as the static one but with piecewise constant values for K as a function of time. We considered three time periods, each with its own value of K ($K_1, K_2$ and $K_3$), motivated by the archaeological and genetic evidence of wolf population changes described in the main text. The first time period was from the present to early Holocene, 0–15,000 years ago. The second time period extended from early Holocene to the Pleistocene and covered the LGM, 15,000–40,000 years ago. Finally, the third time period covered the Late Pleistocene and beyond, that is 40,000 years ago and older.

The population expansion scenarios were based on the static model but with an added population expansion model with founder effects and replacement of local populations (we refer to populations not yet replaced by the expansion as “indigenous”). Starting at time T, the population expanded from the initial deme and replaced its neighbouring populations. The population at the deme of origin was represented as a continuous population through time. After the start of the expansion, the expansion proceeded in fixed steps of $\Delta T$ (in time). At each step, colonized populations replaced neighbouring indigenous populations (if an indigenous deme bordered more than one colonized deme, these demes contributed equally to the colonization of the indigenous deme). In the coalescent framework (which simulates gene genealogies backwards in time) the colonization events corresponds to forced migrations from the indigenous deme to the source deme. If there were more than one source deme, the source of each lineage was chosen randomly with equal probability. Finally, founder effects during the colonization of an indigenous deme were implemented as a local, instantaneous population bottleneck in the deme (after the expansion), with a severity scaled to give a fixed probability x of a coalescent event for each pair of lineages in the deme during the bottleneck (Eriksson & Mehlig, 2004) (x = 1 corresponds to a complete loss of genetic diversity in the bottleneck, and x = 0 corresponds to no reduction in genetic diversity).

Finally, the combined scenario of population expansion and bottlenecks was implemented by making the population size parameter
K in the population expansion model time-dependent as in the pop-
ulation bottleneck model.

**2.8 | Approximate Bayesian Computation analysis**

We used ABC analysis (Beaumont, Zhang, & Balding, 2002) with ABCTOOLBOX (Wegmann, Leuenberger, Neuenschwander, & Excoffier, 2010) to formally test the fit of our different demographic models. This approach allows formal hypothesis testing using likelihood ratios in the cases where the demographic scenarios are too complex for a direct calculation of the likelihoods given the models. We used the most likely tree from BEAST (see Appendix S1 for details) as data, and simulated trees using the coalescent simulations described above.

To match the assumption of random mixing within each
deme in the population genetic model, we removed closely re-
lated sequences if they came from the same geographical loca-
tion and time period, by randomly retaining one of the closely
related sequences to be included in the analysis (Table S1, column
"Samples_used_in_Simulation_Analysis").

To robustly measure differences between simulated and ob-
served trees we use the matrix of the TMRCA for all pairs of samples.
This matrix also captures other allele frequency-based quantities
frequently used as summary statistics with ABC, such as F<sub>ST</sub>, as they
can be calculated from the components of this matrix.

In principle the full matrix could be used, but in practice it is
necessary to use a small number of summary statistics for ABC to
work properly (Wegmann et al., 2010). To this end, we computed the
mean TMRCA between pairs of sequences either within or between
(a) Europe, (b) Middle East, (c) North East Eurasia, Beringia and East
Eurasia combined; and (d) Artic and Continental North America com-
bined. This strategy is based on geographical proximity and genetic
similarity in the data set. We note that this is not the same as mod-
elling the combined demes as a single panmictic deme; structure
between the demes is still modelled explicitly, but the summary sta-
tistics are averaged over multiple demes.

An initial round of fitting the model showed that all scenarios underestimate the deme TMRCA for the Middle East, while the rest of the summary statistics were well captured by the best-fitting demographic scenarios. This could be explained by a scenario where the Middle East was less affected by the reduction in population size during the LGM. However, we currently lack a sufficient number of samples from this area to explicitly test a more complex scenario such as this hypothesis. To avoid outliers biasing the likelihood calculations in ABC (Wegmann et al., 2010) we removed this summary statistic, resulting in nine summary statistics in total.

For each of the 16 scenarios we performed 1 billion simulations
with randomly chosen parameter combinations, chosen from the fol-
lowing parameter intervals for the different scenarios:

* The static scenario: *m* in [0.001,20] and *K* in [0.01,100].

* The bottleneck scenarios: *m* in [0.001,20] and *K*₁, *K*₂, *K*₃ in [0.01,100].

* The expansion scenarios: *m* in [0.001,20], *K* in [0.01,100], *x* in [0,1], *T* in [5,40] and *ΔT* in [0.001,1]. For expansion out of the North American scenario and expansion out of the Arctic North American scenario, glaciation during the LGM in North American and sea-level rise during the deglaciation mean that *T* must be in the range [9,16]

* The combined bottleneck and expansion scenarios: *m* in [0.001,20], *K*₁, *K*₂, *K*₃ in [0.01,100], *x* in [0,1], *T* in [5,40] and *ΔT* in [0.001,1].

The parameter *m* is measured in units of 1/1,000 years, and *T*, *ΔT*, *K*, *K*₁, *K*₂ and *K*₃ are measured in units of 1,000 years. The parameters *x*, *T* and *ΔT* were sampled according to a uniform distribution over the interval, while all other parameters were sampled from a uniform distribution of their log-transformed values. To identify good parameter combinations for ABC, we first calculated the Euclidian square distances between predicted and observed statistics and restricted analysis to parameter combinations within the lowest tenth distance percentile. We then ran the ABCTOOLBOX (Wegmann et al., 2010) on the accepted parameter combinations to estimate posterior distributions of the model parameters, and to calculate the likelihood of each scenario as described in the ABCTOOLBOX manual.

Table S6 provides ABC likelihoods and Bayes factors (BFs) for all demographic scenarios tested. Tables S7 and S8 give posterior probability estimates and Figures S13 and S14 give posterior density distributions for estimated parameters (ΔT, *T*, log₁₀*K*₁, log₁₀*K*₂, log₁₀*K*₃, log₁₀*m*, *x*) in the two most likely models (an expansion out of Beringia with a population size change and an expansion out of East Eurasia with a population size change).

**2.9 | Map plots**

The background map used in Figures 1(a) and 3(a), showing climatic regions on land masses, was generated by downloading the file color_etopo1_ice_low.jpg from ETOPO1 (Amante & Eakins, 2016), a one arc-minute global relief model of the Earth's surface that integrates land topography and ocean bathymetry, and masking out regions where sea depths are greater than 120 m.

**3 | RESULTS**

**3.1 | Population structure of grey wolf across the Northern Hemisphere**

Motivated by the population structure observed in whole genome
studies of modern wolves (Fan et al., 2016), we tested the de-
gree of spatial genetic structure among the modern wolf samples
in our data set, and found a strong pattern of genetic IBD across
Eurasia (ρ = 0.3, *p* < .0001; see Figure S8). Ignoring this population
structure (i.e., modelling wolves as a single panmictic population)
can lead to artefactual results (Mazet, Rodríguez, & Chikhi, 2015;
Mazet, Rodríguez, Grusea, Boitard, & Chikhi, 2016). The use of spa-
tially structured models, in which migration is restricted to adjacent
populations, is a common approach for dealing with such situations (Eriksson et al., 2012; Eriksson & Manica, 2012; Kimura & Weiss, 1964; Wegmann et al., 2010).

To capture the observed geographical structure in our data set, we split the Northern Hemisphere into seven regions, roughly similar in area (Figure 3a). The boundaries of these regions are defined by geographical features, including mountain ranges, seas, and deserts (see Materials and Methods), which are likely to reduce gene flow (Geffen, Anderson, & Wayne, 2004; Lucchini, Galov, & Randi, 2004) and provide an optimal balance between resolution and power given the distribution of samples available for analyses. To quantify how well this scheme represents population structure in modern wolves, we used an AMOVA to separate genetic variance within and between regions. Our regions capture 24.4% of the genetic variation among our modern samples (AMOVA, *p* < .001). This is substantially greater than the ~10% of variance deriving from simple IBD, and supports the hypothesis that the geographical features (major rivers, deserts and mountain chains) define population structure in contemporary wolves across the Northern Hemisphere and therefore constitute obstacles to gene flow (but where the strength of these obstacles may vary).

## 3.2 | Bayesian phylogenetic analysis

All ancient sequences included in the study were subjected to stringent quality criteria with respect to coverage and damage patterns. Of the 45 ancient samples, 38 had well-resolved direct radiocarbon dates. We joined these ancient sequences with 90 modern mitogenome sequences and used BEAST (Drummond et al., 2012) to estimate a wolf mitochondrial mutation rate. By applying the inferred mutation rate we were able to molecularly date the remaining seven ancient sequences (Materials and Methods). We cross-validated this approach through a leave-one-out analysis (Materials and Methods) using all the directly dated ancient sequences and found a very close fit (R² = 0.86) between the radiocarbon and the estimated molecular dates and no systematic biases in our molecularly estimated dates (Figure S9), meriting the inclusion of these sequences and the inferred dates into the spatially explicit analyses.

Our Bayesian phylogenetic analysis suggests that the MRCA of all extant North Eurasian and American wolf mitochondrial sequences dates to ~40,000 years ago, whereas the MRCA for the combined ancient and modern sequences dates to ~90,000 years ago (95% highest posterior density [HPD] interval: 82,000–99,000 years ago) (Figure 2a, see Figures S11 and S12 for node support values and credibility intervals). A divergent clade at the root of this tree consists exclusively of ancient samples from Europe and the Middle East that has not contributed to present-day mitochondrial diversity in our data (see also Thalmann et al., 2013).

The remainder of the tree consists of a monophyletic clade that is made up of ancient and modern samples from across the Northern Hemisphere that shows a pattern of rapid bifurcations of genetic lineages centred on 25,000 years ago. To further quantify this temporal pattern, we made use of a Bayesian skyline analysis (Figure 2b)

that shows a relatively small and stable effective genetic population size between ~20,000 years ago and the present and a decrease in effective population size between ~40,000 and 20,000 years ago. This pattern is consistent with the scenario suggested in whole genome studies (e.g. Fan et al., 2016; Freedman et al., 2014) where wolves had a stable (and probably geographically structured) population across the Northern Hemisphere up to a time point between 20,000 and 30,000 years ago, when the population experienced a bottleneck that severely reduced genetic variation followed by a rapid population expansion.

The samples at the root of this clade are predominantly from Beringia, pointing to a possible expansion out northeast Eurasia or the Americas. However, given the uneven temporal and geographical distribution of our samples, and the stochasticity of a single genetic marker (Nielsen & Beaumont, 2009), it is important to explicitly test the extent to which this pattern can occur by chance under other plausible demographic scenarios.

## 3.3 | Spatiotemporal reconstruction of past grey wolf demography

Having established the phylogenetic relationship between our samples and population structure across the Northern Hemisphere, we tested the ability of different explicit demographic scenarios to explain the observed phylogenetic pattern, while also taking into account the geographical location and age of each sample. To this end, we represented each of the regions in Figure 3(a) as a population in a network of populations connected by gene flow (Figure 3b). We used the coalescent population genetic framework to model genetic evolution in this network, in which each deme constitutes a freely mixing and randomly mating population. The effective population size of demes, as well as movement of individuals between demes, are controlled by parameters covering values that represent different demographic histories.

Using this framework we considered a wide range of different explicit demographic scenarios (illustrated in Figure 3a, see Materials and Methods for details of implementation within the coalescent framework). The first scenario consisted of a constant population size and uniform movement between neighbouring demes. This allowed us to test the null hypothesis that drift within a structured population alone can explain all the patterns observed in the mitochondrial tree. We then considered two additional demographic processes that could explain the observed patterns: (a) a temporal sequence of two population size changes that affected all demes simultaneously (thus allowing for a bottleneck); and (b) an expansion out of one of the seven demes. In the expansion scenarios, the deme of origin had a continuous population through time, while in the remaining demes the indigenous population was sequentially replaced by the expanding population. Scenario 2 was repeated for all seven possible expansion origins, thus allowing us to test continuity as well as replacement hypotheses within each of the seven demes. We considered each demographic event in isolation as well
<figure><img src="image_4.png" /></figure>

**FIGURE 2** (a) Tip calibrated BEAST tree of all samples used in the spatial analyses (diamonds), coloured by geographical region. The circle represents an outgroup (modern Indian wolf, not used in the analyses). (b) The effective population size through time from the BEAST analysis (Bayesian skyline plot). The solid blue line represents the median estimate and the grey lines represent the interquartile range (sold lines) and 95% intervals (dashed lines) [Colour figure can be viewed at wileyonlinelibrary.com]

<figure><img src="image_5.png" /></figure>

**FIGURE 3** (a) Sample locations and geographical regions, with boundaries indicated by dashed lines. The dark blue indicates sea levels shallow enough to be land during the Last Glacial Maximum (sea depth < 120 m). (b) Model network of populations ("demes"), connected by gene flow, corresponding to the regions in (a)

<figure><img src="image_6.png" /></figure>

as their combined effect (resulting in a total of 16 scenarios) and used ABC to calculate the likelihood of each scenario and estimate parameter values (see Materials and Methods for details).

Both the null scenario and the scenario of only population size change in all demes were strongly rejected (BF ≤ 0.1, Figure 4b; Table S6), illustrating the power of combining a large data set of ancient samples with statistical modelling. Scenarios that combined an expansion and replacement with a change in population size (bottleneck) were better supported than the corresponding scenarios (i.e., with the same expansion origin) with constant population size (Figure 4b).

The best-supported scenario (Figure 5) was characterized by the combination of a rapid expansion of wolves out of the Beringian deme ~ 25,000 years ago (95% confidence interval [CI]: 33,000–14,000 years ago) with a population bottleneck between 15,000 and 40,000 years ago, and limited gene flow between neighbouring demes (see Table S7 and Figure S13 for posterior distributions of all model parameters). We also found relatively strong support for a scenario that describes a wolf expansion out of the East Eurasian deme (BF 0.7) with nearly identical parameters to the best-supported scenario (Table S8 and Figure S14). This can be explained by the geographical proximity of East Eurasian and
<figure><img src="image_7.png" /><figcaption>**FIGURE 4** Spatially and temporally explicit analysis. (a) Illustration of the different scenarios, with circles representing one deme each for the seven different geographical regions (see panel b for colour legend and text for full description of the scenarios). Solid lines represent population connectivity. The *static* scenario (far left) shows stable populations through time. The *expansion* scenario (middle left) shows how one deme (here yellow) expands and sequentially replaces the populations in all other demes (from top to bottom). The *population size change* scenario (middle right) illustrates how population size in the demes can change through time (large or small population size shown as large or small circles, respectively). We also show a combined scenario (far right) of both expansion and population size change. (b) Likelihood of each demographic scenario relative to the most likely scenario, shown as Bayes factors, estimated using Approximate Bayesian Computation analyses (see text for details). For expansion scenarios (including the combined expansion and population size changes), we colour code each bar according to the origin of the expansion (see colour legend)</figcaption></figure>

Beringian demes and the genetic similarity of wolves from these
areas.

4 | DISCUSSION

4.1 | Geographical origin of the ancestral wolf population

Recent whole-genome studies (Fan et al., 2016; Freedman et al., 2014; Skoglund et al., 2015) found that modern grey wolves (*Canis lupus*) across Eurasia are descended from a single source popula-
tion. The results of our analyses combining both ancient and mod-
ern grey wolf samples (Figure 1) with a spatially and temporally
explicit modelling framework (Figure 4) suggest that this process
began ~25,000 (95% CI: 33,000–14,000) years ago when a pop-
ulation of wolves from Beringia (or a Northeast Asian region in
close geographical proximity) expanded outwards and replaced
indigenous Pleistocene wolf populations across Eurasia (Figure 5).
This scenario also provides a mechanism explaining the star-like
topology of modern wolf populations observed in whole genome
studies (Fan et al., 2016; Freedman et al., 2014; Skoglund et al.,

2015): the expansion was split up by geographical barriers that
restricted subsequent gene flow between different branches of
the expanding population, which in turn led to the divergence be-
tween different subpopulations observed in contemporary grey
wolves.

In the Americas, the Beringian expansion was delayed due
to the presence of ice sheets extending from Greenland to the
northern Pacific Ocean (Figure 5) (Raghavan et al., 2015). A study
by Koblmüller et al. (2016) suggested that wolf populations that
were extant south of these ice sheets were replaced by Eurasian
wolves crossing the Beringian land bridge. Our data and analyses
support the replacement of North American wolves (following re-
treat of the ice sheets around 16,000 years ago), and our more
extensive ancient DNA sampling, combined with spatially explicit
modelling, has allowed us to narrow down the geographical origin
of this expansion to an area between the Lena River in Russia and
the Mackenzie River in Canada also known as Beringia (Hopkins,
Matthews, & Schweger, 1982). However, due to lack of Pleistocene
wolf samples that pre-date the retreat of the ice sheets in the
area, we are currently not able to resolve the detailed history of
North American wolves. For example, we cannot reject an alter-
native scenario where contemporary North American wolves are
<figure><img src="image_8.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_11.png" /></figure>

**FIGURE 5** The inferred scenario of wolf demography from the Bayesian analysis using our spatially and temporally explicit model (see Figure 4 and the main text). (a) Geographical representation of the expansion scenario (out of Beringia) with median and 95% confidence interval (CI) for the date of the population replacement in each deme given in white boxes next to each deme. (b) Effective population size (thick line, boxes and whiskers show the median, interquartile range and 95% CI, respectively, for each time period). (c) Posterior distribution of migration rate and (d) starting time of expansion [Colour figure can be viewed at wileyonlinelibrary.com]

descendants of a Pleistocene wolf population that was genetically highly similar to the Beringian population but existed south of the ice sheets.

Thus, despite a continuous fossil record through the Late Pleistocene, wolves experienced a complex demographic history involving population bottlenecks and replacements (Figure 5). Our analysis suggests that long-range migration played an important role in the survival of wolves through the wave of megafaunal extinctions at the end of the last glaciation. These results will enable future studies to examine specific local climatic and ecological factors that enabled the Beringian wolf population to survive and expand across the Northern Hemisphere. Furthermore, as the reconstructions in this study are based solely on a maternally inherited genetic marker, our model was thus only able to address a set of simplified demographic scenarios (continuity everywhere, or continuity in one location followed by a replacement expansion from it). Once whole-genome data become available, it will probably be possible to detect contributions from potential refugia at the local scale.

## 4.2 | Implications for the evolution of grey wolf morphology

Morphological analyses of wolf specimens have noted differences between Late Pleistocene and Holocene wolves: Late Pleistocene specimens have been described as craniodentally more robust than present-day grey wolves, as well as having specialized adaptations for carcass and bone processing (Baryshnikov, Mol, & Tikhonov, 2009; Kuzmina & Sablin, 1993; Leonard et al., 2007) associated with megafaunal hunting and scavenging (Fox-Dobbs, Leonard, & Koch, 2008; Germonpré et al., 2017). The early Holocene archaeological record has only yielded a single sample with the Pleistocene wolf morphotype (in Alaska) (Leonard et al., 2007), suggesting that this robust ecomorph had largely disappeared from the Northern Hemisphere by the Pleistocene-Holocene transition. This change in wolf morphology coincides with a shift in wolf isotope composition (Bocherens, 2015), and the disappearance of megafaunal herbivores and other large predators such as cave hyenas and cave lions, suggesting a possible change in the ecological niche of wolves.
To date, it has been unclear whether the morphological change was the result of population replacement (genetic turnover), a plastic response to a dietary shift, or both. Our results suggest that the Pleistocene-Holocene transition was accompanied by a genetic turnover in most of the Northern Hemisphere wolf populations as most indigenous wolf populations experienced a large-scale replacement resulting in the loss of all native Pleistocene genetic lineages (Figure 5). Similar population dynamics of discontinuity and replacement by con specifics have been observed in several other large Pleistocene mam mals in Europe including cave bears, woolly mammoths (Palkopoulou et al., 2013; Stuart, Kosintsev, Higham, & Lister, 2004), giant deer (Stuart et al., 2004) and even humans (Fu et al., 2016; Posth et al., 2016).

The geographical exception to this pattern of widespread re-
placement is Beringia, where we infer demographic continuity be-
tween Late Pleistocene and Holocene wolf populations (Figure 5).
This finding is at odds with a previous suggestion of genetic turnover
in Beringia (Leonard et al., 2007), probably as the result of differences
in both the amount of data available and the analytical methodology
used. Leonard et al. (2007) used a short (427 bases long) segment
of the mitochondrial control region and employed a descriptive phy-
logeographical approach, whereas our conclusions are based on an
expanded data set in terms of both sequence length, sample number,
and geographical and temporal range (Figure 1) and formal hypothe-
sis testing within a Bayesian framework (Figures 4 and 5).

As a consequence, the morphological and dietary shift observed
in Beringian wolves between the Late Pleistocene and Holocene
(Leonard et al., 2007) cannot be explained by population turnover,
but instead requires an alternative explanation such as adaptation
or plastic responses to the substantial environmental and ecologi-
cal changes that took place during this period. Indeed, grey wolves
are a highly adaptable species. Studies of modern grey wolves have
found that differences in habitat, specifically precipitation, tempera-
ture, vegetation and prey specialization, can strongly affect their
craniodental morphology (Flower & Schreve, 2014; Geffen et al.,
2004; Leonard, 2015; O'Keefe, Meachen, Fet, & Brannick, 2013;
Pilot et al., 2006).

The specific causal factors for the replacement of indigenous Eurasian wolves during the LGM by their Beringian conspecifics (and American wolves following the disappearance of the Cordilleran and Laurentide ice sheets) are beyond the scope of this study. However, one possible explanation may be related to the relatively stable climate of Beringia compared to the substantial climatic fluctuations that impacted the rest of Eurasia and Northern America during the Late Pleistocene (Clark et al., 2012). These fluctuations have been associated with dramatic changes in food webs, leading to the loss of most of the large Pleistocene predators in the region (Bocherens, 2015; Hofreiter & Stewart, 2009; Lister & Stuart, 2008; Lorenzen et al., 2011). In addition, the hunting of large Pleistocene predators by late Palaeolithic people (e.g. Cueto, Camarós, Castaños, Ontañón, & Arias, 2016; Germonpré & Hämäläinen, 2007; Münzel & Conard, 2004) may have also negatively impacted large carnivore populations (Fan et al., 2016). An interdisciplinary approach involving morphological, isotopic as well as genetic data is necessary to better understand

the relationship between wolf population dynamics and dietary ad-
aptations in the Late Pleistocene and early Holocene period.

## 4.3 | Implications for the study of wolf domestication

Lastly, the complex demographic history of Eurasian grey wolves re-
ported here (Figure 5) also has significant implications for identifying
the geographical origin(s) of wolf domestication and the subsequent
spread of dogs. For example, our limited understanding of the un-
derlying wolf population structure may explain why previous stud-
ies have produced conflicting geographical and temporal scenarios.
Numerous previous studies have focused on the patterns of genetic
variation in modern domestic dogs, but have failed to consider po-
tential genetic variation present in Late Pleistocene wolf populations,
thereby implicitly assuming a homogeneous wolf population source.
As a result, both the domestication and the subsequent human-medi-
ated movements of dogs were the only processes considered to have
affected the observed genetic patterns in dog populations. However,
both domestication from and admixture with a structured wolf
population will have consequences for patterns of genetic variation
within dogs. In light of the complex demographic history of wolves
(and the resulting population genetic structure) reconstructed by our
analysis, several of the geographical patterns of haplotype distribu-
tion observed in previous studies, including differences in levels of
diversity found within local dog populations (Wang et al., 2016), and
the deep phylogenetic split between Eastern and Western Eurasian
dogs (Frantz et al., 2016), could have resulted from known admixture
between domestic dogs and grey wolves (Fan et al., 2016; Freedman
et al., 2014; Godinho et al., 2011; Verardi, Lucchini, & Randi, 2006).
Future analyses should therefore explicitly include the demographic
history of wolves and demonstrate that the patterns of variation ob-
served within dogs fall outside expectations that take admixture with
geographically structured wolf populations into account.

ACKNOWLEDGEMENTS

We are grateful to Daniel Klingberg Johansson & Kristian Murphy
Gregersen from the Natural History Museum of Denmark; Gabriella
Hürlimann from the Zurich Zoo; Jane Hopper from the Howlett's
& the Port Lympne Wild Animal Parks; Cyrintha Barwise-Joubert
& Paul Vercammen from the Breeding Centre for Endangered
Arabian Wildlife; Link Olson from the University of Alaska Museum
of the North; Joseph Cook & Mariel Campbell from the Museum of
Southwestern Biology; Lindsey Carmichael & David Coltman from
the University of Alberta; North American Fur Auctions; Department
of Environment Nunavut and Environment and Natural Resources
Northwest Territories for DNA samples from the modern wolves. We
are also grateful to the staff at the Danish National High-Throughput
Sequencing Centre for technical assistance in the data generation;
the Qimmeq project, funded by The Velux Foundations and Aage og
Johanne Louis-Hansens Fond, for providing financial support for se-
quencing ancient Siberian wolf samples; the Rock Foundation (New
York, USA) for supporting radiocarbon dating of ancient samples from the Yana site; to Stephan Nylinder from the Swedish Museum of Natural History for advice on phylogenetic analyses and Terry Brown from the University of Manchester for comments on the manuscript. L.L., K.D. and G.L. were supported by the Natural Environment Research Council, UK (grant numbers NE/K005243/1, NE/K003259/1); LL was also supported by the European Research Council grant (339941-ADAPT); A.M. and A.E. were supported by the European Research Council Consolidator grant (grant number 647787-LocalAdaptation); L.F. and G.L. were supported by the European Research Council grant (ERC-2013-StG 337574-UNDEAD); T.G. was supported by a European Research Council Consolidator grant (681396-Extinction Genomics) & Lundbeck Foundation grant (R52-5062); O.T. was supported by the National Science Center, Poland (2015/19/P/NZ7/03971), with funding from EU's Horizon 2020 programme under the Marie Skłodowska-Curie grant agreement (665778) and Synthesis Project (BETAF 3062); V.P., E.P. and P.N. were supported by the Russian Science Foundation grant (N16-18-10265 RNF); A.P. was supported by the Max Planck Society; M.L-G. was supported by a Czech Science Foundation grant (GAČR15-06446S).

## AUTHOR CONTRIBUTIONS

L.L., O.T., M.T.P.G., J.K., G.L., A.E. and A.M. designed the research; O.T., M-H.S.S., V.J.S., K.E.W., M.S.V., I.K.C.L., N.W. and G.S. performed ancient DNA laboratory work with input from J.K., M.T.P.G., H.S., K-H.H., and R.S.M.; M-H.S.S. performed modern DNA laboratory work with input from M.T.P.G; O.T., J.A.S.C. and L.L. performed bioinformatic analyses; L.L., A.E. and A.M. designed the population genetic analyses; L.L. performed phylogenetic analyses; A.E. implemented the spatial analyses framework; L.L and A.E. performed spatial analyses; M.G., J.B., V.V.P., E.Y.P., P.A.N., S.E.F., J.E-L., A.W.K., B.G., H.N., H-P.U. and M.L-G. provided samples; V.V.P., M.G., M.L-G., H.B., H.N., A.W.K., E.Y.P. and P.A.N. provided context for archaeological samples; A.P., M.G., H.B. and K.D. helped setting the results of genetic analyses into an archaeological context; A.M., M.T.P.G., A.J.H., G.L., J.K., E.W. and K.D. secured funding for the project; L.L., O.T. and A.E. wrote the initial draft of the manuscript with input from A.M.; L.L., O.T. and A.E wrote the manuscript and the supplementary information with input from A.P., M.G., H.B., M-H.S.S., M.T.P.G., K.E.W., A.M., G.L and K.D.; V.J.S., L.F., A.W.K., K-H.H., A.J.H., R.S.M., H.S., G.S., V.V.P., E.Y.P., P.A.N. and J.E-L. provided comments on the manuscript and/or on the supplementary information.

## DATA AVAILABILITY STATEMENT

The newly assembled mitochondrial genomes are available from GenBank (accession numbers MK936995-MK937053 [ancient] and MN071185-MN071206 [modern]). The raw sequencing reads used for generating novel ancient mitochondrial genomes can be retrieved from the European Nucleotide Archive under study number PRJEB32023. The code for population genetic simulations of all tested scenarios and scripts for preliminary and output analyses are available on the GitHub repository at https://github.com/LiisaLoog/pleistocene-wolves.

## ORCID

<figure><img src="image_12.png" /></figure>

Liisa Loog https://orcid.org/0000-0002-1770-101X

Mikkel-Holger S. Sinding https://orcid.org/0000-0003-1371-219X

<figure><img src="image_13.png" /></figure>

Gontran Sonet https://orcid.org/0000-0001-7310-9574

<figure><img src="image_14.png" /></figure>

Anders Eriksson https://orcid.org/0000-0003-3436-3726

## REFERENCES

Aggarwal, R. K., Kivisild, T., Ramadevi, J., & Singh, L. (2007). Mitochondrial DNA coding region sequences support the phylogenetic distinction of two Indian wolf species. *Journal of Zoological Systematics and Evolutionary Research*, 45, 163–172. https://doi.org/10.1111/j.1439-0469.2006.00400.x

Amante, C., & Eakins, B. W. (2016). *ETOP01 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis*. NOAA Technical Memorandum NESDIS NGDC-24. Retrieved from https://www.ngdc.noaa.gov/mgg/global/

Barnosky, A. D., Koch, P. L., Feranec, R. S., Wing, S. L., & Shabel, A. B. (2004). Assessing the causes of late Pleistocene extinctions on the continents. *Science*, 306, 70–75. https://doi.org/10.1126/science.1101476

Baryshnikov, G. F., Mol, D., & Tikhonov, A. N. (2009). Finding of the Late Pleistocene carnivores in Taimyr Peninsula (Russia, Siberia) with paleoecological context. *Russian Journal of Theriology*, 8, 107–113. https://doi.org/10.15298/rusjtheriol.08.2.04

Beaumont, M. A., Zhang, W., & Balding, D. J. (2002). Approximate Bayesian Computation in population genetics. *Genetics*, 162, 2025–2035.

Bocherens, H. (2015). Isotopic tracking of large carnivore palaeoecology in the mammoth steppe. *Quaternary Science Reviews*, 117, 42–71. https://doi.org/10.1016/j.quascirev.2015.03.018

Clark, P. U., Shakun, J. D., Baker, P. A., Bartlein, P. J., Brewer, S., Brook, E., ... Williams, J. W. (2012). Global climate evolution during the last deglaciation. *Proceedings of the National Academy of Sciences*, 109, E1134–E1142. https://doi.org/10.1073/pnas.1116619109

Crockford, S. J., & Kuzmin, Y. V. (2012). Comments on Germonpré et al., Journal of Archaeological Science 36, 2009 "Fossil dogs and wolves from Palaeolithic sites in Belgium, the Ukraine and Russia: osteometry, ancient DNA and stable isotopes", and Germonpré, Lázkičková-Galetová, and Sablin, Journal of Archaeological Science 39, 2012 "Palaeolithic dog skulls at the Gravettian Předmostí site, the Czech Republic." *Journal of Archaeological Science*, 39, 2797–2801.

Cueto, M., Camarós, E., Castaños, P., Ontañón, R., & Arias, P. (2016). Under the skin of a lion: Unique evidence of upper Paleolithic exploitation and use of cave lion (*Panthera spelaea*) from the Lower Gallery of La Garma (Spain). *PLoS ONE*, 11, e0163591. https://doi.org/10.1371/journal.pone.0163591

Drake, A. G., Coquerelle, M., & Colombeau, G. (2015). 3D morphometric analysis of fossil canid skulls contradicts the suggested domestication of dogs during the late Paleolithic. *Scientific Reports*, 5, 8299. https://doi.org/10.1038/srep08299

Drummond, A. J., Nicholls, G. K., Rodrigo, A. G., & Solomon, W. (2002). Estimating mutation parameters, population history and genealogy simultaneously from temporally spaced sequence data. *Genetics*, 161, 1307–1320.

Drummond, A. J., Suchard, M. A., Xie, D., & Rambaut, A. (2012). Bayesian phylogenetics with BEAUti and the BEAST 1.7. *Molecular Biology and Evolution*, 29, 1969–1973. https://doi.org/10.1093/molbev/mss075

Druzhkova, A. S., Thalmann, O., Trifonov, V. A., Leonard, J. A., Vorobieva, N. V., Ovodov, N. D., ... Wayne, R. K. (2013). Ancient DNA analysis affirms the canid from Altai as a primitive dog. *PLoS ONE*, 8, e57754. https://doi.org/10.1371/journal.pone.0057754

Eriksson, A., Betti, L., Friend, A. D., Lycett, S. J., Singarayer, J. S., von Cramon-Taubadel, N., ... Manica, A. (2012). Late Pleistocene climate
change and the global expansion of anatomically modern humans. *Proceedings of the National Academy of Sciences*, **109**, 16089–16094.
https://doi.org/10.1073/pnas.1209494109

Eriksson, A., & Manica, A. (2012). Effect of ancient population structure on the degree of polymorphism shared between modern human populations and ancient hominins. *Proceedings of the National Academy of Sciences*, **109**, 13956–13960. https://doi.org/10.1073/pnas.1200567109

Eriksson, A., & Mehlig, B. (2004). Gene-history correlation and population structure. *Physical Biology*, **1**, 220. https://doi.org/10.1088/1478-3967/1/4/004

Excoffier, L., Smouse, P. E., & Quattro, J. M. (1992). Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. *Genetics*, **131**, 479–491.

Excoffier, L., & Yang, Z. (1999). Substitution rate variation among sites in mitochondrial hypervariable region I of humans and chimpanzees. *Molecular Biology and Evolution*, **16**, 1357–1368. https://doi.org/10.1093/oxfordjournals.molbev.a026046

Fan, Z., Silva, P., Gronau, I., Wang, S., Armero, A. S., Schweizer, R. M., ... Wayne, R. K. (2016). Worldwide patterns of genomic variation and admixture in gray wolves. *Genome Research*, **26**, 163–173. https://doi.org/10.1101/gr.197517.115

Flower, L. O. H., & Schreve, D. C. (2014). An investigation of palaeodietary variability in European Pleistocene canids. *Quaternary Science Reviews*, **96**, 188–203. https://doi.org/10.1016/j.quascirev.2014.04.015

Fox-Dobbs, K., Leonard, J. A., & Koch, P. L. (2008). Pleistocene megafauna from eastern Beringia: Paleoeological and paleoenvironmental interpretations of stable carbon and nitrogen isotope and radio-carbon records. *Palaeogeography, Palaeoclimatology, Palaeoecology*, **261**, 30–46. https://doi.org/10.1016/j.palaeo.2007.12.011

Frantz, L. A. F., Mullin, V. E., Pionnier-Capitan, M., Lebrasseur, O., Ollivier, M., Perri, A., ... Larson, G. (2016). Genomic and archaeological evidence suggest a dual origin of domestic dogs. *Science*, **352**, 1228. https://doi.org/10.1126/science.aaf3161

Freedman, A. H., Gronau, I., Schweizer, R. M., Ortega-Del Vecchio, D., Han, E., Silva, P. M., ... Novembre, J. (2014). Genome sequencing highlights the dynamic early history of dogs. *PLoS Genetics*, **10**, e1004016. https://doi.org/10.1371/journal.pgen.1004016

Fu, Q., Posth, C., Hajdinjak, M., Petr, M., Mallick, S., Fernandes, D., ... Reich, D. (2016). The genetic history of Ice Age Europe. *Nature*, **534**, 200–205. https://doi.org/10.1038/nature17993

Geffen, E., Anderson, M. J., & Wayne, R. K. (2004). Climate and habitat barriers to dispersal in the highly mobile grey wolf. *Molecular Ecology*, **13**, 2481–2490. https://doi.org/10.1111/j.1365-294X.2004.02244.x

Germonpré, M., Fedorov, S., Danilov, P., Galeta, P., Jimenez, E.-L., Sablin, M., & Losey, R. J. (2017). Palaeolithic and prehistoric dogs and Pleistocene wolves from Yakutia: Identification of isolated skulls. *Journal of Archaeological Science*, **78**, 1–19. https://doi.org/10.1016/j.jas.2016.11.008

Germonpré, M., & Hämäläinen, R. (2007). Fossil bear bones in the Belgian Upper Paleolithic: The possibility of a proto bear-ceremonialism. *Arctic Anthropology*, **44**, 1–30. https://doi.org/10.1353/arc.2011.0015

Germonpré, M., Lázničková-Galetová, M., Losey, R. J., Räikkönen, J., & Sablin, M. V. (2015). Large canids at the Gravettian Předmostí site, the Czech Republic: The mandible. *Quaternary International*, **359–360**, 261–279. https://doi.org/10.1016/j.quaint.2014.07.012

Germonpré, M., Lázničková-Galetová, M., & Sablin, M. V. (2012). Palaeolithic dog skulls at the Gravettian Předmostí site, the Czech Republic. *Journal of Archaeological Science*, **39**, 184–202. https://doi.org/10.1016/j.jas.2011.09.022

Germonpré, M., Sablin, M. V., Stevens, R. E., Hedges, R. E. M., Hofreiter, M., Stiller, M., & Després, V. R. (2009). Fossil dogs and wolves from

Palaeolithic sites in Belgium, the Ukraine and Russia: Osteometry, ancient DNA and stable isotopes. *Journal of Archaeological Science*, **36**, 473–490. https://doi.org/10.1016/j.jas.2008.09.033

Godinho, R., Llaneza, L., Blanco, J. C., Lopes, S., Álvares, F., García, E. J., ... Ferrand, N. (2011). Genetic evidence for multiple events of hybridization between wolves and domestic dogs in the Iberian Peninsula. *Molecular Ecology*, **20**, 5154–5166. https://doi.org/10.1111/j.1365-294X.2011.05345.x

Gopalakrishnan, S., Sinding, M.-H., Ramos-Madrigal, J., Niemann, J., Samaniego Castruita, J. A., Vieira, F. G., ... Gilbert, M. T. P. (2018). Interspecific gene flow shaped the evolution of the genus *Canis*. *Current Biology*, **28**, 3441–3449.e5. https://doi.org/10.1016/j.cub.2018.08.041

Groucutt, H. S., Petraglia, M. D., Bailey, G., Scerri, E. M. L., Parton, A., Clark-Balzan, L., ... Breeze, P. S. (2015). Rethinking the dispersal of *Homo sapiens* out of Africa. *Evolutionary Anthropology*, **24**, 149–164.

Hofreiter, M., & Stewart, J. (2009). Ecological change, range fluctuations and population dynamics during the Pleistocene. *Current Biology*, **19**, R584–R594.

Hopkins, D., Matthews, J., & Schweger, C. (1982). *Paleoecology of Beringia*, 1st ed. New York, NY: Academic Press.

Kimura, M., & Weiss, G. H. (1964). The stepping stone model of population structure and the decrease of genetic correlation with distance. *Genetics*, **49**, 561–576.

Kingman, J. F. C. (1982). The coalescent. *Stochastic Processes and their Applications*, **13**, 235–248. https://doi.org/10.1016/0304-4149(82)90011-4

Koblmüller, S., Vilà, C., Lorente-Galdos, B., Dabad, M., Ramirez, O., Marques-Bonet, T., ... Leonard, J. A. (2007). Whole mitochondrial genomes illuminate ancient intercontinental dispersals of grey wolves (*Canis lupus*). *Journal of Biogeography*, **43**, 1728–1738.

Kuzmina, I. E., & Sablin, M. V. (1993). Pozdnepleistotsenovyi pesets verhnei Desny. In G. F. Baryshnikov, & I. E. Kuzmina (Eds.), *Materiali po mezoziockoi i kainozoickoi istorii nazemnykh pozvonochnykh* (pp. 93–104). Leningrad, USSR: Trudy 17 Zoologicheskogo Instituta RAN 249.

Lanfear, R., Calcott, B., Ho, S. Y. W., & Guindon, S. (2012). PartitionFinder: Combined selection of partitioning schemes and substitution models for phylogenetic analyses. *Molecular Biology and Evolution*, **29**, 1695–1701. https://doi.org/10.1093/molbev/mss020

Larkin, M. A., Blackshields, G., Brown, N. P., Chenna, R., McGettigan, P. A., McWilliam, H., ... Higgins, D. G. (2007). Clustal W and Clustal X version 2.0. *Bioinformatics*, **23**, 2947–2948. https://doi.org/10.1093/bioinformatics/btm044

Larson, G., Karlsson, E. K., Perri, A., Webster, M. T., Ho, S. Y. W., Peters, J., ... Lindblad-Toh, K. (2012). Rethinking dog domestication by integrating genetics, archeology, and biogeography. *Proceedings of the National Academy of Sciences*, **109**, 8878–8883. https://doi.org/10.1073/pnas.1203005109

Leonard, J. A. (2015). Ecology drives evolution in grey wolves. *Evolutionary Ecology Research*, **16**, 461–473.

Leonard, J. A., Vilà, C., Fox-Dobbs, K., Koch, P. L., Wayne, R. K., & Van Valkenburgh, B. (2007). Megafaunal extinctions and the disappearance of a specialized wolf ecomorph. *Current Biology*, **17**, 1146–1150. https://doi.org/10.1016/j.cub.2007.05.072

Lister, A. M., & Stuart, A. J. (2008). The impact of climate change on large mammal distribution and extinction: Evidence from the last glacial/interglacial transition. *Comptes Rendus Geoscience*, **340**, 615–620. https://doi.org/10.1016/j.crte.2008.04.001

Loog, L., Lahr, M. M., Kovacevic, M., Manica, A., Eriksson, A., & Thomas, M. G. (2017). Estimating mobility using sparse data: Application to human genetic variation. *Proceedings of the National Academy of Sciences*, **114**(46), 12213–12218. https://doi.org/10.1073/pnas.1703642114

Lorenzen, E. D., Nogués-Bravo, D., Orlando, L., Weinstock, J., Binladen, J., Marske, K. A., ... Willerslev, E. (2011). Species-specific responses
of Late Quaternary megafauna to climate and humans. *Nature*, 479, 359-364. https://doi.org/10.1038/nature10574

Lucchini, V., Galov, A., & Randi, E. (2004). Evidence of genetic distinction and long-term population decline in wolves (*Canis lupus*) in the Italian Apennines. *Molecular Ecology*, 13, 523-536. https://doi.org/10.1046/j.1365-294X.2004.02077.x

Mazet, O., Rodríguez, W., & Chikhi, L. (2015). Demographic inference using genetic data from a single individual: Separating population size variation from population structure. *Theoretical Population Biology*, 104, 46-58. https://doi.org/10.1016/j.tpb.2015.06.003

Mazet, O., Rodríguez, W., Grusea, S., Boitard, S., & Chikhi, L. (2016). On the importance of being structured: Instantaneous coalescence rates and human evolution—lessons for ancestral population size inference? *Heredity*, 116, 362-371. https://doi.org/10.1038/hdy.2015.104

Morey, D. F. (2014). In search of Paleolithic dogs: A quest with mixed results. *Journal of Archaeological Science*, 52, 300-307. https://doi.org/10.1016/j.jas.2014.08.015

Münzel, S. C., & Conard, N. J. (2004). Change and continuity in subsistence during the Middle and Upper Palaeolithic in the Ach Valley of Swabia (south-west Germany). *International Journal of Osteoarchaeology*, 14, 225-243. https://doi.org/10.1002/oa.758

Nielsen, R., & Beaumont, M. A. (2009). Statistical inferences in phylogeography. *Molecular Ecology*, 18, 1034-1047. https://doi.org/10.1111/j.1365-294X.2008.04059.x

O'Keefe, F. R., Meachen, J., Fet, E. V., & Brannick, A. (2013). Ecological determinants of clinal morphological variation in the cranium of the North American gray wolf. *Journal of Mammalogy*, 94, 1223-1236. https://doi.org/10.1644/13-MAMM-A-069

Palkopoulou, E., Dalén, L., Lister, A. M., Vartanyan, S., Sablin, M., Sher, A., ... Thomas, J. A. (2013). Holarctic genetic structure and range dynamics in the woolly mammoth. *Proceedings of the Royal Society B: Biological Sciences*, 280, 20131910. https://doi.org/10.1098/rspb.2013.1910

Paradis, E., Claude, J., & Strimmer, K. (2004). APE: Analyses of Phylogenetics and Evolution in R language. *Bioinformatics*, 20, 289-290. https://doi.org/10.1093/bioinformatics/btg412

Perri, A. (2016). A wolf in dog's clothing: Initial dog domestication and Pleistocene wolf variation. *Journal of Archaeological Science*, 68, 1-4. https://doi.org/10.1016/j.jas.2016.02.003

Pilot, M., Jedzejewski, W., Branicki, W., Sidorovich, V. E., Jedzejewska, B., Stachura, K., & Funk, S. M. (2006). Ecological factors influence population genetic structure of European grey wolves. *Molecular Ecology*, 15, 4533-4553. https://doi.org/10.1111/j.1365-294X.2006.03110.x

Posth, C., Renaud, G., Mittnik, A., Drucker, D. G., Rougier, H., Cupillard, C., ... Krause, J. (2016). Pleistocene mitochondrial genomes suggest a single major dispersal of non-Africans and a Late Glacial population turnover in Europe. *Current Biology*, 26, 827-833. https://doi.org/10.1016/j.cub.2016.01.037

Puzachenko, A. Y., & Markova, A. K. (2016). Diversity dynamics of large- and medium-sized mammals in the Late Pleistocene and the Holocene on the East European Plain: Systems approach. *Quaternary International*, 420, 391-401. https://doi.org/10.1016/j.quaint.2015.07.031

Raghavan, M., Steinrucken, M., Harris, K., Schiffels, S., Rasmussen, S., DeGiorgio, M., ... Willerslev, E. (2015). Genomic evidence for the Pleistocene and recent population history of Native Americans. *Science*, 349(6250), aab3884. https://doi.org/10.1126/science.aab3884

Rambaut, A. (2000). Estimating the rate of molecular evolution: Incorporating non-contemporaneous sequences into maximum likelihood phylogenies. *Bioinformatics*, 16, 395-399. https://doi.org/10.1093/bioinformatics/16.4.395

Rieux, A., Eriksson, A., Li, M., Sobkowiak, B., Weinert, L. A., Warmuth, V., ... Balloux, F. (2014). Improved calibration of the human mitochondrial

clock using ancient genomes. *Molecular Biology and Evolution*, 31, 2780-2792. https://doi.org/10.1093/molbev/msu222

Sablin, M. V., & Khlopachev, G. A. (2002). The earliest Ice Age dogs: Evidence from Eliseevichi 1. *Current Anthropology*, 43, 795-799. https://doi.org/10.1086/344372

Sharma, D. K., Maldonado, J. E., Jhala, Y. V., & Fleischer, R. C.. (2004). Ancient wolf lineages in India. *Proceedings of the Royal Society of London. Series B: Biological Sciences*, 271(Suppl_3), S1-S4. https://doi.org/10.1098/rsbl.2003.0071

Sinding, M.-H., Gopalakrishan, S., Vieira, F. G., Samaniego Castruita, J. A., Raundrup, K., Heide Jørgensen, M. P., ... Gilbert, M. T. P. (2018). Population genomics of grey wolves and wolf-like canids in North America. *PLoS Genetics*, 14, e1007745. https://doi.org/10.1371/journal.pgen.1007745

Sinnott, R. (1984). Virtues of the haversine. *Sky and Telescope*, 68, 159.

Skoglund, P., Ersmark, E., Palkopoulou, E., & Dalén, L. (2015). Ancient wolf genome reveals an early divergence of domestic dog ancestors and admixture into high-latitude breeds. *Current Biology*, 25, 1515-1519. https://doi.org/10.1016/j.cub.2015.04.019

Sotnikova, M., & Rook, L. (2010). Dispersal of the Canini (Mammalia, Canidae: Caninae) across Eurasia during the Late Miocene to Early Pleistocene. *Quaternary International*, 212, 86-97. https://doi.org/10.1016/j.quaint.2009.06.008

Stuart, A. J., Kosintsev, P. A., Higham, T. F. G., & Lister, A. M. (2004). Pleistocene to Holocene extinction dynamics in giant deer and woolly mammoth. *Nature*, 431, 684-689. https://doi.org/10.1038/nature02890

Thalmann, O., Shapiro, B., Cui, P., Schuenemann, V. J., Sawyer, S. K., Greenfield, D. L., ... Wayne, R. K. (2013). Complete mitochondrial genomes of ancient canids suggest a European origin of domestic dogs. *Science*, 342, 871-874. https://doi.org/10.1126/science.1243650

Verardi, A., Lucchini, V., & Randi, E. (2006). Detecting introgressive hybridization between free-ranging domestic dogs and wild wolves (*Canis lupus*) by admixture linkage disequilibrium analysis. *Molecular Ecology*, 15, 2845-2855. https://doi.org/10.1111/j.1365-294X.2006.02995.x

Wang, G.-D., Zhai, W., Yang, H.-C., Wang, L. U., Zhong, L. I., Liu, Y.-H., ... Zhang, Y.-P. (2016). Out of southern East Asia: The natural history of domestic dogs across the world. *Cell Research*, 26, 21-33. https://doi.org/10.1038/cr.2015.147

Warmuth, V., Eriksson, A., Bower, M. A., Barker, G., Barrett, E., Hanks, B. K., ... Manica, A. (2012). Reconstructing the origin and spread of horse domestication in the Eurasian steppe. *Proceedings of the National Academy of Sciences*, 109, 8202-8206. https://doi.org/10.1073/pnas.1111122109

Wegmann, D., Leuenberger, C., Neuenschwander, S., & Excoffier, L. (2010). ABCtoolbox: A versatile toolkit for approximate Bayesian computations. *BMC Bioinformatics*, 11, 116. https://doi.org/10.1186/1471-2105-11-116

## SUPPORTING INFORMATION

Additional supporting information may be found online in the Supporting Information section.

**How to cite this article:** Loog L, Thalmann O, Sinding M-HS, et al. Ancient DNA suggests modern wolves trace their origin to a Late Pleistocene expansion from Beringia. *Mol Ecol.* 2020;29:1596-1610. https://doi.org/10.1111/mec.15329