Denisovan Ancestry in Papua New Guinea (PNG)

Fig 3. Levels of population differentiation for the set of cis-regulatory SNPs between western Indonesian and New Guinean populations.

A) Proportion of variants shared across the Indonesian archipelago segregating within the corresponding haplotypes. B) Distribution of the Fst values between western Indonesia and PNG for the three groups of variants.

doi.org/10.1371/journal.pgen.1010470.g003

Discussion
There is significant interest in understanding the functional consequences of archaic introgression. Evidence indicates that both Denisovan and Neanderthal aSNPs, especially those within protein coding and conserved non-coding elements, are mostly deleterious and negatively selected in modern humans [14, 61]. Similar findings have been recently reported for highly pleiotropic enhancers, where aSNPs are depleted likely as a consequence of their potential to perturb gene expression across multiple tissues [12]. Nevertheless, out of the substantial number of archaic variants still segregating within present-day populations, a large fraction falls within genomic regions that show strong evidence of functional activity across a variety of cell types. Indeed, studies conducted primarily on Neanderthal introgressed DNA have suggested a non-negligible contribution to gene expression variation in modern humans [10, 11, 22], with repeated examples of Neanderthal archaic variants falling within regulatory elements or the seed region of mature micro-RNAs predicted to affect transcriptional and post-transcriptional regulatory processes [11, 62].

In this study, we have taken advantage of a recently published dataset [6] to investigate the landscape of archaic introgression in individuals of Papuan genetic ancestry, the functional consequences of which remained poorly understood. This has allowed us to characterise the putative contribution of Denisovan DNA, which is known to account for up to 5% of the genome of present-day Papuans [4, 63], while also comparing it to that of Neanderthal DNA and non-archaic variants that arose following the Out of Africa event. We specifically analysed all of our variants across multiple cell types and functional chromatin states aiming to account for the strong dependency on the cells’ chromatin landscapes [34] when assessing the potential activities of introgressed alleles. While our analyses focus on a small subset of introgressed variants identified by [6], our filtering criteria ensure that we only consider variants with a high likelihood of being truly introgressed. In addition, comparing Denisovan and Neanderthal signals within the same samples serves as an internal control for difference in linkage disequilibrium and background selection between ancestries, hence providing stronger support for any eventual Denisovan- and/or Neanderthal-specific contribution to phenotypic variation in contemporary Papuans.

First, we confirm previous reports that aSNPs mostly occur within non-coding sequences [9, 14, 64], with a large number falling within highly constitutive and/or functionally inert regions, which might be the result of weaker purifying selection acting on these elements. aSNPs are significantly depleted from other inactive chromatin states but over-represented within states putatively containing elements involved in gene regulatory processes. Applying our approach to Neanderthal aSNPs segregating within present-day West Eurasians [31] replicates most of our findings, as well as previously reported signals including a depletion of Neanderthal variants within enhancer-like elements [12]. Previous studies have shown that, while polygenic risk scores [65] and eQTLs [66] have limited cross-population portability, biochemically active functional regulatory elements tend to be under evolutionary constraint within the human lineage [67, 68]. Thus, while our approach is contingent on functional data generated in a population that is not of Papuan ancestry, our focus on functional analyses are likely to be more robust to differences across populations in linkage disequilibrium. Nevertheless, our set of SNPs, and our conclusions, represent a first step that will benefit from functional validation at scale.

Second, we find notable differences between the two archaic ancestries in their patterns of introgression. Such differences are likely consequence of the tissue-specificity nature of CREs, and of enhancers in particular [69]. Indeed we report that variants annotated in these CRE-associated states tend to lie within elements active in a restricted number of cell types, suggesting that the effects of archaic introgression might be highly cell-type-specific. In addition, admixture with Neanderthals might have similar consequences even across distantly related human populations, with Neanderthal aSNPs strongly enriched within elements active within blood & immune T cells and mesenchymal cells both in Papuans and Europeans, replicating previously reported observations [11, 12]. While we report a similar enrichment for Denisovan alleles within these tissues, we further note a Denisovan aSNPs-specific significant enrichment with CREs active within HSC & B cells, suggesting a more pervasive contribution of Denisovan introgression to gene regulatory processes happening within immune-related cells. In line with this, we further show that Denisovan and Neanderthal variants annotated within immune-related CREs might target different sets of genes with limited overlap in the biological processes they are involved in. Indeed, only genes predicted to be regulated by Denisovan aSNPs are strongly involved in active immune responses.

