Worldwide patterns of genomic variation and admixture in gray wolves

Genome data and heterozygosity

In this study, we amassed the full genome sequences of 24 wolves, seven dogs (including the reference genome), and three outgroups (Fig. 1; Supplemental Table S1). Eleven of the individuals were uniquely sequenced in this study using the HiSeq 2000 platform with the remaining sequences obtained from previous studies (Wang et al. 2013; Freedman et al. 2014; Zhang et al. 2014; Supplemental Material).

To quantify genome-wide heterozygosity, we calculated the number of heterozygous SNPs over all sites (Fig. 2; Supplemental Table S2). The Mexican wolf had the lowest autosomal heterozygosity (0.00046), and the two Tibetan wolves also showed very low heterozygosity (0.0007 and 0.00086). Within European wolves, the Portuguese wolf showed the lowest heterozgosity (0.00101). However, the SNP rate was similar in all wolves (Supplemental Table S2). Within dogs, the basenji and dingo had the lowest heterozygosity (<0.001), with dingo having the lowest value of 0.00057 (Supplemental Table S2). We further calculated the heterozygosity of 5 Mb nonoverlapping windows across the 38 autosomal chromosomes (Supplemental Figs. S1–S4). The result confirmed that the low heterozygosity of two Tibetan wolves was evident across their entire genomes (Supplemental Fig. S1). The Mexican wolf exhibited extremely low heterozygosity in about half the chromosomes (Supplemental Fig. S4), and the Portuguese wolf also had very low heterozygosity in more than 15 autosomes (Supplemental Fig. S2). In contrast, the Inner Mongolia wolf 2 had very high heterozygosity across all chromosomes, even higher than other Inner Mongolia wolves, which partly may reflect lower genome coverage and a higher fraction of miscalled sites (Supplemental Figs. S1, S5). We also calculated the heterozygosity of exons and neutral regions (Supplemental Fig. S6).

To avoid the effect of inbreeding in the calculation of heterozygosity, we removed runs of homozygosity (ROH, see below) and recalculated heterozygosity (Fig. 2). The results are similar to that of the full data set. For example, the inbred Mexican wolf still had the lowest heterozygosity within wolves, and two Tibetan wolves had higher values than the Mexican wolf but lower than other wolves.

Genome-wide phylogenetic tree and PCAs

Autosomal SNPs were used to construct a maximum likelihood (ML) tree (Fig. 3). The topology of ML trees with (Fig. 3) and without (Supplemental Fig. S7) the boxer reference genome is consistent with geographical proximity of populations and does not support any specific wolf population as more closely related and possibly ancestral to domestic dogs. Specifically, all the dogs are monophyletic and define a sister taxon with Eurasian gray wolves that excludes a role for New World wolves in dog origins and suggests that the divergence of the modern Eurasian wolf population is nearly coincident with that of domestic dogs. Among gray wolves, East Asian wolves form a single clade, whereas European and Middle Eastern wolves (including Indian wolf) form a separate grouping. The Middle Eastern wolf is aligned with European rather than Asian wolf sequences (Fig. 3). For the New World wolves, two Yellowstone wolves cluster together, and then the Mexican wolf is grouped with them, but the divergence is large, suggesting an ancient separation of the two populations. However, the long branch of the Mexican wolf lineage may reflect the effect of small historic population size as the species went extinct in the wild, which was compounded by an extreme founding bottleneck in the captive population (Fredrickson et al. 2007). Both demographic events would tend to inflate genetic distance values. Nonetheless, this finding supports a previous hypothesis that Mexican wolves represented an early migration into North America (Leonard et al. 2005).

