Stephen Rong,Christopher R. Neil,Anastasia Welch,Chaorui Duan,Samantha Maguire,Ijeoma C. Meremikwu,Malcolm Meyerson,Ben J. Evans,William G. Fairbrother
Large-scale functional screen identifies genetic variants with splicing effects in modern and archaic humans
RNA splicing,modern humans,archaic humans,archaic introgression,massively parallel splicing assay
Significance Our study reveals the impact of genetic differences in pre-mRNA splicing between modern and archaic humans. We used a massively parallel experiment to identify 962 splicing variants that differ between modern and archaic humans, including potentially pathogenic splicing variants unique to Neanderthals and Denisovans, and introgressed splicing variants that may explain differences in modern human traits. Our findings indicate that purifying selection acted against splicing variants in modern human populations, whereas positive selection favored splicing variants in adaptive introgression. By distinguishing causal from linked variants, our study contributes to understanding the functional consequences of genetic variation within extant and extinct hominins and of introgressed variation in modern human populations.
Large-scale functional screen identifies genetic variants with splicing effects in modern and archaic humans

Our study reveals the impact of genetic differences in pre-mRNA splicing between modern and archaic humans. We used a massively parallel experiment to identify 962 splicing variants that differ between modern and archaic humans, including potentially pathogenic splicing variants unique to Neanderthals and Denisovans, and introgressed splicing variants that may explain differences in modern human traits. Our findings indicate that purifying selection acted against splicing variants in modern human populations, whereas positive selection favored splicing variants in adaptive introgression. By distinguishing causal from linked variants, our study contributes to understanding the functional consequences of genetic variation within extant and extinct hominins and of introgressed variation in modern human populations.

Modern humans are the only present-day hominins, but used to coexist with many others, including recently with our closest extinct archaic relatives, the Neanderthals and Denisovans. High-coverage ancient genomes sequenced from Neanderthal and Denisovan fossils have enabled the comparative study of genetic differences between modern and archaic humans (1–3). These genetic differences likely underlie many of the phenotypic differences between modern and archaic fossils, but also phenotypic differences not preserved in the fossil record (4, 5). Remarkably, genomic analyses have revealed multiple admixture events between modern and archaic humans (6–12), notably during modern human expansion out of Africa. These introgression events contributed about 2% Neanderthal ancestry to all non-African populations (13), and up to 5% Denisovan ancestry to populations from Melanesia (11, 12, 14, 15). Because archaic humans had smaller effective population sizes than those of modern humans, they may have accumulated more deleterious alleles that, once introgressed, were subject to purifying selection in modern populations (16–18). Recent methods to detect introgressed segments (7–11, 19, 20) and adaptively introgressed loci (21–23) demonstrate a lasting impact of archaic introgression on several human phenotypes (4, 24), including those involved in local adaptation to novel environments (7, 9, 25), diets (26), and pathogens (27–30). Enrichment analyses of protein-coding and cis-regulatory regions show that introgressed variation has functional consequences in different tissue and cell types (1, 18, 31–33). One common approach is to look for trait-associated introgressed variants in modern human populations, such as those associated with complex traits (4, 24, 34–37), gene expression across tissues (38), or in response to infection (39–41). However, this approach is limited to studying introgressed variants found in well-sampled populations, which are currently biased toward individuals of European ancestry (42), and are ill-suited for studying modern- or archaic-specific variants. Moreover, methods based on association tests are affected by allele frequency (AF) and linkage disequilibrium (LD) and often cannot distinguish causal variants from linked variants. Recent studies have advanced several ways around these limitations. These include prediction models to impute unobserved differences in gene expression and 3D genome organization (43, 44); CRISPR-Cas9 editing of archaic alleles into human-derived stem cells (45, 46); and massively parallel reporter assays (MPRAs) to identify causal variants with cis-regulatory impacts on transcription (47–49). In particular, MPRAs have been applied to study modern-specific and adaptively introgressed variants (50–53). Most of these studies focused on transcriptional regulation even though posttranscriptional processes, such as RNA splicing regulation, are also important to human phenotypes and disease (54–57). Notable exceptions include an adaptively introgressed splicing QTL (sQTL) for isoforms of the OAS1 gene involved in immune response (22, 28, 58, 59), a schizophrenia-associated introgressed variant that creates a novel splice acceptor site in the ADAMTSL3 gene (38), and a survey of Neanderthal introgressed sQTLs involved in immune response in individuals of European ancestry (41). Here, we analyzed a total of 5,169 modern, archaic, and introgressed variants in hominin evolution with a previously developed massively parallel splicing assay (MaPSy) (60). MaPSy is a parallel minigene reporter experimental system for identifying exonic splicing mutations (ESMs) that impact cis-elements involved in exon recognition during the process of pre-mRNA splicing. MaPSy and similar assays have been used to characterize splicing regulatory impacts of many types of variants, including pathogenic and rare variants (61–65). Like MPRAs for transcriptional regulation, MaPSy can distinguish causal from linked variants, and is not limited to variants observed in extant species. It is thus well suited for characterizing the splicing effects of genetic variation within extant and extinct hominins. We use MaPSy-identified ESMs in conjunction with predicted splicing mutations and enrichment of sQTLs to characterize the relative contribution of splicing effects in hominin evolution. We considered the possible effects of natural selection on these variants and provide examples of candidate adaptively introgressed, lineage-specific, and introgressed trait-associated ESMs.

