EPAS1 variants in high altitude Tibetan wolves were selectively introgressed into highland dogs

Samples and genome sequence data

To determine if highland dogs inherited a beneficial allele through admixture with wolves, we obtained publically available genome sequence and SNP genotype data from 21 highland (>3,300 m) dogs (n = 10 indigenous dogs from a putatively admixed population, n = 11 Tibetan mastiffs) in China to compare to 15 reference dogs (n = 13 indigenous Chinese dogs; n = 1 basenji, n = 1 dingo), 22 low elevation Old World (OW) wolves, and six highland wolves from China (Fig. 1, Table S1). For this study, we obtained the genomic sequence and SNP variants for two genomes: one low elevation wolf from Inner Mongolia (im5) and one highland wolf from Tibet (ti3) (Table S1). The whole genome sequencing was performed using an Illumina Hiseq 2,000 at Novogene (China). Two paired-end libraries with insert sizes of 300∼500 bp were generated for each sample. Library preparation and sequencing was performed according to manufacturer’s protocols. Our newly generated sequence data can be found on NCBI’s Sequence Read Archive (SRA) under the accession number SRP096612. The 150 bp pair-end reads of each sample were aligned to the dog genome (CanFam3.1) using Bowtie 2 (Langmead & Salzberg, 2012) under the local alignment algorithm with the very sensitive model. Default options were used for all other parameters. Then, we applied Picard and GATK toolsets (DePristo et al., 2011) to process the alignments to identify SNPs. The pipeline is the same as used in our previous studies (Freedman et al., 2014; Zhang et al., 2014; Fan et al., 2016).

To obtain a set of polymorphic SNP genotypes, we excluded monomorphic sites and applied a strict inclusion threshold of no missing data to obtain 17,709 SNPs on canine chromosome 10. Further, to increase SNP density around the candidate gene EPAS1, we allowed on average 14% missing data, which provided an additional 1,113 SNP genotypes (total number of SNPs = 18,792). To assess if our assumed reference dogs were indeed non-admixed, we used the Ancestry function in Galaxy’s Genomic Diversity toolkit for K = 2–4 (Alexander, Novembre & Lange, 2009; Afgan et al., 2016) and performed a complete population structure analysis for all samples with 17,709 SNP genotypes.

Admixed ancestry analysis

To better explore the hypothesized introgression of EPAS1 in to highland dogs, we utilized ancestry information from a larger reference panel of dogs and gray wolves, and used a different methodological approach than previous work. Specifically, following methods in vonHoldt et al. (2011), we estimated the degree of per-SNP genetic differentiation (F_ST) between the reference dog and Old World wolves and employed a 2 standard deviation (SD) threshold to identify SNPs with significantly high levels of divergence. We then identified SNPs with strongly divergent allele frequencies between the two reference groups of dogs and wolves. Namely, sites with the highest F_ST values, as defined by the 2SD threshold, were classified as ancestry informative markers (AIMs).

We assigned either dog or wolf ancestry to chromosome 10 SNPs from 21 highland dogs that are putatively admixed with OW wolves. Samples were assessed for outlier patterns using principal component analysis (PCA) in Galaxy’s Genomic Diversity toolkit (Afgan et al., 2016) on both the set of 17,709 and 980 AIM SNPs. Further, we surveyed local ancestry proportions for highland dogs using the program Admixture in Galaxy’s Genomic Diversity toolkit from the estimated genotype, two reference populations of dogs and wolves, and a switch penalty value of 10 to estimate per-site ancestry of the putatively admixed highland dogs (Table S1). We then surveyed average wolf ancestry across chromosome 10 in 100-SNP windows with a 50-SNP step in the putatively admixed highland dogs. Outlier regions were identified by a simple outlier threshold of wolf ancestry proportion >2SD above the chromosomal level (average wolf proportion, SD = 0.032, 0.09). We further repeated the PCA and completed a population structure analysis by using the Ancestry function in Galaxy (Alexander, Novembre & Lange, 2009) for SNP genotypes within any putatively introgressed regions. To phase our SNP genotypes using SHAPEIT (Delaneau, Marchini & Zagury, 2012), we constructed a genetic map file as the ratio of inter-SNP distances to the position of the last SNP queried on chromosome 10. Additionally, we phased SNP genotypes using the program PCAdmix and plotted the per-individual loadings for the first two coordinates (Henn et al., 2012). Lastly, we applied a haplotype-based test (Mendez, Watkins & Hammer, 2012a; Mendez, Watkins & Hammer, 2012b) to determine if the shared EPAS1 haplotype in highland dogs and wolves was due to an adaptive introgression event or the alternative hypothesis that selection acted on a haplotype present in the ancestral population. We utilized a recently inferred genetic map (Campbell et al., 2016). We assume a generation time of three years (Freedman et al., 2014; Fan et al., 2016) and that all recombination events are detectable.

Variant effect predictions were conducted with transcripts annotated from the reference dog genome (CanFam 3.1), with SIFT (Sorting Intolerant From Tolerant) scores on non-synonymous substitutions, which range from 0 to 1, used to predict if the non-synonymous variants may be damaging. Specifically, if the SIFT score is ≤0.05 it is considered damaging, or tolerated when the score is >0.05 (McLaren et al., 2010).

Scan for positive selection

We scanned the target genomic region for signals of positive selection in the highland dog genome compared to the reference dogs using cross population extended haplotype homozygosity (XP-EHH) tests (Sabeti et al., 2007). This analysis was based on the phased haplotypes from 836 SNPs that had no missing data within the putatively introgressed region containing PRKCE and EPAS1. Per-SNP F_ST was calculated as described above. We normalized both the F_ST and XP-EHH scores as z-scores with a mean of zero and standard deviation of 1. Their product represents a composite “bivariate percentile score”. We identified outlier loci as those in the 97.5th percentile or greater (z-score > 2) of the bivariate percentile score distribution.

