Evolutionary history of Tibetans inferred from whole-genome sequencing

The indigenous people of the Tibetan Plateau have been the subject of much recent interest because of their unique genetic adaptations to high altitude. Recent studies have demonstrated that the Tibetan EPAS1 haplotype is involved in high altitude-adaptation and originated in an archaic Denisovan-related population. We sequenced the whole-genomes of 27 Tibetans and conducted analyses to infer a detailed history of demography and natural selection of this population. We detected evidence of population structure between the ancestral Han and Tibetan subpopulations as early as 44 to 58 thousand years ago, but with high rates of gene flow until approximately 9 thousand years ago. The CMS test ranked EPAS1 and EGLN1 as the top two positive selection candidates, and in addition identified PTGIS, VDR, and KCTD12 as new candidate genes. The advantageous Tibetan EPAS1 haplotype shared many variants with the Denisovan genome, with an ancient gene tree divergence between the Tibetan and Denisovan haplotypes of about 1 million years ago. With the exception of EPAS1, we observed no evidence of positive selection on Denisovan-like haplotypes.


Introduction
Adaptation to high altitude is among the most notable examples of natural selection in our species. Tibetans exhibit an exceptional capacity to compensate for decreased oxygen availability (hypoxia), UV exposure, cold, and limited food sources as exemplified by continuous habitation at greater than 4000 meters above sea level on the Qinghai-Tibetan Plateau over several millennia [1][2][3][4].
Many candidate genes have been highlighted as likely contributors to Tibetans' trans-generational success at high altitude [5][6][7][8][9]. These findings are largely based upon the overlap of a priori genes, specifically those involved in hypoxia tolerance, and genomic footprints of adaptation identified in patterns within single nucleotide polymorphism (SNP) microarrays [10][11][12] or exome sequencing data [13]. In a few cases, putatively adaptive tagging variants in hypoxia inducible factor (HIF) pathway genes (EPAS1 [10,13], EGLN1 and PPARA [12]) are associated with relatively lower hemoglobin concentration exhibited by many Tibetans at altitude. With the exception of the EGLN1 locus [14,15], however, precise functional variants within these genomic regions are unknown.
Recent advances in whole-genome sequencing (WGS) technology provide progress toward identifying the functional variants that underlie high-altitude adaptation. Rather than focusing on a relatively small portion of the genome using candidate gene approaches, WGS selection analyses often have sufficient power to comprehensively interrogate the genome in an unbiased manner. Moreover, WGS offers the opportunity to perform exhaustive searches for adaptive variation, yielding complete surveys of non-protein coding regulatory, conserved, and structural variants possibly tagged or missed in microarray analyses. WGS data are also important from an evolutionary standpoint, providing unbiased insights into long-standing questions regarding adaptive processes, including the role of adaptive introgression, as recently revealed through the targeted sequence comparisons between Tibetans and the Denisovan genome at the EPAS1 gene locus [16].
To identify regions contributing to high-altitude adaptation and regions with archaic introgressions with high resolution, we performed a comprehensive genomic analysis of WGS data in 27 Tibetans [17]. We estimated the demographic history of Tibetans using MSMC and @a@i, and performed the Composite of Multiple Signals test for fine-mapping loci and variants with recent selective sweep. We also estimated genome-wide levels of archaic admixture, and generated a fine-scale map of Denisovan-like introgression in Tibetans.

