Genomic Approaches Reveal an Endemic Subpopulation of Gray Wolves in Southern China.

Although gray wolves (Canis lupus) are one of the most widely distributed terrestrial mammals, their origins in China are not well understood. We sequenced six specimens from wolf skins, showing that gray wolves from Southern China (SC) derive from a single lineage, distinct from gray wolves from the Tibetan Plateau and Northern China, suggesting that SC gray wolves may form a distinct subpopulation. Of SC gray wolves, one wolf from Zhejiang carries a genetic component from a canid and had gene flow from a population related to or further diverged from wolves than the dhole. This may indicate that interspecific gene flow likely played an important role in shaping the speciation patterns and population structure in the genus Canis. Our study is the first to survey museum gray wolves' genomes from Southern China, highlighting how sequencing the paleogenome from museum specimens can help us to study extinct species.

Introduction

The place of origin for domestic dogs (Canis lupus familiaris) remains a controversial question for the scientific community despite many efforts at studying dog domestication (Botigué et al., 2017, Frantz et al., 2016, Liu et al., 2018, Shannon et al., 2015, Thalmann et al., 2013, Vonholdt et al., 2010, Wang et al., 2016a). Geographic distribution, population structure, and genomic features of wild ancestors are essential factors to determine sources of domestication (Wang et al., 2016b). Gray wolves (Canis lupus) are the closest wild relative of dogs, and they are also one of the most widely distributed terrestrial mammals, originally inhabiting major parts of Eurasia, North America, and North Africa (Loog et al., 2018, Wilson and Reeder, 2005, Young et al., 1944). Previous studies suggested that gray wolves have a complex history (vonHoldt et al., 2011, Wayne et al., 1992), with subpopulation structure related to local niches (Carmichael et al., 2001, Geffen et al., 2004, Pilot et al., 2006, Pilot et al., 2014) and long-term genetic admixture not only with dogs (Freedman et al., 2014, Liu et al., 2018) but also with coyotes (Fan et al., 2016, Monzon et al., 2014, Pilot et al., 2014, vonHoldt et al., 2016, vonHoldt et al., 2011). In China, gray wolves were distributed across the mainland, including most southern regions (Smith and Xie, 2014, Wang et al., 2016b). Genomic approaches using gray wolf specimens from Southern China (SC) may help to shed new light on the demographic history of gray wolves and domestic dogs.

Results and Discussion

From two natural history museums in China, we obtained six historical gray wolf skin samples collected from Mainland China (Figures 1A and S1 and Table 1, detailed description in Wang et al., 2016b). More details for the samples are shown in the methods (Transparent Methods). As skin samples were treated with chemical reagents and underwent special processing for preservation during storage and exhibition in museums, we used a modified ancient DNA (aDNA) protocol (Dabney et al., 2013) to retrieve genetic material from the skin samples. In total, 35 genomic libraries were produced for these six samples using a double-stranded library preparation protocol (Kircher et al., 2012, Meyer and Kircher, 2010), and each was treated with uracil-DNA glycosylase and endonuclease (Endo VIII) to remove characteristic aDNA deamination (Briggs et al., 2007) (Table S1). We sequenced the libraries using 2×150-bp reads on an Illumina HiSeq X platform.

Figure 1

Geographic Locations, Population Structure of 31 Gray Wolves, and Phylogeny of 39 Canids

(A and B) (A) Geographic locations, where the key is shared with the principal-component analysis in (B). In (B), the label “New” represents the six samples sequenced in the study, and the label “Published” represents 25 samples from previous studies.

(C) The maximum likelihood tree of 39 canids, where the Andean fox is used as an outgroup. All the gray wolves are in the red line, and the newly sequenced individuals are marked in bold and italics.

Table 1

Information on Samples Sequenced in This Study

Sample ID	Coverage	SNPs	Sources	Location
W6_Jilin	1.87	6,890,236	National Zoological Museum of China, Beijing, China	Jilin
W7_Zhejiang	0.15	1,866,574		Zhejiang
W9_Heilongjiang	1.61	5,438,606		Heilongjiang
W12_Guizhou	15.33	12,688,043	Kunming Natural History Museum of Zoology, Yunnan, China	Guizhou
W13_Guizhou	12.39	12,968,150		Guizhou
W2_Jiangxi	37.00	13,528,667		Jiangxi