Principal component analysis (PCA) with LD-pruned data excluding the three outgroups and four wolf sequences (Inner Mongolia wolf 2, Eastern Russian wolf, Yellowstone wolf 2, and Yellowstone wolf 3) showed that PC1 (20.2% of variation) divided the samples into three clusters: domestic dogs; highland Chinese wolves; and other gray wolves (Fig. 4A). Further, when two outlier highland wolves were removed (Tibet wolf 1 and Qinghai wolf 1), dogs were more tightly clustered and separated from all wolves on PC1, whereas PC2 distinguished high altitude wolves and the Central Russian wolf from all other wolves (27.7% of variation of both axes combined) (Fig. 4B). PC3 and PC4 of both data sets separate Old and New World wolves, with the Mexican wolf showing the greatest distinction (Fig. 4C,D). The results with only wolves and dogs excluded showed a similar pattern with regard to clusters of wolves (Supplemental Fig. S8). Critically, we found no support for a closer association of Chinese wolves with domestic dogs as suggested by previous studies (Savolainen et al. 2002; Pang et al. 2009; Ding et al. 2012; Wang et al. 2013). We also ran PCAs for different geographical regions (Supplemental Figs. S9–S11). These results are consistent with the tree-based analysis in Figure 3, but do not take into account rate variation between lineages that can bias inferences about the actual amount of divergence. Moreover, PCA is a graphical approach that highlights genetic clusters in the data and should not be used to infer genealogical relationships (Novembre and Stephens 2008).

PSMC

The pairwise sequentially Markovian coalescent model (PSMC) was applied to investigate the timing of population-specific demography. Here, we report only the results for the higher mutation rate (Fig. 5; Supplemental Fig. S12) to be consistent with Freedman et al. (2014), but consider the results from both rates (1.0 × 10⁻⁸ and 0.4 × 10⁻⁸ per generation) in the Discussion, as effective size and the timing of population size changes should differ by a factor of approximately 2.5 (Supplemental Figs. S13, S14). All the wolves exhibited similar demographic trajectories until ∼80 kya; thereafter, the four highland Chinese wolves showed very different trajectories when compared with all other wolves (Fig. 5A). The Tibetan wolves experienced a continuous population decline beginning ∼25 to 55 kya and did not experience further population growth; whereas the Qinghai wolf experienced population growth at the same time as the Tibetan wolf bottleneck (Fig. 5A). However, caution needs to be used in interpreting these results because they might be explained by ancestral population structure or reflect smoothing across time intervals (Freedman et al. 2014). Other wolves experienced population growth or stagnation from 25 to 55 kya, which overlaps the Greatest Lake Period (25–40 kya) (Li and Zhu 2001).

PSMC projections showed that the remaining wolves suffered a worldwide decline of two- to threefold beginning ∼8 to 25 kya (Fig. 5), which is associated with the end of the last glacial period (10.5–25 kya) and megafaunal extinctions. The Chinese wolves showed the most divergent trajectories with the Tibetan wolves, demonstrating a sharp decline beginning >25 kya and followed by a less precipitous decline in Qinghai wolves. In contrast, the lowland Chinese wolf populations do not initiate a decline until ∼10 kya (Fig. 5A). The Middle Eastern wolves (Israeli, Iranian, and Indian wolves) and European wolves exhibited slightly different demographic trajectories between 10 and 80 kya (Fig. 5B). All these wolves show evidence of a population decline beginning 25 kya. Domestic dogs had similar trajectories and experienced a population decline and demographic divergence from wolves beginning ∼50 kya (Fig. 5C). The three Yellowstone wolves had concordant trajectories (Fig. 5D). However, the Mexican wolf experienced a more severe bottleneck, which may reflect both a recent history of decline and demographic smoothing across the last 10,000-yr interval (see Hedrick et al. 1997; Freedman et al. 2014). The Israeli golden jackal had higher N_e than the Kenya golden jackal, and the California coyote exhibited a different trajectory as expected given its status as an independent lineage (Supplemental Fig. S12). In conclusion, a consistent result across all these trajectories is a decline in population sizes during the period of 8 to 25 kya, coincident with the expansion of modern humans worldwide and the development of technology for capturing large game (Van Valkenburgh et al. 2015).

Autozygosity segments

