Denisovan Ancestry in Papua New Guinea (PNG)

new

Admin
Administrator

Posts: 80,516

Denisovan Ancestry in Papua New Guinea (PNG) Dec 10, 2022 0:44:40 GMT

Quote

Post by Admin on Dec 10, 2022 0:44:40 GMT

When modern humans first migrated from Africa to the tropical islands of the southwest Pacific, they encountered unfamiliar people and new pathogens. But their immune systems may have picked up some survival tricks when they mated with the locals—the mysterious Denisovans who gave them immune gene variants that might have protected the newcomers’ offspring from local diseases. Some of these variants still persist in the genomes of people living in Papua New Guinea today, according to a new study.

Researchers have known for a decade that living people in Papua New Guinea and other parts of Melanesia, a subregion of the southwest Pacific Ocean, inherited up to 5% of their DNA from Denisovans, ancient humans closely related to Neanderthals who arrived in Asia about 200,000 years ago. Scientists assume those variants benefited people in the past—perhaps by helping the modern humans better ward off local diseases—but they have wondered how that DNA might still be altering how people look, act, and feel today. It’s been difficult to detect the function of Neanderthal and Denisovan DNA in Melanesians, however, because scientists have analyzed so little genetic data from living humans in Papua New Guinea and other parts of Melanesia.

The new study overcomes that problem by using genetic data from 56 individuals from Papua New Guinea that were recently analyzed for another paper, part of the Indonesian Genome Diversity Project. The researchers, mostly from Australia and New Guinea, compared those genomes with those of Denisovans from Denisova Cave in Siberia, as well as Neanderthals. They found the Papuans had inherited unusually high frequencies of 82,000 genetic variants known as single nucleotide polymorphisms, which arise from differences of a single base or letter in the genetic code—from Denisovans.

The team then looked for those variants in a database that links genes to various functions in different tissues in humans. They homed in on immune-related gene variants that might promote or enhance a nearby gene’s production of proteins, for example, or shut off or dampen its function. These tweaks can help optimize an immune system for the specific pathogens in its environment; too strong an immune response can be as deadly as the infection itself.

In Papuans, the scientists found many Denisovan variants that were located near genes known to impact human immune responses to viruses and other pathogens, such as the flu and chikungunya. Next, they tested the function of eight Denisovan gene variants associated with the expression of proteins produced by two genes in particular, OAS2 and OAS3, “lymphoblastoid”—cell lines of B cells, a type of white blood cell that makes antibodies critical to the body’s immune response. Those cell lines were collected from Papuans by study co-author Christopher Kinipi, a Papuan physician and health services director at the University of Papua New Guinea.

Two of the Denisovan genetic variants found in those Papuan cell lines lowered the transcription or production of proteins that regulate cytokines, part of the immune system’s defense against infections, reducing inflammation. This subdued inflammatory response could have helped Papuans weather a rash of new infections they would have encountered in the region.

“One of the strengths of the study is that they tested the Denisovan variants in Papuan cell lines, which are essentially the cell environment in which they evolved,” says functional genomicist Francesca Luca of Wayne State University, who was not part of the study.

Taken together, these experiments suggest those Denisovan gene variants “might be fine-tuning the immune response” to optimize it to its environment, says human evolutionary geneticist Irene Gallego Romero of the University of Melbourne, lead author of the new study published in PLOS Genetics. “In the tropics where people have high loads of infectious disease, you might want to tone down the immune response a little and not go overboard.”

These findings dovetail with earlier work on the role of Neanderthal variants in living Europeans. Studies of both Neanderthal and Denisovan DNA in different populations are showing how mating with archaic humans—long-adapted to their regions—provided a rapid way for incoming modern humans to pick up beneficial genes, says computational biologist Janet Kelso of the Max Planck Institute for Evolutionary Anthropology. The study shows this sort of gene swap was “an important mechanism for how humans adapted quickly [to new challenges], specifically pathogens,” says human geneticist Luis Barreiro of the University of Chicago.

But he would like to see future work test whether the Denisovan gene variants actually give Papuans a better shot at warding off or surviving specific diseases.

Overall, this study shows “matings which took place tens of thousands of years ago are still influencing the biology of contemporary individuals,” says population geneticist Joshua Akey of Princeton University.

Admin
Administrator

Posts: 80,516

