Whole-genome sequencing of six dog breeds from continuous altitudes reveals adaptation to high-altitude hypoxia.

The hypoxic environment imposes severe selective pressure on species living at high altitude. To understand the genetic bases of adaptation to high altitude in dogs, we performed whole-genome sequencing of 60 dogs including five breeds living at continuous altitudes along the Tibetan Plateau from 800 to 5100 m as well as one European breed. More than 150× sequencing coverage for each breed provides us with a comprehensive assessment of the genetic polymorphisms of the dogs, including Tibetan Mastiffs. Comparison of the breeds from different altitudes reveals strong signals of population differentiation at the locus of hypoxia-related genes including endothelial Per-Arnt-Sim (PAS) domain protein 1 (EPAS1) and beta hemoglobin cluster. Notably, four novel nonsynonymous mutations specific to high-altitude dogs are identified at EPAS1, one of which occurred at a quite conserved site in the PAS domain. The association testing between EPAS1 genotypes and blood-related phenotypes on additional high-altitude dogs reveals that the homozygous mutation is associated with decreased blood flow resistance, which may help to improve hemorheologic fitness. Interestingly, EPAS1 was also identified as a selective target in Tibetan highlanders, though no amino acid changes were found. Thus, our results not only indicate parallel evolution of humans and dogs in adaptation to high-altitude hypoxia, but also provide a new opportunity to study the role of EPAS1 in the adaptive processes.

The mechanisms of organismal adaptation to high-altitude hypoxia are of great interest during recent years. Highland wild animals have a long life history at high altitude, and the whole genomes of yak (descendants of wild yak) (Qiu et al. 2012), Tibetan antelope (Ge et al. 2013), snow leopard (Cho et al. 2013), and wild boar (Li et al. 2013) have been sequenced. In contrast, the human settlement history on highland is rather short, which dates from ∼25,000 yr ago (Zhao et al. 2009). Whole-genome genotyping and re-sequencing have been performed for three typical highland populations including Tibetans (Beall et al. 2010; Bigham et al. 2010; Simonson et al. 2010; Yi et al. 2010; Peng et al. 2011; Xu et al. 2011), Andeans (Bigham et al. 2009, 2010), and Ethiopians (Alkorta-Aranburu et al. 2012; Scheinfeldt et al. 2012).

The increased oxygen uptake and delivery are physiological hallmarks of high-altitude adaptation. On one hand, the capacity of oxygen uptake is determined by hemoglobin concentration and oxygen affinity. For example, the Andean highlanders display a high level of hemoglobin concentration (Beall et al. 2002; Beall 2007). The high oxygen affinity of hemoglobin is found in many highland animals such as yak (Weber et al. 1988), alpaca (Piccinini et al. 1990), deer mice (Storz et al. 2007; Storz et al. 2009), bar-headed goose (Zhang et al. 1996; Liang et al. 2001), and Andean goose (Jessen et al. 1991). On the other hand, the rate of oxygen delivery is determined by blood flow. For example, although Tibetans maintain a nearly normal level of hemoglobin concentration and a low level of oxygen saturation, they display a high level of blood flow, resulting in the increase of oxygen delivery (Beall et al. 2001; Erzurum et al. 2007).

Whole-genome scans revealed that positive selection for human high-altitude adaptation occurred in the hypoxia-inducible factor (HIF) pathway (Bigham et al. 2009; Beall et al. 2010; Bigham et al. 2010; Simonson et al. 2010; Yi et al. 2010; Peng et al. 2011; Xu et al. 2011; Alkorta-Aranburu et al. 2012; Scheinfeldt et al. 2012), which regulates genes associated with blood physiology. In addition, metabolic pathways may also be involved in the adaptive process of yak (Qiu et al. 2012) and Tibetan antelope (Ge et al. 2013).