To assess the history of inbreeding, we quantified genome-wide ROH using PLINK (Fig. 2; Supplemental Fig. S15; Purcell et al. 2007). The Mexican wolf had the longest ROH with a total length of 1,569,600 kb (Fig. 2), which was consistent with a founding bottleneck and subsequent inbreeding (Hedrick et al. 1997; Fredrickson et al. 2007). In fact, the distribution of ROH in the Mexican wolf was distinct from that of all other wolves, and showed the highest fraction of autozygous long segments (Supplemental Fig. S15d), which suggests very recent inbreeding (e.g., Boyko et al. 2010). The two Tibetan wolves had the longest total length of ROH (947,844 kb and 835,018 kb) and the highest fraction of autozygous segments in Old World wolves, especially at small segment size. This result indicates ancient inbreeding in the Tibetan wolf population (Fig. 2; Supplemental Fig. S15a). The Italian wolf had the highest fraction of autozygous segments at smaller ROH sizes among European wolves, whereas the Portuguese wolf had more segments at longer sizes (Supplemental Fig. S15b). This contrasting pattern is consistent with previous genetic analysis, suggesting an ancient population decline in Italian wolves (Lucchini et al. 2004; Pilot et al. 2014) and historical records showing a very recent population decline in Portuguese wolves (Sastre et al. 2011). Within dogs, dingo and basenji had the greatest ROH (dingo: 1,097,810 kb; basenji: 589,502 kb). They also had a higher fraction of autozygous segments, especially in the size range <4 Mb than the Tibetan mastiff and three Chinese indigenous dogs (Fig. 2; Supplemental Fig. S15c), suggesting more ancient inbreeding perhaps in the founding population of dingo that arrived to Australia (>4 kya) and in the origin of the basenji, an ancient breed of domestic dog. These results show that novel demographic insights into population-specific demography are provided by PSMC and ROH analyses, which are consistent with known recent history and past environmental events.

ABBA-BABA

Multiple runs of the ABBA-BABA test were performed to assess gene flow between Old World wolves and dogs (Supplemental Table S3). The results showed that all the European wolves and the Israeli wolf had significant gene flow with basenji and boxer. For the Asian wolves, the two Russian wolves and all the lowland Chinese wolves had significant gene flow with all the Chinese indigenous dogs, Tibetan mastiff, and dingo, whereas the two Tibetan wolves did not show any significant admixture with any dogs (Supplemental Table S3). However, Qinghai wolf 1 showed significant gene flow with two of the three Chinese indigenous dogs, and both Qinghai wolves had gene flow with dingo. The Mexican wolf and Yellowstone wolf did not show any admixture signal with boxer, dingo, or Chinese indigenous dogs (Supplemental Table S3). We note that where admixture is detected from multiple dog samples in one or more wolf populations, it may suggest that gene flow actually occurred from the common ancestor of these dog into a specific wolf population or one that was ancestral to multiple wolf populations.

We estimated the proportion of Chinese indigenous dog ancestry in Asian wolves that had evidence for significant admixture and for which more than one dog defined the comparison pool (Supplemental Table S3, see above; Green et al. 2010; Durand et al. 2011). The proportion of Chinese indigenous dog ancestry in the two Russian wolves varied from 15.3% to 19.52%. The proportion of dog ancestry in the two Xinjiang wolves varied from 9.28% to 11.3%. The average proportion of the dog ancestry in four Inner Mongolia wolves was 10.86%, 12.06%, 13.16%, and 21.59% (Supplemental Table S4). These results suggest substantial dog ancestry in wolf populations worldwide, which is conceivable given the long coexistence of dogs and wild wolf populations (Thalmann et al. 2013; Freedman et al. 2014). The only Old World population not showing any dog ancestry is the Tibetan wolf, which is also the most divergent population in the PCA (Fig. 4), suggesting the dog component of wolf genomes may influence patterns of relationships. The high altitude wolf populations also have a very recent history of exposure to aboriginal dog populations considering that the area was only permanently colonized by humans ∼7 kya (Brantingham et al. 2010; Chen et al. 2015).

Regarding the European and Middle Eastern wolves, we used basenji and boxer to estimate the dog ancestry in these wolf genomes (Supplemental Table S4). The proportion of dog ancestry in Israeli wolf, Western Russian wolf, and Spanish wolf is >20%. Of the others, the Portuguese wolf had the smallest proportion at 7.97%, whereas the Croatian wolf had the largest at 13.76% dog ancestry (Supplemental Table S4). These findings indicate a highly variable but substantial dog ancestry in most all extant wolf populations.

