Denisovan introgression has shaped the immune system of present-day Papuans

Abstract

Modern humans have admixed with multiple archaic hominins. Papuans, in particular, owe up to 5% of their genome to Denisovans, a sister group to Neanderthals whose remains have only been identified in Siberia and Tibet. Unfortunately, the biological and evolutionary significance of these introgression events remain poorly understood. Here we investigate the function of both Denisovan and Neanderthal alleles characterised within a set of 56 genomes from Papuan individuals. By comparing the distribution of archaic and non-archaic variants we assess the consequences of archaic admixture across a multitude of different cell types and functional elements. We observe an enrichment of archaic alleles within cis-regulatory elements and transcribed regions of the genome, with Denisovan variants strongly affecting elements active within immune-related cells. We identify 16,048 and 10,032 high-confidence Denisovan and Neanderthal variants that fall within annotated cis-regulatory elements and with the potential to alter the affinity of multiple transcription factors to their cognate DNA motifs, highlighting a likely mechanism by which introgressed DNA can impact phenotypes. Lastly, we experimentally validate these predictions by testing the regulatory potential of five Denisovan variants segregating within Papuan individuals, and find that two are associated with a significant reduction of transcriptional activity in plasmid reporter assays. Together, these data provide support for a widespread contribution of archaic DNA in shaping the present levels of modern human genetic diversity, with different archaic ancestries potentially affecting multiple phenotypic traits within non-Africans.

Fig 1. Distribution of common-to-high-frequency aSNPs across 15 chromatin states and 111 cell types.

Patterns of enrichment across the 15 chromatin states for the set of common-to-high-frequency Denisovan and Neanderthal aSNPs relatively to the matched background set of naSNPs. Histograms on top indicate the mean number of variants within each chromatin state across all 111 cell types. Asterisks indicate BH-corrected Fisher’s exact test p-values: * = p ≤ 0.05, ** = p ≤ 0.01; *** = p ≤ 0.001; **** = p ≤ 0.0001 (for full statistical results see S5 Table).

Fig 2. Distribution of common-to-high-frequency aSNPs falling within CREs across 18 tissues.

Patterns of enrichment across the 18 different tissue types in Roadmap Epigenomics for the set of common-to-high-frequency aSNPs annotated within chromatin states associated with CREs, relatively to the matched background set of naSNPs. Histograms on top indicate the mean number of variants calculated across all cell types belonging to each tissue. Asterisks indicate BH-corrected Fisher’s exact test p-values: * = p ≤ 0.05, ** = p ≤ 0.01; *** = p ≤ 0.001; **** = p ≤ 0.0001 (for full statistical results see S7 Table).

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 F_st values between western Indonesia and PNG for the three groups of variants.

Fig 4. Biological processes putatively affected by SNPs within immune-related CREs.

For each ancestry the figure shows A) the overlap between the putative sets of target genes for each ancestry; B) the semantic similarity score estimates between the set of significantly enriched GO terms.

Fig 5. Functional validation of the regulatory impact of Denisovan variants near OAS2 and OAS3.

A) The genomic region encompassing the eight Denisovan variants associated with OAS2 and OAS3. The top two tracks display patterns of DNase Hypersensitivity sites in Blood and immune T cells as well as in HSC and B cells [34]. Tested variants are shown along with their calculated Δ PWM. Bottom tracks display the chromatin state information for the same tissues; B) the distribution of the log₂ RNA-seq counts per million in whole blood for OAS2 and OAS3 between the Korowai and the people of Mentawai, from [57]; C) the relative expression changes between Denisovan and non-archaic alleles, or between the alternative and the reference allele for positive control (rs9283753) [45], in two Papuan LCLs. Asterisks mark significant differences from 1, BH-corrected p < 0.05: *** = 0.0001 < p ≤ 0.001; ** = 0.001 < p ≤ 0.01; * = 0.01 < p ≤ 0.05.