Although a lot of studies focused on wildlife and human highlanders, no research was performed on domesticated animals that migrated to the plateau with humans, which represent an adaptation pattern on a short evolutionary time scale of thousands of years. For example, the Tibetan Mastiff is a native dog living in the Tibetan Plateau with an altitude of 3000–6000 m. It is also an ancient dog in the world (Li and Zhang 2012). However, the genetic and physiological mechanisms of its adaptation to high-altitude environments remain elusive.

In this study, we sampled five dog breeds including the Tibetan Mastiff from continuous altitudes along the Ancient Tea Horse Road in southwestern China as well as one European breed. We performed whole-genome sequencing for the dogs and identified candidate genes for high-altitude adaptation using selective sweep mapping. We also measured the hematologic and hemorheologic parameters of the dogs and tested the association between the candidate alleles and blood physiology.

Results

Genetic diversity and population structure

We sampled 60 dogs belonging to six breeds living at continuous altitudes (Table 1). To avoid potential kinships, each dog was sampled from one individual village, which is dispersed along the Ancient Tea Horse Road and far away from each other. The sex ratio was kept as 1:1 for each breed. We classified the breeds into three altitude levels, i.e., high (>3000 m), middle (2000–3000 m), and low (<2000 m). The high-altitude breeds include Tibetan Mastiffs (TM) and Diqing indigenous dogs (DQ). The middle-altitude breed was chosen as Lijiang indigenous dogs (LJ). The low-altitude breeds contain Yingjiang indigenous dogs (YJ), Kunming dogs (KM), and German Shepherds (GS). Grouping multiple breeds living at similar altitudes could contribute to reducing the influence of breed-specific genetic background.

Table 1.

Sample, sequencing, and polymorphism information for the six dog breeds

We performed whole-genome sequencing for each individual dog. The coverage depth after genome alignment is ∼15× for each dog (Supplemental Fig. 1; Supplemental Tables 1, 2), resulting in >150× for each breed for variant calling (Supplemental Tables 3–5). The levels of genome-wide genetic diversity, as measured by Watterson’s θ and pairwise nucleotide diversity π, are relatively higher in the TM and indigenous breeds (DQ, LJ, and YJ) than those in KM and GS (Table 1; Supplemental Fig. 2), consistent with the fact that the former are more ancient breeds than the latter working dogs that have been domesticated recently (Vonholdt et al. 2010; Wang et al. 2013). In agreement with this, the minor allele frequency (MAF) distribution and Tajima’s D indicate a larger proportion of low frequency minor alleles in the ancient breeds than in the modern breeds (Supplemental Fig. 3). Additionally, linkage disequilibrium (LD) and haplotype analyses show that the ancient breeds have a faster decay of pairwise correlation coefficient (r) and lower integrated haplotype homozygosity (iHH) (Voight et al. 2006) than the modern breeds (Supplemental Fig. 4).

To examine genetic relationships among the dogs, we conducted principal component analysis (PCA) based on whole-genome SNPs. In the PCA plot, the LJ, YJ, KM, and GS breeds could be separated by the first and second eigenvectors, while TM and DQ are mixed (Fig. 1A). However, when more eigenvectors are incorporated, most dogs in TM and DQ could still be separated (Fig. 1A). These results suggest close but different genetic backgrounds between TM and DQ. To resolve their phylogenetic relationships, we constructed a neighbor-joining (NJ) tree based on the pairwise genetic distances (Fig. 1B). The tree also splits TM and DQ into distinct branches except one dog in TM. Considering possible admixture among the breeds, we further performed population structure analysis with frappe (Tang et al. 2005), which estimates individual ancestry and admixture proportions assuming K ancestral populations (Fig. 1C). When K = 2, we found a division between the ancient (TM, DQ, LJ, and YJ) and modern breeds (KM and GS). When K = 3, we found a division between the breeds living at high altitude (TM and DQ) and at middle/low altitude (LJ and YJ). When K = 5, TM and DQ are separated, though a slight extent of ingression from DQ happened in the TM breed.