Demographic inference with G-PhoCS

We used the Generalized Phylogenetic Coalescent Sampler method (G-PhoCS) to infer the demographic history of wolves and dogs, including ancestral population sizes, divergence times, and rates of gene flow (Fig. 6). The analysis shows that wolf populations diverged over a relatively short period of time from ∼11,000 to 13,000 yr ago (ya), assuming a per-generation mutation rate of μ = 1.0 × 10⁻⁸ and an average generation time of 3 yr (Fig. 6). If a slower mutation rate of μ = 0.4 × 10⁻⁸ is used as suggested by Skoglund et al. (2015), this period of time is increased by a factor of 2.5 to 27,500–32,500 ya (see Discussion). The divergence of New and Old World wolves is the oldest of these events at 12,500 ya, followed by divergence of Eastern and Western Eurasian wolves at 11,700 ya. The divergence times between sequences from Europe, the Middle East, and Asia fall within a relatively short period of time of ∼1600 yr. New World wolves show an intermediate divergence time of ∼5400 ya. We infer dogs diverged from wolves just before the Eurasian wolf population splits (11,700 ya; CI: 11,100–12,300 ya). This divergence time is only 285–1565 yr (95% CI) more recent than the divergence of New World wolves. The tree implies a considerable preancestry of extant dogs of 1400–2700 yr and a substantial level of divergence among existing dog lineages (dingo, Chinese indigenous dog, and basenji). In contrast, gray wolves have a much more recent common ancestry than expected from their fossil record, and all the population diverged over a narrow time period consistent with a bottleneck followed by a rapid population expansion across Eurasia.

The population ancestral to Old and New World wolves was estimated to have a relatively large effective size of 45,100 individuals (CI: 44,400–45,900), when assuming a per-generation mutation rate of μ = 1.0 × 10⁻⁸. After divergence of Old and New World wolves, both populations experience decline, to 8000 individuals in Old World wolves and 17,300 individuals in New World wolves (Fig. 6; Supplemental Table S5). The Tibetan wolf had the smallest N_e within the Old World wolves (2500 individuals), whereas the lowland Chinese wolves had nearly fourfold larger N_e (Fig. 6; Supplemental Table S5). The ancestral population of New World wolves had a relatively large effective size of 17,300 individuals, implying a fairly modest bottleneck in the founding of the North American population. However, the two sampled populations have much lower inferred sizes of 3500 individuals for the Yellowstone wolf and 600 individuals for the Mexican wolf. The latter likely reflects a history of decline and extinction in the wild (see Discussion).

G-PhoCS models migration bands allowing a test of admixture from D-statistics. We infer relatively high rates of gene flow (5%–21%; aggregated 95% Bayesian credible intervals) from Chinese indigenous dogs to all Asian wolf populations and significant gene flow (2.4%–7.2%) in the opposite direction only for the lowland Chinese wolves (Table 1). Conversely, we find little evidence of admixture with the ancestors of the dingo (migration rates <3% in both directions). For the two highland Chinese wolf populations, G-PhoCS found relatively high rates of gene flow from Chinese indigenous dogs (Tibetan wolf: 5.5%–8.8%; Qinghai wolf: 14.1%–18.9%), although the ABBA-BABA test did not find evidence of admixture. In Western Eurasia, we observe high rates of gene flow between the Israeli wolf and basenji in both directions (4.3%–13.6%) and somewhat lower rates from Croatian wolf to basenji (1.3%–6.4%). Our findings of admixture between Chinese dogs and Asian wolves and Israeli wolf and basenji suggests admixture between the lineages ancestral to these breeds and wolf populations, as their geographic overlap is currently very limited. Consistent with previous results are the high rates of gene flow inferred from the population ancestral to all wolves and dogs into the golden jackal population (11.3%–13.6%) and much lower rates from several sampled wolf populations into the golden jackal population (up to 2.8%) (Freedman et al. 2014). The higher observed value of admixture between golden jackal and the common ancestor of modern wolves suggests an ancient admixture event. Finally, we infer low, but significant, levels of gene flow from the basenji into the Mexican wolf population (1.2%–3.2%), suggesting that like other wolf populations, the New World wolves also experienced admixture with dogs or share a common ancestor with Old World wolves that experienced admixture.

