Evolutionary history of Tibetans inferred from whole-genome sequencing
[['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–4] .', 'Many candidate genes have been highlighted as likely contributors to Tibetans’ trans-generational success at high altitude[5–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–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 .'], ['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 5th 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 .', 'FST 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 FST between Chinese and Japanese populations but lower than the FST 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 Conservation-controlled 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 .', 'FST predicted by the best-fitting demographic model in ∂a∂i produced estimates consistent with observed FST ( 0 . 0147 using Hudson’s FST estimator and 0 . 0148 using the Weir and Cockerham’s estimator ) .', '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 , FST 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 , FST 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–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; p = 1 . 1x10-15 ) .', '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 ( CMSq<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 CMSq<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 Denisovan-like 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 FST 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 3rd 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 FST 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 ( r2 = 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 . 12C>G variant was absent from any other populations in the 1KG dataset and had high FST ( 0 . 70 ) and ΔDAF ( 0 . 70 ) values , indicative of strong positive selection .', 'The c . 380G>C variant has a derived allele frequency of 77 . 8% ( 95% CI: [66 . 7% , 88 . 9%] ) in Tibetans and 44 . 9% ( 95% CI: [40 . 0% , 50 . 0%] ) 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 .', '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 f3 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 f3-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 .'], ['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 gene-flow ( 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 , FST , Δ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 D3 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 hypoxia-inducible 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 lowe
['Introduction', 'Results', 'Discussion', 'Materials and methods']
['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 .']
['The Tibetan population has been residing on high plateau for thousands of years and developed unique adaptation to the local environment .', 'To investigate the demographic history of Tibetans and search for possible adaptive genetic variants , we performed whole-genome sequencing of 27 Tibetan individuals .', 'We found evidence of genetic separation between Han and Tibetans around since 44 and 58 thousand years ago; however , these two populations maintained a high rate of gene flow until 9 thousand years ago .', 'In addition to replicating two previously discovered candidate genes ( EGLN1 and EPAS1 ) for high altitude adaptation , we also found three new candidate genes , including PTGIS , VDR and KCTD12 .', 'We confirmed the high similarity of EPAS1 gene region between Tibetans and Denisovans , but did not detect any evidence of high altitude adaptation from Denisovan gene alleles otherwise .']
['han chinese', 'demography', 'population genetics', 'genetic mapping', 'ethnicities', 'population biology', 'tibetan people', 'introgression', 'comparative genomics', 'people and places', 'haplotypes', 'heredity', 'evolutionary processes', 'genetics', 'population groupings', 'biology and life sciences', 'gene flow', 'genomics', 'evolutionary biology', 'genomics statistics', 'computational biology']