We defined six sets of human evolution–related variants based on comparisons between the 1000 Genomes Project (1KGP) (66); the Genome Aggregation Database (gnomAD) (67); and four high-coverage archaic genomes, the Altai Neanderthal, the Vindija Neanderthal, the Chagyrskaya Neanderthal, and the Altai Denisovan (2, 15, 34, 68) (Fig. 1A). The first four variant sets correspond to fixed or nearly fixed variants with origins on different branches of the human evolutionary tree, which we refer to as X-specific variants for branch X. To polarize reference (REF) or alternate (ALT) alleles as derived (DER) or ancestral alleles (ANC), we used Ensembl-inferred human ancestor alleles (69). “Modern-specific” variants originated on the evolutionary branch leading to all modern humans, defined as having high derived allele frequency (DAF) in modern humans (1KGP and gnomAD global DAF >0.9) (we focused on high-frequency DER variants on the modern human lineage instead of strictly fixed variants to allow for the effect of subsequent archaic introgression) and having being unobserved in archaic humans (0/8 archaic alleles) (to avoid incomplete lineage sorting variants). “Archaic-specific” variants originated on the branch between the split between modern and archaic humans but preceding the split between Neanderthals and Denisovans, defined as having high DAF in archaic humans (≥7 of 8 archaic alleles) and being nearly absent in modern African (AFR) populations (1KGP and gnomAD AFR DAF < 0.01). “Neanderthal-specific” variants originated on the branch leading to Neanderthals after the split from Denisovans (≥5 of 6 Neanderthal alleles, 0 of 2 Denisovan alleles, 1KGP and gnomAD AFR DAF < 0.01). “Denisovan-specific” variants originated on the branch leading to Denisovans after the split from Neanderthals (2 of 2 Denisovan alleles, 0 of 6 Neanderthal alleles, 1KGP and gnomAD AFR DAF < 0.01). (Materials and Methods and SI Appendix, Fig. S3 for a complete flowchart of variant set definitions.) As only one genome was available, Denisovan-specific changes were not used outside the task of identifying example Denisovan-specific ESMs because there was less power to identify fixed or nearly fixed differences. Next, we defined an introgressed variant set, which contains previously published variants linked to archaically introgressed haplotypes, which may be either DER or ANC. These come from either Vernot et al. (8) or Browning et al. (11) and were previously identified in 1KGP populations or Papa New Guinea. Finally, we defined an adaptively introgressed variant set, which contains previously published variants linked to high-frequency adaptively introgressed. These come from either adaptively introgressed haplotypes from Gittelman et al. (22), or introgressed haplotypes from Browning et al. (11) that overlap high-frequency adaptively introgressed loci from Racimo et al. (23) in matching 1KGP populations. Some variants are in multiple variant sets; for example, archaic DER variants can also be introgressed variants due to subsequent introgression into modern humans (SI Appendix, Fig. S4 A and B). A modern-specific variant can also overlap an introgressed variant in the sense that the DER allele became nearly fixed in modern humans, but archaic introgression later reintroduced the ANC into some modern human populations (51). Overall, we identified 1,872,673 variants belonging to one of the human evolution variant sets [available as additional data at https://doi.org/10.5281/zenodo.7158564 (70)], of which all 5,224 variants overlapping internal exons of protein-coding genes that could be incorporated into the MaPSy experiment (limited to testing variants in exons no more than 120 bp in length or exonic variants within 90 bp of a 3′ ss) were synthesized (Fig. 1B and Dataset S1).

MaPSy uses a high-throughput minigene reporter experiment to identify ESMs with cis-regulatory effects on exon inclusion (Fig. 1C) (60). Minigene reporters have a long history of use in determining whether an allelic variant disrupts splicing efficiency (71), and the original MaPSy paper validated the high-throughput approach using endogenous RNA splicing in ENCODE cell lines and patient tissue samples (60). In MaPSy, minigene reporters containing a strong promoter, shared flanking exons, and a variable middle exon containing either DER or ANC alleles are transfected in triplicate into HEK 293T cells, where they undergo transcription followed by RNA splicing, resulting in mRNA transcripts with either exon inclusion or skipping events. MaPSy then estimates a splicing score for each variant as the log2 fold change (FC) between DER/ANC sequencing read counts in the output exon inclusion cDNA/input DNA (SI Appendix, Extended Materials and Methods). We applied MaPSy to 5,224 total human evolution variants, of which 5,169 (98.9%) passed quality control filters for minimum read counts (Dataset S1 and SI Appendix, Table S1). We identified 962 ESMs (19% of assayed variants) that were significant at an FDR threshold of < 0.05 and an MaPSy functional score threshold of |FC| > 1.5. This included 71 (of 322) modern-specific; 51 (of 315) specific to the most recent common ancestor of Neanderthals and Denisovans, hereafter archaic-specific; 169 (of 909) Neanderthal-specific; 229 (of 1,312) Denisovan-specific; 587 (of 3,092) introgressed; and 16 (of 54) adaptively introgressed ESMs (Fig. 1 D–I). We further categorized ESMs as either promoting exon inclusion (positive ESMs) or inhibiting exon inclusion (negative ESMs), and as either strong ESMs (|log2FC| > log2(3)) or weak ESMs (log2(3) > |log2FC| > log2(1.5)). A total of 677 identified ESMs (13%) have weak effects on splicing (weak ESMs) and 285 ESMs (5.5%) have strong effects (strong ESMs).

Many ESMs were found in splice region positions closest to the 3′ ss and 5′ ss, particularly ESMs with the strongest negative effects on splicing (strong negative ESMs) that were expected to result in loss of canonical transcript exon inclusion (72) (SI Appendix, Fig. S5A). The vast majority of ESMs, however, were synonymous and missense mutations found deeper into the exon with relatively moderate positive or negative effects on exon inclusion (weak ESMs) (SI Appendix, Fig. S5B). To better understand the mechanisms by which ESMs affect splicing, we compared MaPSy functional scores to computational predictions of different splicing effects. First, we used SpliceAI, a deep neural network model trained to predict the location of splice sites from RNA sequence, which can be used to identify mutations that disrupt splice sites (73). SpliceAI was trained using the modern human reference genome, but since the core splicing machinery is highly evolutionarily conserved across vertebrates (74, 75), we expect it should be able to predict cis-regulatory splicing variants in archaic genomes as well. Similar models have recently been used to explain splice site differences between humans and chimpanzees, which are many times more diverged than those of modern and archaic humans (76). We used the max of raw SpliceAI scores for predicted donor gain, donor loss, acceptor again, and acceptor loss (SpliceAI max score). SpliceAI max scores tend to increase with MaPSy scores away from zero, and the variants with the highest SpliceAI scores tend to also have the most negative (but not positive) MaPSy scores (Fig. 2A). Next, we used Delta EI score, which is the change in the ability of overlapping hexamers at an SNV to function as exonic splicing enhancers (77). SNVs with positive/negative Delta EI are predicted to increase/decrease exon inclusion. Consistent with this, we find that Delta EI scores increase with MaPSy scores, except when MaPSy scores are extremely negative (Fig. 2A), or when variants are close to a splice site (within 5 bp) (SI Appendix, Fig. S6A). We find similar trends for the change in two alternative hexamer scores of exonic splicing enhancer activity (78) (SI Appendix, Fig. S6 B and C). We cloned and validated using RT-PCR a handful of strong negative (n = 3), weak negative (n = 2), not significant (n = 2), weak positive (n = 2), and strong positive (n = 2) MaPSy variants (SI Appendix, Fig. S8 and Table S2). All the three strong negative ESMs disrupted inclusion of the wild-type exon. Both not-significant variants showed no disruption of splicing. Half of the remaining variants showed evidence of splice disruption consistent with expected direction. In total, we were able to validate the MaPSy ESM category of 8 of the 11 tested variants with RT-PCR. Overall, these results suggest that MaPSy is sensitive at detecting ESMs acting through distinct splicing mechanisms. First, those with strong negative effects (strong negative ESMs) on splicing tend to be found in extended splice regions, affect splicing through the mechanism of splice site disruption, possibly through creation of cryptic splice sites, and are thus more likely to be deleterious. Second, those with weaker effects on splicing (weak ESMs) are enriched in the interior of the exon and affect splicing through changes in exonic splicing regulatory motifs.

Purifying selection is expected to be more efficient in modern humans as compared to Neanderthals and Denisovans owing to the larger effective population size of modern humans. We therefore reasoned that purifying selection might affect the distribution and function of ESMs in different hominin populations (16–18). We tested this hypothesis by comparing proportions of fixed or nearly fixed strong negative MaPSy ESMs and strong SpliceAI scores specific to the modern, archaic, and Neanderthal branches. We focused on these variants because strong negative, but not positive, MaPSy ESMs and strong SpliceAI scores are in general depleted in evolutionarily constrained genes (SI Appendix, Fig. S9) as defined by low loss-of-function observed/expected upper bound fraction (LOEUF) in gnomAD (79). Modern-specific variants were found to be depleted in strong negative ESMs compared to Neanderthal-specific variants (P = 0.036) (Fig. 2B), but not in strong positive ESMs (SI Appendix, Fig. S10A). This is consistent with stronger purifying selection in modern versus archaic humans. They were also nominally enriched in weak ESMs (P = 0.055) (Fig. 2C). Considering a range of SpliceAI scores of increasing severity of splicing effect from weak (0.01 to 0.2), to moderate (0.2 to 0.5), to strong (0.5 to 1.0), modern-specific variants were significantly depleted in strong SpliceAI scores compared to both archaic-specific (P = 0.0088) and Neanderthal-specific variants (P = 0.041) (Fig. 2F), concordant with the MaPSy results. However, all the three branches had similar proportions of weak and moderate SpliceAI scores (Fig. 2 D and E). Together, the MaPSy and SpliceAI results support the hypothesis that modern humans are depleted in ESMs strongly associated with deleterious exon skipping as compared to Neanderthals. This is consistent with expectations based on their relative effective population sizes. Concordant with this, we also observed a similar depletion of missense variants in the modern lineage compared to the archaic (P = 0.012) and Neanderthal (P = 0.0036) lineages, but not of synonymous (as expected) nor nonsense and splice region variants (likely due to low statistical power) (SI Appendix, Fig. S11). No differences in the proportions of weak, moderate, or strong SpliceAI variants were found across modern human populations in the 1KGP (based on variants with MAF ≥ 0.01, SI Appendix, Fig. S12), which is consistent with 1KGP modern Africans and non-Africans having much higher heterozygosity than that of archaic humans (2, 15, 34). We examined whether our results are robust to our choice of AF thresholds for defining modern, archaic, and Neanderthal-specific variants since the sample size for modern humans is orders of magnitude larger than that for archaic humans (SI Appendix, Fig. S13). Because there are only three high-coverage Neanderthal genomes and only four high-coverage archaic genomes in total, we reasoned that even if an allele is observed eight times in eight archaic haplotypes or six times in six Neanderthal haplotypes, it can at best be inferred to be likely at high frequency in the ANC archaic or Neanderthal populations, respectively. Thus, we considered alternate definitions of modern, archaic, and Neanderthal-specific variants on which to subsample the data based on different high-frequency AF thresholds. For archaic-specific, we considered a stricter alternate threshold of AF = 8/8. For Neanderthal-specific, we considered a stricter alternate threshold of AF = 6/6. For modern-specific, we considered a stricter alternate threshold of AF > 0.99 for MaPSy, and varying alternate thresholds from AF > 0.5 to AF > 0.99 for SpliceAI. We still observed statistically significant depletion or nominal depletion of strong negative MaPSy ESMs in modern humans when compared to Neanderthals (SI Appendix, Fig. S13A) and of strong SpliceAI variants in modern humans when compared to both archaic humans and Neanderthals (SI Appendix, Fig. S13F) and that this result did not depend strongly on our choice of AF cutoff when defining lineage-specific variants.

Introgressed variants in modern human populations are known to be depleted in genes, conserved coding regions, and noncoding regulatory regions, which could be due to purifying selection against deleterious introgressed alleles (9, 18, 31–33). However, high-frequency adaptively introgressed variants are known to be enriched in genes and noncoding regulatory regions, which is likely due to positive selection (31, 33). Supporting this, many adaptively introgressed haplotypes have been shown to affect gene expression (22, 28, 29, 39, 40). We tested the hypothesis that adaptively introgressed variants are enriched in splicing effects compared to all introgressed variants using both MaPSy scores and SpliceAI scores. Adaptively introgressed variants were found to be enriched in weak ESMs compared to all introgressed variants (P = 0.014) (Fig. 3A), but not in strong negative or strong positive ESMs (SI Appendix, Fig. S10 B and C). We considered a range of SpliceAI scores from weak (0.01 to 0.2), to moderate (0.2 to 0.5), to strong (0.5 to 1.0). Adaptively introgressed variants were found to be significantly enriched in moderate SpliceAI scores compared to all introgressed variants (P = 0.043) (SI Appendix, Fig. S14B), but not weak or strong SpliceAI scores (SI Appendix, Fig. S14 A and C), which is concordant with the MaPSy results. Together, these findings suggest that adaptively introgressed variants are enriched in splicing variants that may be involved in alternative splicing events compared to all introgressed variants.

The pattern of ESM enrichment detailed above could stem from positive selection favoring beneficial splicing variants or from purifying selection disfavoring deleterious splicing variants. We thus sought to independently evaluate enrichment of splicing variants among the introgressed or adaptively introgressed variant sets relative to matched nonintrogressed control variants. For splicing variants, we used Genotype Tissue-Expression (GTEx) sQTLs and variants with weak (0.01 to 0.2), moderate (0.2 to 0.5), and strong (0.5 to 1) SpliceAI scores. Since archaic introgression itself generates LD, introgressed variants are more likely to tag sQTLs, and thus it is necessary to use an AF- and LD-matched background of nonintrogressed control variants to assess enrichment of sQTLs (SI Appendix, Fig. S15 A and B). Compared to nonintrogressed variants, introgressed variants have overall lower AF and higher LD, whereas adaptively introgressed variants have similar AF but higher LD (Fig. 2 B and C). We thus measured enrichment relative to randomly sampled, nonintrogressed variants in matched AF and LD bins. This approach showed that we could recapitulate known intergenic enrichment and missense depletion among introgressed variants and genic enrichment among adaptively introgressed variants (SI Appendix, Fig. S15 C and D). On analyzing QTL enrichment, we found that introgressed variants were significantly depleted in GTEx-based sQTL variants (those with FDR < 0.05 in at least one tissue) (P < 0.001), adaptively introgressed variants were enriched in sQTLs (P = 0.04), and there was a significant difference between the two sets (P < 0.001) (Fig. 3D). We found a similar pattern for enrichment of overlap with GTEx-based expression QTLs (eQTLs) (Fig. 3E). We also found that introgressed variants were most depleted in strong SpliceAI scores (P = 0.014) (Fig. 3H). Adaptively introgressed variants were enriched in weak (P < 0.001) and moderate SpliceAI scores (P = 0.016), but not strong SpliceAI scores (P = 0.36) (Fig. 3 F and G). Overall, these ESM, sQTL, and SpliceAI enrichment analyses suggest introgressed variants are depleted in variants that result in loss of exon inclusion, which is consistent with the action of purifying selection postintrogression. In contrast, adaptively introgressed variants are enriched in moderate-effect splicing variants that change splicing relative to the hominin ancestor, suggesting the action of positive selection postintrogression. Possible explanations for why we did not observe an enrichment of strong effect splicing variants in adaptive introgression are provided below.

As discussed above, LD often limits the ability to identify an individual variant’s contribution to a QTL, but MaPSy provides a functional test of the splicing phenotype of an individual nucleotide variant. To better understand the contribution of introgressed variants to an sQTL, MaPSy hits were compared to sQTL blocks and predictive measures of splicing effect such as Splice AI prediction. We showed above that introgressed variants overall were depleted in strong negative ESMs, strong SpliceAI scores, and sQTLs, suggesting a role of purifying selection against splicing disruptions. However, the subset of adaptively introgressed variants showed enrichment of moderate splicing effects. One such example is an adaptive Neanderthal-introgressed ESM, rs5743566, in the Toll-like receptor gene TLR1 involved in innate immunity (Fig. 4). The SNV rs5743566 belongs to an adaptive Neanderthal-introgressed haplotype spanning the TLR6/TLR1/TLR10 locus, a hotspot of repeated introgression from both Neanderthals and Denisovans (29). The introgressed haplotypes at this locus have previously shown to have eQTL effects on each of TLR1, TLR6, and TLR10 (22, 29, 38). Moreover, rs5743566 is in near-perfect LD (D′ = 1.0, R2 = 0.9959 across all 1KGP populations) with an intronic variant rs5743593 not predicted by SpliceAI to affect splicing but was previously identified by Rotival et al. (41) to be a Neanderthal introgressed TLR1 sQTL in stimulated monocytes (73, 80). The SNV rs5743566 resides in noncoding exon 2 of TLR1 in the 5′ UTR region 15 bp away from the nearby 5′ ss and has a MaPSy functional score of −0.81, suggesting that the DER G mutation moderately reduces inclusion of the exon compared to the ANC. It has a SpliceAI donor loss score of 0.38 affecting the nearby 5′ ss, which is the highest SpliceAI score of any variant on the adaptively introgressed haplotype (Fig. 4A). An analysis of GTEx sQTLs shows that rs5743566 has strong association with TLR1 splicing phenotype in seven different tissues and resides in a large sQTL block containing dozens of other associated variants in LD (Fig. 4A). Comparisons of genotypes in three Neanderthal and one Denisovan genomes confirm that rs5743566 likely originated in and was nearly fixed in the Neanderthal population (Fig. 4B) and was later adaptively introgressed into modern human populations where it is now found at intermediate frequencies in non-African populations in 1KGP (Fig. 4C). Finally, a motif-based search for RNA-binding protein (RBP) motifs using the ATtRACT database (81) suggests that the C>G mRNA change may modify PTBP1 and ZFP36 binding motifs or create a UG-rich binding motif common to CELF proteins, a family of RBPs known to regulate alternative splicing (82) (Fig. 4D). Moreover, a Delta EI score −0.22 for the C>G mRNA change suggests a moderate loss of exonic splicing enhancer activity. This example evidences the ability of MaPSy to identify causal splicing variants in an adaptively introgressed sQTL block.

We next considered as case studies the strongest Neanderthal-specific, archaic-specific, and modern-specific ESMs that were identified by their MaPSy functional scores (Dataset S2). The strongest Neanderthal-specific ESM that was also introgressed was rs12737091, a synonymous variant in the gene HSPG2. HSPG2 encodes the protein perlecan, a key component of the extracellular matrix that has diverse roles, including contributing to cartilage and bone development and angiogenesis (83). This variant has an MaPSy function score of −4.95 and a SpliceAI donor gain score of 0.98, which is consistent with the C>A mRNA change, creating a GUAAG canonical splice donor motif, and is an sQTL for HSPG2 in thyroid, adipose (subcutaneous), and artery (tibial) tissues (Fig. 5A). The RT-PCR validation confirmed that rs12737091 disrupts inclusion of the wild-type HSPG2 exon in favor of a truncated exon (SI Appendix, Fig. S8). The DER allele is fixed in all the three Neanderthal genomes and is found at low frequencies in non-African modern humans (Fig. 5 B and C). The strongest archaic-specific ESM that was also introgressed was rs17518147 in the splice region of exon 2 of CCDC181, a gene that encodes a microtubule-binding protein in sperm (SI Appendix, Fig. S16) (84). It has an MaPSy functional score of −6.53 and a SpliceAI donor loss score of 0.59, suggesting that the G>A mRNA change negatively disrupts the adjacent splice donor site (SI Appendix, Fig. S16A). The DER allele is fixed across all the four archaic genomes, suggesting an origin prior to the split between Neanderthals and Denisovans (SI Appendix, Fig. S16B), and was later introgressed from Neanderthals into modern humans where it is now found at low-to-intermediate frequencies across non-African populations (SI Appendix, Fig. S16C). It is also identified in GTEx as a CCDC181 sQTL in testis tissue. The strongest modern-specific ESM was rs77724956 in the splice region of exon 12 of FAAH, a gene that encodes a catabolic enzyme for fatty acid amides such as anandamide involved in pain sensation (SI Appendix, Fig. S17) (85). It has an MaPSy functional score of −5.59 and a SpliceAI acceptor loss score of 0.28, suggesting that the G>U mRNA change results in a negative disruption of the adjacent splice acceptor (SI Appendix, Fig. S17A). The DER T allele is not found in any of the four archaic genomes and is nearly fixed in modern humans worldwide (SI Appendix, Fig. S17B). The ancestral G allele is found at a very low frequency in modern humans, with the highest population frequency of 1 to 2% in East and South Asia (SI Appendix, Fig. S17C). Notably, rs77724956 is also an introgressed variant, suggesting that the functional ANC was reintroduced into modern human populations (51). It is nearly absent in European populations, and thus has no associated GTEx sQTL results. However, the eQTLGen consortium, which provides whole-blood eQTLs based on 31,684 individuals across multiple ancestries, identifies rs77724956 as a strongly associated FAAH eQTL (86). The previous case studies of top lineage-specific ESMs were all introgressed. We further sought to identify a subset of high-confidence lineage-specific splicing variants that are private to either modern or archaic humans. We included all variants identified as MaPSy ESMs with a SpliceAI max raw score ≥ 0.1 that are not in the introgressed variant set and are not polymorphic in either 1KGP or gnomAD genomes. This resulted in 36 high-confidence private splicing variants, of which 51 were Denisovan specific, 22 were Neanderthal specific, three were archaic specific, and none were modern specific (Dataset S3). Using more stringent splicing variant criteria of MaPSy strong ESMs with a SpliceAI max raw score ≥0.2 resulted in nine private splicing variants, of which eight were Denisovan specific and one was Neanderthal specific (SI Appendix, Fig. S18). Variants in the more stringent set include i) the Denisovan-specific splice region variant chr8:24342791 G/A in the gene ADAM7, which is important for sperm maturation and predicted by Combined Annotation Dependent Depletion to be deleterious (CADD score = 32) (87, 88); ii) a Denisovan-specific splice region variant chrX:65247366 T/C in the gene VSIG4, which is a macrophage complement receptor and T cell negative regulator (CADD score = 19.8) (89, 90), and which lies in the middle of a previously reported introgression desert (35); and iii) a Neanderthal-specific missense variant chr2:26702482 T/C in the gene OTOF which is involved in nonsyndromic recessive deafness (CADD score = 11.1) (91).