Genetic diversity and relationships of Old and New World wolves

Analysis of complete genome sequence data adds considerable resolution to the evolutionary relationships of gray wolves and domestic dogs. First, the genome-wide phylogenetic tree shows that the earliest split was between New and Old World wolves, which was followed by divergence between Old World wolves and dogs (Fig. 3). This result confirms dogs were domesticated in the Old World. In addition, the finding that no single wolf population is more closely clustered with domestic dogs supports the hypothesis that dogs were derived from a now extinct population of Late Pleistocene wolf (Thalmann et al. 2013; Freedman et al. 2014). However, the divergence time suggested by G-PhoCS (11,700 ya; CI: 11,100–12,300 ya) is more recent than estimates based on ancient DNA analysis of early dogs and wolves (27,000 ya) (Thalmann et al. 2013). These differences might be caused by inflated mutation rates in the neutral regions used in this study, undetected admixture with dogs, or other assumptions of the underlying G-PhoCS model. The existence of dog fossils older than this recent divergence date, and confirmed by mtDNA sequence data, supports a more ancient origination (Thalmann et al. 2013; Skoglund et al. 2015). In fact, if the mutation rates associated with Skoglund et al. (2015) are used, the divergence time increases to ∼29 kya, a value close to their estimate of 27 kya (Fig. 6). Finally, within the Old World clade, wolf and dog represent sister taxa. Therefore, suggestions that the dog or dingo are a separate species (Canis familiaris) (e.g., Crowther et al. 2014) would cause gray wolves to be a polyphetic taxon; and consequently, our results support dogs as a divergent subspecies of the wolf. This result has societal significance as legislation in some countries and regional governments consider wolves and dogs as distinct species restricting the possession, interbreeding, or the use of vaccines and medications in wolves or dog–wolf hybrids if they have only been approved for use in dogs. In this sense, analysis of evolutionary history informs law and veterinary practice, as dog lineages are nearly as distinct from one another as wolves are from dogs, and the justification for treating dogs and wolves differently is questionable.

The evolutionary tree (Fig. 3) and PCA (Fig. 4) show that the Mexican wolf is a divergent form of gray wolf, suggesting it is a remnant of an early invasion into North America (García-Moreno et al. 1996; Leonard et al. 2005; vonHoldt et al. 2011) and contradicting suggestions that it is not a distinct subspecies (Cronin et al. 2015). The ROH and genome-wide heterozygosity results (Fig. 2) also showed that the Mexican wolf is a highly inbred population (vonHoldt et al. 2011). The subspecies had the smallest effective population size of only 600 individuals in the sampled wolf populations (Fig. 6; Supplemental Table S5). Further, the high long-range ROH in the Mexican wolf implies a long-term decline, followed by a small founding population and inbreeding in the captive population (Hedrick et al. 1997; Fredrickson et al. 2007). These results justify immediate conservation actions to protect this endangered and distinct wolf lineage. Further, population numbers should be increased through captive breeding and in situ conservation to prevent additional genetic erosion. Currently, such efforts have been hindered by the lack of an informed management plan (Wayne and Hedrick 2011; Hendricks et al. 2016).