We included 103 canids from previous studies (Auton et al., 2013, Botigué et al., 2017, Decker et al., 2015, Freedman et al., 2014, Marsden et al., 2016, Skoglund et al., 2015, Tang et al., 2019, vonHoldt et al., 2016, Wang et al., 2019a, Wang et al., 2016a, Wang et al., 2013, Wang et al., 2019b, Zhang et al., 2014), including an ancient gray wolf, Taimyr; 25 modern gray wolves chosen from regions overlapping their current ranges; 70 domestic dogs from all over the world; two jackals; two coyotes; one red wolf; one dhole; and one Andean fox (Figure 1, Table S2). We obtained 13.74 million SNPs, including 4.25 million transversions for further analysis (Transparent Methods).

All samples were sequenced to ∼0.15- to 15.3-fold average coverage except for the Jiangxi wolf (W2_Jiangxi), which was sequenced to 37.00-fold coverage (Table 1). Because most samples have low average coverage, we focused on SNPs previously found for modern canids (Marsden et al., 2016, Wang et al., 2016a), for which we called haploid alleles using randomly chosen sequence reads for every sample except W2_Jiangxi. We applied a filter where we ignored fragments with length less than 30, ignored the first and last two base pairs of each fragment, and required a base pair quality higher than 20 and mapping quality of at least 30. For the high-coverage W2_Jiangxi (Table 1) as well as the two medium-coverage Guizhou wolves (W12_Guizhou and W13_Guizhou), we called heterozygotes using the software GATK with the Unified Genotyper parameter to determine diploid calls (DePristo et al., 2011). Results using diploid calls for W12_Guizhou and W13_Guizhou are similar to results using random calls without heterozygous sites.

Phylogeny and Population Structure

To investigate the relationship of the newly sampled individuals to wolf and dog populations, we calculated pairwise allele-sharing distances among all pairs of wolf and dog populations. We applied a principal-components analysis (PCA) to the resulting pairwise distance matrix using SMARTPCA (version: 13050) (Patterson et al., 2006). The first principal component distinguishes between gray wolf and dog populations, whereas the second principal component distinguishes between East Asian and European dogs (Figure S2), consistent with previous studies (Frantz et al., 2016, Vonholdt et al., 2010, Wang et al., 2016a). To obtain increased resolution, we redid the PCA excluding dog populations (Figure 1B). The resulting PCA shows that the two new Guizhou wolves (W12_Guizhou and W13_Guizhou) cluster with a Chinese wolf from the San Diego Zoo, whose origin is not recorded (labeled as “China_X” in this study) (Freedman et al., 2014), and two gray wolves from Shanxi, China, sampled in 1988 near the border of SC (Wang et al., 2016a). We also find that the new Jilin and Heilongjiang wolves (W6_Jilin and W9_Heilongjiang) cluster with gray wolves from Inner Mongolia and Liaoning. The new Jiangxi wolf (W2_Jiangxi) is closest to the Qinghai wolf, whereas the Zhejiang wolf (W7_Zhejiang) is closest to the cluster containing the gray wolves from Inner Mongolia and Liaoning (Figure S1).

We constructed a maximum likelihood (ML) tree (Figure 1C) and a neighbor-joining tree with only wolves (Figure S3), from which we further find that the Zhejiang wolf forms a clade with the Guizhou, Shanxi, and China_X wolves (Figure 1C). We define these gray wolves as the gray wolves from SC and find that they are most closely related to the Qinghai, Tibet, and Jiangxi wolves, where the Tibetan gray wolf (Canis lupus chanco) is a gray wolf subspecies that occupies habitats on the Qinghai-Tibet Plateau (Pocock, 1941) and possesses adaptations to high-altitude environments (Zhang et al., 2014). Other gray wolves from Northern Asia (NA, e.g., W6_Jilin, W9_Heilongjiang, and Liaoning) form a clade with the SC and Tibetan gray wolves relative to other NA, Middle Eastern, and European gray wolves that form a distinct clade with each other. The 35,000-year-old wolf from the Taimyr Peninsula in Northern Siberia joins at the base of the Eurasian wolf phylogeny, and American wolves separate the earliest from other wolves, consistent with Freedman et al. (2014).

We used TreeMix (v. 1.13) (Pickrell and Pritchard, 2012) to investigate the genetic relationship between historical and present-day wolves. TreeMix determines population structure using ML trees and allows for both population splits and potential gene flow by using genome-wide allele frequency data and a Gaussian approximation of genetic drift. The ML tree (Figure S4) without admixture (m = 0) is consistent with previous patterns, wherein gray wolves from East Asia form three groups: Guizhou and Zhejiang wolves form a clade with Shanxi and China_X wolves, the Jiangxi wolf forms a clade with Tibet and Qinghai wolves, and Jilin and Heilongjiang wolves, like other NA wolves, are outgroup to SC and Tibetan wolves. All three of these groups form a clade relative to non-Asian wolves.

