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

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.


| INTRODUC TI ON
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 geographical 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 radiation 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 understanding 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 history of dog domestication.
Reconstructing past demographic events solely from modern genomes 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 eliminate 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 samples 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 nuclear 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 allows 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 contractions, that involved space as well as time are particularly challenging 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 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 spatiotemporally explicit models. We conclude with a discussion of how our findings relate to earlier studies and implications for future research.

| 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 ancient 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 spectrometry [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 criteria 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 sequences that contained more than 0.1% of singletons showing signs of deamination damage typical for ancient DNA (C to T or A to G singletons). 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 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 form 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).  (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. F I G U R E 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] Sample age [1,000 years BP]  CGG16  TU14  TH14  TU4  TU8  TU5  CGG20  CGG21  TU13  TU11  TH11  TH15  TH7  TH5  CGG33  CGG19  TH10  TU15  TH1  TH8  CGG17  TH3  TH6  CGG34  TU10  TU9  CGG26  CGG22  CGG25  CGG28  TU3  CGG14  CGG27  CGG15  TU2  TH12  SK1  TU1  TU6  CGG18  TU7  TH4  CGG12 CGG32 CGG29 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.

| 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 neighbour- 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 radiocarbon-dated samples) or as the mean of the age distribution from the beast analyses (for molecularly dated samples).

| 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: where G is the distance in kilometres between individuals i and j; φ i and φ j are the latitude coordinates of individuals i and j, respectively; λ i and λ 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).

| Geographical deme definitions
We represented the wolf geographical range as seven demes, defined by major geographical barriers through time. (1)

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

| 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 −6 ) with 24.4% of the variance within the data set explained by the chosen demes.

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

| 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
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 related sequences if they came from the same geographical location 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 observed 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 ST , 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 bined. This strategy is based on geographical proximity and genetic similarity in the data set. We note that this is not the same as modelling the combined demes as a single panmictic deme; structure between the demes is still modelled explicitly, but the summary statistics are averaged over multiple demes. The parameter m is measured in units of 1/1,000 years, and T, ΔT, K, K 1 , K 2 and K 3 are measured in units of 1,000 years. The parameters x, T and ΔT were sampled according to a uniform distribution over   Figures S13 and S14 give posterior density distributions for estimated parameters (ΔT, T, log 10 K 1 , log 10 K 2 , log 10 K 3 , log 10 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).

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

| 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 degree 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 spatially structured models, in which migration is restricted to adjacent populations, is a common approach for dealing with such situations 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).

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

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

| 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 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 retreat 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 alternative scenario where contemporary North American wolves are

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

F I G U R E 5
The inferred scenario of wolf demography from the Bayesian analysis using our spatially and temporally explicit model (see Figure  Effective population size

Prob. density
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).
The geographical exception to this pattern of widespread replacement is Beringia, where we infer demographic continuity between 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 phylogeographical 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 hypothesis 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 ecological 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, temperature, 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).  (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 adaptations in the Late Pleistocene and early Holocene period.

| Implications for the study of wolf domestication
Lastly, the complex demographic history of Eurasian grey wolves reported 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 underlying wolf population structure may explain why previous studies 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 potential 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-mediated 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 distribution observed in previous studies, including differences in levels of diversity found within local dog populations , 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 observed within dogs fall outside expectations that take admixture with geographically structured wolf populations into account.

ACK N OWLED G EM ENTS
We

DATA AVA I L A B I L I T Y S TAT E M E N T
The newly assembled mitochondrial genomes are available from