Whole-genome sequencing of native sheep provides insights into rapid adaptations to extreme environments

Whole-Genome Sequencing of Native Sheep Provides Insights into Rapid Adaptations to Extreme Environments Ji Yang, †,1 Wen-Rong Li, †,2 Feng-Hua Lv, †,1 San-Gang He, †,2 Shi-Lin Tian, †,3 Wei-Feng Peng, 1,4 Ya-Wei Sun, 2,5 Yong-Xin Zhao, 1,4 Xiao-Long Tu, 3 Min Zhang, 1,6 Xing-Long Xie, 1,4 Yu-Tao Wang, 7 Jin-Quan Li, 8 Yong-Gang Liu, 9 Zhi-Qiang Shen, 10 Feng Wang, 11 Guang-Jian Liu, 3 Hong-Feng Lu, 3 Juha Kantanen, 12,13 Jian-Lin Han, 14,15 Meng-Hua Li,* ,1 and Ming-Jun Liu* ,2 1 CAS Key Laboratory of Animal Ecology and Conservation Biology, Institute of Zoology, Chinese Academy of Sciences (CAS), Beijing, China 2 Animal Biotechnology Research Institute, Xinjiang Academy of Animal Science, Urumqi, China 3 Novogene Bioinformatics Institute, Beijing, China 4 University of Chinese Academy of Sciences (UCAS), Beijing, China 5 College of Animal Science and Technology, Shihezi University, Shihezi, China 6 School of Life Sciences, University of Science and Technology of China, Hefei, China 7 College of Biological and Geographic Sciences, Kashgar University, Kashgar, China 8 College of Animal Science, Inner Mongolia Agricultural University, Hohhot, China 9 College of Animal Science and Technology, Yunnan Agricultural University, Kunming, China 10 Shandong Binzhou Academy of Animal Science and Veterinary Medicine, Binzhou, China 11 Institute of Sheep and Goat Science, Nanjing Agricultural University, Nanjing, China 12 Green Technology, Natural Resources Institute Finland (Luke), Jokioinen, Finland l3 Department of Environmental and Biological Sciences, University of Eastern Finland, Kuopio, Finland 14 CAAS-ILRI Joint Laboratory on Livestock and Forage Genetic Resources, Institute of Animal Science, Chinese Academy of Agricultural Sciences (CAAS), Beijing, China 15 International Livestock Research Institute (ILRI), Nairobi, Kenya † These authors contributed equally to this work. *Corresponding author: E-mail: [email protected]n; [email protected]. Associate editor: Yuseob Kim Abstract Global climate change has a significant effect on extreme environments and a profound influence on species survival. However, little is known of the genome-wide pattern of livestock adaptations to extreme environments over a short time frame following domestication. Sheep (Ovis aries) have become well adapted to a diverse range of agroecological zones, including certain extreme environments (e.g., plateaus and deserts), during their post-domestication (approximately 8–9 kya) migration and differentiation. Here, we generated whole-genome sequences from 77 native sheep, with an average effective sequencing depth of 5for 75 samples and 42for 2 samples. Comparative genomic analyses among sheep in contrasting environments, that is, plateau (>4,000 m above sea level) versus lowland (<100 m), high-altitude region (>1500 m) versus low-altitude region (<1300 m), desert (<10 mm average annual precipitation) versus highly humid region (>600 mm), and arid zone (<400 mm) versus humid zone (>400 mm), detected a novel set of candidate genes as well as pathways and GO categories that are putatively associated with hypoxia responses at high altitudes and water reabsorption in arid environments. In addition, candidate genes and GO terms functionally related to energy metabolism and body size variations were identified. This study offers novel insights into rapid genomic adaptations to extreme environments in sheep and other animals, and provides a valuable resource for future research on livestock breeding in response to climate change. Key words: extreme environment, rapid adaptation, Ovis aries,whole-genome sequencing,climate change. Introduction Globalclimatechangehasaconsiderableimpactonorganisms, including livestock (Brown and Funk 2008;Hoffmann 2010). In particular, the influence of climate change on living organisms is exacerbated at the interface of extreme environments, such as on plateaus and in desert regions (Easterling et al. 2000). Within this context, it is important to understand the genetic basis of well-adapted local livestock breeds in extreme environments to develop appropriate breeding Article ßThe Author 2016. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] Open Access 2576 Mol. Biol. Evol. 33(10):2576–2592 doi:10.1093/molbev/msw129 Advance Access publication July 8, 2016 programs under scenarios of future climate change (Howden et al. 2007;Lobell et al. 2008). After domestication in the Fertile Crescent 8,000–9,000 years ago (8–9 kya; Legge 1996), sheep (Ovis aries) spread and became adapted to a wide range of agroecological conditions, especially those distributed on plateaus or in desert regions, which are sensitive toclimatechange(Easterling et al. 2000;Lv et al. 2014). Thus, these animals provide an excellent model to gain novel insights into genetic mechanisms underlying the rapid adaptations of livestock to extreme environments within a short period of time. In recent years, to characterize adaptive genetic variations, whole-genome sequencing studies have been performed on a wide range of organisms that live in harsh or extreme environments (Simonson et al. 2010;Mackelprang et al. 2011; Scheinfeldt et al. 2012;Cho et al. 2013;Ge et al. 2013;Ma et al. 2013;Qu et al. 2013;Liu et al. 2014). Studies conducted on livestock are limited, although they include work on adaptations to high altitudes in yak (Bos grunniens;Qiu et al. 2012), Tibetan mastiff (Canis lupus familiaris;Gou et al. 2014)andTibetanchicken(Wang et al. 2015), hot and arid environments in goat and sheep (Kim et al. 2016), severe desert conditions in domestic Bactrian camel (Jirimutu et al. 2012;Wu et al. 2014) and subarctic cold environments in Yakutian horse (Librado et al. 2015). However, to our knowledge, no study has characterized the rapid genetic adaptations of livestock to various extreme environments based on whole-genome sequences. Out of the dispersal center in the Mongolian region, Chinese native sheep breeds only have diverged for several thousands of years (e.g., Lv et al. 2015). Previous investigations showed significant differences in hematologic parameters (e.g., erythrocyte count, hematocrit, and hemoglobin; Jiang et al. 1991;Liu and Chen 2006)andbodysize(Du 2011) between Tibetan sheep on the Qinghai–Tibetan Plateau and those at relatively low altitudes as well as robust morphology (i.e., good body constitution; Du 2011)insheeporiginating from the Taklimakan Desert region. In this study, we sequenced the whole genomes of 77 sheep (O. aries) including those from habitats in extreme (orharsh)environments:Tibetan areas on the Qinghai– Tibetan Plateau (defined as ‘plateau’, altitude >4,000 m above sea level), high-altitude region (altitude >1,500 m), Taklimakan Desert region (defined as ‘desert’, average annual precipitation <10 mm), and arid zone (<400 mm, representing arid and semi-arid regions; Piao et al. 2010;fig. 1Aand B, supplementary tables S1 and S2,Supplementary Material online). The set of samples represented 21 native breeds of different genetic and geographic origins in China (supplementary table S1,Supplementary Material online). Through comparisons of the genomes of sheep from extreme environments with those from contrasting environments, we aimed to identify the candidate genes, functional Gene Ontology (GO) categories and signaling pathways responsible for the rapid adaptations (i.e., over thousands of years) of sheep to plateau and desert environments. In addition, to elucidate the evolutionary history of Chinese native sheep, a comprehensive analysis of the genomic diversity, population structure and demographic history of these animals was performed based on genomic data. Results and Discussion Genomic Variants We generated the genomic sequences of 77 native sheep and, for genomic comparison and as phylogenetic outgroups, three wild species of the subfamily Caprinae,vizO. a. musimon,Ovis ammon polii and Capra ibex (one animal from each species). This yielded 1,150 Gb of aligned highquality genomic data at a total of 438effective sequence depth for the subsequent analyses. An average of 94.66% (92.14–95.30%) of the SNPs identified in the 77 native sheep were validated in the sheep dbSNP database (supplementary table S3,Supplementary Material online), indicating the high reliability of the called SNP variations in this study (e.g., Park et al. 2015). For the 3,720,550 SNPs not confirmed by dbSNP database, 1,582,868 SNPs were further validated by their presence in 2–77 individuals of the samples. Of the common SNPs for our sequencing data and the Ovine 50K Beadchip array in the 33 native sheep (7,200– 15,016 SNPs in each individual), an average 95.51% (89.46–96. 21%) of the SNP genotypes were consistent (supplementary table S4 and fig. S1, Supplementary Material online), thereby demonstrating the high quality of our SNP calling. The 4.49% discrepancy was resulted from genotyping errors and missing data in the whole genome sequencing process (supplementary table S5,Supplementary Material online). Moreover, a large number of shared SNPs (15.58–17.48 million) between different sheep groups from different regions (i.e., Qinghai–Tibetan Plateau, Yunnan–Kweichow Plateau, Northern and Eastern China; supplementary table S6, Supplementary Material online) also support for the validity of the SNP calling. We obtained 21.26 million high-quality SNPs for 77 native sheep, 2.77 million for O.a.musimon, 4.71 million for O.a. polii, and 10.10 million for C.ibex (supplementary table S7, Supplementary Material online). There were 15.81 million high-quality SNPs unique to the native sheep (supplementary table S7,Supplementary Material online). Based on one individual per species, the overall numbers of SNPs per bp in the exonic, intronic, intergenic regions and chromosomes of domestic sheep (averaged over 77 individuals) were higher than those of O.a.musimon, but lower than those of O.a.polii and C.ibex (supplementary table S8,Supplementary Material online). Of the three identified genetic groups of native sheep, Northern and Eastern Chinese breeds possessed the largest amounts of the total (19,227,904) and unique (784,796) SNPs (supplementary table S6 and fig. S2A, Supplementary Material online), indicating their overall relatively high genomic diversity and early population history (Lv et al. 2015). Most of the high-quality SNPs in the individual genome of sheep were located in intergenic regions (2.678 million, 71.83%), with only 0.76% (0.028 million) located in exonic regions (supplementary table S8,Supplementary Material online). A total of 11,167 nonsynonymous SNPs and 16,879 synonymous SNPs were located within exons, which resulted in a Rapid Adaptation of Sheep to Extreme Environments .doi:10.1093/molbev/msw129 MBE 2577 AB CD E G F FIG.1.Geographic distribution and population genetics analyses of 21 native sheep breeds. (A) Sampling sites in this study. A total of 77 Chinese sheep representing 21 native breeds were included. The elevation (m) of the study area is also visualized. (B) Geographic variation of the annual mean precipitation (mm). Precipitation data during the 21 years from 1981 to 2001 were retrieved from the Abdus Salam International Center for Theoretical Physics, Italy (ICTP, https://www.ictp.it, last accessed January 5, 2015). The 200, 400, and 800 mm average annual precipitation lines are also visualized. (C) Diagram of the average shared SNPs between O.aries individuals as well as between O.aries individuals and the individuals of the three wild species (O.a.musimon,O.a.polii, and C.ibex). (D) Decay of linkage disequilibrium (LD) in the Chinese native sheep breeds, with one line per breed. (E) NJ tree constructed using p-distances between individuals. The three wild species are used as outgroups. (F) Principal components 1 and 2 for the 77 native sheep. (G) Population genetic structure of the 77 native sheep inferred from the program FRAPPE v1.1. The length of each colored segment represents the proportion of the individual genome inferred from ancestral populations (K¼2–4). See supplementary table S1, Supplementary Material online for the abbreviations of the breeds and individuals. Yang et al. .doi:10.1093/molbev/msw129 MBE 2578 nonsynonymous/synonymous ratio of 0.662. This ratio was within the range of ratios (0.40–1.42) for humans and other animals (e.g., mouflon, argali, ibex, and goat; supplementary table S9,Supplementary Material online). There were 15.98% (10.40–27.56%) of SNPs shared between individuals within breeds, 14.20% (8.86–26.41%) between breeds, 11.08% (9.13–13.05%) for O.aries–O.a.musimon, 7.87% (6.28–10.13%) for O.aries–O.a.polii,and1.60% (1.32–2.07%) for O.aries–O.a.musimon–O.a.polii (fig. 1C). For the three wild species, the shared proportions of SNPs were 6.93%, 2.65%, 4.28%, and 0.64% for O.a.musimon–O.a. polii,O.a.musimon–C.ibex,O.a.polii–C.ibex and O.a. musimon–O.a.polii-C.ibex, respectively. These statistics showed the proportions of SNPs with a given taxonomic range (i.e., Ovine species, O.aries and breeds) and thus of different origins (i.e., from predating the divergence of ovine species to predomestic, and then to postdating breed divergence). Moreover, O.aries–O.aries showed the largest average number of shared SNPs per pair of individuals, followed by O.aries–O.a.musimon,O.aries–O.a.polii,andO.aries–C.ibex (supplementary table S7,Supplementary Material online). The average numbers of private SNPs per pair of individuals ranged from 2.11 million for O.aries–O.a.musimon to 9.77 million for O.a.musimon–C.ibex (supplementary table S10, Supplementary Material online). In addition, we identified an average of 264,871 indels (<100 bp) per individual for the native sheep, 192,315 for O.a.musimon, 354,619 for O.a.polii and 755,867 for C.ibex (supplementary table S11,Supplementary Material online). For the two high-coverage samples (ZNQ24 and LOP41), we detected 177,371 copy number variations (CNVs, 200 bp–5 Mb; Freeman et al. 2006)(supplementary table S12,Supplementary Material online) and 156,134 structural variations (SVs, 100 bp–chromosome level; Alkan et al. 2011) which included insertions, deletions, inversions, intrachromosomal translocations and inter-chromosomal translocations (supplementary table S13,Supplementary Material online). Most of the indels, CNVs and SVs were located in intergenic regions. Patterns of Genomic Variation and Linkage Disequilibrium The genome-wide average h p valueforthenativesheepwas 2.28 10 3 (1.9–2.5 10 3 ), which was within the ranges of the parameters for other animals and humans (supplementary table S14,Supplementary Material online). The genomic variation parameters [average h p value, average number of regions of homozygosity (ROHs) and mean ROH size] for the three sheep groups showed congruent patterns. In other words, the Northern and Eastern Chinese breeds showed the greatest genomic diversity (h p ) and exhibited the fewest ROH number (879) and the smallest average ROH size (171. 34 kb), whereas the lowest genomic variations were observed in the breeds from the Yunnan–Kweichow Plateau (supplementary tables S15 and S16 and fig. S2B–G, Supplementary Material online). Inbred individuals were not observed in the 77 native sheep samples according to the IBS score (IBS <0.9). The breeds from the Qinghai– Tibetan Plateau (e.g., ZLZ and ZCD) and the Yunnan– Kweichow Plateau (e.g., WNS and SPS) showed an overall slow decay rate and a high level of LD, whereas the breeds from Northern and Eastern China exhibited a rapid decay rate and a low level of LD (fig. 1D). The differences in the genome-wide LD patterns between breeds likely reflected larger effective population sizes (N e ) for the Northern and Eastern Chinese breeds compared with the other breeds. Overall, the patterns of genomic diversity and LD strongly suggested a northern origin for Chinese sheep, from which Qinghai–Tibetan and Yunnan–Kweichow breeds were subsequently derived. This result is consistent with the migration routes of eastern Eurasian sheep, being from northern China to southern China, which were deduced in our recent mitogenomic study (Lv et al. 2015). Population Genetic Structure The neighbor-joining (NJ) tree revealed strong clustering of the native sheep into three genetic groups that presented a low level of genetic differentiation (mean F ST ¼0.014–0.040; supplementary table S17,Supplementary Material online): Qinghai–Tibetan breeds (n¼29,ZCD,ZNQ,ZRK,ZLZ, GDS, GZS, and MXS), Yunnan–Kweichow breeds (n¼20, TCS, WNS, SPS, and DQS) and Northern and Eastern Chinese breeds (n¼28, BRK, LOP, NMS, HUS, WDS, ALS, BYK, KAZ, HTS, and TSK; fig. 1E,supplementary fig. S3, Supplementary Material online). A principal component analysis (PCA) and a Bayesian model-based clustering analysis provided additional corroborating evidence for these groupings (fig. 1Fand G,supplementary figs. S4 and S5, Supplementary Material online). In the clustering analysis, when K¼2, the native sheep were genetically divided into the Yunnan–Kweichow breeds and the remaining breeds; when K¼4, the Northern and Eastern Chinese breeds, the Qinghai–Tibetan breeds and population ZLZ of the Tibetan sheep formed three additional genetic clusters (fig. 1G). We observed a clear signature of genetic admixture between the sheep breeds from the Qinghai– Tibetan Plateau and Northern and Eastern China as well as genetic introgression from Northern and Eastern Chinese breeds into the Yunnan–Kweichow breeds (fig. 1G,supplementary fig. S5,Supplementary Material online). Compared with the traditional view that Chinese native sheep were derived from three old breeds (i.e., Kazakh, Mongolian and Tibetan sheep; Du 2011), our results provided evidence that breeds from the Qinghai–Tibetan Plateau, Yunnan–Kweichow Plateau, and Northern and Eastern China are genetically distinct (fig. 1G). This finding is consistent with a recent study on 10 representative Chinese sheep breeds genotyped using the 50K SNP Beadchip (Wei et al. 2015). Demographic History and Population Admixture The two high-coverage samples (ZNQ24 and LOP41) exhibited concordant demographic trajectories, with two apparent expansions and two severe bottlenecks that mirrored the glacial cycles and sea level fluctuations that affected the areas (fig. 2A). The first bottleneck occurred 0.6 Mya Rapid Adaptation of Sheep to Extreme Environments .doi:10.1093/molbev/msw129 MBE 2579 (million years ago; fig. 2A) and is consistent with the Naynayxungla glaciation (0.78–0.50 Mya), which was the most extensive glaciation during the Quaternary Period (Zheng et al. 2002;Lehmkuhl and Owen 2005;Ehlers and Gibbard 2007). At that time, the sea level frequently fluctuated and was maintained at an overall low level (fig. 2A). After the retreat of the Naynayxungla glaciation, the ancestral sheep population expanded and reached a pinnacle (0.18 Mya; fig. 2A) during the Penultimate glaciation (0.30–0.13 Mya; Zheng et al. 2002;Lehmkuhl and Owen 2005;Ehlers and Gibbard 2007). The cold-climate interval and rising sea level at this stage could have contributed to a population expansion because an increase in grassland was likely under such environmental conditions (Lorenzen et al. 2011). The second bottleneck (70 kya; fig. 2A)wasdetected towards the end of the interglacial period (0.13–0.07 Mya; Zheng et al. 2002;Lehmkuhl and Owen 2005;Ehlers and Gibbard 2007), which presented environmental conditions similar to that of the present (Orlando et al. 2013). The ancestral sheep population size then peaked again at 15 kya during the Last Glacial Maximum (LGM) because the glaciations were less extensive and the grasslands had expanded (Zheng et al. 2002;Lehmkuhl and Owen 2005;Ehlers and Gibbard 2007;fig. 2A). Subsequently, the population size gradually decreased to a small number (fig. 2A), which was most likely the result of climate warming and associated grassland FIG.2.Demographic history of the sheep population. (A) PSMC analysis results for the representative individuals sequenced at a high read coverage (42) exhibit inferred variations in N e over the last 10 6 years. The N e of mouflon, argali, ibex, dog, horse, and wild boar are also rescaled. The glaciation periods, atmospheric surface air temperature (C) and global relative sea level data over the last 10 6 years are included. (B)@ a @ i analysis showing the demographic history of Chinese native sheep from 4,000 years ago to the present. Two population divergences occurred at 3,272 (t1) and 2,134 (t2) years ago. The average number of migrants per year between groups is shown between the black arrows. (C) Phylogenetic network of the inferred relationships among the 21 native breeds with three inter-group migration edges. The colored regions in the phylogenetic tree represent three inferred genetic groups. Arrows indicate migration events, and a spectrum of heat colors indicate the migration weights of the migration events. The scale bar shows 10the average standard error of the entries in the sample covariance matrix. NEC, Northern and Eastern Chinese breeds; QT, Qinghai–Tibetan breeds; YK, Yunnan–Kweichow breeds. For the abbreviations of the breeds, see supplementary table S1, Supplementary Material online. Yang et al. .doi:10.1093/molbev/msw129 MBE 2580 contractions after the LGM (Lorenzen et al. 2011). Despite relatively low estimates of N e , the pairwise sequentially Markovian coalescent (PSMC) results for the low-coverage samples (supplementary fig. S6,Supplementary Material online) showed demographic trends that were consistent with those inferred from the two high-coverage samples (fig. 2A). In general, our findings are consistent with previous results, indicating historical population expansion events for Chinese sheep as inferred from mtDNA sequence variability (Wang et al. 2007;Zhao et al. 2013). Compared with other animals (e.g., horse, Orlando et al. 2013; wild boar, Li et al. 2013;dog, Freedman et al. 2014; mouflon, argali, and ibex in this study), the demographic trend of sheep is most similar to horse (e.g., close timings of peaks and bottlenecks; fig. 2A)andmouflon (supplementary fig. S6,Supplementary Material online). These similar demographic trends might have been attributed to the similar habitat requirements for these herbivores. The differences in demographic shapes such as the extent of the changes in effective population sizes and timing of expansion/contraction were observed between sheep and the two wild Caprinae species (i.e., argali and ibex; supplemen tary fig. S6,Supplementary Material online). This may result from smaller N e and higher susceptibility to environmental changes (e.g., habitat fragmentation) of the wild herbivores compared with sheep and various environmental fluctuations in geographic regions where they inhabited. Among the five diffusion approximation for demographic inference (@a@i) models tested, model 3 achieved a maximum log-likelihood value of 16643 and was selected as the optimal model (supplementary table S18, Supplementary Material online). This model showed that Northern and Eastern Chinese breeds diverged from the remaining breeds 3,200 years ago, with an estimated N e of 24,220 for Northern and Eastern Chinese breeds and 498 for the other breeds (fig. 2B,supplementary table S19, Supplementary Material online). Approximately 2,100 years ago, breeds in areas other than Northern and Eastern China further diverged into the Qinghai–Tibetan breeds and Yunnan–Kweichow breeds, both of which have gradually increased and exhibited estimated N e values from 2,881 to 5,122 and 1,707 to 2,667, respectively (fig. 2B,supplementary table S19,Supplementary Material online). The current N e value for Northern and Eastern Chinese breeds is 4,873, which has obviously decreased since their initial appearance (24,220; fig. 2B, supplementary table S19,Supplementary Material online). An overall low to moderate level of gene flow (0.082–11. 791 migrants per year) was detected between the sheep groups (fig. 2B,supplementary table S19,Supplementary Material online). The maximum-likelihood (ML) tree without migration events inferred from the TreeMix analysis divided the 77 native sheep into three clusters (fig. 2C), which are similar to the population structuring patterns identified from the phylogenetic tree, PCA and genetic structure analysis (fig. 1E–G). When potential migration edges (i.e., migration events between the branches) were added to the ML tree, strong migration events were detected among the three clusters. Up to 98.77% of the variance between breeds was explained by a model with five migration events (supplementary fig. S7, Supplementary Material online). In the model, we observed three migration edges among clusters from Northern and Eastern China to the Yunnan–Kweichow Plateau, from the Yunnan–Kweichow Plateau to the Qinghai–Tibetan Plateau and from the Qinghai–Tibetan Plateau to Northern and Eastern China (fig. 2C). This result is consistent with the inter-group gene flow found in the @a@ianalysis and the genetic admixture between sheep groups observed in the genetic structure analysis. Genome-Wide Selective Sweep Test Using the top 5% of F ST values and h p ratio cutoffs (Z(F ST )>1.834, 1.859, 1.731, and 1.823 and log 2 (h p ratio) >0.352, 0.277, 0.247 and 0.198 for the Control group (HUS and WDS from Eastern China; see supplementary table S1,Supplementary Material online)/Tibetan group, Control group/Taklimakan Desert group, Low-altitude group/High-altitude group and Humid group/Arid group, respectively; figs. 3Aand 4A), we identified 731 candidate genes associated with plateau adaptations (>4,000 m; supplementary table S20,Supplementary Material online), 452 candidate genes associated with desert adaptations (average annual precipitation <10 mm; supplementary table S21,Supplementary Material online), 465 candidate genes involved in highaltitude adaptations (>1,500 m; supplementary table S22, Supplementary Material online) and 603 candidate genes involved in arid adaptations (average annual precipitation <400 mm; supplementary table S23,Supplementary Material online) for the native sheep. In addition, 261 (top 5% outliers, XP-EHH value >1.076) and 145 genes (top 5% outliers, XP-EHH value >0.897) were positively selected in XPEHH analysis, and 147 (top 1% outliers, log 10 (P)>0.433) and 95 genes (top 1% outliers, log 10 (P)>0.456) were positively selected in LFMM analysis regarding the plateau and desert environments, respectively. Overall, 112 and 75 genes were positively selected across all the four methods under the plateau and desert environments, respectively (figs. 5 and 6, supplementary tables S20 and S21,Supplementary Material online). The phylogenetic tree reconstructed based on the SNPs of the candidate genes of Tibetan sheep showed that most of the 77 native sheep (71/77, 92.21%) were correctly assigned into separated clades between the high-altitude (i.e., the Qinghai– Tibetan Plateau and Yunnan–Kweichow Plateau) and lowaltitude regions (supplementary fig. S8A,Supplementary Material online). Similarly, the phylogenetic tree that included the SNPs of the candidate genes of breeds from the Taklimakan Desert region largely divided the 77 sheep (74/77, 96.10%) into two clades between arid (i.e., average annual precipitation <200mm)andnon-aridregions(supplementary fig. S8B,Supplementary Material online). This topological pattern was different from that inferred from the genome-wide SNPs (fig. 1E), which supported the credibility of the candidate genes under selection. Regarding the large-effect SNP analyses of the two categories of candidate genes, we observed significantly (P<0.05) higher frequency differences between the contrasting groups than between the other Rapid Adaptation of Sheep to Extreme Environments .doi:10.1093/molbev/msw129 MBE 2581 pairwise groups (supplementary fig. S9,Supplementary Material online). Moreover, the Z(F ST )values(supplementary fig. S10,Supplementary Material online) and log 2 (h p ratio) values (supplementary fig. S11,Supplementary Material online) were much higher for the selected genomic regions of the sheep from the plateau region, desert region, high-altitude area and arid zone than the values at the whole-genome level. Compared with the overlap expected from chance, there was a significant excess of overlapping selective signals shared between the candidate genes identified here and the predefined gene panel (i.e., the previously published candidate genes in other mammalian species under similar extreme -4 -2 02 20 40 60 80 100 Cumulative (%) 0 20406080100 Cumulative (%) 2 4 6 8 10 Frequency (%) 0.2 0.4 0510 4 log2(θπ·control/θπ·Tibetan)Frequency (%) 0 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 -1.0 0.0 3.0 FST log2(θπ·control/θπ·Tibetan) Tajimas’D FST log2(θπ·control/θπ·Tibetan) Tibetan sheep Control sheep Ovis aries-M Ovis aries-W Bos taurus Tursiops truncatus Pan troglodytes Gorilla gorilla gorilla Homo sapiens Callithrix jacchus Microcebus murinus Pteropus vampyrus Loxodonta africana Cavia porcellus Oryctolagus cuniculu s Mus musculus 0.000.10 CCCTGAGGGAA CCCTGCGGGAA CCCTGCGGGAA CCCTGCGGGAA CCCTGCGGGAG CCCTGCGGGAG CCCTGCGGGAG CCCTGCGAGAG CCCTGCGGGAG CCTTGCGGGAA CCCTGCGGGAG CCCTGCGCGAG CTCTGCGCGAG CCCTGCGCGAG CTGATTGACTA CTGATCGACTA CTGATCGACTA CTGATCGACTA CTGATCGACTA CTGATCGACTA CTGATCGACTA CTGATCGACTA CTGATCGACTA CTGATCGACTA CTGATCGACTA CTGATTGACTA CTGATCGACTA CTGATTGACTA CTCTATACATC CTCTACACATC CTCTACACATC CTCTACACATC CTCTACACGTC CTCTACACGTC CTCTACACGTC CTCTACACGTC CTCTACACGTC CTCTACACATC CTCTACACGTC CTGTACACGTC CTCTACACATC CTGTATACATC SOCS2 Peyer's patch abomasum abomasum mucosa adrenal gland alveolar macrophage biceps brain brain stem caecum cardiac ventricle cerebellum cerebrum cervix colon corpus luteum duodenum epididymis heart hypothalamus kidney liver longissimus dorsi lung mammary gland mesenteric lymph node omentum ovarian follicles ovary pituitary gland placenta (including membranes) prescapular lymph node rectum renal cortex renal medulla rumen skin skin from back skin from side spleen testis thyroid gland tonsil uterus ventricle white adipose whole embryo                                             01.9 log10(FPKM) Roslin - female, adult Roslin - female, juvenile Roslin - female, adult Jiang et al. (2014) - Texel SOCS2                                                                                                                                                                                         2945 410 275 49 (6) 4 (0.3) 67 (19) 49 (13) Candidate Genes Tibetan Candidate Genes High Altitude Candidate Genes Predefined gene panel Z(FST) A B C D E F 2.0 1.0 SOCS2 FIG.3.Genomic regions with strong selective signals in Tibetan sheep. (A) Distribution of log 2 (h p ratios) and Z(F ST ) values calculated in 100-kb sliding windows with 50-kb increments between Tibetan group (including breeds ZNQ, ZCD, and ZRK from the plateau environment) and control group (including breeds HUS and WDS from East China). The data points in red (corresponding to the top 5% of the empirical log 2 (h p ratios) ratio distribution with values >0.35 and the top 5% of the empirical Z(F ST ) distribution with values >1.83) are genomic regions under selection in Tibetan sheep. (B) Comparison between the overlap of candidate genes and the overlap expected by chance. Numbers in the intersection regions are the observed overlapping genes among the candidate genes in Tibetan sheep, sheep breeds from the high-altitude regions and the predefined gene panel (i.e., the previously published candidate genes in other mammalian species under the high-altitude environment, including human, dog, wolf, yak, pig, and Tibetan antelope). Numbers in parentheses show the number of genes expected by chance. The total numbers of genes for the sheep and the gene panel involved in the test are 18,013 and 65,029, respectively. (C) log 2 (h p ratios) and F ST values around the SOCS2 locus. The black and red lines represent the log 2 (h p ratios) and F ST values, respectively. (D) Tajima’s Dvalues around the SOCS2 locus. The blue and purple lines represent the Tibetan sheep and control sheep, respectively. (E) Evolutionary analysis of the SOCS2 gene. The inter-species NJ tree is derived from the 12 vertebrate orthologous sequences, and the mutations are marked in red. O. aries M., namely O. arise mutant, represents Tibetan sheep in which the mutations were observed. O. aries W., namely O. aries wild, refers to other sheep breeds in which the nucleotides were conserved. (F) Gene expression of SOCS2 in different sheep tissues is based on four different experiments deposited in the EBI Gene Expression Atlas database. The FPKM (fragments per kilobase of transcript per million mapped reads) value is used to measure the expression level. Yang et al. .doi:10.1093/molbev/msw129 MBE 2582 environments; figs. 3Band 4B). Similarly, the number of overlapping candidate genes among the F ST ,h p ratio, XP-EHH, and LFMM analyses significantly exceeded those expected by chance (figs. 5Cand 6C). These results provided further support that the candidate genes were reliable. Adaptive Mechanisms in Plateau Environments Among the candidate genes for the plateau adaptations, five (IFNGR2,MAPK4,NOX4,SLC2A4,andPDK1) were located in the classical HIF-1 (hypoxia-induced factors) pathway, which plays a central role in regulating cellular responses to hypoxia -4 -2 0 2 4 6 8 -4 -2 0 2 4 3 5 8 10 13 0 0.2 0.4 20 40 60 80 100 20 40 60 80 100 Cumulative (%) Cumulative (%) Frequency (%) Frequency (%) Z(FST) log2(θπ·control/θπ·Taklimakan Desert) GPX3 Peyer's patch abomasum abomasum mucosa adrenal gland alveolar macrophage biceps brain brain stem caecum cardiac ventricle cerebellum cerebrum cervix colon corpus luteum duodenum epididymis heart hypothalamus kidney liver longissimus dorsi lung mammary gland mesenteric lymph node omentum ovarian follicles ovary pituitary gland placenta (including membranes) prescapular lymph node rectum renal cortex renal medulla rumen skin skin from back skin from side spleen testis thyroid gland tonsil uterus ventricle white adipose whole embryo Roslin - female, adult Roslin - female, juvenile Roslin - female, adult Jiang et al. (2014) - Texel 4.4 log10(FPKM) GPX3 Candidate Genes Taklimakan Desert Candidate Genes Arid and semi-arid Candidate Genes Bactrian camel 287 377 250 7 (3.6) 27 (5.3) A B -1.0 0.0 1.0 2.0 3.0 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 FST log2(θπ·control/θπ·Taklimakan Desert) Tajimas’D -1.0 0.0 1.0 2.0 3.0 4.0 Taklimakan Desert sheep Control sheep FST log2(θπ·control/θπ·Taklimakan Desert) GPX3 C D E F 0%50%100% 0% 50%100% 0% 50%100% 0% 50%100% 0% 50%100% 0% 50%100% HUS WDS TCS DQS SPS GZS ZNQ ZCD ZRK ZLZ GDS NMS LOP BRK A G G C C T C T T C C T 0 FIG.4.Genomic regions with strong selective signals in sheep breeds from the Taklimakan Desert region. (A) Distribution of log 2 (h p ratios) and Z(F ST ) values calculated in 100-kb sliding windows with 50-kb increments between Taklimakan Desert group (including breeds LOP and BRK from the desert environment) and control group (including breeds HUS and WDS from East China). The data points in red (corresponding to the top 5% of empirical log 2 (h p ratios) ratio distribution with values >0.28 and the top 5% of empirical Z(F ST ) distribution with values >1.86) are genomic regions under selection in the sheep breeds from the Taklimakan Desert region. (B) Comparison between the overlap of candidate genes and the overlap expected by chance. Numbers in the intersection regions are the observed overlapping genes among the candidate genes in sheep breeds from the Taklimakan Desert region, sheep breeds from the arid regions and the predefined gene panel (i.e., the previously published candidate genes in other mammalian species under the arid environment, including Bactrian camel). Numbers in parentheses show the number of genes expected by chance. The total numbers of genes for sheep and Bactrian camel involved in the test are 18,013 and 20,251, respectively. (C) log 2 (h p ratios) and F ST values around the GPX3 locus. The black and red lines represent log 2 (h p ratios) and F ST values, respectively. (D) Tajima’s Dvalues around the GPX3 locus. The blue and purple lines represent the Taklimakan Desert sheep and control sheep, respectively. (E) Allele frequencies of six SNPs within the GPX3 gene across Chinese native sheep breeds. The desert breeds include BRK and LOP, whereas the non-desert breeds comprise all other breeds. (F) Gene expression of GPX3 in different sheep tissues is based on four different experiments deposited in the EBI Gene Expression Atlas database. The FPKM (fragments per kilobase of transcript per million mapped reads) value is used to measure the expression level. For the abbreviations of the breeds, see supplementary table S1,Supplementary Material online. Rapid Adaptation of Sheep to Extreme Environments .doi:10.1093/molbev/msw129 MBE 2583 (Frede and Fandrey 2013); four were found in the corresponding downstream vascular endothelial growth factor (VEGF; RRAS and MAPK4) and glycolysis/gluconeogenesis (KIF2A and KHSRP) pathways; and five (ARHGEF19,STK17A, MAPK4,PPP1R1A and PRKG1)werefoundinthevascular smooth muscle contraction (VSMC) pathway, which adjusts the diameter of blood vessels and the delivery of blood oxygen (fig. 7A). We also found nine positively selected genes (PLIN2,TWIST1,KIF2A,PTPN9,SOCS2,NCOA3,STK17A, PDGFD,andLONP1) regulating specific genes in the HIF-1 pathway, two genes (NF1 and HAND2) affecting the VEGF signaling pathway, one gene (LONP1) influencing the glycolysis/gluconeogenesis pathway and two genes (NOSIP and SPSB4) regulating the nitric oxide (NO) compound in the VSMC pathway (fig. 7A). The network of these relevant pathways indicated that hypoxia-induced factors, angiogenesis, vasodilatation and glycolysis metabolism are the most important factors that allow sheep to manage extreme hypoxic environmental pressure. Furthermore, the functions of the proposed pathway genes are related to the regulation of the cardiovascular system and energy metabolism under hypoxic conditions in human and murine cells (e.g., Mora et al. 2003;Tang et al. 2003;Kim et al. 2006;Jian et al. 2010; Sawada et al. 2012;Pinti et al. 2015). For a detailed description of the relevant functions of the pathway candidate genes, please see supplementary information,Supplementary Material online. In addition, six candidate genes (FSTL1,PVR,EXT2,ALT4, SOX6,andHAND2) ranking within the top 20 Z(F ST )values and two additional candidate genes (PDGFD and BMPR2; supplementary table S20,Supplementary Material online) were functionally involved in bone, muscle, craniofacial, limb, skin, and embryonic development in vertebrates based on the NCBI annotations and Ensembl databases as well as previous functional studies in humans, mice and zebrafish (e. g., Norton et al. 2005;Kayserili et al. 2009;Geng et al. 2011; Long et al. 2015;supplementary information,Supplementary Material online). These findings indicated that these genes may have played a dominant role in regulating body morphology in Tibetan sheep, suggesting that selection on morphology is an important adaptive mechanism under the plateau environment, which is consistent with the small A B C FIG.5.Analysis of the signatures of positive selection in the genome of Tibetan sheep. Genomic landscape of the (A) XP-EHH values and (B) P-values in the LFMM analysis in the genome of Tibetan sheep. The genes visualized in (A) and (B) are the candidate genes from the signaling pathways of Tibetan sheep in fig. 7A.(C) Number of candidate genes identified in Tibetan sheep by the four methods listed in each of the Venn diagram components. Numbers in parentheses show the number of genes expected by chance. Yang et al. .doi:10.1093/molbev/msw129 MBE 2584 References Alkan C, Coe BP, Eichler EE. 2011. Genome structural variation discovery and genotyping. Nat Rev Genet. 12:363–376. Barrett JC, Fry B, Maller J, Daly MJ. 2005. Haploview: analysis and visualization of LD and haplotype maps. Bioinformatics 21:263–265. Breyer MD, Breyer RM. 2000. Prostaglandin receptors: their role in regulating renal function. Curr Opin Nephrol Hypertens. 9:23–29. Brown ME, Funk CC. 2008. Food security under climate change. Science 319:580–581. Carroll MA, Balazy M, Margiotta P, Huang DD, Falck JR, McGiff JC. 1996. Cytochrome P-450-dependent HETEs: profile of biological activity and stimulation by vasoactive peptides. Am J Physiol. 271:R863–R869. Cho YS, Hu L, Hou H, Lee H, Xu J, Kwon S, Oh S, Kim H-M, Jho S, Kim S, et al. 2013. The tiger genome and comparative analysis with lion and snow leopard genomes. Nat Commun. 4:2433. de Magalh~ aes JP, Costa J. 2009. A database of vertebrate longevity records and their relation to other life-history traits. JEvolBiol. 22:1770–1774. Du L-X. 2011. Animal genetic resources in China. Beijing: China Agriculture Press. Easterling DR, Meehl GA, Parmesan C, Changnon SA, Karl TR, Mearns LO. 2000. Climate extremes: observations, modeling, and impacts. Science 289:2068–2074. Ehlers J, Gibbard PL. 2007. The extent and chronology of Cenozoic global glaciation. Quatern Int. 164–165:6–20. Esworthy RS, Chu FF, Paxton RJ, Akman S, Doroshow JH. 1991. Characterization and partial amino acid sequence of human plasma glutathione peroxidase. Arch Biochem Biophys. 286:330–336. Frede S, Fandrey J. 2013. Cellular and molecular defenses against hypoxia. In: Swenson ER, B€ artsch P, editors. High altitude: human adaptation to hypoxia. New York: Springer Press. p. 26. Freedman AH, Gronau I, Schweizer RM, Ortega-Del Vecchyo D, Han E, Silva PM, Galaverni M, Fan Z, Marx P, Lorente-Galdos B, et al. 2014. Genome sequencing highlights the dynamic early history of dogs. PLoS Genet 10:e1004016. Freeman JL, Perry GH, Feuk L, Redon R, McCarroll SA, Altshuler DM, AburataniH,JonesKW,Tyler-SmithC,HurlesME.2006.Copynumber variation: new insights in genome diversity. Genome Res. 16:949–961. Frichot E, Schoville SD, Bouchard G, Franc¸ois O. 2013. Testing for associations between loci and environmental gradients using latent factor mixed models. Mol Biol Evol. 30:1687–1699. Gatalica Z, Lilleberg SL, Koul MS, Vanecek T, Hes O, Wang B, Michal M. 2008. COX-2 gene polymorphisms and protein expression in renomedullary interstitial cell tumors. Hum Pathol. 39:1495–1504. GeR-L,CaiQ,ShenY-Y,SanA,MaL,ZhangY,YiX,ChenY,YangL, Huang Y, et al. 2013. Draft genome sequence of the Tibetan antelope. Nat Commun. 4:1858. Geng Y, Dong Y, Yu M, Zhang L, Yan X, Sun J, Qiao L, Geng H, Nakajima M, Furuichi T, et al. 2011. Follistatin-like 1 (Fstl1) is a bone morphogenetic protein (BMP) 4 signaling antagonist in controlling mouse lung development. Proc Natl Acad Sci USA. 108:7058–7063. GouX,WangZ,LiN,QiuF,XuZ,YanD,YangS,JiaJ,KongX,WeiZ,etal. 2014. Whole-genome sequencing of six dog breeds from continuous altitudes reveals adaptation to high-altitude hypoxia. Genome Res. 24:1308–1315. Gutenkunst RN, Hernandez RD, Williamson SH, Bustamante CD. 2009. Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genet. 5:e1000695. Hancock AM, Clark VJ, Qian Y, Di Rienzo A. 2011. Population genetic analysis of the uncoupling proteins supports a role for UCP3 in human cold resistance. Mol Biol Evol. 28:610–614. Hoffmann I. 2010. Climate change in context: implications for livestock production and diversity. In: Odongo NE, Garcia M, Viljoen GJ, editors. Sustainable improvement of animal production and health. Vienna: IAEA-FAO. p. 33–44. Howden SM, Soussana JF, Tubiello FN, Chhetri N, Dunlop M, Meinke H. 2007. Adapting agriculture to climate change. Proc Natl Acad Sci USA. 104:19691–19696. Jian B, Wang D, Chen D, Voss J, Chaudry I, Raju R. 2010. Hypoxia-induced alteration of mitochondrial genes in cardiomyocytes—role of Bnip3 and Pdk1. Shock 34:169–175. Jiang J-C, He M-L, Karma R, Penpa T. 1991. A comparative study on several hematologic values of Tibet Plateau sheep at different altitudes. Chinese J Anim Sci. 27:11–14. JiangY,XieM,ChenW,TalbotR,MaddoxJF,FarautT,WuC,Muzny DM, Li Y, Zhang W, et al. 2014. The sheep genome illuminates biology of the rumen and lipid metabolism. Science 344:1168–1173. JirimutuWZ,DingG,ChenG,SunY,SunZ,ZhangH,WangL,HasiS, Zhang Y, Li J, et al. 2012. Genome sequences of wild and domestic Bactrian camels. Nat Commun. 3:1202. KayseriliH,UzE,NiessenC,VargelI,AlanayY,TuncbilekG,YigitG, UygunerO,CandanS,OkurH,etal.2009.ALX4 dysfunction disrupts craniofacial and epidermal development. Hum Mol Genet. 18:4357–4366. Kijas JW, Lenstra JA, Hayes B, Boitard S, Neto LRP, San Cristobal M, Servin B, McCulloch R, Whan V, Gietzen K, et al. 2012. Genome-wide analysis of the world’s sheep breeds reveals high levels of historic mixture and strong recent selection. PLoS Biol. 10:e1001258. Kim E-S, Elbeltagy AR, Aboul-Naga AM, Rischkowsky B, Sayre B, Mwacharo JM, Rothschild MR. 2016. Multiple genomic signatures of selection in goats and sheep indigenous to a hot arid environment. Heredity 116:255–264. Kim JW, Tchernyshyov I, Semenza GL, Dang CV. 2006. HIF-1-mediated expression of pyruvate dehydrogenase kinase: a metabolic switch required for cellular adaptation to hypoxia. Cell Metab. 3:177–185. Korneliussen TS, Albrechtsen A, Nielsen R. 2014. ANGSD: analysis of next generation sequencing data. BMC Bioinformatics 15:356. Legge T. 1996. The beginnings of caprine domestication. In: Harris DR, editor. The origins and spread of agriculture and pastoralism in Eurasia. New York: Smithsonian Institution Press. p. 238–262. Lehmkuhl F, Owen LA. 2005. Late Quaternary glaciation of Tibet and the bordering mountains: a review. Boreas 34:87–100. Li H, Durbin R. 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25:1754–1760. Li H, Durbin R. 2011. Inference of human population history from individual whole-genome sequences. Nature 475:493–496. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. 2009. The sequence alignment/map format and SAMtools. Bioinformatics 25:2078–2079. Li M, Tian S, Jin L, Zhou G, Li Y, Zhang Y, Wang T, Yeung CK, Chen L, Ma J, et al. 2013. Genomic analyses identify distinct patterns of selection in domesticated pigs and Tibetan wild boars. Nat Genet. 45:1431–1438. Librado P, Der Sarkissian C, Ermini L, Schubert M, J onsson H, Albrechtsen A, Fumagalli M, Yang MA, Gamba C, Sequin-Orlando A, et al. 2015. Tracking the origins of Yakutian horses and the genetic basis for their fast adaptation to subarctic environments. Proc Natl Acad Sci USA. 112:E6889–E6897. Liu F-Y, Chen Q-H. 2006. Effect of acute hypoxia on the parameters of blood gas in Tibetan Plateau sheep. Chinese J Zool. 41:48–52. LiuS,LorenzenED,FumagalliM,LiB,HarrisK,XiongZ,ZhouL, Korneliussen TS, Somel M, Babbitt C, et al. 2014. Population genomics reveal recent speciation and rapid evolutionary adaptation in polar bears. Cell 157:785–794. Lobell DB, Burke MB, Tebaldi C, Mastrandrea MD, Falcon WP, Naylor RL. 2008. Prioritizing climate change adaptation needs for food security in 2030. Science 319:607–610. LongL,OrmistonML,YangX,SouthwoodM,Gr € af S, Machado RD, Mueller M, Kinzel B, Yung LM, Wilkinson JM, et al. 2015. Selective enhancement of endothelial BMPR-II with BMP9 reverses pulmonary arterial hypertension. Nat Med. 21:777–785. Lorenzen ED, Nogue´s-Bravo D, Orlando L, Weinstock J, Binladen J, Marske KA, Ugan A, Borregaard MK, Gilbert MT, Nielsen R, et al. 2011. Species-specific responses of Late Quaternary megafauna to climate and humans. Nature 479:359–364. Lv F-H, Agha S, Kantanen J, Colli L, Stucki S, Kijas JW, Joost S, Li M-H, Ajmone Marsan P. 2014. Adaptations to climate-mediated selective pressures in sheep. Mol Biol Evol. 31:3324–3343. Rapid Adaptation of Sheep to Extreme Environments .doi:10.1093/molbev/msw129 MBE 2591 Lv F-H, Peng W-F, Yang J, Zhao Y-X, Li W-R, Liu M-J, Ma Y-H, Zhao Q-J, Yang G-L, Wang F, et al. 2015. Mitogenomic meta-analysis identifies two phases of migration in the history of eastern Eurasian sheep. Mol Biol Evol. 32:2515–2533. MaT,WangJ,ZhouG,YueZ,HuQ,ChenY,LiuB,QiuQ,WangZ, Zhang J, et al. 2013. Genomic insights into salt adaptation in a desert poplar. Nat Commun. 4:2797. Mackelprang R, Waldrop MP, DeAngelis KM, David MM, Chavarria KL, Blazewicz SJ, Rubin EM, Jansson JK. 2011. Metagenomic analysis of a permafrost microbial community reveals a rapid response to thaw. Nature 480:368–371. Maquire JJ, Johnson CM, Mockridge JW, Davenport AP. 1997. Endothelin converting enzyme (ECE) activity in human vascular smooth muscle. Br J Pharmacol. 122:1647–1654. Mi H, Muruganujan A, Casagrande JT, Thomas PD. 2013. Large-scale gene function analysis with the PANTHER classification system. Nat Protoc. 8:1551–1566. Mi H, Poudel S, Muruganujan A, Casagrande JT, Thomas PD. 2016. PANTHER version 10: expanded protein families and functions, and analysis tools. Nucleic Acids Res. 44:D336–D342. MoraA,DaviesAM,BertrandL,SharifI,BudasGR,Jovanovi  cS,Mouton V, Kahn CR, Lucocg JM, Gray GA, et al. 2003. Deficiency of PDK1 in cardiac muscle results in heart failure and increased sensitivity to hypoxia. Embo J. 22:4666–4676. Mysterud A, Stenseth NC, Yoccoz NG, Langvatn R, Steinheim G. 2001. Nonlinear effects of large-scale climatic variability on wild and domestic herbivores. Nature 410:1096–1099. Mysterud A, Yoccoz NG, Langvatn R, Pettorelli N, Stenseth NC. 2008. Hierarchical path analysis of deer responses to direct and indirect effects of climate in northern forest. Philos Trans R Soc Lond B Biol Sci. 363:2359–2368. Nielsen R, Korneliussen T, Albrechtsen A, Li Y, Wang J. 2012. SNP calling, genotype calling, and sample allele frequency estimation from newgeneration sequencing data. PLoS One 7:e37558. Norton WHJ, Ledin J, Grandel H, Neumann CJ. 2005. HSPG synthesis by zebrafish Ext2 and Extl3 is required for Fgf10 signalling during limb development. Development 132:4963–4973. Oakeshott JG, Wilson SR, Knibb WR. 1988. Selection affecting enzyme polymorphisms in enclosed Drosophila populations maintained in a natural environment. Proc Natl Acad Sci USA. 85:293–297. Oleksyk TK, Smith MW, O’Brien SJ. 2010. Genome-wide scans for footprints of natural selection. Philos Trans R Soc Lond B Biol Sci. 365:185–205. Orlando L, Ginolhac A, Zhang G, Froese D, Albrechtsen A, Stiller M, Schubert M, Cappellini E, Petersen B, Moltke I, et al. 2013. Recalibrating Equus evolution using the genome sequence of an early Middle Pleistocene horse. Nature 499:74–78. Park SD, Magee DA, McGettigan PA, Teasdale MD, Edwards CJ, Lohan AJ, Murphy A, Braud M, Donoghue MT, Liu Y, et al. 2015. Genome sequencing of the extinct Eurasian wild aurochs, Bos primigenius, illuminates the phylogeography and evolution of cattle. Genome Biol 16:234. Parker KL, Robbins CT. 1985. Thermoregulation in ungulates. In: Hudson RJ, White RG, editors. Bioenergetics of wild herbivores. Boca Raton, FL: CRC Press, Inc. p. 161–182. Patterson N, Price AL, Reich D. 2006. Population structure and eigenanalysis. PLoS Genet. 2:e190. Petryszak R, Burdett T, Fiorelli B, Fonseca NA, Gonzalez-Porta M, Hastings E, Huber W, Jupp S, Keays M, Kryvych N, et al. 2014. Expression Atlas update—a database of gene and transcript expression from microarrayand sequencing-based functional genomics experiments. Nucleic Acids Res. 42:D926–D932. Piao S, Ciais P, Huang Y, Shen Z, Peng S, Li J, Zhou L, Liu H, Ma Y, Ding Y, et al. 2010. The impacts of climate change on water resources and agriculture in China. Nature 467:43–51. Pickrell JK, Pritchard JK. 2012. Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genet. 8:e1002967. Pinti M, Gibellini L, Liu Y, Xu S, Lu B, Cossarizza A. 2015. Mitochondrial Lon protease at the crossroads of oxidative stress, aging and cancer. Cell Mol Life Sci. 72:4807–4824. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, et al. 2007. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 81:559–575. Qiu Q, Zhang G, Ma T, Qian W, Wang J, Ye Z, Cao C, Hu Q, Kim J, Larkin DM, et al. 2012. The yak genome and adaptation to life at high altitude. Nat Genet. 44:946–949. Qu Y, Zhao H, Han N, Zhou G, Song G, Gao B, Tian S, Zhang J, Zhang R, Meng X, et al. 2013. Ground tit genome reveals avian adaptation to living at high altitudes in the Tibetan plateau. Nat Commun. 4:2071. Rico-Bautista E, Flores-Morales A, Fern andez-Pe´rez L. 2006. Suppressor of cytokine signaling (SOCS) 2, a protein with multiple functions. Cytokine Growth Factor Rev. 17:431–439. Sabeti PC, Varilly P, Fry B, Lohmueller J, Hostetter E, Cotsapas C, Xie X, Byrne EH, McCarroll SA, Gaudet R, et al. 2007. Genome-wide detection and characterization of positive selection in human populations. Nature 449:913–918. SawadaJ,UrakamiT,LiF,UrakamiA,ZhuW,FukudaM,LiDY,Ruoslahti E, Komatsu M. 2012. Small GTPase R-Ras regulates integrity and functionality of tumor blood vessels. Cancer Cell 22:235–249. Scheinfeldt LB, Soi S, Thompson S, Ranciaro A, Woldemeskel D, Beggs W, Lambert C, Jarvis JP, Abate D, Belay G, et al. 2012. Genetic adaptation to high altitude in the Ethiopian highlands. Genome Biol. 13:R1. Sezgin E, Duvernell DD, Matzkin LM, Duan Y, Zhu CT, Verrelli BC, Eanes WF. 2004. Single locus latitudinal clines and their relationship to temperate adaptation in metabolic genes and derived alleles in Drosophila melanogaster. Genetics 168:923–931. SimonsonTS,YangY,HuffCD,YunH,QinG,WitherspoonDJ,BaiZ, Lorenzo FR, Xing J, Jorde LB, et al. 2010. Genetic evidence for highaltitude adaptation in Tibet. Science 329:72–75. Stopka T, Zivny JH, Stopkova P, Prchal JF, Prchal JT. 1998. Human hematopoietic progenitors express erythropoietin. Blood 91:3766–3772. Tajima F. 1989. Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123:585–595. Tang H, Peng J, Wang P, Risch NJ. 2005. Estimation of individual admixture: analytical and study design considerations. Genet Epidemiol. 28:289–301. TangM,WangGR,LuP,KarasRH,AronovitzM,HeximerSP, Kaltenbronn KM, Blumer KJ, Siderovski DP, Zhu Y, et al. 2003. Regulator of G-protein signaling-2 mediates vascular smooth muscle relaxation and blood pressure. Nat Med. 9:1506–1512. Turnley AM. 2005. Role of SOCS2 in growth hormone actions. Trends Endocrinol Metab. 16:53–58. WangM-S,LiY,PengM-S,ZhongL,WangZ-J,LiQ-Y,TuX-L,DongY, Zhu C-L, Wang L, et al. 2015. Genomic analyses reveal potential independent adaptation to high altitude in Tibetan chickens. Mol Biol Evol. 32:1880–1889. Wang X, Ma YH, Chen H, Guan WJ. 2007. Genetic and phylogenetic studies of Chinese native sheep breeds (Ovis aries) based on mtDNA D-loop sequences. Small Ruminant Res. 72:232–236. WeiC,WangH,LiuG,WuM,CaoJ,LiuZ,LiuR,ZhaoF,ZhangL,LuJ, et al. 2015. Genome-wide analysis reveals population structure and selection in Chinese indigenous sheep breeds. BMC Genomics 16:194. Weir BS, Cockerham CC. 1984. Estimating F-statistics for the analysis of population structure. Evolution 38:1358–1370. WuH,GuangX,Al-FageehMB,CaoJ,PanS,ZhouH,ZhangL, Abutarboush MH, Xing Y, Xie Z, et al. 2014. Camelid genomes reveal evolution and adaptation to desert environments. Nat Commun. 5:5188. ZhaoZ,XuD,WangL,HaoJ,WangJ,ZhouX,WangW,QiuQ,HuangX, Zhou J, et al. 2016. Convergent evolution of rumen microbiomes in high-altitude mammals. Curr Biol. 26(14):1873–1879. Zhao E, Yu Q, Zhang N, Kong D, Zhao Y. 2013. Mitochondrial DNA diversity and the origin of Chinese indigenous sheep. Trop Anim Health Pro. 45:1715–1722. Zheng B, Xu Q, Shen Y. 2002. The relationship between climate change and Quaternary glacial cycles on the Qinghai-Tibetan Plateau: review and speculation. Quatern Int. 97–98:93–101. Yang et al. .doi:10.1093/molbev/msw129 MBE 2592