We sought to identify potential phenotypic impacts of causal splicing variants by exploring whether they were associated with trait variation. Most variants identified by genome-wide association studies are not causal due to LD. Statistical fine-mapping methods identify likely causal variants within a given phenotype-associated genomic region (92, 93). We intersected MaPSy ESMs with variants with posterior inclusion probability (PIP) ≥ 0.1 from a recent fine-mapping study across 94 traits in the UK Biobank (UKB) and 79 traits in BioBank Japan (BBJ) (94). This identified seven variants across and eight phenotypes with evidence of both splicing and phenotypic impacts (Dataset S4). However, six of these ESMs were also missense variants, suggesting alternative mechanisms of influencing phenotype. We identified one synonymous introgressed ESM, rs11045818 in the gene SLCO1B1, which has a moderately high PIP = 0.626 for total bilirubin (SI Appendix, Fig. S19A). SLOC1B1 is primarily expressed in the liver and encodes the protein organic anion transporting polypeptide 1B1 (OATP1B1), which transports bilirubin and other compounds into the liver (95). The variant has an MaPSy functional score of −2.05 and a SpliceAI acceptor loss score of 0.09. The DER allele is fixed across all the three Neanderthal genomes, and it narrowly misses our criteria for being identified as a Neanderthal-specific variant (1KGP and gnomAD AFR AFs are slightly above 1%) (SI Appendix, Fig. S19 B and C). A search for RBP-binding motifs using the ATtRACT database (81) suggests that the G>A mRNA mutation may introduce potential binding sites for splicing factors SRSF2 and NOVA1/NOVA2 (SI Appendix, Fig. S19D), though only SRSF2 is highly expressed in liver tissue (96). The variant resides in a credible set with two intronic variants, rs17329885 (PIP = 0.177) and rs11045820 (PIP = 0.197), on the same introgressed haplotype. While rs17329885 is not predicted by SpliceAI to affect splicing, rs11045820 is predicted to have a potential small effect on splicing (acceptor loss score of 0.08), and thus may be a potential candidate causal variant in the credible set. Alternative splicing of SLCO1B1 is not reported in GTEx V8 adult liver tissue (96), but has been reported in pediatric liver tissue (97) and in GTEx V9 long-read sequenced novel transcripts (98). We also intersected variants with SpliceAI max raw score ≥ 0.1 in the human evolution variant sets with fine-mapped variants with PIP ≥ 0.1. This identified an additional nine intronic introgressed variants, including one splice region variant, across eight distinct genes and 12 phenotypes (SI Appendix, Fig. S20 and Dataset S5). This included one adaptively introgressed intronic variant, rs17730088 in the keratin gene KRT71 fine-mapped to balding type four (PIP = 0.147), and two high PIP (≥0.5) introgressed intronic variants, rs2272783 in FECH fine-mapped to mean corpus hemoglobin (PIP = 0.989) and rs8192327 in SFTPC fine-mapped to FEV/FVC ratio (PIP = 0.875).

