Neanderthals and humans interbred in Europe for 5,400 years

Figure 1
Schematic Describing the Method of Detecting Allele-Specific Expression Neanderthal-Introgressed Haplotypes

The sequencing of the Neanderthal genome revealed that approximately 2% of the ancestry of each non-African modern human traces to ancient gene flow from Neanderthals (Green et al., 2010, Prüfer et al., 2014). Recent studies have extended this observation by developing approaches to identify Neanderthal-inherited haplotypes in individual modern human genomes from globally diverse populations (Vernot and Akey, 2014, Sankararaman et al., 2014, Vernot et al., 2016). Many of these sequences encompass protein-coding genes and other functional genomic features, raising the intriguing possibility that they may influence variation in modern human traits. Lending support to this hypothesis, several introgressed haplotypes overlap annotated hits from published genome-wide association studies (Sankararaman et al., 2014). Furthermore, a recent analysis of electronic medical records presented evidence that Neanderthal alleles are associated with a range of clinical traits, including depression, actinic keratosis, hypercoagulation, and tobacco use (Simonti et al., 2016). While these studies suggest phenotypes that may be impacted by past hybridization, the functional mechanisms by which these associations arise remain poorly characterized. Regulatory variation influencing gene expression is a key source of phenotypic variation within and between species (King and Wilson, 1975). We thus sought to characterize the functional legacy of ancient gene flow by systematically investigating the contribution of Neanderthal-introgressed sequence to the landscape of modern human cis-regulatory variation.

A powerful approach for detecting and quantifying the impacts of cis-regulatory variation is to test for allelic differences in transcript abundance for individuals heterozygous at transcribed polymorphisms (Yan et al., 2002, Skelly et al., 2009; Figure 1A). By contrasting counts of reads supporting each allele within heterozygous individuals, this approach is less susceptible than expression quantitative trait locus (eQTL) mapping to batch effects and other confounding variables. While sophisticated methods have been developed to detect allele-specific expression (ASE) in a single sample (Skelly et al., 2009), across individuals (van de Geijn et al., 2015), or across tissues (Pirinen et al., 2015), no generalized approach has been developed to combine all information simultaneously. Integrating this information is complicated, however, by the complex structure of large-scale RNA sequencing (RNA-seq) datasets such as Genotype-Tissue Expression (GTEx) Project (GTEx Consortium, 2015), wherein many individuals are sampled across many tissues, alleles vary in frequency, and genes vary in expression level across individuals and tissues. These challenges require a flexible method that capitalizes on the wealth of information contained in such complex data structures. To this end, we developed a Bayesian generalized linear mixed model (GLMM) approach (Figure 1A; STAR Methods) to combine expression data to estimate allele-specific effects, augmenting statistical power by integrating information across multiple individuals and tissues (Figure 1B). Applying this method on a genome-wide scale to the GTEx dataset (214 individuals, 52 tissues) revealed abundant cis-regulatory effects of Neanderthal-introgressed sequences and evidence of tissue-specific variation in regulatory divergence.


Figure 2
Allelic Effect Estimates for Neanderthal-Introgressed Variants

Applying our statistical framework to the GTEx data identified 1,236 introgressed SNPs (24.5%) in 767 genes that showed significant ASE at a false discovery rate (FDR) of 10% (Figures 2 and S2). Introgressed SNPs showing significant ASE were significantly enriched for directionally concordant single-tissue eQTL identified in previous GTEx analyses (Fisher’s exact test: odds ratio [OR] = 2.51, 95% credible interval [CI] [2.06, 3.06], p < 1 × 10−10), as were non-introgressed SNPs showing significant ASE (Fisher’s exact test: OR = 2.40, 95% CI [2.36, 2.44], p < 1 × 10−10). The magnitude of enrichment was greatest at low minor allele frequencies, reflecting the fact that rare variants must have large effects to be called as significant eQTL and are thus more likely to show concordant ASE (Figure S3). Of introgressed SNPs showing concordant ASE and eQTL effects, 80% fall within an eGene whose lowest p value SNP has an alternative allele that matches an Altai Neanderthal allele, supporting the ability of our analysis to tag potential causal regulatory variants of Neanderthal-introgressed origin. Most of the remaining 20% are also expected to be Neanderthal in origin, with mismatches attributable to diversity within and between the introgressing Neanderthal and Altai Neanderthal populations.

We observed that the Neanderthal allele was upregulated for 49.8% and downregulated for 50.2% of SNPs showing significant ASE, indicating no overall directional bias (binomial test: 95% CI [0.469, 0.526], p = 0.887). Notable examples of variants exhibiting extreme ASE include linked SNPs (r2 = 1.0) rs73236617, rs3924112, rs5743557, and rs5743556 in TLR1, an innate immunity gene previously suggested to be a target of adaptive introgression (Dannemann et al., 2016, Gittelman et al., 2016, Quach et al., 2016, Nédélec et al., 2016, Deschamps et al., 2016). The Neanderthal haplotype is associated with significantly increased expression of TLR1 (rs5743557, binomial GLMM: β = 1.122, 95% CI [0.708, 1.563], p = 5.50 × 10−9, MAFEUR = 0.221), consistent with recent eQTL results from Dannemann et al. (2016)).


Figure 3
Introgression Introduced Functional Variants with Complex Regulatory Effects

Several additional introgressed variants exhibiting significant ASE were previously linked to human disease traits in published genome-wide association studies (GWAS), suggesting potential cis-regulatory mechanisms underlying these associations. In total, we identified eight Neanderthal regulatory variants associated with nine distinct phenotypes, including rs3765107 (binomial GLMM: β = −0.428, 95% CI [–0.475, −0.384], p < 1 × 10−10, MAFEUR = 0.111), which lies within the lysosomal transporter-encoding gene SLC15A4 and is associated with systemic lupus erythematosus (Table 1; Han et al., 2009). SLC15A4 is required for endosomal Toll-like receptor (TLR) signaling and secretion of proinflammatory cytokines by plasmacytoid dendritic cells (Blasius et al., 2010). This example thus adds to growing evidence that Neanderthal introgression contributed to risk of autoimmune disorders (Sankararaman et al., 2014) and innate immune response (Quach et al., 2016, Nédélec et al., 2016).

Disease-associated Neanderthal regulatory variants furthermore reveal how hybridization contributed to the genomic complexity of modern humans. These include rs950169, a SNP in extracellular matrix protein ADAMTSL3 that is significantly associated with both height (Table 1; Wood et al., 2014) and schizophrenia risk (Table 1; Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014). In the GTEx data, rs950169 shows tissue-wide downregulation of the Neanderthal-introgressed allele (binomial GLMM: β = −0.413, 95% CI [–0.463, −0.362], p < 1 × 10−10, MAFEUR = 0.273), though a SNP on the same haplotype (rs2135551, r2 = 1.0) shows nearly balanced expression of the two alleles (binomial GLMM: β = −0.020, 95% CI [−0.030, 0.073], p = 0.0151, MAFEUR = 0.273). This pattern is consistent with a model of splicing regulation proposed by Need et al. (2009)), which provides a detailed mechanistic explanation for the observed ASE. Specifically, the Neanderthal allele introduces a splice acceptor site in exon 30 (Figure 3), resulting in alternative splicing and truncation of the protease and lacunin (PLAC) domain encompassing rs950169.