The Tibetan wolf was found to be the most highly divergent Old World wolf, given its distinct position in the phylogenetic tree (Fig. 3) and the PCA plot (Fig. 4). It also exhibited extremely low heterozygosity (Fig. 2; Supplemental Figs. S1–S4), suggesting that it experienced a historical bottleneck, and only recently recolonized much of the Tibetan Plateau. Indeed, PSMC revealed that the Tibetan wolf suffered a substantial population bottleneck that began ∼55 kya (mutation rate 1.0 × 10⁻⁸) or >100 kya, assuming a slower mutation rate, and then declined to the present day (Fig. 5A). Notably, all other wolves showed evidence of growth during the Greatest Lake Period from ∼25 to 55 kya (Fig. 5). The severe habitat loss during glaciations probably contributed to the dramatic population decline of the Tibetan wolf between 10 and 55 kya (Xu and Shen 1995; Yi et al. 2005; Clark et al. 2009; Chevalier et al. 2011; Heyman 2014; Zhang et al. 2014). In addition, both archaeological and genetic analysis suggest that the first colonization might be as early as 30 kya (Aldenderfer 2011), and little evidence exists for permanent human occupation before 7 kya (Brantingham et al. 2010; Chen et al. 2015). Therefore, the appearance of human settlements may have contributed to the decline of Tibetan wolf population but did not initiate the population bottleneck more than 50 kya. Finally, Tibetan wolves had the longest total length of ROHs of the Old World wolves (Fig. 2), and a large proportion of their ROHs are in relatively short segments (Supplemental Fig. S15), which suggests that it experienced ancient inbreeding. Moreover, the ABBA-BABA test did not detect substantial gene flow between Tibetan wolf and dogs, suggesting dog admixture did not contribute to ROH (see discussion below). In summary, we suggest that the unique high altitude environment and history of the Tibetan Plateau made wolves there more susceptible to habitat loss, genetic isolation, and allowed for local adaptation. Consequently, these conditions resulted in the evolution of the most distinct wolf population in the Old World.

Geographical structure is evident within Old World wolves (Fig. 3). Previously, analysis of Eurasian wolves with mtDNA control region sequences did not reveal any distinct genetic partitions and suggested modern wolves originated over 250 kya (Vilà et al. 1999). However, European and Middle Eastern partitions were apparent in genome-wide SNP data (vonHoldt et al. 2011). The one-million-year divergence time between wolves and coyotes used previously was based on fossil occurrence data; and given the dynamics of morphological turnover in the coyote lineage (Meachen and Samuels 2012), first occurrence of coyote-like specimens may not accurately reveal the ancestry of modern forms. Our results and those from previous studies (Freedman et al. 2014; Koepfli et al. 2015; Skoglund et al. 2015) suggest the one-million-year divergence time may be inflated by a factor of 20 or more, and modern Eurasian wolves coalesce ∼13 kya or ∼32.5 kya, the latter using the slower mutation rate from Skoglund et al. (2015). Importantly, the slower rate leads to divergence dates more consistent with the presence of ancient dog fossils well before 15,000 yr ago. Nonetheless, the Skoglund et al. rate needs additional confirmation because it is based on a single fossil specimen with only onefold sequencing depth and used only a subset of DNA sequence to calculate the rate (Skoglund et al. 2015). Thus, until more direct measurements of mutation rates become available, fossil calibration will remain the main source of uncertainty in the timing of key events in canid evolutionary history.

The PSMC results revealed that all wolves shared a similar trajectory before ∼100–125 kya (mutation rate 0.4 × 10⁻⁸) or ∼30–50 kya (mutation rate 1.0 × 10⁻⁸). In combination, the dating and PSMC results suggested that over the last million years, numerous wolf-like forms existed but that turnover was high, and modern wolves were not the lineal ancestors of dogs (Leonard et al. 2005; Thalmann et al. 2013; Freedman et al. 2014). Indeed, the population size of the Croatian wolf reduced about 10-fold compared to the wolf ancestor, and Yellowstone wolf and Mexican wolf also reduced five- to 28-fold (Fig. 6; Supplemental Table S5). This pattern of population reduction and turnover also is supported by recent mtDNA sequence analysis of modern and ancient wolves from the Last Glacial Maximum (Leonard et al. 2007; Pilot et al. 2010; Thalmann et al. 2013) and the dynamic pattern of turnover in other large carnivores such as brown and polar bears, hyenas, and lions as inferred from genetic data (Miller et al. 2012; Cho et al. 2013).