Key to our understanding of differences between modern and archaic humans is an appraisal of the functional effects of noncoding variants that existed in the evolutionary past, but which are no longer extant in present-day populations (5, 6, 41, 99, 100). Most genomic approaches developed so far have focused on transcriptional (31–33, 39, 40, 43) rather than posttranscriptional variant effects (but see refs. 41 and 32), are limited to studying introgressed variants associated with phenotypes in well-sampled modern human populations (4, 24, 35, 36, 38), and do not distinguish causal variants from linked variants. Recently, MPRAs have been used to address the latter two limitations, as shown in recent applications to detecting modern-specific and adaptively introgressed variants that affect cis-regulatory elements (50–53). In this study, we applied a MaPSy (60) to detect 962 ESMs out of a total of 5,169 tested variants, including nearly all modern-specific, archaic-specific, Neanderthal-specific, archaically introgressed, and adaptively introgressed variants that could fit in the assay (Fig. 1). We found MaPSy to be sensitive at detecting ESMs affecting both splice sites and exonic splicing regulatory elements (Fig. 2A and SI Appendix, Figs. S5 and S6). To better understand the evolutionary landscape of splicing effects across the different variant sets in human evolution, we used a combination of insights from MaPSy functional scores, SpliceAI predicted splicing effects, and enrichment of GTEx sQTLs. We found that modern humans are depleted in loss of exon variants compared to the Neanderthal lineage (Fig. 2 B–F), which is consistent with stronger purifying selection acting in the modern lineage due to their larger long-term effective population sizes (2, 16, 17). We next found that adaptively introgressed variants were enriched in moderate effect splicing variants relative to all introgressed variants (Fig. 3A). Analysis of enrichment or depletion of GTEx sQTLs and SpliceAI scores (relative to nonintrogressed control variants with matching AF and LD) showed that this was likely due to both increased purifying selection acting against introgressed variants that result in loss of exon inclusion and positive selection in favor of moderate-effect, but not weak or strong effect, splicing variants that altered exon inclusion frequencies in modern humans (both positively and negatively) relative to the hominin ancestor (Fig. 3 B–H and SI Appendix, Fig. S14). This was somewhat unexpected and suggests the hypothesis that adaptive introgression may favor functional alleles with intermediate effect sizes. Weak-to-moderate splicing variants are more likely to contribute to tissue-specific alternative splicing (73). Selection coefficients may be more likely to be positive if variants only disrupt splicing in some tissues, and more likely to be negative if variants strongly disrupt splicing throughout the body. Finally, we provide evidence that MaPSy can identify candidate causal splicing variants. As an illustrative example, we characterized an adaptively introgressed ESM, rs5743566, in the innate immunity gene TLR1 (Fig. 4). The TLR6/TLR1/TLR10 locus is a hotspot for repeated archaic introgression due to an important role in recognizing pathogens (29). Introgressed haplotypes at this locus have been shown to affect both TLR1 gene expression and splicing (22, 29, 38–41). The MaPSy-identified variant, rs5743566, lies in an internal exon in the 5′ UTR of TLR1 away from the splice regions. It is thus a seemingly unlikely candidate for splicing effect annotation. However, MaPSy functional scores, SpliceAI scores, predicted changes in RBP binding motifs, and predicted changes in exonic splicing enhancer activity, all supported the conclusion that rs5743566 is a causal splicing variant with a moderate negative effect on TLR1 exon two inclusion (Fig. 3). Moreover, GTEx identifies rs5743566 as a significant sQTL in seven out of 49 tissues, which could be explained by tissue-specific alternative splicing (96). Interestingly, Jagoda et al. (53) recently found another putatively causal variant, rs73236616, at this adaptively introgressed locus in a cis-regulatory element that impacts gene expression not in TLR1, TLR6, or TLR10, but the neighboring FAM114A1 gene. Thus, multiple causal variants affecting different genes may be found on the same adaptively introgressed haplotype, any of which could be a target of positive selection. We further showed that MaPSy can identify causal splicing variants that are fixed or nearly fixed on the modern, archaic, and Neanderthal lineages (Fig. 5 and SI Appendix, Figs. S14–S16), such as the DER Neanderthal-specific ESM, rs12737091, in the gene HSPG2 that encodes the ubiquitous extracellular matrix protein perlecan, and the DER modern-specific ESM, rs77724956, in the gene FAAH that is also an example of Neanderthal introgression reintroducing an ANC functional allele back into modern populations at low frequencies (SI Appendix, Fig. S17) (51). This also includes a set of high-confidence splicing variants private to Neanderthals or Denisovans that are not found in 1KGP or gnomAD genomes and would be unlikely to have been discovered even with deep resampling of transcriptomes across diverse modern human populations (SI Appendix, Fig. S18). These include variants that strongly disrupt splicing and are predicted to be highly deleterious by CADD scores in genes related to sperm maturation, ADAM7, and immune response, VSIG4. Finally, we identified MaPSy- and SpliceAI-identified candidate splicing variants that are also candidate causal variants for trait associations in the UKB and BBJ (94), including variants underlying associations with total bilirubin, balding type 4, mean corpus hemoglobin, and FEV/FVC ratio (SI Appendix, Figs. S19 and S20). There are notable limitations of MaPSy that apply to our study. While our findings are suggestive of evolutionary phenomena in ancient hominids, low statistical power of assayable lineage-specific variants may have hindered more robust evaluation. We validated a half exon approach for testing exonic variants near the 3′ ss of longer exons (SI Appendix, Fig. S2), but there were nearly 5,500 exonic variants that we were unable to test using MaPSy. Intronic variants, which include those that disrupt canonical splice signals, such as the 3′ ss, 5′ ss, polypyrimidine tract, and branchpoint signal, and over 11,000 proximal intronic variants (within 50 bp of a splice site), were also not tested using MaPSy. These could be tested using newer splicing assay designs based on barcoding or fluorescence-activated cell sorting (63, 64). MaPSy also shares limitations with other MPRA-like approaches. The minigenes used by MaPSy only include a small part of the endogenous sequence and thus lack important sequence and epigenetic information present in vivo that are further away and that also help determine splicing (73), such as potential competing splice sites of neighboring exons (101). Another shared limitation is that MaPSy is limited to testing a single cell type at a time, and here we only transfected HEK 293T. This is less of a problem for MaPSy than for MPRAs that test cis-regulatory elements affecting gene expression, because splicing variation, especially at splice sites, is shared across tissues more than expression variation (96, 102). The notable exception is brain-related tissues, which have a distinct splicing program compared to other tissues (102), and express tissue-specific splicing factors that regulate alternative splicing, such as NOVA1 (103). Neuronal cell lines and primary cells would thus be a high priority for future MaPSy experiments. As well, MaPSy only evaluates individual variants with cis-effects on splicing. Trans-effects on splicing are also important, and notably, all modern humans carry a fixed DER mutation in NOVA1 that has been shown to alter transcriptomes and neurodevelopment compared to the ANC found in archaic humans (46). It is also conceivable that multiple variants in cis could have synergistic effects on splicing (61, 104); this aspect of splicing is also not interrogated by MaPSy, and around 3% of the MaPSy tested exons contained multiple variants in our study. Our use of MaPSy provides unique insights into the role of natural selection acting on splicing variation in hominin evolution and enabled the identification of causal splicing variants, not limited to those in well-sampled present-day populations. Approaches such as MaPSy thus open the door to functionally characterizing variants found in diverse, understudied, and extinct populations. MaPSy is but one of the many recently developed multiplexed assays of variant effect (105). As these assays and other emergent methods for interpreting variant-to-function improve, such as gene editing and machine learning models, we will increasingly be able to apply them to answer questions about the origins of our species’ unique molecular and complex phenotypes.