Figure 1.

Genetic relationships and population structure of the 60 dogs based on all autosomal SNPs. (A) Principal component plot. The first (PC1) and second component (PC2) are shown in the main figure, while the third (PC3) and fourth component (PC4) are shown in the inset. The percentages indicate the proportion of variance explained by each component. (B) Unrooted neighbor-joining tree. The evolutionary distance is measured by the number of net nucleotide substitutions between individuals. (C) Population structures with the number of ancestral clusters K from 2 to 6. Each color represents one ancestral cluster and each vertical bar represents one dog. The length of colored segments represents corresponding ancestry attributions.

Signature of selection associated with high altitude

To detect the signature of natural selection associated with high altitude, we searched the dog genome for regions with a high coefficient of nucleotide differentiation FST among the breeds living at high (TM, DQ), middle (LJ), and low altitude (YJ, KM, GS) (Fig. 2A). The approach of population differentiation has been widely used to detect selective sweeps based on whole-genome SNPs, such as artificial selection during dog domestication (Axelsson et al. 2013) and diversification (Akey et al. 2010; Vonholdt et al. 2010). We scanned the autosomes with a nonoverlapping 100-kb window and calculated the FST value for each window. As previous studies (Axelsson et al. 2013), we focused on the regions with extremely high Z-transformed FST values [Z(FST) > 5] in the genome-wide empirical distribution. In total, 28 unique autosomal regions containing 141 candidate genes were identified (Supplemental Table 6), and five genes of them including EPAS1, MSRB3, HBB, CDK2, and GNB1 belong to the gene ontology (GO) categories “response to oxygen levels” and “response to oxidative stress” (Fig. 2A). To date the sweeps during the migration of dogs from low to high altitude, we also calculated the pairwise FST and Z(FST) between the three altitude levels (Supplemental Fig. 5; Supplemental Table 7).

Figure 2.

Selective sweep analysis of the breeds from different altitude levels. (A) Manhattan plot of FST among the high- (TM, DQ), middle- (LJ), and low- (YJ, KM, GS) altitude breeds. The FST was calculated for each 100-kb autosomal window. The dashed line denotes a threshold of Z(FST) = 5. (B) Diversity π and iHH around the EPAS1 loci. The region with Z(FST) > 5 is shaded. π and iHH were calculated for each 10-kb window and Z-transformed on the genome scale. The values were smoothed by lowess regression. (C) π and iHH around the beta hemoglobin loci.

The region with the strongest differentiation signal [Z(FST) = 10.83] (Fig. 2A) among the three altitude levels incorporates endothelial Per-Arnt-Sim (PAS) domain protein 1 (EPAS1), a gene encoding the HIF 2A (Patel and Simon 2008). The region is also highly differentiated when the high-altitude breeds were compared with the low- and middle-altitude ones [Z(FST) = 10.97 and 8.02, respectively] (Supplemental Fig. 5), suggesting that it is involved in a selective sweep during the move of dogs from middle to high altitude (Supplemental Table 7). Interestingly, EPAS1 was also identified as a selective target in Tibetan people (Beall et al. 2010; Yi et al. 2010; Peng et al. 2011; Xu et al. 2011), which is associated with their hemoglobin concentration (Beall et al. 2010; Yi et al. 2010). Network analysis indicated that the other candidate hypoxia-response genes we identified would all be regulated by HIF signaling pathway (Supplemental Fig. 6), suggesting an essential role of EPAS1 in the adaptation of high-altitude dogs.

