Population structure analyses performed on the extended dataset including Tibetan, Sherpa, and lowland East-Asian individuals.

(A) Admixture analysis showed the best predictive accuracy when seven (K=7) population clusters were tested. Populations included in the dataset are labelled according to population names and acronyms reported in Supplement table 1. (B) Map showing geographic location and admixture proportions at K=7 of the high-altitude groups included in the extended dataset. The label Tibetans_WG indicates whole genome sequence data for individuals of Tibetan ancestry analysed in the present study. Additional information about the considered samples (e.g., number of individuals per group, reference study, and used abbreviations) are reported in Supplement table 1. (C) PCA plot considering PC1 vs PC2 and summarizing genomic divergence between high-altitude Tibetan/Sherpa people and the cline of variation observable for lowland East-Asian populations. The enlarged square displays clustering between Tibetan samples sequenced for the whole genome (i.e., blue dots) and Tibetan samples characterized by genome-wide data (i.e., light-blue squares).

Distribution of VolcanoFinder statistics suggestive of putative adaptive introgrossed loci across the TBC1D1 and PRKAG2 genomic regions.

On the x-axis are reported genomic positions of each SNV, while on the y axis are displayed the related statistics obtained. Pink background indicates the chromosomal interval occupied by the considered genes, while the grey background identifies those genes (i.e., PGM2 in the TBC1D1 downstream genomic region and the RHEB gene in the upstream PRKAG2 region) possibly involved in regulatory transcription mechanisms. The dashed red line identifies the threshold set to filter for significant LR values (i.e., top 5% of LR values). For both these genomic regions the distribution of LR and -logα are concordant with those observed at the EPAS1 positive control for AI. (A) A total of 50 significant LR values (red stars) and -logα (grey diamonds) values resulted collectively elevated in both the TBC1D1 gene and its downstream genomic regions. A remarkable concentration of significant LR values characterizing 19 SNVs was especially observable in the first portion of the gene. (B) The entire PRKAG2 genomic region was found to comprise 46 SNVs showing significant LR values, with the greatest peaks being located in the downstream region associated to such gene. Peaks detected for the LR statistic are accompanied by peaks of -logα values.

Significant gene networks including Denisovan-like derived alleles according to the Signet analysis.

(A) Schematic representation of the activation of the RAS/MAPK(ERK) axis after interaction of the bradykinin receptors with their ligands (e.g., ANG II) within the framework of the Pathways in Cancer network. Genes supported by both Signet (i.e. belonging to the significant network associated to Pathways in cancer) and VolcanoFinder (i.e., including at least a SNV showing LR value within top 5% of the obtained results) analyses as potentially introgressed loci, are highlighted in red and present solid outline. Grey circles with dotted-dashed contour instead indicate genes supported only by Signet, while loci marked with stars are those including genomic windows showing LASSI T statistic within top 5% of the related distribution.

After the interaction between ANGII (active enzyme angiotensin II) and bradykinin receptors, activation of the Ras protein encoded by KRAS mediated by RAS-GTPases (e.g., RASGRF2) comports a series of phosphorylation reactions that eventually promotes angiogenesis (Kranenburg et al., 2004). In detail, phosphorylation of the MAPK1 protein and prevention of MAPK1-DAPK-1-dependent apoptosis leads to increased MAPK1 activity (Kanehisa & Goto, 2000; Stevens et al., 2007) that causes improved FOS mRNA expression (Monje et al., 2005). FOS together with other proteins (e.g., Jun) forms the AP-1 transcription factor, which bounds to the VEGF promoter region upregulating its expression in endothelial cells (Catar et al., 2013) and sustaining angiogenesis when the HIF-1 signalling cascade is inhibited (Lorenzo et al., 2014). (B) Gene network built by setting co-expression as force function and by displaying the entire set of genes identified by the Signet algorithm as belonging to significant pathways including Denisovan-like derived alleles. Genes whose variation pattern was supported by both VolcanoFinder and Signet analyses (e.g., TBC1D1) as shaped by archaic introgression are displayed with a solid black outline. The EPAS1 positive control locus that has been previously proved to have mediated adaptive introgression in Tibetan populations was represented as light-blue octagonal. Genes included in pathways involved in angiogenesis (e.g., RASGRF2) and/or activated in hypoxic conditions (e.g., PRKAG2) are reported as dark red and light-blue circles respectively, while the remining fraction of significant genes are represented as light-grey circles. The closeness or the distance between all nodes reflects the tendency to be co-expressed with each other and all the connections inferred are characterized by a confident score ≥ 0.7.

Representation of genetic distances between modern and archaic haplotypes and barplots showing haplotype frequency spectra for TBC1D1 and RASGRF2 candidate AI genes.

Haplotypes are reported in rows, while derived (i.e., black square) and ancestral (i.e., white square) alleles are displayed in columns. Haplotypes are ranked from top to bottom according to their number of pairwise differences with respect to the Denisovan sequence. (A) Heatmap displaying divergence between Tibetan, CHB and YRI TBC1D1 haplotypes with respect to the Denisovan genome. A total of 33 TBC1D1 haplotypes (i.e., 61% of the overall haplotypes inferred for such a region) belonging to individuals with Tibetan ancestry are plotted in the upper part of the heatmap thus presenting the smallest number of pairwise differences with respect to the Denisovan sequence. (B) Heatmap displaying divergence between Tibetan, CHB and YRI RASGRF2 haplotypes with respect to the Denisovan genome. A total of 16 Tibetan haplotypes in the RASGRF2 genomic region present no differences with respect to the Denisovan sequence.

As regards barplots, on the x axis are reported the haplotypes detected in the considered genomic windows, while on the y axis is indicated the frequency for each haplotype. The black and dark-grey bars indicate the more frequent haplotypes (i.e., the putative adaptive haplotypes inferred by the LASSI method), while red stars mark those haplotypes carrying Denisovan-like derived alleles. (C) TBC1D1 haplotype frequency spectrum. The TBC1D1 gene presents a haplotype pattern qualitatively comparable to that observed at EPAS1 (Figure 4 – figure supplement 3A), with a predominant haplotype carrying archaic derived alleles and reaching elevated frequencies in Tibetan populations. In line with this observation, such a pattern was inferred by LASSI as conformed with a non-neutral evolutionary scenario, even if it seems to be characterized by a soft rather than a hard selective sweep due to the occurrence of three sweeping haplotypes. (D) RASGRF2 haplotype frequency spectrum. A soft selective sweep was inferred also for the considered RASGRF2 genomic window, although frequencies reached by the sweeping haplotypes turned out to be more similar with each other. The second most represented haplotype was that carrying the archaic derived alleles and, reached a frequency of 29% in the Tibetan group.