Individual variant sets were described in the main text in the section “Defining Human Evolution Variant Sets Used in the MaPSy Experiment” and also in detail in SI Appendix, Fig. S3. All analyses used the hg19 reference. For lineage-specific variants, we consider the cases where the human reference (REF) allele is either the ANC or DER allele separately, based on inferred ANC alleles from Ensembl (69) from https://ftp.ensembl.org/pub/release75/fasta/ancestral_alleles/homo_sapiens_ancestor_GRCh37_e71.tar.bz2, and ignoring sites with ambiguous ANC alleles. These were then filtered based on 1KGP phase 3 global or African AF from http://hgdownload.cse.ucsc.edu/gbdb/hg19/1000Genomes/phase3/ and gnomAD (v2.1.1 genomes) global or African AF from https://gnomad.broadinstitute.org/downloads (67). Finally, these were then filtered by archaic allele counts in the four high-coverage archaic genomes from http://cdna.eva.mpg.de/neandertal/, and by the recommended filters for each archaic genome. 1KGP and gnomAD AF information could be missing due to either a true absence of the ALT allele in the population or because they were not genotyped. If AF is missing in one of 1KGP or gnomAD but not the other, we assumed that it was not genotyped in the dataset where it is missing. If AF is missing in both of 1KGP or gnomAD, we assumed that it was not genotyped or truly absent depending on whether it was found at low or high AF across the four high-coverage archaic genomes. Introgressed variants were downloaded from https://akeylab.princeton.edu/downloads.html (8) and https://data.mendeley.com/datasets/y7hyt83vxr/1 (11). Adaptively introgressed variants were downloaded from https://doi.org/10.1016/j.cub.2016.10.041 (SI Appendix, Table S1 of ref. 22) and Racimo https://doi.org/10.1093/molbev/msw216 (SI Appendix, Table S3 of ref. 23). This resulted in 1,872,673 total variants across variant sets available as additional data at https://doi.org/10.5281/zenodo.7158564 (70). Variant set intersections were plotted using UpSetR (106).