The region with the second largest FST value [Z(FST) = 9.75] (Fig. 2A) among the three altitude levels is overlapped with methionine sulfoxide reductase B3 (MSRB3). Both the high- and middle-altitude breeds show high differentiation to the low-altitude ones in this region [Z(FST) = 9.77 and 6.75, respectively] (Supplemental Fig. 5), implying that the selective sweep occurred on the migration of dogs from low to middle altitude (Supplemental Table 7). MSRB3 has been identified as a gene associated with ear types in dogs (Boyko et al. 2010; Vaysse et al. 2011), but its association with altitude has not been reported previously. The beta hemoglobin (HBB) cluster region is also among the top three outliers [Z(FST) = 8.51] (Fig. 2A) and has a large differentiation between the high- and low-altitude breeds [Z(FST) = 8.84] (Supplemental Fig. 5), indicating a selective sweep on the dogs’ settlement from low to high altitude (Supplemental Table 7). The increases in oxygen affinity of hemoglobin have been found in diverse animal species to mediate the adaptive response to high altitude hypoxia (Storz and Moriyama 2008).

In addition to the increase of population differentiation, selective sweeps could result in the decrease of genetic diversity and increase of LD around the selective target (Nielsen 2005). Therefore, we compared π and iHH of the six breeds around the EPAS1 (Fig. 2B) and HBB loci (Fig. 2C), respectively. The π value of TM is decreased in the EPAS1 region and is the lowest among the six breeds. The iHH value of TM is increased in the EPAS1 region and is the highest among the six breeds (Fig. 2B). Both of the results supported that a selective sweep occurred in the EPAS1 region in TM. The DQ breed also shows an increased iHH value in the EPAS1 region as TM. However, the π value is increased in DQ, which may be caused by genetic admixture (Fig. 2B). The diversity of LJ in the EPAS1 region is the highest among the six breeds, suggesting extensive admixture and weak directional selection in the middle altitude dogs (Fig. 2B). When we inspected the beta hemoglobin cluster region, we observed a decrease of π value and increase of iHH value in both TM and DQ (Fig. 2C), confirming the FST screening that a selective sweep occurred in the loci of high-altitude dogs.

Mutations in the EPAS1 gene and physiological association

We examined the mutations in the EPAS1 gene in detail and found four nonsynonymous mutations in our sample, which cause amino acid variants Gly305-Ser (G305S), Asp494-Glu (D494E), Val500-Met (V500M), and Pro750-Ser (P750S) in the translated protein (based on Ensembl annotation) (Fig. 3A; Supplemental Table 8; Flicek et al. 2014). The four mutant alleles are in complete LD (pairwise r = 1) in all breeds, which are part of an 8-kb LD block (Supplemental Fig. 7). Haplotype tree (Fig. 3B) and network (Supplemental Fig. 7) of the LD block show that the haplotypes comprising the four mutant alleles all belong to the high-altitude clade, which have a long evolutionary distance from the haplotypes of the low-altitude clade. Also, the phylogenetic tree of the six breeds reconstructed from the LD block show separated clades between high and middle/low altitude (Supplemental Fig. 8), which is quite distinct from the genome-wide tree (ancient-modern separation). Interestingly, the frequency of the four mutant alleles is increased with the elevated altitude. It keeps rare in low-altitude breeds (KM 0%, GS 0%, YJ 5%), moderate in middle-altitude breeds (LJ 30%) but reaches high frequency in high-altitude breeds (DQ 80%, TM 100%) (Fig. 3C). To confirm this observation based on high-throughput sequencing, we performed genotyping of the four alleles on additional samples with Sanger sequencing and the result was similar: KM 0% (n = 11), GS 0% (n = 10), YJ 0% (n = 20), LJ 20% (n = 32), DQ 73% (n = 30), TM 94% (n = 35) (Fig. 3C). The high divergence of the haplotypes and alleles among the breeds indicates that the nonsynonymous mutations may be responsible for the selective sweep at EPAS1.

Figure 3.