Third, we have characterised the impact of introgressed variants on transcription factor binding sites. A substantial proportion of aSNPs within CREs are predicted to modify the interactions between TFs and their cognate binding sites, although we do not detect an overall higher impact of aSNPs relative to naSNPs. Our data suggest that, despite aSNPs possibly disrupting individual DNA motifs, at a genome-wide level admixture with archaic hominins is unlikely to have resulted in large rewiring of transcription factor regulatory networks. While this might contrast the findings of previous studies [70], in our study we focused on families of DNA motifs (as clustered by [41]) rather than consider enrichment within individual PWMs, an approach that we believe is better suited to controlling for the typical redundancy of PWMs. This, however does not exclude that aSNPs might not disrupt specific DNA motifs. Indeed, our approach also identifies Neanderthal aSNPs disrupting the PWM of TFs such as CREB1, USF1, MYB and STAT6, similarly to what has been reported by [70]. We also note that there is a significant (hypergeometric p < 2.2 ⋅ 10−16) overlap between genes targeted by N-eQTLs in [70] and our own set of possible target genes (S9 Fig).

Finally, we performed a plasmid reporter assay experiment to quantify the molecular impact of Denisovan TFBS-disrupting variants annotated within CREs active in immune-related cells and predicted to affect the expression level of OAS2 and OAS3. These genes belong to a family of pattern-recognition receptors involved in innate immune responses against both RNA and DNA viruses, with OAS3 considered to be essential in reducing viral titer during Chikungunya, Sindbis, influenza or vaccinia viral infections [71]. At least two previous studies have shown the presence of both Neanderthal and Denisovan archaic haplotypes segregating at this locus respectively within European [52] and Papuan [56] individuals. Sams et al. found two variants (rs10774671, rs1557866) within these Neanderthal haplotypes which are respectively associated with the codification of different OAS1 splicing isoforms and with a reduction in OAS3 expression levels, the latter only upon viral infection [52]. Our previous work has found that both OAS2 and OAS3 are differentially expressed between western Indonesians and Papuans [57]. Here we report a set of eight SNPs (rs368816473, rs372433785, rs139804868, rs146859513, rs143462183, rs370655920, rs375463218 and rs372139279), located roughly 41 kb upstream of OAS2 and OAS3, all of which are predicted to strongly alter the ability of different transcription factors, including IRF4, NFKB2 and TAL1, to bind to their underlying DNA motifs. In all cases, the reference allele is fixed within western Indonesian populations, whereas the archaic alleles segregate at frequencies between 0.2 and 0.4 in Papuans. We show that at least 5 of these variants lie within sequences that can regulate expression in reporter gene plasmid assays, and that in two of these variants (rs139804868 and rs146859513) the Denisovan allele is associated with significantly lower transcriptional activity compared to the non-archaic allele in immune cells from two different Papuan donors, again suggesting that these SNPs are of biological importance.

Recent work by [72] has shown that variants which arose in the human lineage following its split from the Neanderthal-Denisovan common ancestor, and which have reached (near) fixation in modern humans since then can exhibit significant differences in regulatory activity relative to the ancestral state. Similarly, Jagoda et al. has shown large impact on gene expression associated with the presence of Neanderthal aSNPs in vitro [73]. Our results suggest that Denisovan alleles segregating within modern human populations are also likely to actively participate in gene regulatory processes, especially those that take place within immune-related cells. This agrees with recent findings from a study that analysed the genome of present-day people of Taiwan, the Philippines, the Solomon Islands and Vanuatu [50]. While further experimental validation of our observations is necessary in order to characterise the genome-wide impact of archaic introgression, the results presented in this study argue for a potential contribution of Denisovan variants to immune-related phenotypes amongst early modern humans in the region, potentially favouring adaptation to the local environment [15, 22].