GENCODE v32 basic transcripts with liftOver to hg19 were downloaded from https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/gencode.v32.basic.annotation.gff3.gz. A single canonical transcript for each gene was retained based on APPRIS annotations (selection criteria in descending order of priority: appris_principal_1 through 5, appris_alternate_1 through 2, and then longest transcript) (107). Variants in the human evolution variant sets were then overlapped with internal exons in the transcript set. For short exons (10 ≤ length ≤ 120 bp), all exonic variants were retained. For long exons (121 ≤ length ≤ 500), only exonic variants within the first 90 bp from the 3′ ss were retained. This resulted in an initial set of 6,369 variant–exon pairs for the MaPSy experiment (Dataset S1). This includes additional variant–exon pairs found only in two earlier studies (1, 2) of lineage-specific variants that were later dropped and not retained in the final analysis. Removing these, this left 5,224 variants belonging to at least one of the human evolution variant sets.

See additional text in SI Appendix, Extended Materials and Methods.

Precomputed genome-wide raw SpliceAI scores (v1.3) were downloaded from https://basespace.illumina.com/s/otSPW8hnhaZR (73). SpliceAI max scores were defined as the max of raw SpliceAI donor gain, donor loss, acceptor gain, and acceptor loss scores. Variants that overlap multiple genes were assigned to the one with highest SpliceAI max raw score. Variants were binned by SpliceAI max scores into no splicing (score of 0), weak (between 0.01 and 0.19), moderate (between 0.2 and 0.49), and strong (0.5 to 1) effect categories. Pair-wise differences in proportions were tested using Fisher’s exact test. Enrichment analyses for each SpliceAI effect category were based on AF-matched nonintrogressed controls. SpliceAI scores are available in additional data at https://doi.org/10.5281/zenodo.7158564 (70).