EPAS1 mutations in the coding region. (A) Structural and evolutionary analysis of the four amino acid variants. The protein coordinate is based on Ensembl ID ENSCAFP00000003819. The coordinate of NCBI RefSeq XP_531807.2 is slightly different but the variants remain the same. The upper panel shows the Pfam domains of the protein. The orthologous protein sequences from 17 vertebrates are aligned with the mutant residues shown in red. The NJ tree derived from the multiple alignment is shown in the left panel. The red stars indicate two hydroxylation sites (Pro405 and Pro531), which are essential for oxygen sensing (Patel and Simon 2008). (PAS) Per-Arnt-Sim; (HIF) hypoxia-inducible factor; (CTAD) C-terminal transactivation domain. (B) NJ tree of the 8-kb haplotypes comprising the four nonsynonymous mutations. The evolutionary distance is defined as the number of different nucleotides between haplotypes. The haplotypes from high- and low-altitude breeds were grouped into two distinct clades. (C) Percentages of homozygotes and heterozygotes in the six breeds. The mutant and reference alleles are represented by “+” and “−”, respectively. Besides the 60 dogs sequenced with the Illumina technology, additional samples were genotyped with Sanger sequencing. The sample size is shown beside the bars. (D) Association between EPAS1 genotypes (mutant allele, A; reference allele, G) and blood flow resistance in DQ. The blood flow resistance was measured at different shear rates. In the boxplot, the line within the box defines the median; the ends of the boxes define the 25th and 75th percentiles; and the error bars define the 10th and 90th percentiles. The ANOVA F-test was performed, taking the age as covariant.

Among the four variants, only one (G305S) occurred in a well-defined protein domain (PAS domain, Fig. 3A). This domain facilitates heterodimerization with the aryl hydrocarbon receptor nuclear translocator (ARNT), resulting in a functional transcriptional complex (Patel and Simon 2008). Taking the crystal structure of PAS domain in human EPAS1 as template, our homology modeling suggested that the mutation G305S occurred in a beta sheet, which may affect the thermodynamic stability of the domain (Supplemental Fig. 9). To further evaluate the functional impact of the variants, we aligned the mutant EPAS1 protein with its ortholog proteins in diverse vertebrates (Fig. 3A). The comparison reveals that G305S is also a quite conserved amino acid mutation, which is invariant among all the other animals we examined. In contrast, the remaining three variants could all be found in other species. Prediction of functional effects of the variants supports that only G305S is deleterious, while the others are tolerated (Supplemental Table 8). All of these results imply that G305S is the most likely causal mutation for the EPAS1 sweep in high-altitude dogs.

To explore the physiology of the dogs’ adaptation to high-altitude hypoxia, we measured the hematologic and hemorheologic parameters for indigenous dogs living in a range from 500 to 3000 m (Supplemental Table 9). There are evident tendencies of the blood-related traits with the increase of altitude. Especially, the adult hemoglobin concentrations among middle- and high-altitude dogs reach 170 g/L as compared with 130 g/L among low-altitude dogs, which is a clear advantage in the rise of arterial oxygen content. However, with the increase of hemoglobin concentration, the blood viscosity and vascular resistance are also increased (Supplemental Table 9), which brings about potential risks for cardiovascular systems (Somer and Meiselman 1993; Allen and Patterson 1995). As the association between SNPs at EPAS1 and blood-related traits has been well established in humans (van Patot and Gassmann 2011), we conducted association testing for the variant G305S in DQ (Supplemental Table 10), the high-altitude breed where enough homozygotes (n = 40) and heterozygotes (n = 29) could be collected. Although no evident relationship with hemoglobin concentration was found, we discovered a statistically significant association between the genotypes and vascular resistance (Fig. 3D). The homozygotes with two mutant alleles show decreased vascular resistance as compared with the heterozygotes. The results suggest that the G305S variant may play important roles in the improvement of blood fluidity, which confers advantages for oxygen delivery and cardiovascular functions.

Discussion