Global genomic analysis and demographic history estimate
We conducted principal components analysis (PCA) on all SNVs with observed non-reference allele frequency greater than 5% among the combined samples of Tibetans and five 1000 Genomes Project [18] populations (Peruvian, Punjabi, Chinese, Yoruba and European) (Supplementary Methods). Tibetans and Han Chinese were closely clustered in the first 4 principal components (S1A and S1B Fig) similar to our previous observation [12], but were clearly differentiated along 5 th principal component (S1C Fig). Similarly, ADMIXTURE analysis [19] (Supplementary Methods) on the same six populations shows that when the number of ancestral populations (K) was set to less than 6, Tibetans and Han Chinese were genetically similar. However, when K was set to 6 or greater, these two populations exhibited distinct genetic profiles (S2 Fig). The estimated proportion of Tibetan genetic ancestry was greater than 99% for 19 out of the 27 Tibetan individuals. For the remaining Tibetan individuals, the average estimated proportion of Tibetan genetic ancestry was 83.1%. For the demographic history analysis, we excluded Tibetan individuals with less than 99% Tibetan genetic ancestry. F ST between Tibetans and Chinese based on our SNV data was 0.0148 using the Hudson estimator [20] and 0.0149 using the Weir and Cockerham estimator [21], similar to our previous analyses based on SNP microarray data [22], which is slightly higher than the F ST between Chinese and Japanese populations but lower than the F ST between pairs of HapMap populations from different continents [23].
We next annotated the functional impact of SNVs that are frequent in Tibetans but uncommon in other populations. We selected all germline SNVs with an alternative allele frequency above 10% in Tibetans but below 1% among Yoruba, Han Chinese, and Europeans from the 1000 Genomes Project. In total, we identified 10,702 such SNVs, among which there were 65 nonsynonymous variants that occurred within 58 genes (S1 Table). We predicted the impact of these protein-altering SNVs on genome functions using PolyPhen-2 and the Conservationcontrolled Amino Acid Substitution Matrix (CASM) in VAAST 2.0 [24]. Of the 65 nonsynonymous variants, Polyphen-2 predicted 12 to be functional and CASM predicted a distinct set of 3 variants to be functional. One of the variants, Chr2:233244223 A->C (p.T104P; all genomic coordinates in this manuscript are relative to GRCh37 reference sequence), is in the gene ALPP with an observed frequency of 0.14% in dbSNP and 14.81% in Tibetans. ALPP encodes an alkaline phosphatase, which is mainly expressed in placental and endometrial tissues. Rare nonsynonymous variants in ALPP have previously been associated with decreased risk of spontaneous abortion and in vitro fertilization failure [25].
We inferred the demographic history of Tibetans using the multiple sequential Markovian Coalescent (MSMC) [26] (Fig 1) and @a@i [27] (Fig 2; Table 1) (Supplementary Methods). MSMC estimated that the relative cross coalescence rate (a measurement of the amount of gene flow between two populations) between Han and Tibetans fell below 80% around 7 kya (bootstrap 95% CI: 3kya to 10kya). We then derived a more detailed demographic model using @a@i, which predicted that the initial divergence between Han and Tibetan was much earlier, at 54kya (bootstrap 95% C.I 44 kya to 58 kya). However, for the first 45ky, the two populations maintained substantial gene flow (6.8x10 -4 and 9.0x10 -4 per generation per chromosome). After 9.4 kya (bootstrap 95% C.I 8.6 kya to 11.2 kya), the gene flow rate dramatically dropped (1.3x10 -11 and 4x10 -7 per generation per chromosome), which is consistent with the estimate from MSMC. F ST predicted by the best-fitting demographic model in @a@i produced estimates consistent with observed F ST (0.0147 using Hudson's F ST estimator and 0.0148 using the Weir and Cockerham's estimator).