Finally, even assuming a slower mutation rate, our results imply a remarkably recent coalescence of extant wolves several hundred thousand years after the appearance of wolf-like canids (Wayne and Ostrander 2007). Both slow and fast mutation rate estimates are consistent with the possibility that modern humans impacted the demography of gray wolves as they colonized Eurasia, encountered wolves, domesticated some, and possibly caused the decline of others. Humans are the most effective competitor of large carnivores and could have readily removed them from ecosystems as they do today. Additionally, the presence of large domestic dogs may have accelerated the rate of decline of carnivores that competed with humans (e.g., Shipman 2015). Our results imply that the effect of humans on large predators may have preceded the megafaunal extinctions ∼10 kya and may represent one of the earliest anthropogenic causes of decline in animal populations.

Admixture and relationships to domestic dogs

None of our wolf sequences cluster exclusively with domestic dogs, supporting the hypothesis based on only three wolf genomes (Freedman et al. 2014) and ancient DNA (Thalmann et al. 2013) that the immediate gray wolf ancestor of dogs is now extinct. Nonetheless, modern gray wolves have likely influenced the recent history of domestic dogs through admixture. Both the ABBA-BABA tests and G-PhoCS support the notion of extensive admixture between dogs and wolves, with up to 20% of the genome of East Asian wolves showing signs of dog ancestry. We also detected that the genomes of European and Middle Eastern wolves had ∼7%–25% dog ancestry. Most of the observed gene flow events have not been reported previously. Interestingly, the two highland wolf populations of the Tibetan Plateau showed no evidence of admixture in the ABBA-BABA tests, but G-PhoCS did infer elevated migration rates from Chinese indigenous dogs into these populations (Tibetan wolf: 5.5%–8.8%; Qinghai wolf: 14.1%–18.9%). Conceivably, this finding may be a result of gene flow from dogs into the population ancestral to all modern wolves, which influenced the distribution of coalescent times but cannot be detected using D-statistics because it similarly affected all wolf populations.

For comparison, using ABBA-BABA tests, it was found that modern humans admixed with Neanderthals over 40 kya, but no more than 5% of the modern human genome could be attributed to admixture (Green et al. 2010), suggesting that wolves and dogs have more extensive and regionally based admixture. As in Neanderthals, admixture may have enhanced adaptation in wolves. For example, admixture of pre-Columbian dogs and wolves in North America transferred the black coat color locus to wolves, conferring greater longevity and resulting in a continent-wide selective sweep (Anderson et al. 2009; Coulson et al. 2011). However, in the North American wolves sampled, no apparent trace of admixture remains elsewhere in the genome (Anderson et al. 2009). The persistent admixture between dogs and wolves suggests a unique mode of evolution in which mutations that occur independently in dogs and wolves, under dramatically different selective regimes, can be shared and potentially accelerate the process of evolution. Such coupled evolutionary histories may exist in other vertebrate species as well, such as in wild and domestic pigs and in brown bears and polar bears (Groenen et al. 2012; Miller et al. 2012).

Finally, admixture can have a confounding effect on inferences about dog domestication history. Specifically, past inferences about dog origins based on private SNPs shared with dogs (vonHoldt et al. 2010), greater genome-wide similarity between Chinese wolves and dogs (Wang et al. 2013), or lower LD (Shannon et al. 2015) may reflect regional admixture with wolves and gene flow among dog populations rather than the geographic origin of domestication. Similarly, highly divergent breeds may have more admixture and wolf ancestry retained in their genome. Potentially, this bias might be removed by applying analytical approaches that excise dog segments from wolf genomes. However, direct tracking of genetic changes in wolves and dogs through ancient DNA analysis may be a more robust approach (Grimm 2015).

Samples and sequencing

We sequenced genomes of dogs, wolves, and other wild canids from Africa, Asia, Europe, the Middle East, and North America. Together with published canid genomes, we generated a final data set with 34 full genome sequences at 9–28× coverage with an average coverage of 29.8× (Fig. 1; see Supplemental Material).

Mapping short reads and genotyping