In this study, we performed whole-genome sequencing of 60 dogs including five breeds from continuous altitudes along the Ancient Tea Horse Road in the Tibetan Plateau as well as one breed of GS. We catalogued millions of SNPs of each breed for evolutionary and genetic researches. Specifically, this is the first study to characterize the genetic polymorphisms of TM, which are known as “oriental deiform dogs” in Tibet and one of the most ancient and ferocious dogs in the world (Li and Zhang 2012). Our population structure analyses revealed that TM have a close but different genetic background with other indigenous dogs living at high altitude (Fig. 1), though their origin and relationship with plateau wolves need further studies. We also found the existence of admixture among TM and indigenous dogs (Fig. 1), which raised the demands for genetic testing in the identification and breeding of TM.

We conducted selective sweep mapping for the breeds from different altitudes and identified several candidate regions with a high extent of differentiation on the genome scale (Fig. 2A). Two of the genes contained in the top regions, EPAS1 and HBB, were previously reported as targets of natural selection for hypoxia adaptation in humans and other animals (Storz and Moriyama 2008; van Patot and Gassmann 2011). The two loci also displayed reduced genetic diversity and/or increased LD in high-altitude dogs especially TM (Fig. 2B,C). These results suggest that, driven by the same hypoxic environment, natural selection could act on a similar group of genes in the adaptive process of different species. Interestingly, the parallel evolution on the molecular level may extensively exist between humans and dogs in metabolic and neurological processes, probably due to their similar environments in the recent past (Wang et al. 2013).

We discovered and validated four novel nonsynonymous mutations at EPAS1 in high-altitude dogs, of which G305S occurred in a well-defined domain and quite conserved site (Fig. 3A). The mutant alleles and corresponding haplotypes have a large divergence between low- and high-altitude dogs (Fig. 3B,C), implying that the amino acid variations may be responsible for the signal of selective sweep at the EPAS1 locus.

The transcription factor encoded by EPAS1 plays a key role in transcriptional response to hypoxia in various physiological and pathologic conditions, including up-regulating the expression of erythropoietin (EPO) (Patel and Simon 2008). Especially, mutations at EPAS1 in humans are strongly associated with hematologic phenotypes (van Patot and Gassmann 2011). On one hand, clinical case studies found that several gain-of-function mutations occurred in close vicinity to the primary hydroxylation site (Pro531) at EPAS1 cause erythrocytosis (Percy et al. 2008a,b). On the other hand, the EPAS1 polymorphisms unique to the native Tibetan people were found to be associated with their lower hemoglobin concentrations (Beall et al. 2010; Yi et al. 2010), suggesting a loss-of-function role of EPAS1 in high-altitude adaptation. As erythrocytosis is a common symptom of chronic mountain sickness which will lead to high blood viscosity and cardiovascular disorders, the decrease in hemoglobin level may provide a protective mechanism for Tibetan people (Beall et al. 2010; van Patot and Gassmann 2011). In fact, comparing with Andeans and other highland migrants with shorter settlement history, Tibetans do not show markedly elevated hemoglobin concentration (Beall 2007), probably due to their unique EPAS1 genotypes (van Patot and Gassmann 2011). Instead, Tibetans display a higher blood flow for oxygen delivery (Erzurum et al. 2007), which could compensate the effect of lower hemoglobin concentration.

Nonetheless, it is worth noting that all EPAS1 variations detected previously in the Tibetan population are at introns (Beall et al. 2010; Yi et al. 2010), inhibiting further study of its molecular mechanism in the adaptive process. Fortunately, the key amino acid mutation we identified at EPAS1 in dogs, G305S, is predicted to be damaging, which is also likely to cause the loss of function of EPAS1. Moreover, we demonstrated that the variation is associated with lower blood flow resistance in high-altitude dogs (Fig. 3D). Therefore, we propose that similar mechanisms may be adopted by both humans and dogs in adaptation to high-altitude hypoxia. Although the physiological responses need further comparison between humans and dogs, these novel variations of EPAS1 could shed new insights into the study of the gene function.

Methods

Sample preparation and sequencing