Whole-genome positive selection scans in the Tibetan genome
We used the Composite of Multiple Signals (CMS) test [28] to map regions under positive selection in the Tibetan genome (Fig 3). At each genomic variant, CMS evaluates the likelihood of five test statistics (iHS, XP-EHH, ΔiHH, F ST and ΔDAF) under two models: a null  model that assumes neutral evolution and an alternative model that incorporates a variety of scenarios involving recent strong positive selection (see S1 Supplementary Methods). CMS then combines information from the five test statistics in a Bayesian framework to create an aggregated CMS score. CMS and its components (iHS, F ST and XP-EHH) have been successfully applied to identify adaptive sweeps originating both from the same population [12,28,29] and from archaic introgression [30][31][32].
We first calculated the CMS score at each SNV using the 27 Tibetan genomes, with 62 Han Chinese genomes from the 1000 Genomes (1KG) Project [18,33] included as a comparison population (Supplementary Methods). In order to evaluate the statistical significance of our top candidates, we performed coalescence simulations with cosi [34] under a neutral model, using the demographic model that we estimated with @a@i [27] (Supplementary methods). Comparing the CMS scores from our observed data and the null distribution from coalescence simulations, we identified 377 candidate SNVs with FDR < 0.3; 100 of these SNVs are within 350kb of EGLN1 and 199 are within 350kb of EPAS1. We listed all the SNVs with FDR<0.3 in S2 Table. Because CMS scores between nearby SNVs are often correlated due to linkage disequilibrium, we also performed a region-based CMS test by dividing the genome into 200kb consecutive regions [12], with the test statistic equal to the highest CMS score within each region. Table 2 lists the top 10 regions with the highest CMS scores. The top regions are in the proximity of EPAS1 and EGLN1 genes, which were previously reported to be responsible for high-altitude adaptation among Tibetans [12,13]. One of the regions overlaps the gene VDR, which encodes the vitamin D3 receptor. Mutations in this gene can cause an abnormality in vitamin D metabolism, which may in turn lead to rickets [35] and preterm birth [36]. We examined the evidence of positive selection on the three SNVs with the best-CMS scores in this 200kb region. Based on the Complete Genomics (CG) 1000 Genome panel data [18,33], the alternative alleles of top 2 SNVs (Chr12:48363253 G->C and Chr12:48359527) are infrequent in Han Chinese (MAFs: 1.6% and 3.2%), but the MAFs were comparable between Tibetans (31.5% and 31.5%) and Europeans (32.1% and 31.4%). The third ranking SNV (Chr12: 48337328 C->T) has a minor allele frequency of 2.4% in Han Chinese, 0.0% in Europeans, 17.0% in Yoruba and 25.9% in Tibetans. The iHS value at this location is 4.22 standard deviations higher than the genome-wide mean, which is about 13,000 times more likely under the alternative model than the null. Fig 4a displays the haplotype map in Tibetans surrounding this region. The focal SNV is 394 bp upstream of the first exon of the transcript AK309587, overlapping the binding site of three transcription factors: TCF7L2, MAFK and MAFF [37]. Both a DNaseI hypersensitivity assay and H3K27me3 histone marker signature in the K562 cell line indicates this is a potential noncoding regulatory region. Another top 10 ranking region overlapped the KCTD12 gene (potassium channel tetramerization domain containing 12), a component of the G-protein-coupled receptors for γ-aminobutyric acid (GABA) receptor [38]. This region contains 3 SNVs with FDR < 0.05 from our single-variant CMS test. The strongest signal (Chr13: 77399462) among these three SNVs has   an observed minor allele frequency of 18.5% in Tibetans, 0.06% in dbSNP (build 142) and 0% among Yoruba, Europeans, and Han Chinese in the CG 1000 Genomes panel. The normalized integrated haplotype score (iHS) value at this location was 4.90 standard deviations higher than the genome-wide mean, which is~48,000 times more likely to occur under the alternative model than under the null model. This indicates that the haplotype homozygosity around the derived allele is much higher than around the ancestral allele, a scenario consistent with strong positive selection. The three SNVs are located in genomic regions with potential regulatory function. Chr13: 77399462 T->C and Chr13:77398842 G->A overlap with histone mark H3K27Ac annotation and are within enhancer regions, while Chr13:77405193 G->A is at a DNase I hypersensitive site and has a strong signature of H3K27me3 in the K562 (myelogenous leukemia) cell line [37,39]. We next examined Gene Ontology (GO) [40] annotations associated with genes in the top 0.2% of 200kb windows. We evaluated the statistical significance for over-representation of any GO term by repeated sampling without replacement from all 200kb regions in our dataset ( Table 3). The following GO terms related to hypoxia responses were over-represented: GO:0071456 (cellular response to hypoxia), GO:0036294 (cellular response to decreased oxygen levels), and GO:0071453(cellular response to oxygen levels). In addition to two genomic regions surrounding EPAS1 and EGLN1, we identified one additional gene related to hypoxia response in GO:0071456: prostaglandin I2 synthase (PTGIS). PTGIS converts prostaglandin H2 to prostacyclin, an effective vasodilator and inhibitor of platelet activation [41]. The SNV with the highest CMS score in this region (Chr20: 48175598 C->T) has a normalized XP-EHH value of 3.67, which was 51 times more likely under the alternative model than the null. The derived allele frequency was 85.2% in Tibetan, 42.7% in Han Chinese, 32.1% in Europeans, and 8.9% in Yoruba. The SNV is in the first intron of PTGIS, overlapping DNase I hypersensitivity regions, potential ZNF217 transcription factor binding sites and H3K36me3 histone mark signatures in the NT2-D1 (pluripotent embryonic carcinoma) cell line [37].
Because of the lower variant-calling qualities at genomic insertions and deletions (indels), we could not obtain reliable results from haplotype-based tests (XP-EHH, iHS and ΔiHH) on indels. Instead, we used the Population Branch Statistic (PBS) [13] to search for signatures of positive selection. To calculate the significance levels of the PBS scores, we estimated the null distribution on the same simulation dataset as in the CMS test. We found that indels with the best PBS signals tended to fall into regions having the highest CMS scores. Specifically, 6 out of the top 10 indels are also within the top 10 200kb regions with the highest CMS scores, all in the EPAS1 and EGLN1 regions (S3 Table;

Mapping of variants under positive selection in the EPAS1 and EGLN1 regions
We fine-mapped SNVs under positive selection in the EPAS1 region by combining the CMS scores with the population allele frequencies in major continental population groups (Europeans, Han, Native Americans and Yoruba) (Fig 4b). We discovered 199 SNVs with FDR < 0.3 from the CMS test (CMS q<0. 3 ) in the EPAS1 region (S4 Table), which spans 321 kb from the first intron of EPAS1 to 257 kb downstream of the last exon of EPAS1. Of these 199 SNVs, 64 were in the intronic region and the remaining 135 were in the 3' downstream region. These SNVs are enriched for alleles present in the Denisovan genome and rare in non-Tibetan populations (MAF<0.05 among Yoruba, Europeans, Native Americans and Asians). Specifically, 16.6% of EPAS1 CMS q<0. 3 SNVs meet this criterion, as compared to 0.3% genome-wide (p<2x10 -16 ).
A previous study identified a 32.7kb region (Chr2: 46567916-46600661) in EPAS1 that was highly differentiated between Han and Tibetans and also has a high proportion of Denisovanlike variants [16]. In our CMS test result, 2 out of the top 3 SNVs with the highest CMS score in the EPAS1 region (46597756 and 46598025 on Chromosome 2) were in this 32.7kb region. After excluding rare variants (SNVs with MAF <0.05 in both Han and Tibetans), the average F ST between Han and Tibetans in this 32.7kb region was 0.28. However, the haplotype extends far beyond the previously reported 32.7 kb region and contains genomic regions at least 150 kb downstream of EPAS1. Fig 4b illustrates the Tibetan haplotypes in EPAS1 and its 3' region. The linkage disequilibrium between SNVs (both Denisovan-like and not) was high among Tibetans (S5 Table), suggesting a slow decay of haplotype homozygosity. The small indel with the highest PBS score was a 4bp deletion located in the 2nd intron of EPAS1 (Chr2:46577800-46577803). Three out of the four nucleotides affected by this deletion are highly conserved in mammals. This indel overlaps an activating H3K27Ac mark in seven cell lines and the binding sites of three transcription factors: the polymerase subunit POLR2A, the HIF co-activator EP300 [42], and GATA2, which is key in the control of erythroid differentiation [43]. The observed frequency of the derived allele (deletion) was 62.9% in Tibetans, 0.8% in Han Chinese and 0% in Europeans, corresponding to a PBS score of 0.600, which was ranked 3 rd among all the genomic variant (indels and SNVs) in the EPAS1 gene.
In addition to the 32.7 kb block, we also identified a second 39.0kb candidate block downstream of EPAS1 coding region on the haplotype (Chr2: 46675505-46714553); this block was highly differentiated between Han and Tibetans and enriched for Denisovan variants (S3 Fig). Within this block, 88.2% (60 out of 68) of common Tibetan SNVs (derived allele frequency >50%) are present in the Denisovan genome, as compared to 81.4% (48 out of 59) in the previously reported 32.7kb region. The proportion of SNVs shared with the Neanderthal genome is 26.5% (18 out of 68) for the new 39.0kb block as compared to 25.4% (15 out of 59) for the previous 32.7kb block. In this block, the average F ST between Han and Tibetans for SNVs common in at least one of the two populations is 0.304 (as compared to 0.279 for the previously reported region). Intriguingly, this region also fully contains a 3.4kb Tibetan-enriched deletion region that showed footprints of positive selection in a previous study [44] (Fig 4b). Among our Tibetan samples, this deletion has an allele frequency of 59.2% and is in high LD (r 2 = 0.86) with the SNV with the highest CMS in the EPAS1 region (Chr2: 46597756; S6 Table).
We estimated the origin age of the selective sweep on the adaptive EPAS1 haplotype to be 12 kya (95% CI: 7 kya to 28 kya) using a maximum likelihood approach based on coalescence simulations of the demographic model described in Table 1 (see S1 Supplementary Methods for details). Our estimate is consistent with a previous estimate of 12 kya to 15 kya [44]. Using this point estimate of selection start time and assuming that the adaptive haplotype indeed originated from Denisovans or a closely related population, we further inferred that the introgression of this haplotype most likely occurred between 32 and 12 kya. We also estimated the gene-tree divergence time between the adaptive haplotype in Tibetans and the available Denisovan DNA sequence to be 997 kya (95% CI: [868, 1,139]). From the gene-tree divergence and calibrated models of Denisovan demographic history [45], we estimated a population divergence time of 868 kya (95%CI: 952 kya to 238 kya) between the Densivan population and the archaic population that contributed the adaptive EPAS1 haplotype.
Previously, we identified two non-synonymous variants in the EGLN1 gene (c.12C>G, p. D4E; c.380G>C, p.C127S) which together act as a co-adapted gene complex contributing to high-altitude adaptation in Tibetans [14]. In the current dataset, the c. in Han Chinese, consistent with previous findings [14]. The CMS score is -60.43 for c.12C>G and -72.08 for c.380G>C, both ranked among the top 0.01% variants genome-wide. We also identified a novel rare EGLN1 variant (c.C358T, p.P120L) within one Tibetan individual, who also carries the c.380G>C but not c.12C>G variant. In our dataset this variant has a frequency of 1.85% (95% CI: [0.05%, 9.89%]) in Tibetans, 0.81% (95% CI: [0.02%, 4.41%]) in Han Chinese and is absent from other populations in the 1KG dataset.

Mapping regions with high Denisovan ancestry in the Tibetan genome
To evaluate whether any genomic regions with Denisovan introgression other than EPAS1 may contribute to Tibetan high-altitude adaptation, we first calculated D-statistics [46] to evaluate the overall level of Denisovan admixture in the Tibetan genome. When comparing Tibetans to Yoruba, we observed a D of 0.011 (95% CI: 0.0068 to 0.0163), indicating a significant excess of Denisovan admixture among Tibetans. However, we observed no significant difference in D between Tibetans and Han Chinese (D = -0.0007, 95% CI -0.0027 to 0.0021), demonstrating that the genome-wide level of Denisovan admixture does not differ substantially between Tibetans and Han Chinese, and is consistent with a previous report [47]. Using the Q statistic of Rogers and Bohlender [48], we estimated that the proportion of Denisovan admixture in Tibetans is 0.4% (95% CI: 0.2% to 0.6%, see Methods and S4 Fig). Briefly, Q is an f 3 ratio estimator of admixture [48]. Given a sample from three populations X, Y and N, so that X and Y are more recently diverged from each other than either is from N, the calculation of Q uses a ratio of expectations to calculate the total admixture from N into Y, in a similar way as other f 3 -ratio estimators [46]. Q has the advantage of requiring only a single archaic individual, and uses external information about ghost admixture to correct bias due to admixture from an unsampled archaic population (Neanderthal).
To investigate the variation in Denisovan admixture across the genome among Tibetans, we first used an LD-based statistic S Ã [32] to identify introgressed regions originated from Denisovans. S Ã searches for haplotypes that share a substantial amount of similarity with the archaic genome and are present in the population of interest but not in the comparison population (see S1 Supplementary Methods). In our study, we used Han Chinese as the comparison population and Denisovan alleles as the archaic genome to search for Denisovan haplotypes absent among the Han Chinese genomes and present in at least one Tibetan genome (Supplementary Methods). In total, S Ã identified 660 50kb candidate regions using this procedure (S5 Fig; S7 Table). As an alternative approach, we also developed a new statistic, D Ã , which normalizes the D-statistics to enable the identification of admixed regions with calibrated Type I error (Methods). We used D Ã to compare Tibetans to Han Chinese for all 200kb genomic regions (S5 Fig). With this procedure, we identified 11 regions with significant proportions of Denisovan admixture after multiple-testing correction (FDR < 0.30). The region containing EPAS1 ranked second, with D Ã = 6.2 (p = 6.1 x 10 −6 ). Six of 11 regions identified by D Ã overlapped with the 660 regions identified by S Ã , including the region containing EPAS1 (Table 4).
Next, to test for evidence of adaptive Denisovan introgression, we compared the location of the Denisovan DNA regions identified by either S Ã or D Ã to the regions with the highest CMS scores. The top candidate identified by this analysis is EPAS1. Since we were interested in identifying novel adaptive Denisovan introgression, we removed the EPAS1 region and evaluated the statistical significance of the observed number of overlaps. To calculate the p-value, we permuted the CMS scores of all 200kb genomic regions to generate the distribution of overlaps under the null hypothesis that no additional Denisovan-like DNA was positively selected in Tibetans. When we compared the Denisovan-like genomic regions to the top 0.2% of CMS regions, we found no overlapping region with D Ã and only one with S Ã (Chr2:42200001-42400000). In the top 1% of CMS regions, we identified no overlapping regions with D Ã and three overlapping regions (Chr2: 42200001-42400000; Chr1: 79200001-79400000; Chr1:1200001-1400000) with S Ã ; however, the p-values from all tests were non-significant. Therefore, we find no evidence that additional Denisovan-like DNA has contributed to Tibetan high-altitude adaptation. We summarized the CMS scores of Denisovan variants that are present in Tibetans but uncommon among major population groups in Fig 3b.

Discussion
Our demographic history estimates from @a@i suggest that population structure existed between the Tibetan and Han ancestral subpopulations as early as 44 to 58 kya, with high rates of admixture maintained between the two subpopulations until around 9 kya. This agrees well with the findings of Lu et al. [49] that the early colonization of the Tibetan plateau occurred between 62kya and 38kya, and the post-Last Glacial Maximum (LGM) arrival at the Tibetan plateau could be between 15kya and 9kya. Moreover, archaeological evidence [1,50] suggests that the Tibetan plateau was occupied during the Late Pleistocene, roughly coinciding with the time of Han-Tibetan separation in our demographic model. Interestingly, MSMC did not detect a decrease in relative cross coalescence rate between Han and Tibetans until 3 to 9 kya. We hypothesized that this discrepancy was due to an inability of MSMC to differentiate between panmixia and high rates of gene flow between Tibetan and Han Chinese populations between 9kya and 54kya. To test this hypothesis, we conducted coalescent simulations to evaluate the ability of MSMC to accurately estimate the Tibetan-Han divergence time under a model where the two populations diverged at 54 kya but experienced a large amount of geneflow (6.8x10 -4 and 9.0x10 -4 per generation per chromosome) until 9 kya. As shown in S6 Fig,  MSMC was indeed unable to detect the ancestral population split at 54 kya. The observed pattern of estimated relative cross-coalescence rate was consistent between the observed data and simulations based on the best-fitting @a@i model, with no decrease in the relative cross-coalescence rate prior to 10 kya. Conversely, when we conducted coalescence simulation using the simplified demographic model predicted by MSMC (i.e., Han and Tibetan separated 7,000 years ago) and ran @a@i, we found that @a@i accurately recovered the simulated divergence date (S8 Table), providing support for the @a@i model with an earlier divergence date but a high migration rate between the two populations until approximately 9 kya. Our admixture analysis in the EPAS1 region confirmed findings from previous studies that a common Tibetan haplotype in this region contains excessive Denisovan-like DNA [16], but also showed that on a genome-wide level, the amount of Denisovan admixture in Tibetans (0.4%) is similar to that of Han Chinese. This suggests that at least one major introgression occurred from the archaic population to the common ancestor of Han and Tibetans. The high-altitude adaptive haplotype in the EPAS1 region may have been acquired prior to Han-Tibetan divergence, but was either lost or present at low frequencies in Han due to the lack of selection. This explanation agreed with the hypothesis proposed by Huerta-Sanchez et al [16], who posit that the EPAS1 haplotype may have been introduced prior to the separation of Han and Tibetans, considering the presence of the Tibetan EPAS1 haplotype in a single Han Chinese individual. Alternatively, it is possible that the low frequency of the Tibetan EPAS1 haplotype in China is the result of relatively recent migration from Tibet [16] and that the haplotype was originally introduced into Tibet after the Han-Tibetan divergence. This alternative explanation is supported by our 95% C.I. for the introduction of the EPAS1 haplotype into the ancestral Tibetan population (32 to 12 kya) which is later than the previous estimates for the date of Denisovan admixture into modern humans (44.0 to 54.0 kya) [51]. Other than the EPAS1 haplotype, we observed no evidence of other Denisovan-like DNA segments contributing to high-altitude adaptation in Tibetans, given that we did not detect significant overlaps between introgressed regions and the top 1% of CMS regions. However, we cannot rule out the possibility that with larger sample sizes, subtle adaptive introgression signals may be detected.
In addition to the previously reported 32.7kb haplotype block in EPAS1 [16], we identified a much larger haplotype that contains both the EPAS1 genic region and at least 150kb 3' of EPAS1. However, the extensive LD within the EPAS1 and its 3' region prevented us from pinpointing the location of the adaptive mutation. Interestingly, the selected haplotype (Fig 4b; S3  Fig) contains many Denisovan-like and non-Denisovan-like variants, and many top candidate SNVs in our analysis, as well as the previously reported 3.4kb Tibetan enriched deletion, were absent from the high-coverage Denisovan genome. This observation is reflected in the estimated population divergence time of 952 to 238 kya between the population represented by the Denisovan reference genome and the archaic population that admixed with Tibetans (Supplementary Methods). This estimate is broadly consistent with previous evidence that the archaic Denisovan-like population that admixed with modern human populations separated from the population represented by the Denisovan reference genome between 276 kya and 403 kya [52]. Thus, the advantageous EPAS1 haplotype shows far greater similarity to Denisovans than would be expected in the absence of archaic admixture, but substantial mismatch should be expected given the apparent diversity of archaic populations outside of Africa [53].
Our CMS test combines five statistical tests of positive selection (iHS, XP-EHH, F ST , ΔiHH and ΔDAF), offering a more robust performance across a wide variety of scenarios compared to any single test [28]; this allowed us to identify novel candidate genes contributing to adaptations on the Tibetan plateau. One of the candidate gene for adaptive selection is VDR, a gene encoding nuclear hormone receptor for 1, 25 dihydroxyvitamin D 3 and functions in Vitamin D metabolism. Vitamin D is a secosteroid nutrient that plays an important role in calcium homeostasis and bone mineralization, and mainly comes from two sources: 1) through skin exposure to sunlight and 2) food, such as fish, milk and egg yolk. Deficiency in vitamin D leads to impairment of bone mineralization and skeletal deformities, known as rickets in children and osteomalacia and osteoporosis in adults [54]. Low level of vitamin D has also been associated with cancers, diabetes, autoimmune diseases, hypertension, and infectious disease [35,55]. In a study involving vitamin D status in a cohort of 63 Tibetans [54], the proportion of nomadic Tibetans with vitamin D deficiency (25(OH)D <75 nmol/L) was 100% with 80% of people having severe deficiency (25(OH)D < 30nmol/L); the proportion of non-nomadic Tibetans with vitamin D deficiency ranges from 40% to 83%. Consistently, in another study, 61% of Tibetan children suffer from rickets and 51% have stunted growth [56]. Such a high prevalence of vitamin D deficiency may be explained by the traditional Tibetan diet consisting of barley, yak meat and butter tea which are poor sources of vitamin D, and clothing habits in cold temperatures which allow for minimal skin exposure to the sunlight [54]. Therefore, we hypothesize that VDR gene is positively selected to compensate for the lack of vitamin D, the mechanism of which remains to be determined. Previously, it has been shown that genomic regions bound by VDR were under adaptive selection in the human genome [57].
A second promising candidate identified in our CMS test is PTGIS, a gene associated with three hypoxia-response GO terms (0071456, 0036294,0071453) and regulated by hypoxiainducible factor 1 (HIF-1). PTGIS encodes prostaglandin I2 synthase, which participates in the synthesis of prostanoid. Previous studies have shown that the expression of PTGIS is induced by hypoxic conditions in human lung fibroblast cells and cancer cell lines [58], and can activate vascular endothelial growth factor. Therefore, it is possible that a Tibetan-unique genomic variant may induce vasodilation and angiogenesis in response to hypoxia by altering the expression of PTGIS.
The CMS test also identified KCTD12 as a new high-altitude adaptation candidate. Previously, it has been shown in B16 murine melanoma cells that KCTD12 is up-regulated in hypoxic conditions [59], and thus it may play a role in the hypoxic responses. This gene is predominantly expressed in fetal organs, suggesting a potentially important role in early development, but dramatically lower levels in adult tissue including brain and lung [60]. KCTD12 is differentially expressed in human pulmonary endothelial cells upon 48 hours of hypoxia exposure [59] and noted as one of 40 CpG sites with the greatest difference in methylation levels between highland and lowland Oromo Ethiopians [61].
Our functional annotation of Tibetan genomic variants identified a frequent protein-coding SNV in ALPP, a gene associated with pregnancy losses [25]. Since Tibetans tend to have higher birth weight in high altitude environments compared to lowland native populations living at similar altitudes [62], it is possible that the nonsynonymous variant in ALPP may affect the birth weight, development, and pregnancy outcomes in Tibetans. For every 1000 meters of altitude, birth weight decreases an average of 100 grams due to restricted fetal growth [63][64][65][66], with a three-fold increase in the number of infants born small for gestational age (SGA) at high altitude [67]. Weights of infants born at high altitudes to mothers of highland origin, however, are greater than those of lowland origin. This is shown specifically in infants born to Bolivian versus European women in the Andes, whose birth weights are about 300g higher in the former group, while no difference is reported by ethnicity at low altitude [67]. Whether these putatively adaptive markers play roles in this process remains to be determined.
In summary, we performed a comprehensive genomic analysis on whole-genome sequence data from 27 Tibetan individuals. Our analyses detected evidence of population structure between the ancestral Han and Tibetan subpopulations beginning between 44 and 58 kya, although admixture rates between the two subpopulations remained high until around 9 kya. The Denisovan EPAS1 haplotype introgressed into the Tibetan population between 12 and 32 kya, and positive adaptive pressure on this haplotype began between 7 and 28 kya. We summarized the dates of important demographic events in Fig 5. The Our CMS test identified novel candidate genes for high-altitude adaptation including KCTD12, VDR and PTGIS, and also generated a list of candidate variants within the EPAS1 gene region. We estimated that 0.4% of the Tibetan genome are introgressed DNA from Denisovans, although EPAS1 is probably the only introgressed locus that was influenced by strong positive selection in Tibetans. Our study provided a rich genomic resource of the Tibetan population and generated hypotheses for future positive selection tests.