We measured the shared genetic drift between each newly sequenced individual (X) and other dogs and wolves (Y) since their separation from an outgroup, Dhole, using f(X, Y; Dhole) (Raghavan et al., 2014, Wang et al., 2019a) and found a similar pattern as above, where the Zhejiang wolf shares the most genetic drift with gray wolves from Guizhou (Figure S5). We then used D-statistics (Patterson et al., 2012) of the form D(Fox, Test; X, Y), where X and Y are previously published wolves, to formally test the relationship these historical wolves have with different wolf populations. For Guizhou and Zhejiang wolves, we find that they share more alleles with the SC gray wolves (Shanxi and China_X) than with all other wolves, as D(Fox, W12_Guizhou/W13_Guizhou/W7_Zhejiang; X, SC) > 0 (9.8 < Z < 38.0, Table S3). The Jiangxi wolf shares more allele with the Tibet and Qinghai gray wolves than with other wolves, as D(Fox, W2_Jiangxi; X, Tibet/Qinghai) > 0 (11.9 < Z < 27.6, Table S3), whereas Jilin and Heilongjiang wolves share the most alleles with NA gray wolves. These results support our above-mentioned analyses, again grouping the Guizhou and Zhejiang wolves with the SC wolves, the Jiangxi wolf with the Tibetan wolves, and the Jilin and Heilongjiang wolves with the NA wolves.

In summary, our results revealed that the lowland Chinese wolves (Fan et al., 2016) consisted of two major populations: SC and NA wolves. Gray wolves from Zhejiang and Guizhou group most closely with and share the most genetic drift with SC wolves, which includes present-day populations in Shanxi and China_X. The Jilin and Heilongjiang gray wolves share the most genetic similarity to the NA gray wolves, which is the other clade in Northern China and Eastern Russia.

Testing for Admixture in Gray Wolves

Using D(Fox, X; Test, Y), we find that all gray wolves (X) share more alleles with the Guizhou, Jilin, and Heilongjiang wolves (Test) than with the gray wolves from Tibet and Qinghai, i.e., D(Fox, X; Test, Qinghai/Tibet) < 0 (Table S5). The Jiangxi wolf shares more alleles with the Tibet and Qinghai wolves than with the other wolves (Table S3), and here, we find that the Tibet and Qinghai wolves share more alleles with the Jiangxi wolf, i.e., D(Fox, Qinghai/Tibet; W2_Jiangxi, X) < 0 (−34.5 < Z < −15.4, Table S5), emphasizing that the Jiangxi, Tibet, and Qinghai wolves form a clade. However, whereas D(Fox, X; W2_Jiangxi, Qinghai)∼0 (−2.3 < Z < 1.8, Table S5), indicating a symmetric relationship as expected for the Jiangxi and Qinghai wolves forming a clade, we observe D(Fox, X; W12_Jiangxi, Tibet) < 0 (Table S5), suggesting that the Jiangxi wolf has a connection to other wolves relative to the Tibetan wolf. Thus, we observe that the Tibet and Qinghai gray wolves act as an outgroup to most SC and NA gray wolf populations, and although the Jiangxi wolf forms a clade with the Tibet and Qinghai gray wolves, the Jiangxi wolf also shows connections to non-Tibetan wolf populations. To test for admixture between the SC, NA, and Tibetan wolves, we used f3(Test; X, Y), where a significantly negative value (Z < −3) suggests that the Test population is a mixture of ancestry related to X and Y, two other wolf populations. Testing all combinations as both a source population and the admixed population, we found that f3(W12_Jiangxi/Qinghai; Tibetan, SC/NA) < 0 (−19.101 < Z < −9.141, Table S4), suggesting that both the Jiangxi and Qinghai wolves show evidence of ancestry from populations related to both Tibetan and SC gray wolves, explaining why the Jiangxi gray wolf shares a connection to non-Tibetan gray wolf populations.