Two regions are enriched for wolf ancestry in highland dogs

To determine if admixture mapping is useful in searching for adaptively introgressed alleles in Chinese highland dogs, we first conducted a survey of genetic structure for 64 Old World canids along chromosome 10 (SNPs n = 17, 709) (Fig. 1; Table S1). We found that the reference wolves represent two genetic clusters, high and low elevation individuals, with little evidence for admixture with dogs (Figs. 2, 3A). We further identified the reference dogs to be largely non-admixed with wolves, with the only exception being the basenji (bas; Fig. 2). Although the highland and reference dogs were similar in the PCA analysis (Fig. 3A), these two populations were distinct at K = 4 and displayed varying degrees of admixture with each other (Fig. 2). Further, highland dogs contained low chromosome-wide proportions of wolf ancestry on chromosome 10 (max = 0.027) (Table 1).

We surveyed ancestry across multiple dog and wolf populations based on 18,794 SNP genotypes on canine chromosome 10 to provide a higher density survey across the candidate gene, EPAS1. Using the 2SD threshold, we classified 980 F_ST outlier SNPs as AIMs that were informative for dog and wolf ancestry (all SNPs: F_ST average = 0.04, SD = 0.10, +2SD = 0.23; AIMs: F_ST average = 0.38, SD = 0.13), with an average density of 1 AIM per 70Kb across chromosome 10 (median = 26 Kb). Though non-model based, a PCA of 980 AIMs clearly identified four clusters, partitioned by species and elevation (Fig. 3B; Fig. S1). A sliding window analysis of AIMs along chromosome 10 revealed an enrichment of wolf ancestry in two regions (chr10: 8147438–8288850 containing the MSRB3 gene, average wolf proportion = 0.278; chr10: 48151501–48761234: PRKCE and EPAS1 genes, average wolf proportion = 0.355), thus identifying them as candidate outliers (Fig. 4). Genetic analysis of each outlier region varies in the degree of structure (Figs. 3C and 3D; Fig. S1) and observed outliers are found in three putatively introgressed genes. Following previous work and with the objective to explore the introgressed history of EPAS1 in high elevation adaptation, we focused on the region containing PRKCE and EPAS1 in subsequent analyses. We conducted a PCA of phased haplotypes from this genomic region that contains both an excess of wolf ancestry in dogs and two functionally-relevant candidate genes. Our objective was to identify the potential source population of this putatively introgressed region in highland dogs. Low elevation wolves are spatially discrete compared to those from the highland, where 10 of 12 haplotypes are nested within the putatively admixed highland dogs (Fig. 5). This result suggests that the putatively introgressed region in the highland dog genome originated from highland wolves in Qinghai and Tibet, China. From the phased haplotypes of the introgressed region containing PRKCE and EPAS1, we defined a 396 Kb subregion (n_SNPs = 467; chr10: 48305073–48700082) that had had a high proportion of wolf ancestry (wolf prop. > 0.45) across all highland dogs implying the origin of this region derives from admixture with highland gray wolves.

To better understand functional variation in this introgressed region, we determined the presence of damaging variation. We restricted our assignment of functional changes to SNPs within this candidate subregion and found one missense damaging non-synonymous variants in each of the PRKCE and EPAS1 genes (Table 2). The PRKCE variant was only observed in a single highland wolf (qh4), and the EPAS1 variant found in a putatively admixed highland Tibetan mastiff (DQZA81) and one low elevation indigenous reference dog (dc5).

Using the recent canine genetic map, we inferred that the 60 Kb length block containing the EPAS1 gene has a genetic map length of 0.032 cM. Consequently, the time until one of the two identical-by-descent copies of a EPAS1 haplotype are broken by a recombination event is exponentially distributed with a parameter of (0.00032 *  2) per generation. We found that one of the two identical-by-descent copies of the EPAS1 haplotype has a 50%, 95% and 99% probability of at least one recombination event in 3,261, 14,092 and 21,663 years, respectively. On the other hand, estimates of the divergence time of dogs and wolves range from 11,700 or 29,000 years ago depending on the mutation rate used (Freedman et al., 2014; Skoglund et al., 2015; Fan et al., 2016). These divergence time estimates imply that recombination would have broken one of the two identical by descent copies of the EPAS1 haplotype at least once (with a 95% probability) if the ancestral haplotype was present in the population ancestral to dogs and wolves. Contrary to this expectation, highland dogs and wolves share the full 60 Kb EPAS1 haplotype, which is evidence in favor of an adaptive introgression event rather than selection on an ancestral haplotype in common ancestor of dogs and wolves.

Signal of selection in EPAS1

To confirm the signal of selection on EPAS1 as detected in previous studies, we conducted a cross-population analysis of extended haplotype homozygosity and genetic differentiation between the highland and low elevation dog genomes. This analysis supports the previous finding that EPAS1 is under strong positive selection (chr10: 48565127–48625694). However, this signal is discrete and does not encompass the adjacent gene, PRKCE (Fig. 6A), further supporting the specific role of EPAS1 in high elevation adaptation in domestic dogs. Additionally, phasing the 363 SNPs in EPAS1 also confirmed a shared haplotype between the highland dogs and wolves from Tibet (samples ti1, ti2, ti3) and Qinghai (samples qh1, qh2, qh4), China (Fig. 6B).