Exome sequencing reveals genetic differentiation due to high-altitude adaptation in the Tibetan cashmere goat (Capra hircus)

Abstract

Background: The Tibetan cashmere goat (Capra hircus), one of the most ancient breeds in China, has historically been a critical source of meat and cashmere production for local farmers. To adapt to the high-altitude area, extremely harsh climate, and hypoxic environment that the Tibetan cashmere goat lives in, this goat has developed distinct phenotypic traits compared to lowland breeds. However, the genetic components underlying this phenotypic adaptation remain largely unknown.

Results: We obtained 118,700 autosomal SNPs through exome sequencing of 330 cashmere goats located at a wide geographic range, including the Tibetan Plateau and low-altitude regions in China. The great majority of SNPs showed low genetic differentiation among populations; however, approximately 2-3% of the loci showed more genetic differentiation than expected under a selectively neutral model. Together with a combined analysis of high- and low-altitude breeds, we revealed 339 genes potentially under high-altitude selection. Genes associated with cardiovascular system development were significantly enriched in our study. Among these genes, the most evident one was endothelial PAS domain protein 1 (EPAS1), which has been previously reported to be involved in complex oxygen sensing and significantly associated with high-altitude adaptation of human, dog, and grey wolf. The missense mutation Q579L that we identified in EPAS1, which occurs next to the Hypoxia-Inducible Factor-1 (HIF-1) domain, was exclusively enriched in the high-altitude populations.

Conclusions: Our study provides insights concerning the population variation in six different cashmere goat populations in China. The variants in cardiovascular system-related genes may explain the observed phenotypic adaptation of the Tibetan cashmere goat.

Fig. 1

Phylogenetic analysis and Fst simulation test. a Phylogenetic tree analysis based on all markers. b Cumulative distribution of observed and simulated (assuming neutrality) Fst values. c Histogram of Fst values in the simulated and observed datasets (note the truncated y-axis). d Phylogenetic tree based on those global SNPs showing significant genetic differentiation

Fig. 2

Overlap of SNPs from the global SNP dataset and high-lowland dataset

Fig. 3

EPAS1 mutation in the coding regions. a EPAS1 protein sequence analysis. The protein coordinate is based on NCBI RefSeq XP_005686651.1. The upper panel shows the Pfam domains of the protein. The double arrows represent domains of goat EPAS1. The orthologous protein sequences from 17 vertebrates are aligned with the mutant residues shown in the box. Sheep, ENSOARP00000006140; cattle, ENSBTAP00000004836; pig, ENSSSCP00000009011; human, ENSP00000406137; dog, ENSCAFP00000003819; elephant, ENSLAFP00000010336; mouse, ENSMUSP00000024954; opossum, ENSMODP00000001136; zebra finch, ENSTGUP00000004086; anole lizard, ENSACAP00000004025; turkey, ENSGACP00000015093; xenopus, ENSXETP00000031612; zebrafish, ENSDARP00000074832; lamprey, ENSPMAP00000000148; stickleback, ENSGACP00000015093. (HLH) helix-loop-helix domain; (PAS) Per-Arnt-Sim; (HIF) hypoxia-inducible factor; (CTAD) C-terminal transactivation domain; (DAG1) Dystroglycan (Dystrophin-associated glycoprotein 1), the blue shade represent the HIF-1 domain. b Percentages of the reference allele and variant allele in larger population samples that were genotyped with Sanger sequencing technology. The number on the bar represents the sample size. The altitude information is shown at the right side of the bars. c Association analysis between EPAS1 genotypes (mutant allele, T; reference allele, A) and the MCHC in BG population. The ANOVA F-test was performed, and we found a significant association between the genotypes and MCHC