Hexamer-based scores for exonic splicing regulatory based on high-throughput minigene experiments were downloaded from two previous studies. Enrichment index (EI) hexamer scores were downloaded from https://genome.cshlp.org/content/21/8/1360/suppl/DC1 (SI Appendix, Table S1 of ref. 77). Exonic alternative 3′ ss (exonic A3SS) and 5′ ss (exonic A5SS) hexamer scores were downloaded from https://github.com/Alex-Rosenberg/cell-2015/blob/master/ipython.notebooks/Cell2015_N4_Motif_Effect_Sizes.ipynb (78). Delta EI scores for a variant in a gene were calculated as the mean of six overlapping hexamer EI scores on the coding strand for the DER sequence minus that of the ANC sequence. Delta exonic A3SS and Delta exonic A5SS scores were defined similarly. Hexamer score annotations for MaPSy variants are available in Dataset S1.

The Ensembl VEP command line tool (v101.0) (108) was used to annotate variants with VEP consequences using the parameters: --cache --pick --assembly GRCh37. Analysis of VEP proportions by MaPSy functional score bins was restricted to exonic variants with splice region, stop gain, missense, synonymous, 5′ UTR, and 3′ UTR VEP annotations. Enrichment analyses of VEP annotations were based on AF-matched nonintrogressed controls. VEP annotations are available in Dataset S1.