Denisovan Ancestry in Papua New Guinea (PNG) Dec 15, 2022 19:14:16 GMT

Quote

Post by Admin on Dec 15, 2022 19:14:16 GMT

Denisovan introgression has shaped the immune system of present-day Papuans
Davide M. Vespasiani,Guy S. Jacobs,Laura E. Cook,Nicolas Brucato,Matthew Leavesley,Christopher Kinipi,François-Xavier Ricaut,Murray P. Cox,Irene Gallego Romero
Published: December 8, 2022
doi.org/10.1371/journal.pgen.1010470

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.

Author summary
Humans of Papuan ancestry owe roughly 5% of their genome to Denisovans, a poorly characterised archaic hominin. While introgressed DNA segments can be readily identified, understanding their biological consequences remains challenging. By examining the distribution of introgressed DNA against existing functional genomics datasets, it is possible to predict the phenotypes they impact. In Papuans, Denisovan DNA, but not Neanderthal, strongly and consistently affects immune cells and immune-related processes of potential evolutionary relevance. In vitro testing of introgressed variants confirms these predictions, suggesting Denisovan variants can impact gene regulation in vivo. Variation in gene expression might be key to understanding the consequences of admixture between modern humans and archaic hominins, as has been observed with Neanderthal DNA in other human populations.

Introduction
Modern humans are known to have interbred with Neanderthals [1], Denisovans [2] and possibly other archaic hominins [3]. While genetically similar populations of Neanderthals are thought to have contributed approximately 2% to non-African genomes, Denisovan introgression has been observed to be more variable [4]. Particularly, Denisovan ancestry accounts for up to 5% of the genomes of the Indigenous peoples of Island Southeast Asia and Australia [4, 5]. In addition, these components exhibit a deep divergence from the reference Altai Denisovan genome, providing strong evidence for the occurrence of multiple Denisovan introgression events across time and space [6, 7].

At the genomic level these introgressed archaic alleles are mostly observed outside protein-coding sequences, distributed over non-functional and regulatory regions [8, 9]. Enhancers, in particular, are amongst the top targeted elements of archaic hominin introgression [10, 11]. Here, archaic alleles are thought to drive phenotypic differences by altering gene pre- and post-transcriptional regulatory processes [11]. Furthermore, both Neanderthal and Denisovan variants seem to preferentially affect enhancers in a tissue-specific manner, with highly pleiotropic elements being depleted of archaic variation [12]. However, beyond general agreement that introgressed archaic DNA has mainly been deleterious and actively removed from coding sequences and conserved non-coding elements [8, 12–14], the actual phenotypic consequences of this variation are not well understood. Several lines of evidence highlight associations between archaic DNA and risk for disease traits, including autoimmune diseases [8, 15, 16], or with traits of possible evolutionary advantage for early non-Africans [17–20]. For example, Neanderthal variants within immune genes and immune-related cis-regulatory elements (CREs) have been associated with differential responses to viral infections among present-day Europeans [21, 22].

Unfortunately, two main factors limit our interpretation of the biology of these alleles. First, we still lack a detailed characterisation of global levels of human genetic diversity, vital to identifying differential archaic hominin contributions across modern populations [23, 24]. Because Denisovan DNA is not present in the genomes of modern Europeans, which make up the vast diversity of biomedical genetics cohorts, it has been impossible to link them to phenotype, unlike segregating Neanderthal variants [16, 25]. Second, and more challenging, most of these alleles lie within non-coding sequences where, despite their acknowledged contributions to human evolutionary history [26, 27], an understanding of their actual biological functions remains elusive.

To gain insights into the consequences of archaic introgression in Papuans, we have analysed a previously published dataset of 56 present-day genomes sampled across the island of New Guinea [6]. By comparing the distribution of archaic single nucleotide polymorphisms (aSNPs) and non-archaic SNPs (naSNPs) segregating within these populations, across multiple genomic elements and cell types, we find that aSNPs are enriched within functional cis-regulatory elements (CREs) and transcribed regions, particularly those active within immune-related cells. We also observe that the presence of archaic alleles within these elements can lead to substantial disruption of transcription factor binding sites (TFBSs), with a consistent Denisovan-specific signal affecting immune-related processes. We finally validate these results through experimental testing of multiple Denisovan aSNPs, and find they are associated with significant transcriptional changes in reporter gene experiments.