Sequencing samples and variants calling
Our study subjects included 27 Tibetan individuals recruited from Tibetans living in San Diego, California (n = 3), Salt Lake City, Utah (n = 5), and the United Kingdom (n = 19). The sample type and place of origin information are listed in S9 Table. Whole-genome sequencing and variant-calling was performed by CG (v2.0 for the first 17 samples and v2.5 for the remaining 10 samples). Variants with genotyping quality less than 30 (i.e., more than 1 calling error among 1000 variants) were converted to missing genotype calls, and sites with more than 5% missing genotype rate were removed. To identify additional variants enriched for higher genotyping error rates, we compared 62 genomes that had been sequenced in both CG public genome data [33,68] and 1KG Project [18] Phase I data. We identified 652,899 SNVs with more than 5% discordant genotype calls between these two datasets, and removed them from all further analyses. We then performed haplotype phasing and imputation using Shapeit2 [69]. After previous steps, 10,405,415 SNVs remained for subsequent analyses.

Whole-genome positive selection scan using Composite of Multiple Signals (CMS) test
The CMS test combines the signals from five different tests (ΔiHH, iHS, XP-EHH, F ST and ΔDAF) to create a single test statistic [28]. On each variant, CMS calculates the posterior probability of the variant being selected in a naïve Bayes framework by taking the product of the posterior probabilities from each of the five tests. A higher CMS score is consistent with a stronger signal of positive selection. Our implementation of the CMS test closely follows the recommendations given in the previous report [28]. Briefly, we first performed coalescence simulations using cosi [34] based on the demographic model estimated by @a@i. To simulate genetic data under various scenarios of selective sweeps, we used the following selection coefficients (s): 0.02, 0.03 and 0.04; the following start times of selection: 0, 100, 200, 300 and 400 generations ago (the last one represents the oldest selection start time that cosi may simulate within our model); and the following non-reference allele frequencies at the conclusion of the sweep: 0.2, 0.4, 0.6 and 0.8. In addition, we also simulated the scenario where no selective sweep had occurred, which corresponded to the null model in the CMS test. In each situation, we performed 1000 simulations. We calculated the CMS scores for all SNVs that passed our quality-control procedures. To calculate their p-values, we also ran the CMS test over the null model that we simulated above. The p-value for a certain CMS score was calculated as the proportion of SNVs in the null model that has the same or higher CMS scores.