The 100-bp pair-end (PE) short reads of each sample were aligned to the dog genome (CanFam3.1) using Bowtie 2 (Langmead and Salzberg 2012) under the local alignment algorithm with very sensitive model and proper insert sizes of each sample. Default options were used for other parameters. Then, we applied Picard and GATK toolsets (DePristo et al. 2011) to process the alignments to SNP calls. The whole pipeline converted the short reads to BAM format alignment files, and then generated genotype calls in Variant Call Format (VCF). The pipeline is the same as used in our previous studies (Fan et al. 2014; Freedman et al. 2014; Zhang et al. 2014). We applied a series of data quality filters to improve the quality of genotype calls (see Supplemental Material).

Phylogenetic tree and PCA

A ML tree from whole-genome SNP data was constructed using SNPhylo (Lee et al. 2014). SNPhylo transforms genotype data into a structured data array (Bioconductor gdsfmt) and then generates and aligns SNP sequences and constructs the phylogenetic trees. The program was run with 100 bootstrap repetitions, and only one outgroup was used (Israeli golden jackal) due to the software's internal limitations.

PCA was performed using the pairwise allele-sharing genetic distance. Following vonHoldt et al. (2011), sites exhibiting apparent strong local linkage disequilibrium (R² > 0.5) were filtered using the --indep option in PLINK (--indep 50 5 0.2) (Purcell et al. 2007). To improve resolution among wolves, the two golden jackals and coyote were removed from the PCA because they were too divergent from dogs and wolves, and their inclusion compressed the scatter among wolves on the first few PCs. The lower coverage genomes (<10-fold; Inner Mongolia wolf 2, Eastern Russian wolf, and Yellowstone wolf 3) were also removed due to their potential high genotype error (Supplemental Fig. S5). Additional PCA was also performed excluding one Tibetan wolf and one Qinghai wolf based on the observation that highland Chinese wolves were similar to one another but were highly divergent from all other wolves (Zhang et al. 2014). Finally, additional PCAs were performed with samples from specific geographic regions, such as Asia or Europe, and including only gray wolves. In both tree analyses and PCAs, we excluded the Yellowstone wolf 2 because it is the offspring of Yellowstone wolf 1 (mother) and Yellowstone wolf 3 (father).

Inference of population size changes through time with PSMC

We used PSMC (Li and Durbin 2011) to infer demographic history. The following parameters were used: numbers of iterations = 25; time interval = 64 × 1; and generation time = 3 (Freedman et al. 2014; Zhang et al. 2014). Our previous studies used a mutation rate of 1.0 × 10⁻⁸ per generation, a commonly applied value. However, one recent study based on ancient wolf genome sequences estimated that the mutation rate was only 0.4 × 10⁻⁸ per generation (Skoglund et al. 2015). Therefore, we used both mutation rates in our study to bracket estimates of divergence time and effective population size.

Runs of homozygosity analysis

ROHs were calculated with PLINK (Purcell et al. 2007). The input data here is the same data in PCA. We then searched for ROHs spanning at least 500 homozygous SNPs in 1 Mb nonoverlapping windows, allowing for a maximum of five heterozygous sites and 30 missing genotypes per window.

Detection of gene flow using the D-statistic

To test whether there was gene flow between wolves and dogs in each geographical region after their divergence, we applied the ABBA-BABA test (D-statistic) between closely related populations by detecting differences in allele sharing between two lineages (P1 and P2) with a third lineage (P3) (Durand et al. 2011; see Supplemental Material). Under the assumption of one gene flow event that is recent compared to the divergence of dogs and wolves, we further used the Durand et al. (2011) equation to estimate the proportion of dog ancestry in the wolf genomes (see Supplemental Material).

Demographic inference with G-PhoCS

The G-PhoCS method (Gronau et al. 2011) was used to estimate divergence times, population sizes, and migration rates. Considering the computational resources required, our analysis focused on a subset of 22 of the 33 genomes that represent all 11 wolf populations, four dog populations, and the Israeli golden jackal as outgroup (Supplemental Table S6). Alignments of the 22 genomes were done over the 13,647 neutral loci designed by Freedman et al. (2014) for use in demographic inference. We ran several analyses (see Supplemental Material) using the standard settings and assumed standard priors for model parameters, as described by Gronau et al. (2011).