Admin
Administrator

Posts: 80,516

Denisovan Ancestry in Papua New Guinea (PNG) Dec 17, 2022 5:08:10 GMT

Quote

Post by Admin on Dec 17, 2022 5:08:10 GMT

Materials and methods
Ethics statement
All human genome data in this study was taken from Jacobs et al. [6]. All collections followed protocols for the protection of human subjects established by institutional review boards at the Eijkman Institute (EIREC #90) and Nanyang Technological University (IRB-2014–12-011); the analyses in this publication were additionally approved by University of Melbourne’s Human Ethics Advisory Group (1851585.1). Permission to conduct research in Indonesia was granted by the State Ministry of Research and Technology (RISTEK). All individuals gave their full informed written consent to participate in the study.

The two lymphoblastoid cell lines (LCLs) used in the plasmid reporter experiments were established from donors sampled at the University of Papua New Guinea, in Port Moresby, Papua New Guinea. Collection of material from healthy adult donors was coordinated by Dr Christopher Kinipi (Director of Health Services at the University of Papua New Guinea) and approved by the Medical Research Advisory Committee of the National Department of Health of the Government of Papua New Guinea (permit number MRAC 16.21), by the University of Melbourne (ID: 1851585.1), and by the French Ethics Committees (Committees of Protection of Persons 25/21_3, n°SI:21.01.21.42754). Permission to conduct research in Papua New Guinea was granted by the National Research Institute of Papua New Guinea (permit 99902292358), with full support from the School of Humanities and Social Sciences, University of Papua New Guinea. These approvals commit our team to following all ethical guidelines mandated by the government of Papua New Guinea. All individuals gave their full informed written consent to participate in the study.

Variant filtering and selection
A list of 7,791,042 SNPs segregating within non-archaic, Denisovan and Neanderthal haplotypes identified in 56 genetically Papuan individuals was taken from [6] (for detailed information about samples, genotyping and inference of introgressed state, see [6]). To define a high confidence set of introgressed SNPs and a matched set of SNPs found outside archaic haplotypes to use as controls, we applied a series of stringent filtering steps to the genome-wide data, as follows. First, all singletons were removed, together with variants lacking ancestral state information in the 6 non-human primate EPO alignment [28]. Variants were then labelled as possibly archaic if they fell within archaic haplotypes defined by [6], and further classified as Denisovan or Neanderthal aSNPs on the basis of the inferred ancestry of the haplotype. Due to incomplete lineage sorting and the relatively short divergence time between Denisovan and Neanderthal lineages [2], a single archaic SNP can segregate in both Neanderthal and Denisovan haplotypes. We observed 18,518 aSNPs that segregated in both Denisovan and Neanderthal haplotypes. These variants were assigned to one of the two archaic ancestries first by considering in which archaic haplotype they segregated at higher frequencies and, then, by comparing the state of the main introgressed allele in Papuans to the Altai Neanderthal [29] or the Altai Denisovan [2] reference genomes. This resulted in 6,098 and 9,744 variants that were respectively assigned to Denisovan and Neanderthal ancestries. The 2,676 variants that could not be disambiguated in this manner were labelled as ambiguous and excluded from downstream analyses.

To further control for incomplete lineage sorting between humans, Neanderthal and Denisovans, and possible errors in haplotype calling, allele frequencies were calculated by considering the number of archaic or non-archaic haplotypes carrying the allele of interest relative to the total number of chromosomes in the sample. Specifically, for naSNPs we calculated the derived allele frequency (DAF) as the number of derived alleles (relative to the EPO alignment) segregating within non-archaic haplotypes divided by the total number of chromosomes genotyped at that site. Likewise, for each aSNP, we calculated the main introgressed allele frequency (MIAF) by counting the number of observations of the most abundant allele found within archaic haplotypes divided by the total number of chromosomes genotyped at that site. Variants with MIAF or DAF < 0.05 or ≥ 0.05 were respectively labelled as low-frequency or common to high-frequency. We also calculated the frequency of the main introgressed allele within non-archaic haplotypes for all aSNPs, and then examined the differences in frequency across modern and archaic haplotype backgrounds. This showed that in a large number of cases, alleles that segregate within archaic haplotypes also segregate within individuals that carry a non-archaic haplotype at the same site. To address this confounder, we calculate allele frequencies in archaic and non-archaic haplotypes, and removed all SNPs where the difference in frequency between archaic and non-archaic haplotypes was < 0.25 (S1 Fig). This step does also exclude some aSNPs segregating at overall low frequencies in our samples and thus reduces the overall size of our aSNP sets, but leads to a higher-confidence dataset. Finally, we removed all aSNPs and naSNPs for which the main introgressed allele or derived allele, respectively, segregated at frequencies ≥ 0.005 within the sub-Saharan African populations in the 1000 Genomes Project [30].

Following these steps, we retained a set of 140,916, 88,625 and 1,108,747 variants of putative Denisovan, Neanderthal and non-archaic ancestry, respectively. A detailed overview of the number of variants retained after each filtering step can be found at S1 Table. Inspection of the site frequency spectrum of these three sets (SFS) revealed a different distribution of the allele frequencies between aSNPs and naSNPs (S2(A) and S2(B) Fig). To account for this, we randomly sub-sampled a total of 229,541 naSNPs (corresponding to the combined number of aSNPs) from the set of non-archaic variants to match the SFS of both Denisovan and Neanderthal aSNPs and used this set for all downstream comparisons (S2(C) Fig). The final sets of Denisovan, Neanderthal and non-archaic variants are listed in S2–S4 Tables.

Definition of a control set of Neanderthal variants segregating within West Eurasians
We additionally defined a control set of Neanderthal variants following the same principles as above. Briefly, Neanderthal haplotype inferences were taken from 75 West Eurasian individuals from the Simons Genome Diversity Project [31] that were analysed in [6], retrieving the Neanderthal haplotypes inferred and extracting allele counts for each SNP in Neanderthal and non-Neanderthal inferred haplotypes. Filtering steps to obtain a high-confidence set of Neanderthal variants segregating in West Eurasians were performed as above.

Quantifying the levels of background selection
To ensure the observed differences between aSNPs and naSNPs were not caused by different levels of background selection we downloaded from github.com/gmcvicker/bkgd the bed file containing the list of pre-computed values for the B-statistic in hg19 coordinates [32]. We then intersected these regions with the list of aSNPs and naSNPs to test whether the distribution of the corresponding values were comparable across the three ancestries. As a further control, we also calculated the distribution of the B-values for all know human protein-coding sequences (S3(A) and S3(B) Fig).

Functional annotation
To functionally annotate the refined set of SNPs we first used the annotatr R package [33] to determine the distribution of aSNPs across multiple genomic elements and whether they were significantly depleted from protein-coding sequences. To investigate the function of aSNPs across chromatin states and human cell types, we downloaded the 15-state hg19 mnemonics BED files containing the segmented chromatin state predictions for each of the 111 cell types included in the Roadmap Epigenomics project [34], egg2.wustl.edu/roadmap/web_portal/index.html. Variants were intersected with chromatin state files, and the 111 cell types were grouped into the the 18 tissues as in [34]. We counted the number of cell types across which each SNP was annotated within the same chromatin state as a measure of its potential pleiotropic activity. Finally, we downloaded the list of V8 cis-eQTLs and their eGenes from the GTEx web portal (https://www.gtexportal.org/home/) [35], retaining only the positional information for the list of significant cis-eQTLs (i.e., qval < 0.05). Variants associated with multiple genes were counted only once. The genomic coordinates for the resulting list of cis-eQTLs were then lifted over from hg38 to hg19 using liftOver from the rtracklayer R package [36]. We used this dataset to identify shared instances first only with our set of common-to-high-frequency variants annotated within putative CREs and then with our entire set of 140,916 Denisovan, 88,625 Neanderthal and 229,541 non-archaic variants. In both analyses, we retained only those variants where the reference and alternative alleles matched between the datasets.

Computing Fst values between western Indonesia and New Guinea populations
Fst values were calculated between western Indonesian populations and New Guinean populations using the BCF files from [6] and the VCFtools –weir-fst-pop command. The western Indonesian population set was defined as in [6].

Enrichment of aSNPs across chromatin states and tissues
To quantify the impact of aSNPs within each chromatin state and/or across each chromatin state/tissue combination, we counted the number of Denisovan, Neanderthal and non-archaic variants annotated in each combination of interest. We then computed the odds ratio (OR) between the number of aSNPs and naSNPs for Denisovan and Neanderthal variants. Statistically significant deviations from an OR = 1 were tested using Fisher’s exact test. Unless otherwise noted, resulting p-values were corrected for multiple hypothesis testing using the Benjamini-Hochberg (BH) method [37]. Enrichment scores together with their corresponding p-values for each comparison of interest can be found at S5–S7 Tables.

Impact of aSNPs on transcription factor binding sites
To quantify the impact of archaic alleles on transcription factor binding sites (TFBSs), we focused on variants annotated to TssA, TssAFlnk, TxFlnk, Enh and EnhG in at least one cell type. We then used the motifbreakR R package [38] to identify putative TFBS-disrupting variants. Position weight matrices (PWMs) were retrieved from Jaspar 2018 [39] and HOCOMOCO v.11 [40]. Motifs were matched to the hg19 genome using p-value threshold of 1−5, and the allelic impact on each PWM (PWM score) was quantified via the built-in sum of log-probabilities algorithm specifying a background probability of 0.3 and 0.2 for A/T and G/C, respectively. Differences in PWM score between the two alleles at a site, defined as Δ PWM, were polarised by the introgressed or derived state of aSNPs and naSNPs, respectively. Differences in PWM score between the two alleles at a site, defined as Δ PWM, were polarised by the introgressed or derived state of aSNPs and naSNPs, respectively. When variants were predicted to disrupt a motif associated with the same TF across both databases we retained only the SNP-PWM association with the highest absolute Δ PWM score.

To assess the genome-wide impact of aSNPs on TFBSs while avoiding redundancy between TFs due to similar binding motif preferences, we grouped all TFs into families based on motif clustering information from [41]. We removed 67 HOCOMOCO and 62 Jaspar PWMs that were either not associated with any cluster or that have mismatches in the PWM versions between MotifDb [42] and [41]. We then tested whether any motif cluster contained an excess of TFBS-disrupting aSNPs compared to naSNPs by performing a χ2-test. Differences in Δ PWM score distributions across ancestries were tested via a Wilcoxon test.

Gene ontology enrichment analysis
To understand the biological processes potentially affected by Denisovan and Neanderthal introgression, we retrieved all CRE-associated variants annotated in at least one haematopoetic stem cell (HSC) & B and/or blood & T related cell type. We then used the rGREAT R package package [43] to assign each variant to the regulatory domain of the nearest gene located within 1 Mb of distance and to perform a gene ontology (GO) enrichment analysis [44], using as a background set the union of all the TFBS-disrupting variants associated genes (both archaic and non-archaic). Importantly, by linking variants to the regulatory domain of its nearest gene, this approach does not yield multiple SNP-genes associations for a single SNP which could occur when in presence of gene clusters. Significantly enriched GO terms were defined as those with BH-corrected p-value < 0.01 and can be found in S8 Table.

Quantifying the impact of Denisovan alleles on gene expression using a plasmid reporter assay
For each of the 5 selected SNPs (see Results section), we designed 2 DNA oligonucleotides of 170 nucleotides (nt) centred on the two alleles of interest (84 bp on the 5’ and 85 bp on the 3’ end). To ensure the effectiveness and reproducibility of our design, we also included a positive control variant (rs9283753) from Tewhey et al. 2016 [45]. All variants were synthesised on the positive strand, regardless of variant orientation with respect to its predicted target gene, as enhancer activity in reporter assays has been found to be largely independent of orientation [46]. In addition, all oligos contained two common 15 bp adapter sequences flanking each regulatory element, a 44 bp sequence centred on a BsiWI restriction site, a oligo-specific 15 bp barcode sequence and a universal 20 bp buffer sequence, resulting in a final oligo length of 284 bp. Oligos were synthesised as a pool by IDT. A list with the full sequences and DNA primers used to amplify the oligo library (see below) is available at S9 Table. The pMPRA1 (Addgene, 49349) was chosen as plasmid backbone for the molecular cloning.

Synthesised oligos underwent two consecutive runs of Q5 (NEB, M0491) low-cycle PCR using two different sets of primers, first to make them double-stranded and then to add a SfiI restriction site at each end of the sequence. After each PCR amplification, oligos were purified with AMPure XP beads (Beckman Coulter, A63881) following the manufacturer’s protocol. Amplified oligos and pMPRA1 vector were both SfiI digested (NEB, R0123) for an hour. Digestion solutions were then ran on an agarose gel, extracted using the gel extraction kit (NEB, T1020) and ligated using T4 DNA Ligase (NEB, M0202), again following the corresponding manufacturer’s protocol. Ligated products were then transformed into E.coli competent cells (Promega, JM109), which were grown overnight in 25 mL of LB media supplemented with 100 μg/mL of ampicillin (Thermo Fisher Scientific, FSBBP1760–5) prior to plasmid purification (Zymo, D4209). To create the final library, assembled plasmids were linearised using the BsiWI-HF (NEB, R3553) restriction enzyme. An amplicon containing a minimal promoter, GFP open reading frame (ORF), synthesised as a gblock by IDT (S9 Table), was then inserted by Gibson assembly (NEB, E2621). The resulting ligated products were again transformed into bacterial cells, which were overnight grown before subsequent plasmid extraction (Zymo, D4200). For each step, the transformation efficiency was estimated to be > 106 cfu. We ensured correct assembly of each of the DNA oligos by Sanger sequencing (Macrogen).

The two LCLs used for experimental validation (PNG15 and PNG21) were grown in RPMI (Lonza, 12–702Q) supplemented with 1% Glutamax, 1% NEAA and 10% FBS. Cells were maintained at a cell density of 5 ⋅ 105 cells/mL. Each cell line was electroporated in triplicate using the Neon transfection system (Thermo Fisher Scientific, MPK5000). For each transfection, 1 ⋅ 106 cells were centrifuged at 700x g and resuspended in 100 μl of R buffer with 1 μg of plasmid library before applying 3 pulses of 1200 V for 20 ms each [45]. After transfection, cells were allowed to recover in 2 mL of culture media for 24 hours then collected by centrifugation and frozen at -80°C. For each transfection replicate, total DNA, including plasmid DNA (pDNA), and RNA were extracted from cells using Qiagen AllPrep kit (Qiagen, 80204) following the manufacturer’s protocol, including on-column DNase digestion of the RNA fraction. pDNA and total RNA samples were quantified using the Qbit broad range dsDNA (Thermo Fisher Scientific, Q32850) and RNA kits (Thermo Fisher Scientific, Q10211), respectively. First-strand cDNA was synthesised from 1 μg of mRNA for each sample using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, 4368814) and an universal custom reverse primers binding downstream to the oligo-specific barcode sequence (see S9 Table).

The regulatory activity of each oligo was then quantified via qPCR using the NEB Luna Universal qPCR Master Mix (NEB,M3003). We first assessed relative primer binding efficiencies through a series of 4 (1:10) serial dilutions of an aliquot of cDNA as in [47]. These efficiencies were then incorporated in the calculation of each oligo’s relative expression (see below). For each transfection replicate we quantified the relative abundance of pDNA and cDNA for each oligo using a common forward primer binding within the GFP ORF and 12 different reverse primers, each targeting the oligo-specific barcode sequence. All reactions were performed in triplicate. Relative expression of each oligo for each transfection replicate was calculated with the ΔCt method, normalising it by the relative amount of pDNA. Expression differences between archaic and non-archaic alleles were then computed using the relative expression ratio, using a t-test to quantify any significant deviation from an expected ratio of 1 [48], BH-correcting the resulting p-values.

Admin
Administrator

Posts: 80,516

Denisovan Ancestry in Papua New Guinea (PNG) Dec 19, 2022 5:09:39 GMT

Quote

Post by Admin on Dec 19, 2022 5:09:39 GMT

Results
Building a high-confidence set of archaic variants
To curate a list of high-confidence archaic SNPs to analyse for possible functional contributions to present- day Papuans, we took advantage of a recently described dataset of both archaic and non-archaic haplotypes segregating within 56 individuals of Papuan ancestry [6]. Due to incomplete lineage sorting and recombination, a substantial fraction of variants segregating in these archaic haplotypes is not expected to be of archaic origin [9]. Thus, to generate a high-confidence set of SNPs of archaic origin (aSNPs) we applied a series of filtering steps to our starting list of 7,791,042 genomic SNPs (see Methods and S1 Table). This allowed us to generate two high-confidence sets of 140,916 Denisovan and 88,625 Neanderthal aSNPs (S2 and S3 Tables, respectively). In parallel we defined a set of 229,541 (the sum of all high-confidence aSNPs) variants segregating in present-day Papuans outside archaic haplotypes and which likely arose after the Out of Africa event that matched the allele frequency distribution of our list of Denisovan and Neanderthal aSNPs (see Methods and S2(C) Fig) to act as a control set in all analyses. We refer to these variants as non-archaic SNPs (naSNPs, S4 Table).

Following calculation of the main introgressed allele frequency (MIAF) and of the derived allele frequency (DAF) for aSNPs and naSNPs respectively (see Methods), we found that 82,013 Denisova, 54,415 Neanderthal and 136,428 non-archaic variants are segregating at common-to-high-frequencies in PNG (MIAF or DAF > 0.05, see S1 Table). Importantly, we found that Denisovan and Neanderthal and non-archaic variants, whether segregating at low or common-to-high-frequencies, occur within genomic regions subjected to similar levels of background selection which are also higher (corresponding to less conserved regions) than those observed for protein-coding sequences (S3(A) and S3(B) Fig). We thus reasoned that, despite retaining only 10.5%, 9.2% and 4% of the original set of Denisovan, Neanderthal and non-archaic variants (S1 Table), our filtering criteria yielded a high-confidence set of archaic variants and an adequate corresponding background set of non-archaic alleles. Throughout the text and unless otherwise noted, to prioritise variants of potential evolutionary interest we focused on the above defined set of 82,013 Denisova, 54,415 Neanderthal and 136,428 non-archaic variants common-to-high-frequency variants.

Similar patterns of introgression for Denisovan and Neanderthal aSNPs across chromatin states
To gain insights into the potential consequences of aSNPs in present-day Papuans, we began by annotating the set of aSNPs and naSNPs across multiple genomic elements (see Methods). We found that SNPs lie within non-coding elements, regardless whether considering the entire set of variants or only those instances segregating at MIAF or DAF ≥ 0.05 in New Guinea Island (S4(A) and S4(B) Fig). As previously reported [8, 9, 14], we also noted a significant depletion of both Denisovan and Neanderthal aSNPs from protein coding sequences, relatively to naSNPs (n = 3,431 and 1,919, p = 1.41 ⋅ 10−11 and 6.61 ⋅ 10−19 for Denisovan and Neanderthal aSNPs, respectively).

To investigate the patterns of archaic introgression at a higher granularity, we intersected the set of common-to-high-frequency variants against epigenomic maps of 15 different chromatin states across 111 cell types generated by the Roadmap Epigenomics Project [34], and computed odds ratios (OR) between aSNPs and naSNPs to identify states containing an excess of Denisovan and/or Neanderthal variants. aSNPs are generally depleted from inactive chromatin states (states not associated with gene expression as defined by [34]) as well as promoters and transcription start sites, but enriched within transcribed and subsets of enhancer elements (Fig 1; full results are available at S5 Table). In particular, we observed a significant depletion of Denisovan and Neanderthal aSNPs from PolyComb repressed (ReprPC pDeni = 7.7 ⋅ 10−65), bivalent enhancers, promoters and transcription start sites (EnhBiv pDeni = 5.0 ⋅ 10−69, BivFlnk pDeni = 1.74 ⋅ 10−20 and TssBiv pDeni = 4.95 ⋅ 10−20, respectively). We also observed a Neanderthal-specific significant enrichment for heterochromatin (Het pNean = 2.9 ⋅ 10−18), whereas Denisovan aSNPs are significantly depleted from this state (pDeni = 4.6 ⋅ 10−6). An exception to this trend was an excess of aSNPs from both archaic ancestries that fell within quiescent chromatin (Quies pDeni = 1.5 ⋅ 10−15), as previously reported by [12]. This observation is not surprising given that quiescent chromatin is the predominant state in each epigenome [34], and this finding might indicate that a substantial fraction of aSNPs are functionally inert in most contexts, although silencers have also been found enriched within this state [49].

Admin
Administrator

Posts: 80,516

Denisovan Ancestry in Papua New Guinea (PNG) Dec 20, 2022 1:43:14 GMT

Quote

Post by Admin on Dec 20, 2022 1:43:14 GMT

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).

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