Testing for Denisovan admixture using the D statistic
The D statistic tests whether the amount of archaic admixture in one population exceeds that of another population by examining variant allele frequency spectrum in the case population, control population, an archaic population and an outgroup population [70]. In our study, these were selected as Tibetans, Han Chinese/Yoruba, Denisovan and chimpanzee. Letp i 1 ,p i 2 ,p i 3 , andp i 4 , be the allele frequency of the i-th variant in the case, control, archaic and outgroup populations, respectively. D is defined as: We estimated the confidence interval of D using bootstrap sampling with 200 replicates. Specifically, we divided the whole genome into 1MB blocks, and then randomly sampled with replacement the same number of 1MB blocks as in the actual data, each time calculating a sample D statistic. The 2.5% and 97.5% quantiles of the 200 D statistics values were used as the two end points of the 95% confidence interval.

Identifying regions of Denisovan admixture using the D* statistic
We developed the D Ã statistic as a complement to S Ã to identify local genomic regions with excess archaic admixture in a population. We first intended to use the D statistic to identify 200kb genomic regions with archaic admixture in the case population, but noticed that the variance of the D-statistic was high if the number of archaic variants in a target region was small. This sensitivity of the D-statistic to regions of low archaic gene flow has been noted previously, and alternatives to the D-statistic have been suggested [71]. To illustrate this sensitivity, consider the extreme scenario where only one archaic variant exists in a 200kb genomic region, and this variant was present in only one case genome and absent from all control genomes. In this situation, D will be the largest possible value Eq (1) despite the fact that the evidence supporting an archaic admixture in this region is poor. This property of D makes it unsuitable for prioritizing candidate introgression regions, since the majority of regions with high D values will be those with few archaic variants in modern human genomes.
Previous work has suggested using direct estimators of the quantity of gene flow as an alternative to a method like D-statistics due to the sensitivity of the statistic to population history and other exogenous parameters [71]. However, D-statistics and the related f-statistics vary only in magnitude as a result of these biases, and will only deviate from zero, in expectation, as a result of gene flow [46]. Normalization of D-statistics preserves their desirable characteristics, while eliminating a source of bias.
To normalize D-statistics, we calculated the following statistic (U) for each 200kb genomic region: As can be seen, U is the denominator of the D in Eq (1) and is negatively correlated with the standard deviation associated with D ( S7 Fig). To create a normalized D-statistic (D Ã ) that can be used to compare admixture signals across regions, we sort all 200-kb genomics regions according to their U values and then divide the regions into 20 equal-sized groups. Then D Ã of a 200kb region is equal to its D value divided by the within-group standard deviation. We estimated the significance of a given D Ã by comparing it to the null distribution of D Ã generated by coalescence simulation.

Estimating the proportion of Denisovan admixture
Because only one Denisovan genome is available, we use the statistic Q of Rogers and Bohlender (Equation 11 of [48]), interpreting the primary source of introgression as Denisovan rather than Neanderthal. The expectation of Q depends not only on the fraction, m D , of Denisovan admixture but also on m N , the fraction of ghost admixture from another archaic such as Neanderthal. Equating observed and expected values defines m D as an implicit function of m N . To evaluate this function, we assumed the parameter values in table 3 of [48]. The result is shown as a solid black line in S4

Availability of data and material
The database of Genotypes and Phenotypes (dbGAP) accession number for the whole-genome sequencing data reported in this paper is phs001338.