In contrast to all other gray wolves from China, the Zhejiang wolf shows a markedly different pattern, where all other wolf populations, including the Tibetan and Qinghai gray wolves and more distantly related wolves from further west, form a clade with each other relative to the Zhejiang wolf. That is, we observe that D(Fox, X; W7_Zhejiang, Y) > 0 (8.0 < Z < 81.3, Table S5), where X and Y are all other gray wolves, including the Taimyr. Earlier, we found that D(Fox, W7_Zhejiang; X, Shanxi/China_X) > 0 (Table S3), indicating that the Zhejiang wolf shows connections to the wolves from Shanxi and China_X. These results suggest that the Zhejiang wolf shows a close relationship to gray wolves from Shanxi and China_X, but that this wolf also possesses an ancestral component that is older than the common ancestral population of the Taimyr and all other gray wolves. The error rate for the Zhejiang wolf (0.4%) is higher than that estimated for other wolves sequenced in this study (0.1%–0.2%), likely because of its low coverage (Table 1). After simulating an error rate similar to that observed for the Zhejiang wolf in these other wolves, we find that our results remain consistent with our previous results. That is, the Zhejiang wolf shows a distinct pattern from that observed in other wolves for both lower and higher error rates.

We use the genomic data from canids typically outgroup to all wolves and dogs, the Dhole, Jackal, Coyote, and Red wolf, to understand how deeply the old component found in the Zhejiang wolf separated from other canid populations. Other canids separated from wolf populations very early, with the Dhole diverging earliest, followed by the Jackal and most recently the Coyote (Lindblad-Toh et al., 2005, Macdonald and Sillero-Zubiri, 2004, Wayne et al., 1997). The Red wolf is genetically very similar to the Coyote and shows substantial gene flow from gray wolves (vonHoldt et al., 2016, vonHoldt et al., 2011). First, comparing the Jackal to wolves (X) and the Coyote, we find that for all wolves but the Zhejiang wolf, D(Fox, Jackal; X, Coyote) ranges from −0.042 to −0.017 (−16.5 < Z < −6.4, Table S6), indicating a connection between the Jackal and gray wolves. We find the reverse for the Zhejiang wolf, however, where the Jackal shares more alleles with the Coyote than with the Zhejiang wolf, i.e., D(Fox, Jackal; W7_Zhejiang, Coyote)=0.134 (Z = 34.1, Table S6). We find similar results replacing the Coyote with the Red wolf. The large contrast between the results for the Zhejiang wolf compared with other gray wolves suggests that the old component came from a population that diverged deeply in the past, who separated before the common ancestor of jackals and coyotes.

We also observe that for all gray wolves, we find D(Fox, Dhole; X, Jackal) > 0 (Table S6), suggesting that gray wolves share a deep lineage older than the separation of the Jackal and Dhole or that there is a direct genetic connection between the Jackal and Dhole. However, whereas D ranges from 0.004 to 0.028 (4.5 < Z < 7.9, Table S6) for most gray wolves, using the Zhejiang wolf greatly increases the D value to 0.117 (Z = 20.2, Table S6). We find that D(Fox, Dhole; W7_Zhejiang, Jackal) remains significantly positive (Z = 12.2) using only transversions, suggesting that the result for the Zhejiang wolf is not related to aDNA damage and likely reflects an unusual admixture history. If the Zhejiang wolf was no different from other gray wolves, especially the two Guizhou individuals to which they share the closest relationship (Table S5), we would expect to find that D(Fox, Dhole; X, Zhejiang)∼0, which would indicate that the Zhejiang wolf and other gray wolves are similarly related to the Dhole. However, we observe that D(Fox, Dhole; X, W7_Zhejiang) < 0 (−22.4 < Z < −11.7, Table S6), indicating that the Dhole shares more alleles with other gray wolves than with the Zhejiang wolf. These patterns suggest that the ancestral component found in the Zhejiang wolf came from a population that diverged earlier than the common ancestor of the Jackal and Dhole, which is older than the separation of the Jackal and Dhole from wolves, suggesting that the Zhejiang wolf possesses very deep ancestry not found in other gray wolves.