When considering active chromatin states, we noticed a significant depletion of Denisovan and Neanderthal ancestries from transcription start sites (TssA pDeni = 3.4 ⋅ 10−4), promoters (TssAFlnk pDeni = 0.026) and states marked by both promoter and enhancer epigenetic marks (TxFlnk pDeni = 0.0093; Fig 1) and of Neanderthal aSNPs from enhancers (Enh, pNean = 1.3 ⋅ 10−8). On the other hand, we observed a significant enrichment for Denisovan and Neanderthal aSNPs within genic enhancers (EnhG pDeni = 1.6 ⋅ 10−4) transcribed (Tx pDeni = 3.9 ⋅ 10−30) and within zinc-finger protein-genes and repetitive sequence (ZNF/Rpts pDeni = 0.033) states, this latter particularly for Neanderthal variants (Fig 1). Finally, we found a significant enrichment of Denisovan aSNPs within weakly transcribed regions (TxWk pDeni = 2.3 ⋅ 10−10) (Fig 1). Importantly, we did not find any substantial differences in these patterns when considering introgressed SNPs regardless of allele frequency (S5 Fig and S6 Table).

To validate this annotation approach we repeated our analyses using a set of 11,422,574 Neanderthal variants identified within 75 individuals of West Eurasian ancestry [31], which we processed as above. We found similar patterns of enrichment across chromatin states in the Papuan and West Eurasian datasets, with Neanderthal variants being significantly enriched within Quies, Het, Tx and TxWk states (S6(A) Fig). Similarly to [12], we also observed a significant depletion of Neanderthal variants from Enh and EnhBiv, although we detected a significant enrichment within EnhG (S6(A) Fig); however, the three states were collapsed into a single category in [12] and fewer SNPs were associated with this state than with Enh. These results not only validate our SNP filtering criteria; they also suggest that, on a genome-wide scale, different instances of admixture with archaic hominins might broadly contribute to modern human populations primarily through gene regulatory change.