GTEx (V8) sQTLs (eQTLs) significant at FDR < 0.05 for 49 tissues defined on GRCh38 were downloaded from https://gtexportal.org/home/datasets (96). These were then converted to hg19 coordinates using CrossMap (v0.5.2) (109). Variants involved in an sQTL (eQTL) association in at least one tissue were said to overlap an sQTL (eQTL) available as additional data at https://doi.org/10.5281/zenodo.7158564 (70).

Enrichment analyses were conducted using the subset of 1000 Genomes Project (1KGP) phase three autosomal single-nucleotide polymorphisms (SNPs) found in the European superpopulation with a minor allele frequency (MAF) of at least 0.01. This was because the vast majority of individuals in GTEx are of European descent, and GTEx QTL results are only available for SNPs with MAF > 0.01. 1KGP SNPs on chromosomes 1 to 22 in hg19 were downloaded from http://hgdownload.cse.ucsc.edu/gbdb/hg19/1000Genomes/phase3/ (66). 1KGP SNPs were initially binned by their European allele frequency (EUR AF) into 50 AF bins (allele counts of 1 through 10 plus 40 equal sized bins for allele counts >10), and by a univariate summary of LD in Europeans (EUR LD, defined for a given SNP as the number of SNPs with European minor allele frequency (MAF) >0.01 within 1 Mb with r2 > 0.2) into 20 equal sized LD bins. LD between pairs of 1KGP SNPs was calculated using vcftools (v0.1.16) (110) using the parameters: --hap-r2 --maf 0.01 --ld-window-bp 1000000 --min-r2 0.2 --max-missing 1. 1KGP SNPs overlapping regions of low recombination (<0.1 cM/Mb) [based on the sex-averaged recombination rate map from (111)] and outside the 1KGP phase 3 accessibility mask (downloaded from the UCSC Table Browser) were excluded from the analysis (66). Finally, the set of nonintrogressed controls was defined as the subset of 1KGP autosomal SNPs with European MAF > 0.01 not already found in the introgressed variant set available as additional data at https://doi.org/10.5281/zenodo.7158564 (70).

We describe the enrichment analyses for GTEx sQTLs among introgressed variants, but all enrichment analyses are conducted similarly. For introgressed variants in an AFxLD bin, we sample with replacement an equal number of nonintrogressed variants belonging to the same AFxLD bin. We then calculate the log2 odds ratio that the introgressed variants overlap an sQTL in at least one tissue compared to the set of matched nonintrogressed controls. We do this for 1,000 randomly sampled sets of matched nonintrogressed controls to quantify the variation in log2 odds ratio estimates. The P-value for the two-tailed test of the null hypothesis that the log2 odds ratio equals 0 is calculated as 2*min (fraction of estimates above 0, fraction of estimates below 0). Finally, to test for a difference between the log2 odds ratio estimates for introgressed versus adaptively introgressed variants, we calculate a two-tailed P-value in the same way but for the 1,000 differences between their log2 odds ratios. VEP and SpliceAI annotations do not depend on LD, whereas GTEx QTLs are highly dependent on LD. Enrichment analyses for VEP annotations were done using 100 randomly sets of AF-matched nonintrogressed controls, for GTEx eQTLs and sQTLs using 1,000 randomly sampled sets of AF- and LD-matched nonintrogressed controls, and for SpliceAI categories using 1,000 randomly sampled sets of AF-matched nonintrogressed controls. Because the background sets used are the set of exonic variants included in the MaPSy experiment, and the set of variants that could be scored by SpliceAI within 50 bp of an exon, this analysis controls for the known enrichment of genes overlapping adaptive introgression.

Genome tracks for MaPSy functional scores, SpliceAI scores, and GTEx sQTL P-values were visualized using Gviz (v1.32.0) (112). SpliceAI scores were only shown for other variants that share a variant set with the target variant. GTEx results were shown only for other variants associated with the same gene as the target variant. 1KGP population allele frequencies were downloaded from the Geography of Genetic Variants Browser (https://popgen.uchicago.edu/ggv/) (113). The ATtRACT RBP motif database was used to search the DER or ANC 11-mer sequence window centered on a target variant for RBP-binding motifs of at least 6 bp in length (81). eQTLGen Consortium whole-blood eQTLs are from https://www.eqtlgen.org/cis-eqtls.html (86).

Fine-mapping results for 94 complex traits in the UKB and BBJ were downloaded from https://www.finucanelab.org/data and https://pheweb.jp/ (94). Fine-mapped variants with PIP ≥ 0.1 were overlapped with either MaPSy ESMs or variants with SpliceAI max raw score ≥ 0.1 taking the higher of FINEMAP or SuSiE PIP scores (Datasets S4 and S5).

Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file. Click here for additional data file.