Using the tree model with no admixture from TreeMix (Figure S6), we visualized the matrix of residuals (Figure S7) to determine how the estimated genetic relationship between each pair of canids fit the model. A high residual indicates that the pair does not fit the tree model and may be candidates for an admixture event. We find two candidate admixture events, the first between the Andean fox and Zhejiang wolf (Figure 2) and the second between gray wolves from Portugal and Iberia. In an ML tree allowing two admixture events, admixture from the lineage leading to the Andean fox to the lineage leading to the Zhejiang wolf is included, whereas gray wolves from Portugal and Iberia are grouped into the same cluster (Figure S8). Admixture between the outgroup Andean fox and the Zhejiang wolf supports our conclusions from the D-statistic analysis (Tables S5 and S6), in which the Zhejiang wolf possesses an ancestral component that came from a population that diverged earlier than the Jackal or Dhole did from wolves. The estimated value of the migration event in the Zhejiang individual is 12.3% ± 0.4% (p = 2.2 × 10−308). In the TreeMix analysis, we used the Andean fox as an outgroup, whose distance from the included canids would result in weak phylogenetic constraints. In addition, we also used the F4-ratio test to estimate the proportion of this deep ancestry, and as it is older than the separation of the Jackal and Dhole, we used an unrooted phylogeny where the Fox is used as a proxy for the source of the deep ancestry (Figure S9). Thus, we estimate the proportion of ancestry related to the Fox, which is given by:where X is each gray wolf in turn.

Figure 2

The Maximum Likelihood Tree Based on TreeMix with m = 1

The scale bar shows 10 times the average standard error of the entries in the sample covariance matrix. We have used the prefix “Gray_Wolf_” for highlighting gray wolves.

Using this method, we found that the estimated admixture proportion of the deep ancestry for the Zhejiang wolf is 11.7% ± 0.5%, whereas all other gray wolves have an estimated admixture proportion close to zero. Similarly, in the qpGraph analysis, we estimated an admixture proportion for the deep ancestry in W7_Zhejiang of 14% (Figures S10–S14). Thus the TreeMix, F4-ratio and admixture graph analyses all support the presence of gene flow from an ancient canid population into the ancestors of the Zhejiang wolf.

Conclusion

The distribution of gray wolves in East Asia is controversial because some studies have claimed that gray wolves never existed (Callaway, 2013, Larson and Fuller, 2014) or are now extinct in SC (Lau et al., 2010), whereas others sources, especially those based on Chinese literature, have stated that they are present across all of mainland China (Wang et al., 2016b). In this study, we provide the first comparative genomic analysis of gray wolves from East Asia, focusing particularly on wolves from SC, where some believed no gray wolves were distributed (Callan et al., 2013, Larson and Fuller, 2014). Previously, Asian wolves could be divided into two populations: Tibetan gray wolves (Canis lupus chanco) and Chinese lowland wolves (Fan et al., 2016, Zhang et al., 2014). Here, using ancient genome-wide data, we reconstruct the phylogeny and evaluate the population structure and shared genetic drift between East Asian gray wolves to show that they form three major groups, which we call Southern China (SC), Northern Asia (NA), and Tibetan, based on their geographic distribution. Interestingly, specimens from SC gray wolves were all collected from 1963 to 1988. Our results highlight that the population in SC is endemic, and with the fast growing economic development of China, it is paramount to protect and restore their ecological habitat. Through our study, we also emphasize the value of using paleogenomic approaches to study the numerous museum specimens available (Min-Shan Ko et al., 2018, Roy et al., 1994), and we address the importance of using population genomics to determine current or future conservation efforts.

Finally, our analyses show that admixture played a large role between the different Asian wolves, and we highlight two instances here. First, we find that a wolf as far southeast as Jiangxi province shows evidence of being a mixture of Tibetan-related wolves and other wolves in China. Second, we traced an unusual admixture event in the Zhejiang wolf. In many analyses, this wolf behaved similarly to other gray wolves in China, particularly those from SC. However, tests of admixture suggest that the Zhejiang wolf shows gene flow from a canid that diverged earlier than the separation of wolves and jackals. D-statistic analysis suggests that this may be from the Dhole, a species distributed in SC and Southeast Asia (Iyengar et al., 2005), or another canid that separated earlier than the divergence between wolves and the Dhole. Whatever the source of this ancestry, estimates of the admixture proportion from this deeply diverging population are estimated to be ∼12%–14%. Our results, taken together with previous research (Gopalakrishnan et al., 2018), reveal that canids are an ideal system in which to study how gene flow can shape speciation in a genus, and highlight the need for greater study of ancient gray wolf populations.

Limitations of Study

Although our results showed support for a SC gray wolf population and strong evidence for interspecific gene flow to the W7_Zhejiang wolf, the source of the interspecific gene flow into W7_Zhejiang is unclear. The availability of few specimens makes sampling in appropriate regions difficult, and where specimens are available, obtaining high coverage genomic data is not easy owing to natural degradation processes and loss of genetic material from museum storage methods. More samples are needed in future studies to further understand the admixture dynamics in canids.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.