Finally, to confirm whether these signals were caused by SNPs occurring within highly pleiotropic or cell-specific functional elements, we counted the number of different cell types across which aSNP and naSNP -containing elements were annotated with any given chromatin state. This yielded an estimate of the potential pleiotropic activity of each functional element. We observed similar degrees of pleiotropy across all the three ancestries, with variants within Quies, TssA, Tx and TxWk exhibiting higher pleiotropic activities than those within TssAFlnk, TxFlnk, Enh and EnhG (S7 Fig). Given that 64,944 Denisovan, 41,667 Neanderthal and 108,282 non-archaic variants fell within chromatin states associated with transcribed regions and CREs in at least one cell type (Fig 1 and S5 Fig), these results indicate that a non-trivial fraction of introgressed variants might have functional consequences, particularly for gene regulatory processes occurring in a limited number of cell types.

enisovan and Neanderthal aSNPs might have different effects on gene regulatory processes across tissues
Given the observed enrichment of aSNPs within chromatin states associated with CREs, we sought to investigate the effects of Denisovan and Neanderthal aSNPs on gene regulation. Hence, we restricted our analyses to all common-to-high-frequency SNPs within TssA, TssaFlnk, TxFlnk, Enh and EnhG chromatin states, which resulted in three sets of 39,636, 24,653 and 67,098 Denisovan, Neanderthal and non-archaic variants. Next, to identify tissues containing an over-representation of Denisovan and/or Neanderthal variants, for each tissue we again computed the OR between the number of aSNPs and naSNPs (see Methods). Despite roughly similar numbers of aSNPs falling within CREs across all tissue types, patterns of enrichment vary substantially by tissue type (Fig 2; for full results see S7 Table). For instance, Denisovan and Neanderthal aSNPs are significantly depleted from iPSCs (pDeni = 1.2 ⋅ 10−7, pNean = 9.6 ⋅ 10−5) and Neurospheres (pDeni = 0.02, pNean = 0.0016), but significantly enriched within mesenchymal (pDeni = 1.7 ⋅ 10−7, pNean = 3.3 ⋅ 10−5) and blood & immune T cells (pDeni = 3.5 ⋅ 10−15, pNean = 3.3 ⋅ 10−6). In addition, enrichment patterns vary by aSNP ancestry, with Neanderthal-specific enrichment within smooth muscle (pNean = 8.2 ⋅ 10−6) and Denisovan aSNPs significantly overrepresented within haematopoietic stem cells (HSCs) & B cells (pDeni = 5.0 ⋅ 10−4).