For each dog, the genomic DNA was extracted from 200 μL peripheral venous blood with the QIAamp DNA Blood Mini kit (Qiagen). The quality and integrity of DNA was controlled by A260/280 ratio and agarose gel electrophoresis. For sequencing library preparation, the genomic DNA was sheared to fragments of 300–500 bp, which were then end-repaired, “A”-tailed, and ligated to Illumina sequencing adapters. The ligated products with sizes of 400–500 bp were selected on 2% agarose gels and then amplified by LM-PCR. The libraries were sequenced on Illumina HiSeq 2000 with 2 × 100 bp paired-end mode, which was controlled by Illumina HiSeq Control Software.

Reads processing and variant calling

The raw reads were processed by two rounds of quality control (QC) (Supplemental Table 1). Firstly, low quality reads marked by Illumina sequencers in the FASTQ files (Casava v1.8) were filtered. Secondly, for each read, the low quality end with base quality scores <20 were trimmed. The read was removed if its length was <35 after trimming. Only paired reads were preserved during QC. Reads after QC were mapped to the dog reference genome assembly canFam3 (http://hgdownload.soe.ucsc.edu/goldenPath/canFam3/bigZips/canFam3.fa.gz) using BWA-MEM (v0.7.4) (Li and Durbin 2009) with default parameters on a per-individual basis. We removed duplicated reads and calculated mapping statistics including the coverage of depth from the .bam files by SAMtools (v0.1.19) (Li et al. 2009; Supplemental Table 2).

The SNPs and small indels were called using the SAMtools pipeline (Li 2011) on a per-breed (10 individuals) basis (Supplemental Table 3). Firstly, .mpileup files were generated by SAMtools mpileup with the parameters “-u -C50 -DS -q20”. Secondly, .vcf files containing raw variants were generated by bctools view with the parameters “-evcgN”. Finally, the raw variants were filtered by vcfutils.pl varFilter with the minimum depth “-d 20” and maximum depth “-D 300”. The cutoffs were set according to the sequencing depth distribution per individual (Supplemental Fig. 1). Other parameters for filtering were default, including P-value of Hardy–Weinberg equilibrium (0.0001). Given the position of a variant, the genotype with the maximum posterior probability of each individual was given in the .vcf file. The variants after filtering were annotated by ANNOVAR (v2013-06-21) (Wang et al. 2010) according to Ensembl (Flicek et al. 2014) gene annotation (http://hgdownload.soe.ucsc.edu/goldenPath/canFam3/database/ensGene.txt.gz) (Supplemental Tables 4, 5).

Population genetics analysis

The .vcf files for the six breeds were merged by vcf-merge in VCFtools (v0.1.11) (Danecek et al. 2011). When a genotype was missing in one breed, it was set to the reference allele with the option “-R 0/0”. Only biallelic SNPs in autosomes were preserved for subsequent population structure and polymorphism analysis. The merged .vcf file was converted to tab-delimited text file by vcf-to-tab in VCFtools, and to PLINK format files (.ped and .map) by VCFtools and PLINK (v1.07) (Purcell et al. 2007) when appropriate. The PCA was performed by GCTA (v1.13) (Fig. 1A; Yang et al. 2011), which first generated the genetic relationship matrix (GRM) with the “–make-grm” option and then estimated the first four principal components with the “–pca 4” option. The pairwise genetic distances between individuals were measured by the number of net nucleotide substitutions (formula 12.67 in Nei and Kumar 2000), which were calculated for each SNP sites and then summed up. Based on the distance matrix, the NJ tree was constructed and displayed by MEGA (v5.2) (Fig. 1B; Tamura et al. 2011). The population structure was inferred by frappe (v1.1) (Tang et al. 2005) with a maximum likelihood method (Fig. 1C). The number of ancestral clusters K ranged from 2 to 6, with 10,000 iterations each run. The graphs of the population structure were displayed by distruct (v1.1) (Rosenberg 2004) on a per-individual basis.

Standard population genetic statistics, including Watterson’s θ, pairwise nucleotide diversity π (Supplemental Fig. 2), and Tajima’s D (Supplemental Fig. 3) were calculated for each breed by the Bio∷PopGen∷Statistics package in BioPerl (v1.6.1) (Stajich et al. 2002). The coefficient of nucleotide differentiation FST among the breeds was calculated by the Bio∷PopGen∷PopStats package in BioPerl (Supplemental Table 6). The statistics were estimated over a sliding window with a fixed size along the genome. The Manhattan plot of genome-wide FST values was generated using the mhtplot function in the R package gap (Fig. 2A; Supplemental Fig. 5). The level of LD measured by the correlation coefficient (r) between unphased markers (Supplemental Fig. 4A) was calculated by Haploview (v4.2) (Barrett et al. 2005) with the option “-dprime -maxdistance 50 -minMAF 0.1 -hwcutoff 0.05”. We also phased the genotypes of all individuals into haplotypes by fastPHASE (v1.4.0) (Scheet and Stephens 2006) with the parameters “-T10 -K8 -u”. The iHH (Voight et al. 2006) for each breed was calculated by XP-EHH (Pickrell et al. 2009) based on the phased data (Supplemental Fig. 4B). The iHH values were log-transformed and averaged over a sliding window with fixed size along the genome.

EPAS1 mutation analysis

The protein domains of EPAS1 were predicted by the Pfam web service (Fig. 3A; Punta et al. 2012). The ortholog proteins of EPAS1 in 17 representative vertebrates were retrieved from Ensembl (Flicek et al. 2014) (release 72) and only one-to-one orthologs were included (Fig. 3A). The multiple sequence alignment was performed by ClustalX (v2.1) (Larkin et al. 2007) and displayed by Jalview (v2.8) (Waterhouse et al. 2009). The LD block containing the nonsynonymous mutations was detected and visualized by Haploview (v4.2) (Supplemental Fig. 7A; Barrett et al. 2005). The haplotype NJ tree was constructed by MEGA (v5.2) (Tamura et al. 2011) with the default parameters (Fig. 3B). The haplotype network was constructed and plotted by Network (v4.611, Fluxus Engineering) (Supplemental Fig. 7B). It was calculated with the median-joining method (Bandelt et al. 1999) (ɛ = 0), followed by reduction with the MP option (Polzin and Daneshmand 2003). The sequence of dog EPAS1 was searched against the PDB database (Rose et al. 2013) using BLAST-P (Altschul et al. 1997) to identify suitable templates for homology modeling (Supplemental Fig. 9). I-TASSER (v2.1) (Roy et al. 2010) was used for structure modeling by threading the sequence onto the selected templates. The three-dimensional model structures were visualized with PyMOL (v1.6, Schrödinger, LLC). Structure-based energy calculation was performed by the Eris server (Yin et al. 2007) using the flexible backbone method upon single amino acid substitutions. The functional consequences of amino acid mutations were predicted by SIFT (Kumar et al. 2009), which was invoked through Ensembl’s (Flicek et al. 2014) Variant Effect Predictor (Supplemental Table 8).

Physiological measurement and association test

The foreleg venous blood was collected to measure the physiological parameters in situ. The hematologic and hemorheologic parameters were measured by BC-2800Vet Auto Hematology Analyzer (Mindray Co., Ltd.) and ZL1000 Auto Blood Rheology Analyzer (Zonci Co., Ltd.), respectively (Supplemental Tables 9, 10). The nonsynonymous mutations at EPAS1 were genotyped by traditional Sanger sequence technology. The association testing was performed with each physiological parameter as the response variable, genotype as the predictive variable, and age (<1 yr vs. >1 yr) as the covariate. We also modeled the gender as a covariate but no difference was found in this factor. The F test in the analysis of variance (ANOVA) was used to evaluate the significance of genotypes (Fig. 3D).

Data access

All sequencing data from this study have been submitted to the NCBI Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra) under accession number SRP035294.