Archaic adaptive introgression in modern human reproductive genes

We found 47 archaic segments collectively in global, modern human populations that overlap genes associated with reproduction. These segments represent 37.88 Mb (million base pairs/megabases) of sequence data with at least one archaic variant in each segment surpassing a frequency of 40% (Supplementary Data  1 ), which is 20 times higher than the typical introgressing archaic DNA in modern humans 23 . Regarding regional-specific segments suggestive of adaptive introgression, we found 26 segments in American, 17 in East Asian, 6 in European, 1 in Middle Eastern, and 6 in Oceanic populations, respectively (Supplementary Data  1 ). Each of these 47 independent segments intersects at least three previously published datasets describing archaic segments recovered from modern humans (Supplementary Text); therefore, we are confident our segments represent authentic signatures. We extracted the 47 global segments from each of our populations and grouped them by region to visualise the admixture relationships between modern humans and the archaic samples within these specific regions. These relationships are shown in Fig.  1 where we present the percentage of recovered variants from SPrime 10 per geographic region as labelled by map_arch 41 . Fig. 1 Variants recovered in adaptively introgressed reproductive genes. Putative archaic variants recovered from segments overlapping adaptively introgressed reproduction-associated genes. The colours denote map_arch labelling 41 and represent the percentage of total recovered variants for the A Altai Neanderthal, B Chagyrskaya Neanderthal, C Vindija Neanderthal, and D Denisovan archaic samples that were recovered in each geographic region by SPrime 10 . The population codes refer to: African (AFR), American (AMR), East Asian (EAS), European (EUR), Middle Eastern (EAS), Oceanic (OCE), and South Asian (SAS). The plots were generated using the ggplot2 package 140 in R 120 . E SNP density plot of variants overlapping our 47 independent segments (Supplementary Data  1 ), which was generated using rMVP 141 . The chromosome length is truncated at the last SNP in our analysis per chromosome. High-density regions are shown in red and low-density regions in green. Putative archaic variants recovered from segments overlapping adaptively introgressed reproduction-associated genes. The colours denote map_arch labelling 41 and represent the percentage of total recovered variants for the A Altai Neanderthal, B Chagyrskaya Neanderthal, C Vindija Neanderthal, and D Denisovan archaic samples that were recovered in each geographic region by SPrime 10 . The population codes refer to: African (AFR), American (AMR), East Asian (EAS), European (EUR), Middle Eastern (EAS), Oceanic (OCE), and South Asian (SAS). The plots were generated using the ggplot2 package 140 in R 120 . E SNP density plot of variants overlapping our 47 independent segments (Supplementary Data  1 ), which was generated using rMVP 141 . The chromosome length is truncated at the last SNP in our analysis per chromosome. High-density regions are shown in red and low-density regions in green. Even with our specific focus on only genes associated with reproduction, we were able to recreate many of the previously established introgression patterns in modern human groups from archaic populations. These include the Chagyrskaya and Vindija Neanderthals sharing more alleles with the introgressed Neanderthal sequence in modern humans 38 , 39 , and Oceanic populations having a higher proportion of Denisovan ancestry than other global populations 5 , 11 , 17 (Fig.  1A–D ). As expected, African populations have considerably less archaic segments relative to non-Africans 13 (Fig.  1A–D ), where less than 1% of the genome-wide variants in Africans were direct matches to the archaic allele from any of our four archaic samples, as highlighted by SPrime 10 and map_arch 41 . The list of archaic alleles with frequencies of at least 40% in modern humans can be seen in Supplementary Data  2 . We applied a broad range of selection tests to both the general introgressed segments identified by SPrime 10 , and the smaller core haplotypes within each segment, which only include regions where the maximum archaic allele frequency variant overlaps one of the reproductive genes of interest in our analysis. The results presented here will focus only on core haplotypes. Some of the general introgressed segments identified by SPrime can be quite large, nearing 1 Mb, overlapping a substantial number of genes. Thus, it may be possible that putative selective signatures identified in these large segments may be due to selection acting on genes that are relatively far from our reproductive genes of interest. After filtering our data to identify core haplotypes (see “Materials and Methods”), we were able to generate 11 regions that overlap 15 genes, which we believe are the most probable segments of adaptive introgression in our analysis, where these regions are presented in Supplementary Data  3 . The PNO1-ENSG00000273275-PPP3R1 region (chr2:68342443-68503920) in the Chinese Dai in Xishuangbanna, China (CDX), AHRR segment (chr5:416118-450777) in the Finnish in Finland (FIN), and FLT1 region (chr13:28962942-28997886) in the Peruvian in Lima Peru (PEL) had extended haplotype homozygosity (EHH) 42 , F ST 43 , and Relate 44 selection test results suggesting positive selection within the core haplotype (Supplementary Data  3 A– C, F ). AHRR showed the highest number of variants at the top 1% of the genome-wide distribution for Relate’s statistic, with 10 variants surpassing this threshold (Fig. 2A ), while the other two core haplotypes only had one variant each that reached the cut-off. For this reason, we believe AHRR to be the strongest candidate of adaptive introgression in our data set. Fig. 2 Evidence of selection in the AHRR core haplotype. The A AHRR core haplotype region on chromosome 5, where 10 variants within the core haplotype (green dots), and 80 additional variants in the SPrime 10 located archaic segment (blue dots), all show evidence of selection by Relate 44 . Each of these variants was at the top 1% of the genome-wide distribution (red dashed) within the FIN. The AHRR main core haplotype allele, chr5:421193 (rs67898037) is shown (red diamond), however, it did not reach the top 1% threshold. Despite this, 90 variants in the introgressed segment do show evidence of selection, with some having the highest selection scores in the entire region. The plot was generated using the locuszoomr package 142 in R 120 . B EHH 42 decay curve for the archaic allele (red) and non-archaic allele (blue) for a 300,000 bp window around rs67898037. The archaic, derived allele has a larger homozygous region compared to the non-archaic allele, suggesting selection may be favouring the archaic allele. C Rapid increase in allele frequency of the main core haplotype allele, rs67898037, in the FIN, and D chr5:416118 (rs56279338) and E chr5:443236 (rs115554641) ancestral recombination graphs showing a derived mutation (red branches) expanding in the FIN relative to the Yoruba in Ibadan, Nigeria (YRI), all indicating positive selection within the core haplotype region. EHH and allele frequency graphs were generated using ggplot2 140 in R 120 , while the ancestral recombination graphs were created using Relate 44 . The A AHRR core haplotype region on chromosome 5, where 10 variants within the core haplotype (green dots), and 80 additional variants in the SPrime 10 located archaic segment (blue dots), all show evidence of selection by Relate 44 . Each of these variants was at the top 1% of the genome-wide distribution (red dashed) within the FIN. The AHRR main core haplotype allele, chr5:421193 (rs67898037) is shown (red diamond), however, it did not reach the top 1% threshold. Despite this, 90 variants in the introgressed segment do show evidence of selection, with some having the highest selection scores in the entire region. The plot was generated using the locuszoomr package 142 in R 120 . B EHH 42 decay curve for the archaic allele (red) and non-archaic allele (blue) for a 300,000 bp window around rs67898037. The archaic, derived allele has a larger homozygous region compared to the non-archaic allele, suggesting selection may be favouring the archaic allele. C Rapid increase in allele frequency of the main core haplotype allele, rs67898037, in the FIN, and D chr5:416118 (rs56279338) and E chr5:443236 (rs115554641) ancestral recombination graphs showing a derived mutation (red branches) expanding in the FIN relative to the Yoruba in Ibadan, Nigeria (YRI), all indicating positive selection within the core haplotype region. EHH and allele frequency graphs were generated using ggplot2 140 in R 120 , while the ancestral recombination graphs were created using Relate 44 . We created a haplotype network (Fig.  3A ) and generated ancestral recombination graphs (Figs.  2D, E and 3B ) for AHRR due to it likely being the best candidate for both adaptive introgression and positive selection in our data set. The introgressed segments containing AHRR and PNO1-ENSG00000273275-PPP3R1 both have the most similarity with the Chagyrskaya Neanderthal, while FLT1 has the highest match/mismatch ratio with the Vindija Neanderthal (Supplementary Data  3A ). A list of all putative archaic donors is listed in Supplementary Data  3 . We also generated a contour plot 41 for AHRR to visualise how closely related the AHRR core haplotype region is to all four of the archaic samples in our analysis (Fig.  3C ), where again the Chagyrskaya Neanderthal has the highest match/mismatch ratio to this region in the FIN. Fig. 3 Haplotype relationships for AHRR. The haplotype network plot A was generated using PopArt 138 for AHRR , the core haplotype with the largest amount of evidence of positive selection within its core region. The archaic haplotypes are labelled with arrows, and the number of mutations along each edge between nodes is shown in brackets. B Ancestral recombination graph for chr5:423595 (rs56346498), one of the archaic alleles in the AHRR core haplotype. This variant displays an archaic-like haplotype signature with a derived branch (red) originating before 1,000,000 years ago, followed by a long, non-recombining branch that, once introduced to modern human populations, sees rapid expansion in a short time span. This variant shows high frequency in the FIN population but is absent in the YRI reference. The plot was generated using Relate 44 . C A contour plot illustrating the match/mismatch ratio of the AHRR core region (white crosshair) to the archaic samples in our analysis. The heatmap displays the number of introgressed segments genome-wide at each represented match/mismatch ratio, with the Chagyrskaya Neanderthal having the highest match rate for this particular segment. The plot was created using the kde2d function in the MASS package 125 from R 120 using modified code given in ref. 41 . The haplotype network plot A was generated using PopArt 138 for AHRR , the core haplotype with the largest amount of evidence of positive selection within its core region. The archaic haplotypes are labelled with arrows, and the number of mutations along each edge between nodes is shown in brackets. B Ancestral recombination graph for chr5:423595 (rs56346498), one of the archaic alleles in the AHRR core haplotype. This variant displays an archaic-like haplotype signature with a derived branch (red) originating before 1,000,000 years ago, followed by a long, non-recombining branch that, once introduced to modern human populations, sees rapid expansion in a short time span. This variant shows high frequency in the FIN population but is absent in the YRI reference. The plot was generated using Relate 44 . C A contour plot illustrating the match/mismatch ratio of the AHRR core region (white crosshair) to the archaic samples in our analysis. The heatmap displays the number of introgressed segments genome-wide at each represented match/mismatch ratio, with the Chagyrskaya Neanderthal having the highest match rate for this particular segment. The plot was created using the kde2d function in the MASS package 125 from R 120 using modified code given in ref. 41 . Our identified segments (Supplementary Data  1 ) overlap 880 variants within 118 genes and intergenic segments (Supplementary Data  2 ). Of these variants, 327 of them, located within 47 genes, are found at genome-wide significant thresholds. These results are listed in Supplementary Data  4 . From this larger list of variants, 114 were found to overlap 10 of our core haplotypes. Next, we explored potential functional consequences of variants by running SNPnexus 45 , 46 , which we provide in Supplementary Data  5 . We wanted to see if any of our genome-wide significant variants (Supplementary Data  4 ) were associated with gene expression levels (i.e. were expression quantitative trait loci (eQTLs)) and used the SNP2GENE function within FUMA GWAS 47 , 48 to analyse these variants. In total, 308 variants were discovered to be eQTLs regulating 176 genes, while 113 eQTLs regulating 44 genes were found within our core haplotype regions. Supplementary Data  6 lists all of the eQTLs, while Supplementary Data  7 lists only the eQTLs from the core haplotypes. Approximately 94% of the archaic variants showing genome-wide significant associations with complex traits and diseases are also eQTLs (Supplementary Data  6 ). Over 74% of these eQTLs modulate gene expression in ovaries, prostate, testes, uterus, and vagina. When filtering for eQTLs overlapping with core haplotype segments, this figure is increased to 81% of the eQTL set, save the uterus where no core haplotype eQTLs regulated gene expression in this tissue (Supplementary Data  7 ). We noted that only 63% of the eQTLs in our results are related to several major organs, including the brain, liver, lung, spleen, and stomach. Interestingly, 73% of the eQTLs are found expressed in the heart and musculoskeletal tissues, which is similar to the number of eQTLs found in reproductive tissues, but when looking at core haplotypes specifically, the number of eQTLs expressed in heart and musculoskeletal tissues dropped to 69.9% (Supplementary Data  7 ). Therefore, archaic eQTLs in our analysis, particularly in the core haplotypes, appear to be important in helping regulate genes expressed in the ovaries, prostate, testes, uterus, and vagina with a higher number of eQTLs regulating gene expression in these tissues compared to others. For all genes in our dataset ( n  = 118), we used BioMart 49 to explore gene ontology (GO) associations generally, which are included in Supplementary Data  8 . Expanding on this analysis, we further examined our genes to see if they were significantly enriched in any pathways or ontologies. We did this in four separate runs, by first including all genes from Supplementary Data  2 with genome-wide significant variants (set 1; n  = 47), core haplotype genes with genome-wide significant variants (set 2; n  = 10), all genes that are regulated by archaic eQTLs (set 3; n  = 176), and genes that are regulated by archaic eQTLs found within a core haplotype region (set 4; n  = 44). We used ShinyGO 50 and Enrichr 51 – 53 to examine each of the four datasets above. Within set 1, we found no significant ontologies; however, nine pathways were significantly enriched in our ShinyGO results, five were significantly enriched according to the Reactome 2022 54 database, and two were significantly enriched in the KEGG 2021 human 55 pathways identified by Enrichr. The ShinyGO results for set 1 can be seen in Supplementary Data  9 , while the Enrichr results can be seen in Supplementary Data  10 . For set 2, ShinyGO and Enrichr were both unable to identify any significant pathways. However, Enrichr located 58 significant GO biological processes and 10 significant GO molecular functions 56 , 57 . The results for set 2 are listed in Supplementary Data  11 . In set 3, there were no significant ontologies or pathways found according to Enrichr 51 – 53 . Our ShinyGO 50 analysis for set 3, alternatively, was able to identify three significant pathways, which included carbon metabolism, central carbon metabolism in cancer, and DNA replication. These results can be viewed in Supplementary Data  12 . Lastly, for set 4, ShinyGO and Enrichr both indicated no significant pathways for these genes. Enrichr also did not find any significant GO biological processes 56 , 57 , but did highlight four significant GO molecular processes, which can be seen in Supplementary Data  13 .

Our study identified 47 independent segments within modern humans harbouring high-frequency archaic variants from Neanderthals and Denisovans, suggesting adaptive introgression. These segments are located within genes that have been previously described as being related to reproduction and fertility in either humans or mice 33 . In total, these segments overlap 880 variants within 118 genes or intergenic segments, where the variants have archaic allele frequencies ≥40%. Further filtering of these segments revealed 11 core haplotypes, where three of these core haplotypes displayed evidence of positive selection directly in their identified core region. We were able to identify that variants overlapping genes previously described as being related to reproduction within humans and mice show pleiotropic expressions, with genome-wide significant markers exhibiting a wide array of traits. As discussed above, other studies have documented evidence of purifying selection across the genome against archaic sequences 28 – 30 , and this may explain, to some extent, the absence of significant GWAS results for reproductive traits for the archaic variants identified in our analysis. Nonetheless, we were able to identify 15 high-frequency archaic variants across four genes that reduce the risk of prostate cancer when considering the archaic allele; however, they were not within core haplotypes. Additional evidence for archaic alleles potentially being important in cancer comes from our enrichment analysis, where three of our gene sets have significant enrichment for pathways previously implicated in either cancer development or suppression. We reported that adaptively introgressed variants are significantly enriched in gene pathways related to the development of muscle and heart tissues and with the Notch signalling pathway, which has been previously described to be associated with proper foetal organ development. Lastly, we identified 308 archaic eQTLs that regulate gene expression levels across a variety of tissues. Over 74% of these archaic eQTLs regulate genes with expression in the ovaries, prostate, testes, uterus, and vagina, with only eQTLs expressed in the testis and vagina being seen significantly more than expected. This may suggest that archaic haplotypes may play a currently unidentified role, or assist in connection with other processes, in developmental pathways crucial in embryonic development and later assist in other important pathways throughout a person’s lifetime. Our analysis contributes to the study of archaic introgression and its impact on traits within the modern human genome. To our knowledge, this is the first assessment of positive selection and adaptive introgression in segments of the modern human genome known to be introgressed from archaic hominins that also overlap genes associated with reproduction. Our results describe an intriguing relationship for archaic variants within a modern human genetic background. These variants are interesting due to their strong associations with developmental-related pathways, cancer incidence, and potential roles in helping indirectly balance important endogenous levels that have been previously tied to fertility outcomes. Despite these results, future work is necessary to fully understand the connections between these variants and their roles within the human genome. In conclusion, our work has clearly outlined the adaptive value of a number of variants and genes that are the direct result of admixture with archaic hominins and the roles they may have regarding human health, disease, and regulating normal development.

EHH 42 was the most sensitive statistic in our analysis, showing higher levels of haplotype homozygosity for the archaic allele compared to the non-archaic allele for all of the core haplotypes, with the exception of AMPH , as expected under an adaptive introgression scenario (Supplementary Data  3B ). Following this, F ST 43 highlighted seven core haplotypes that were very highly genetically differentiated at the main core haplotype allele relative to the mean F ST value genome-wide (Supplementary Data  3C ). This differentiation is particularly unusual because the main core haplotype allele had average F ST values 7.5 times higher relative to the rest of the genome, and the core haplotype mean F ST averaged over 5.5 times higher than the genome-wide mean F ST values. As discussed above, Relate 44 identified three core haplotypes where the core haplotype had at least one variant within the top 1% of the population-specific empirical distribution for their selection statistic (Supplementary Data  3F ), with AHRR showing 10 variants that passed this threshold (Fig.  2A ). Additionally, it also showed the main core allele, chr5:421193 (rs67898037) quickly increasing in frequency (Fig.  2C ), and evidence of several archaic variants within the core rapidly expanding through the FIN (Fig.  2D, E ). In short, we are confident that our analysis was able to highlight archaic segments with evidence of selection. The remaining selection tests, number of segregating sites by length ( nS L ) 58 , RAiSD 59 , saltiLASSI 60 , Tajima’s D 61 , and cross-population number of segregating sites by length (XP-nSL) 62 showed no evidence of selection within any core haplotype region (Supplementary Data  3A ). However, the normalised windowed results of nS L and XP-nSL, and the saltiLASSI Λ statistic highlighted selection in the AHRR general segment outside of our core haplotype. The nS L windowed results scored the chr5:800116-900116 region in the top 5% of windows genome-wide, and XP-nSL scored it in the top 1% of windows genome-wide for the FIN (Supplementary Data  3D, I ). The saltiLASSI Λ statistic highlighted the chr5:905191-914835 region as being at the top 1% of the Λ empirical distribution for the FIN as well (Supplementary Data  3G ). Prior research has attempted to determine if archaic alleles were immediately adaptive within humans (immediate selection) or became adaptive at some point later (standing selection), where immediate selection would be seen fairly equally in Eurasian populations and standing selection would be isolated to geographic regions or even individual populations 63 . Briefly, evidence for both immediate and standing selection was discovered across Eurasia 63 , which was also supported in a later analysis 64 . We noted a general trend of region-specific selection in our results by way of elevated archaic allele frequencies within geographically similar groups (Supplementary Data  2 ). Additionally, over 18% of our core haplotype segments were discovered in American and East Asian populations and over 54% of them were discovered in Oceanic populations (Supplementary Data  3A ). This contrasts with only 9% of the core haplotypes being found in European populations and none being found in Middle Eastern or South Asian populations (Supplementary Data  3A ), which is counter to the immediate selection hypothesis, as all Eurasian groups should have relatively equivalent allele frequencies for our adaptively introgressed variants under this scenario. To investigate this further, we used Relate 44 to determine the timing of selection acting on these regions by tracking allele frequencies of the main core haplotype variants. Recently, it has been shown that almost all of the gene flow between modern humans and Neanderthals occurred in a very short window between 50,500 and 43,500 years ago 65 , with Denisovan gene flow events likely occurring over several different pulses 10 , 11 . Current estimates suggest that non-African populations have no more than 2% Neanderthal ancestry 13 and less than 1% Denisovan ancestry 14 , 15 , except in Oceania and some populations in the Philippines, which have substantially higher levels of Denisovan ancestry than the current global average 5 , 11 , 16 , 17 . Due to this, we calculated when allele frequencies rose above 10%, as we believe this to be a reasonable signal of selection on an archaic SNP since it is at least five times higher than the typical introgressing archaic allele frequencies seen globally. Our saltiLASSI 60 results clearly indicated that all of our core haplotypes are part of incomplete, soft sweeps (Supplementary Data  3G ), and as a result, we can see evidence of standing selection using Relate. For instance, the main core allele for AHRR , chr5:421193 (rs67898037), retains allele frequencies below 10% up until ~10,000 years ago, where it begins a rapid increase in frequency to its present-day levels of ~44% (Fig.  2C ; Supplementary Data  14 ). The same pattern can also be seen in FLT1 for its main core allele, chr13:28997886 (rs11620315), where it too begins to expand after 10,000 years ago to its current frequency at ~42% with the PEL (Supplementary Data  14 ). For PNO1-ENSG00000273275-PPP3R1 , the allele frequency trajectories within the core haplotype show that selection may have begun as early as ~26,800 years ago and as late as ~13,900 years ago, depending on the tested variant, up to their current levels of just under 41% in the CDX (Supplementary Data  14 ). Therefore, these variants remained at relatively low levels in these populations for at least 23,000 years before reaching appreciable frequencies. Taken together, our results suggest recent positive selection acting on a regional level, favouring archaic haplotype segments, rather than selection favouring archaic variants immediately after admixture. Our haplotype network (Fig.  3A ) analysis determined that a number of European haplotypes cluster only two to three mutations away from the Chagyrskaya Neanderthal haplotype, providing further evidence that this is the highest probability donor of this core haplotype. The Vindija Neanderthal also plots closely to the Chagyrskaya Neanderthal haplotype, while the Altai Neanderthal and Denisovan haplotypes plot substantially farther away. Evidence further supporting this can be seen in our contour plot with closer relationships to the Chagyrskaya and Vindija Neanderthals at the AHRR core haplotype segment relative to the Altai Neanderthal and Denisova 3 (Fig.  3C ). The archaic haplotypes also generally have a large number of mutations between them and African haplotypes, suggesting isolated evolutionary trajectories within this core haplotype region. To visualise the AHRR core haplotype in more detail, we generated a haplotype heatmap. Our results show that several FIN samples are the most closely related haplotypes to the Chagyrskaya Neanderthal haplotypes (Supplementary Text; Supplementary Figs.  1 , 2 ). Additionally, our ancestral recombination graphs also support our beliefs that the AHRR core haplotype is of archaic origin. There are several variants within the core haplotype displaying derived mutations absent in Africans that originated at least 1,000,000 years ago and have long, non-recombining branches where the mutation was isolated in archaic populations until it was introduced to modern humans, where it then spread quickly throughout these populations (Figs.  2D, E and 3B ). In summary, multiple lines of evidence all suggest the AHRR core haplotype in the FIN is derived from an archaic population similar to the Chagyrskaya Neanderthal, and further, that this haplotype was not shared with African populations. Our analysis identified 11 core haplotypes, three of which likely have evidence of positive selection within their core regions, and include AHRR , FLT1 , and the PNO1-ENSG00000273275-PPP3R1 region (Supplementary Data  3A ). AHRR was included in Greer and colleagues 33 because of its association with micropenis in humans 66 . AHRR has also been associated with infertility in Asian men when exposed to persistent organic pollutants 67 . Intriguingly, one variant overlapping the AHRR core haplotype but absent in our data, chr5:422955 (rs2292596), has also been implicated with sperm apoptotic behaviour in European men who were exposed to the same type of pollutants 68 . In women, AHRR was previously described to increase endometriosis risk in Japan 69 ; however, other studies did not come to the same conclusions 70 or were inconclusive 71 . Further, studies of AHRR have highlighted its potential role as a tumour suppressant for several forms of cancer 72 , 73 and how it may have a role in inflammatory disease 74 . FLT1 also showed evidence of selection (Supplementary Data  3A ), and was included in Greer et al. 33 due to its association with post-implantation embryo development. Prior research has established a causal marker near FLT1 linked with preeclampsia 75 – 77 and the gene has been speculated as able to negatively regulate angiogenesis in mice, ultimately leading to embryonic death if divergent genotypes are found 78 . Additional phenotypes for FLT1 include pregnancy loss, low birth weights, and placental inflammation 79 . Our data are consistent with previous research that described signatures of positive selection within FLT1 generally 80 and more specifically in the PEL population as a result of archaic admixture 81 . The last core haplotype in our analysis displaying putative selection is the PNO1-ENSG00000273275-PPP3R1 region on chromosome 2, where PNO1 and PPP3R1 were included in the Greer et al. 33 dataset for their roles in pre-implantation embryo development and foetal development in mice, respectively. In knockout mice experiments, strains with no PNO1 stopped developing during the embryo stage and ultimately did not survive this phase of development, while heterozygous and transgenic mice had no apparent anomalies and were fertile 82 . In humans, PNO1 was shown in males to be associated with differentiating spermatogonia 83 in the Male Fertility Gene Atlas 84 and has also exacerbated progression of certain cancers 85 – 87 . No studies to our knowledge have described the role of PPP3R1 regarding its impact on human reproduction; however, previous research has linked this gene to Alzheimer’s disease pathogenesis 88 . Unfortunately, there were no genome-wide significant markers overlapping any of our results that were associated with any sort of reproductive trait (Supplementary Data  4 ), making it difficult to understand the impact archaic introgression may have had on those phenotypes. Nonetheless, our genome-wide association study (GWAS) analysis (Supplementary Data  4 ) identified 29 archaic alleles overlapping AHRR that had reached genome-wide significance thresholds in the IEU OpenGWAS Project database 89 , 90 for 10 blood-related traits. Our results show that the archaic alleles are significantly positively associated with haemoglobin, haematocrit, and red blood cell levels in Europeans 91 . Two of the archaic variants, chr5:434722 (rs34453673) and chr5:421455 (rs11741954), had coding mutations with the former being a synonymous and the latter a non-synonymous change for the archaic allele, respectively (Supplementary Data  5 ). The IEU OpenGWAS database shows that one archaic allele, chr13:29069039 (rs55927955) overlapping FLT1 , has a genome-wide significant negative association with lung function 92 (Supplementary Data  4 ), but it was not associated with any protein changes (Supplementary Data  5 ). Within the PNO1-ENSG00000273275-PPP3R1 core haplotype there are 31 genome-wide significant archaic alleles listed in the IEU OpenGWAS data, none of which had protein coding changes for the archaic allele (Supplementary Data  5 ). The 31 archaic alleles were shown to be negatively associated with various vein-related traits in both the FinnGenn 93 and UK Biobank datasets 94 , but are also positively associated with height 91 , 92 and waist-to-hip ratio 92 , 95 . Outside of the core haplotypes with evidence of selection, we identified 15 variants in four genes that have significant associations with prostate cancer risk (Supplementary Data  4 ) within the chr2:173199359-173589239 segment recovered within the Tujia population (Supplementary Data  2 ). These genes, ENSG00000225205 , ENSG00000226963 , ITGA6 , and ITGA6-AS1 , have archaic alleles that reduce the risk of prostate cancer (Supplementary Data  4 ) reported in large GWAS cohorts from Asia and Europe 91 , 94 , 96 , 97 . Previous research has shown that one variant in our analysis, chr2:173311553 (rs12621278), was also found to be protective for prostate cancer onset in Asian populations 98 . Our enrichment analyses also show that several of the genes with high-frequency archaic segments are important in cancer pathways, where ShinyGO 50 highlighted four significantly enriched cancer pathways within our set 1 genes (Supplementary Data  9 ) and one significantly enriched pathway, central carbon metabolism in cancer, within set 3 (Supplementary Data  12 ). Further evidence for our genes being associated with prostate cancer comes from set 2, where one of our core haplotypes, PPP3R1 , has been shown to be associated with pathways related to prostate cancer previously 99 and is also enriched in several GO biological processes and molecular functions 56 , 57 including calcium ion binding, calcium-mediated signalling, and calcineurin-mediated signalling (Supplementary Data  11 ). Calcium ions and calcineurin have both been shown to have roles regarding prostate cancer development and growth 100 , 101 . Lastly, some genes regulated by core haplotype eQTLs are significantly enriched in palmitoyltransferase activity, which previous research implicated in cancer initiation and modulation 102 . Our results from four genes in set 1 show ties with the Notch signalling pathway (Supplementary Data  9 and 10 ). The Notch signalling pathway was described before as being important for controlling proper foetal development of several vital organs 103 , 104 . Some of these genes have reproductive associations within mouse databases 105 . For example, mutant forms of JAG1 were associated with sterility and premature death 106 . In humans, KAT2B , a hormone receptor, was differentially expressed in polycystic ovary syndrome oocytes 107 . Additionally, our core haplotype segments (set 2) identified 68 GO biological processes 56 , 57 , 18 of which are related to cardiac, muscle, or skeletal biological processes (Supplementary Data  11 ). One gene, TTN , has displayed various skeletal and cardiac anomalies using mouse models 105 . This includes an embryonic shrunken head phenotype associated with an inability to generate proper circulation, eventually leading to death of the foetus 108 and increased left ventricular diastolic stiffness, causing a reduction in activity tolerance 109 . In modern humans, disease models show that TTN is tied to dilated cardiomyopathy and two forms of muscular dystrophy 105 . Our GWAS results for TTN also support some of the associations addressed here, as 21 SNPs overlapping this gene have an archaic allele associated with lower pulse rates 110 (Supplementary Data  4 ), which may signal weakening heart abilities. Also for set 2, significant enrichment for calcium ion binding was identified (Supplementary Data  11 ). Calcium is important in mammalian reproduction, having roles in sperm development, embryogenesis, and healthy development of the foetus 111 . In our set 4 genes, we have significant GO molecular functions related to palmitoyltransferase activity, which have functions related to neuritogenesis and signalling in early development 112 . Taken together, three of the gene sets are significantly enriched in ontologies and pathways that are associated with the healthy development of the foetus across several important organs. Additional genes within these sets linked to reproduction suggest that these gene pathways may be important in modulating the proper development and function of several important processes across an organism’s lifespan. As described above, we noted a high number of our genome-wide significant archaic alleles being eQTLs for reproductive-associated tissues. To test if these results were significant, we created a randomised set of 5534 genome-wide eQTLs to use as a comparison to our results and ran these through FUMA GWAS 47 , 48 , finding that only 56.7% of these ( n  = 3138) were expressed in reproductive tissues (Supplementary Text). The archaic eQTLs in our results tend to be associated with reproductive tissues in a Chi-square analysis, however, only testis were observed significantly more than expected, X 2 (1, N  = 308) = 19.207, p  = 0.0000117, while archaic eQTLs expressed in the uterus were observed significantly less than expected, X 2 (1, N  = 308) = 33.408, p  = 0.0000000747 (Supplementary Text; Supplementary Data  15 ). Within the core haplotype archaic eQTL set we observed that archaic eQTLs expressed in the ovaries, prostate, and uterus were seen less than expected, with ovaries, X 2 (1, N  = 113) = 13.339, p  = 0.00026, and uterus, X 2 (1, N  = 113) = 13.939, p  = 0.000189, being significant (Supplementary Text; Supplementary Data  15 ). We also found that core haplotype archaic eQTLs were expressed significantly more than expected in the testis, X 2 (1, N  = 113) = 32.466, p  = 0.0000000121, and vagina, X 2 (1, N  = 113) = 4.6175, p  = 0.031646 (Supplementary Text; Supplementary Data  15 ). However, all significant associations described in this section are relatively weak, regardless of the eQTL set (Supplementary Text; Supplementary Data  15 ). Our research is not without limitations. Several of the Human Genome Diversity Project (HGDP) target populations we analysed here have sample sizes below 15, a recommended minimum number of included individuals for SPrime to maintain high levels of accuracy in identifying introgressed segments 10 . However, due to sampling bias present in many genomic studies 113 , research at this time must rely on populations with limited sample sizes at this time to explore patterns in populations outside of Eurasia. A possible consequence of this is that some of our segments listed in our results may represent false positives, however, the likelihood of this remains small as these segments did pass our filtering parameters to be considered an authentic segment (see “Materials and Methods”), and further, many of these segments are also shared with populations with adequate sample sizes (i.e. those from the 1000 Genomes Project (1KGP)) suggesting they are attributed accurately as introgressed. A second consequence of this is elevated allele frequencies in the HGDP samples, which may not be found in populations with more than 15 samples. We acknowledge that using strictly allele frequency as a measure of adaptive introgression may lead to artificially elevated detection of adaptively introgressed segments in our results. Therefore, we suggest caution in interpreting some of the results from the HGDP due to their low sample sizes and concomitant elevated allele frequencies unless they are paired with higher-than-expected archaic allele frequencies in other populations with adequate sample sizes. Another possible drawback is our use of the SPrime software, which has been shown to mask some instances of admixture if the segments are too similar to those of the human outgroup 114 . It is possible that using other software that is reference-population-free to repeat our study, such as IBDMix 114 , may uncover additional interesting results. It is also important to note that we have used a stringent archaic allele frequency of 40% in our filtering pipeline. We chose a relatively high threshold because our main goal was to identify putative instances of adaptive introgression increasing the frequencies of archaic haplotypes in modern human populations, but this high threshold means that we may have removed other archaic segments with elevated frequencies relative to typical introgression levels worth exploring. Lastly, we also used GWAS summary statistics to infer the direction of effect of the archaic variants. The interpretation of some of these effects is not always straightforward, given that phenotypes reported in several studies may not be clearly defined. Additionally, although a number of archaic alleles display genome-wide significant effects for a wide range of traits, this does not necessarily prove they are the causal variants, and the associations may be simply due to the archaic allele being in linkage disequilibrium with the actual causal variants. Because identifying causal variants was not the aim of this project, future research should explore fine-mapping and functional validation to uncover whether the variants we present here are causal in these trait effects or simply a byproduct of linkage disequilibrium. Exploration of polygenic adaptation signals in modern humans, individually, or through regionally shared signatures, may also be useful to see if the archaic alleles highlighted in this study contribute to broader polygenic selection signals associated with reproduction, an analysis that was outside of the scope of our current work.

Anatomically modern humans ( Homo sapiens ), Denisovans, and Neanderthals share an evolutionary history, with a common ancestor 765,000–550,000 years ago 1 . Modern humans evolved in Africa about 300,000 years ago 2 and began the first of many expansions out of Africa by at least 85,000 years ago 3 . Upon dispersal, they encountered and interbred with archaic Homo , the Neanderthals and Denisovans, episodically in diverse geographic regions 4 – 13 . As a result of these admixture events in the past, estimates show that, on average, most non-African populations have just less than 2% of their DNA composed of Neanderthal sequence 13 . Estimates of Denisovan contributions to some human populations are less than 1% on average 14 , 15 , save for some Oceanic populations that have been shown to have nearly 5% of their genome shared with Denisovans 5 , 11 , 16 . However, the highest level of Denisovan ancestry recovered to date occurs in the Philippine Ayta 17 . The consequences of these admixture events within the context of modern human evolution have been mixed. Some archaic DNA segments have been shown to be adaptively introgressed into modern human populations, where these sequences, conferring some fitness benefit, saw rapid expansion within some groups to high frequencies in relatively short evolutionary timespans. These include, but are not limited to, several immunity genes 18 – 21 , keratin genes related to skin and hair phenotypes 7 , 22 – 24 , and high-altitude adaptation in Tibetan populations 25 . Contrary to this, research over the past decade has shown that introgression deserts exist within the modern human genome, suggestive of the removal of archaic segments 15 , 22 , 26 , 27 . Many of these loci were thought to be weakly deleterious and purged via purifying selection in modern human populations 28 , 29 . Direct evidence of this reduction in archaic sequences has been documented, where examination of 51 ancient modern human samples with ages between 45,000 and 7000 years old showed a steady decline of Neanderthal ancestry towards the younger samples, consistent with purifying selection 30 . Some research has explored the potential phenotypic expression of this purifying selection. Among these, one compelling discovery was a greater-than-expected number of archaic loci being removed from the modern human X chromosome relative to the other chromosomes, which may be evidence for male sterility in modern human–archaic hybrids 15 , 22 . Further evidence comes from research showing that a Neanderthal Y chromosome haplotype has not been documented in modern humans so far, and mutations in this region have been associated with spermatogenesis, also suggesting male sterility of hybrids may have occurred 31 . Support for this can be found in Haldane’s Rule, stating that in instances of hybridisation, chromosomal pair fusion difficulties would manifest in a sterility phenotype in the heterogametic sex 32 . However, more recent research has explored reproductive phenotypes, and through modelling, they found that reproductive incompatibilities between modern humans and archaics likely did not occur, but archaics and modern humans did share some of the same deleterious variants in reproductive genes 33 . Lastly, research focusing specifically on the impact of loci associated with reproduction derived from Neanderthals into modern humans has found mixed results. For example, Li et al. 34 documented a missense variant, rs1042838, within the PGR locus associated with preterm birth that derived from Neanderthals in some European populations at frequencies as high as 18%. Later analysis showed that a Neanderthal haplotype within the PGR gene was linked with a reduction in the number of miscarriages and decreased bleeding during pregnancy, which suggests that within a modern human genetic background, this haplotype may increase fertility 35 . While prior analysis has demonstrated both beneficial and deleterious consequences as a result of admixture, to our knowledge, there has been no widespread analysis of the role introgression has played within genes related to reproduction within modern human populations. Archaic alleles within these genes could have a significant influence on sperm or oocyte development or function, embryogenesis, and other downstream functions, such as foetal development. Therefore, research of this nature warrants analysis. Prior research targeting specific genes or phenotypes has uncovered associations with archaic ancestry, such as a high-frequency Neanderthal haplotype conveying severe COVID-19 response in South Asian populations 36 , and we wanted to explore if similar mechanisms are acting on genes associated with reproduction. We speculated that, like the evidence from PGR 35 , other archaic-derived loci within genes related to reproduction may show functional benefits in modern human populations, being brought to higher frequencies suggestive of adaptive introgression within specific populations. The use of large scale genomic datasets is becoming commonplace to help better understand admixture events in the past between modern humans and archaic hominins that may have shaped phenotypes or adaptations of modern humans (for example, McArthur et al. 24 and Li et al. 37 ). We expand on these recent developments by exploring for evidence of adaptive introgression within genes related to reproduction in both mice and modern humans listed in Greer et al. 33 . We examined for evidence of archaic introgression within 1692 autosomal genes associated with reproduction using the previously published high coverage Altai Neanderthal 6 , Chagyrskaya Neanderthal 38 , Vindija Neanderthal 39 , and Denisova 3 5 , 40 samples in 76 worldwide, modern human populations.

Our main workflow steps are visualised in Fig.  4 . Detailed descriptions of our materials and methods are provided below. Fig. 4 Methods flowchart. The main workflow steps in our analysis. Input data for the archaic and modern human samples are shown as diamonds. Links to these resources are in the data availability statement, and descriptions of these data are provided in the Materials and Methods. Boxes represent main workflow steps, and octagons demonstrate end points of the analyses corresponding to sections outlined in detail in the Materials and Methods. ARGs refer to ancestral recombination graphs. The main workflow steps in our analysis. Input data for the archaic and modern human samples are shown as diamonds. Links to these resources are in the data availability statement, and descriptions of these data are provided in the Materials and Methods. Boxes represent main workflow steps, and octagons demonstrate end points of the analyses corresponding to sections outlined in detail in the Materials and Methods. ARGs refer to ancestral recombination graphs. We used SPrime 10 with the standard settings to identify regions likely of archaic origin within modern human populations. To accurately identify these segments, SPrime uses a scoring algorithm to flag putatively introgressed segments, where it is recommended to filter out segments that score below 150,000. Our modern human data came from the recently released phased gnomAD 1KGP + HGDP callset 115 , a high-resolution data panel combining the HGDP ( n  = 51 populations) and 1KGP ( n  = 25 populations) mapped to GRCh38 (hg38) coordinates. To filter our input data, we followed the SPrime protocol described by Zhou and Browning 41 . In order to effectively identify putative archaic segments, SPrime requires a non-admixed reference panel to compare the target panel against, which it uses to mask potentially ancestral alleles from being recognised as of archaic ancestry. Following the SPrime protocol, we used the YRI population from the 1KGP ( n  = 106 samples), as our reference panel and merged them with each of the other target populations within our analysis set. We updated the variant IDs for each SNP using the dbSNP database annotation files 116 for known variants so we could match variants in downstream analyses. We kept only biallelic markers and removed duplicated variants after filtering using BCFtools v1.13 117 . Our modern human data is mapped to hg38 coordinates; therefore, we used the UCSC LiftOver Linux executable 118 to convert our results to the GRCh37 (hg19) genome build, as that matches the published Neanderthal and Denisovan genomes. In case LiftOver had mapped any variants incorrectly, we removed variants that had moved from their original chromosomes. We used the Denisova 3 5 , 40 , Altai Neanderthal 6 , Vindija Neanderthal 39 , and Chagyrskaya Neanderthal 38 as our archaic samples to explore adaptive introgression within human populations. The archaic genomes and their mask files were downloaded from their hosting sites and are provided in the data availability statement. Next, we used map_arch 41 to match our results to the archaic samples. This software takes as input the SPrime 10 output, an archaic VCF, and the associated mask file to list, for each variant in the SPrime output, whether the modern human marker is a match, mismatch, or is not comparable relative to an archaic sample. To generate allele frequencies for the modern human samples, we used BCFtools 117 and subsequently merged these allele frequencies with our initial output using the dplyr package 119 in R v4.1.2 120 . Using the output file, we were able to determine the archaic allele and generated the allele frequency of this allele within our population of interest (written as archaic allele frequency). Our research focused on identifying whether any archaically derived segments overlapped genes known to be involved in reproductive traits. To do this, we downloaded the list of 1692 genes found on autosomes, which were examined in ref. 33 , from Ensembl’s BioMart GRCh37 Release 112 49 and intersected the whole gene coordinates with our SPrime 10 results. These genes potentially impact 11 different biological processes related to reproduction in either mice or modern humans and include: epididymis development or functionality, fertilisation, foetal development, oogenesis, ovary development or functionality, pre- and post-implantation embryo development, sex determination, spermatogenesis, testis development or functionality, and uterus development or functionality 33 . We focused on autosomal loci because our archaic samples are all female 5 , 6 , 38 – 40 , therefore lacking a Y chromosome, and because previous analyses have already identified a dearth of archaic segments on the modern human X chromosome 15 , 22 , 27 . We aimed to normalise our results by first reporting on variants rather than general segments, and secondly, describing our results using hg19 coordinates to match those of the archaic samples. We used the dplyr package 119 in R 120 to extract variants overlapping our genes of interest in our results by matching the variant rsIDs. We selected segments within our results that have at least 30 variants per segment that can be compared to at least one archaic sample and have match/mismatch ratios (the ratio of the number of matches divided by the total number of matches and mismatches found by map_arch) of >50% for the Neanderthals and >40% for the Denisovan, following the recommended SPrime filtering 10 . After removing segments that failed this step, we annotated our VCF files with SnpEff v5.2a 121 , allowing us to merge gene and consequence data with our population files based on matching the rsIDs. We used BEDTools v2.30.0 122 to intersect the results from each population against every other population to search for potential regional signatures within our results. Additionally, because it is possible that one archaic sample may harbour a signal absent in the other archaic populations, we generated intersection data for every modern human population relative to the four archaic samples in our analysis. To locate regions of contiguous blocks of introgression within our results, we extracted the putatively introgressed segments identified in our analysis from each population that had at least one variant overlapping a gene listed by Greer et al. 33 . Following this, we combined these segments together with populations from the same region according to the gnomAD sample metatable, and reduced them to the minimum number of non-overlapping segments using the GenomicRanges package v3.19 123 in R 120 . We also repeated this analysis to generate a global introgression map regarding our results, where we combined all of the archaic segments regardless of geographic origin in our analysis and reduced them to non-overlapping segments. Match/mismatch ratios have previously been used to gauge the similarity of genomic segments to one archaic sample compared to the others 10 , 41 , 124 . Introgressed segments with 30 or more variants with match/mismatch ratios over 60% to a Neanderthal sample with a concomitant match/mismatch ratio below 40% to the Denisovan are likely of Neanderthal affinity 10 . Similarly, segments most likely of Denisovan origin will showcase match/mismatch ratios greater than 40% to the Denisovan and less than 30% to the Neanderthals. To explore this, we assessed which segments are clearly introgressed from one archaic population relative to the others. When a segment passed these donor thresholds with a match/mismatch ratio 5% higher than the other three archaic samples, we suggest the segment is derived from that archaic population. If the match/mismatch ratio is less than 5% higher than the other archaic samples, we consider this segment to be general archaic ancestry. We applied these calculations to a series of core haplotypes and visualised them using contour plots based on the scripts provided by the SPrime Protocol 41 using the kde2d function from MASS 125 in R 120 . Next, we filtered out any variants with archaic allele frequencies below 40% in the SPrime 10 identified segment. For the remaining variants left in each segment that passed our segment and variant filters, we located variants with allele frequencies of at least 40% that directly matched the archaic allele of at least one archaic sample and were also the highest archaic allele frequency variant within their segment. Once we located these variants, we selected only the maximum archaic allele frequency variants that intersected at least one of the genes listed in Greer et al. 33 , described as being associated with reproductive phenotypes. After identifying these core alleles, we then preferentially selected archaic variants directly adjacent to these core alleles, which had allele frequencies no more than 5% below the maximum archaic allele frequency variant. The remaining windows around the maximum frequency archaic allele are what we consider to be core haplotypes within our analysis. We focused this analysis on populations from the 1KGP because of the small sample sizes of the HGDP. One exception is the Melanesian and Papuan population samples of the HGDP, which were included because they are not represented in the 1KGP reference dataset. By including these, we could investigate both potential Neanderthal and Denisovan contributions. Three of these core haplotypes were found at the edge of their larger introgressing segment (Supplementary Data  3A ), however, no high frequency variants ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ge$$\end{document} ≥ 40% allele frequency and no more than 5% below the maximum archaic allele frequency of the segment) that we used to define core haplotypes extended outside of their respective segments into adjacent segments ( Supplementary Text ). These findings illustrate the uniqueness of these particular segments compared to the other archaic segments identified throughout the rest of the genome. We wanted to identify core haplotypes with evidence of adaptive introgression; therefore, we wanted to see if any of these core regions showed evidence of positive selection using 10 different methods to detect selection. To begin, we calculated F ST using the Weir and Cockerham method 43 (--weir-fst-pop) implemented in VCFtools v0.1.16 126 in 100,000 base pair (bp) sliding windows (--fst-window-size 100,000) using a 10,000 bp stepping pattern (--fst-window-step 10,000). F ST is used to evaluate allele frequency differences between pairs of populations. Previous research has indicated that F ST can be influenced by allele frequency 127 and sample size 128 ; therefore, we tried to minimise these effects by size-matching the reference population to the target population. To do this, we used the YRI for the 1KGP populations and the Yoruba from the HGDP for the Oceanic populations. We additionally calculated the per-site F ST values with no additional parameters (--weir-fst-pop) using VCFtools. Following this, we then calculated Tajima’s D statistic 61 using VCFtools within 100,000 bp non-overlapping windows (--TajimaD 100,000) for the whole genome. Next, we used several methods capable of detecting selection on haplotypes, which identify selection based on atypical levels of homozygosity surrounding a specified region. We used Selscan v2.0.3 129 to examine population-specific statistics that measure this phenomenon, implementing EHH42 (--ehh) and nS L 58 (--nsl). The former has been shown to have comparable accuracy to newer methods when detecting selection 130 , and the latter is robust to a variety of differing population structures and selection scenarios 58 . We calculated EHH for our main core haplotype alleles and also included a 300,000 bp window (--ehh-win 300,000) surrounding the main allele to visualise the levels of homozygosity for both the archaic and non-archaic alleles at each position. When the archaic allele haplotype was larger relative to the non-archaic allele, measured when the EHH score fell below 0.25 in both directions from the main core position (Supplementary Data  3B ), we deemed this to have evidence of selection based on EHH. Per-site ancestral and derived allele nS L scores (--nsl) were obtained for the core haplotypes in their respective target population. Further, we also used XP-nSL 62 (--xpnsl), a method that does cross-population comparisons of haplotypes, using the target population and the YRI population as reference. We normalised both the nS L and XP-nSL outputs using selscan’s norm function (norm --nsl --bp-win --winsize 100000 --qbins 10 --min-snps 10; norm --xpnsl --bp-win --winsize 100000 --qbins 10 --min-snps 10) per site and additionally in a windowed method using 100,000 bp windows where each window has to have at least 10 variants. The windows were then placed into 10 quantile bins for normalisation. The pipeline was described previously by Szpiech et al. 62 . Following this, we applied two recent composite likelihood methods, RAiSD v2.9 59 and saltiLASSI v1.2.1 60 . RAiSD identifies selection in regions based on three criteria, namely, a concentrated reduction in polymorphic loci, increased linkage disequilibrium around a particular mutation, and changes to the derived variant site frequency spectrum, which are then compiled into the μ statistic. In contrast, SaltiLASSI examines the haplotype frequency spectrum for atypical patterns to identify selection, and consists of three statistics. The number of high-frequency haplotypes is reported by m . When m  = 1, a hard sweep is occurring, and when m  > 1, it suggests a soft sweep. The second statistic, A , attempts to measure the size of the sweep to infer the strength and timing of a selection event, while the final statistic, Λ, scores changes in the haplotype frequency spectrum relative to other nearby haplotypes. Lastly, we used Relate v1.2.1 44 , a coalescent tree-based approach that measures the speed a mutation proliferates throughout a population compared to other lineages within the same group as evidence of selection. We applied an outlier approach to gauge whether any of our core haplotypes showcased evidence of selection. For each statistic, we used the genome-wide output from each population and identified values within the top 1% of the empirical distribution. In the case of two-sided tests, we calculated both the top and bottom 1% values from the genome-wide empirical distribution. Tajima’s D is a two-sided statistic where values greater than zero indicate increased heterozygosity, or population expansion, and values below zero may suggest increased homozygosity, or population contraction, within a specific region 61 . Both nS L 58 and XP-nSL 62 are two-sided statistics. Large, positive values indicate selection in the direction of the derived allele when using nS L and in the direction of the target population when using XP-nSL. Conversely, large, negative values suggest selection for the ancestral allele when applying nS L and for the reference population when applying XP-nSL. All variants and windows that were within each test’s population-specific 1% threshold for both one and two-sided tests were taken as putative evidence of selection. Finally, when normalising outputs in selscan 129 , the software will report windows within the top 5 and 1% of the empirical distribution of windows genome-wide, both values which we applied to our analysis to determine selection. We generated a haplotype network for all of our modern human populations and the archaic samples for the AHRR core haplotype. This core haplotype tied for the highest number of tests showing evidence of selection with FLT1 and PNO1-ENSG00000273275-PPP3R1 (Supplementary Data  3 ); however, these latter two core regions only had one variant at the top 1% of the empirical distribution for Relate’s 44 selection statistic, while AHRR had 10 that met this criteria (see Results). For this reason, we believe AHRR is the strongest signal of adaptive introgression and subsequent positive selection in our dataset. Because haplotypes require phased data, and the archaic samples are not phased, we took advantage of the fact that our Neanderthal and Denisovan samples have very long runs of homozygosity 38 , 39 , 131 , 132 . Therefore, removing heterozygous sites from these samples will artificially phase the genomes because both alleles will be the same. Further, prior analyses have made use of a similar approach and were able to generate effective networks showing haplotype similarities between Denisovans and Tibetan populations 25 . After filtering out heterozygous sites using BCFtools 117 on chromosome 5, we removed 22,976 sites (0.0001926% of the total sites on chromosome 5) from the Altai Neanderthal, 17,938 sites (0.0001504% of the total sites on chromosome 5) from the Chagyrskaya Neanderthal, 26,736 sites (0.0002242% of the total sites on chromosome 5) from Denisova3, and 20,533 sites (0.0001721% of the total sites on chromosome 5) from the Vindija Neanderthal. As a consequence, from the AHRR core haplotype we removed four sites from the Altai Neanderthal (0.000144% of the total variants), 8 sites (0.000288% of the total variants) from the Chagyrskaya Neanderthal, three sites (0.000108% of the total variants) from Denisova3, and 10 sites (0.000361% of the total variants) from the Vindija Neanderthal, respectively. In total, we were left with 939 variants over a 34,659 bp region for analysis. We generated a FASTA format file for AHRR from the merged VCF using PGDSpider v2.1.1.5 133 followed by using MUSCLE v3.8.31 134 to align the file. We converted the aligned FASTA file to NEXUS format 135 using MEGA11 v11.0.13 136 , and then imported the file into DnaSP6 v6.12.03 137 , where we assigned population sequence sets and generated a haplotype file. PopART v1.7 138 was used to build the haplotype network using the median-joining method. We combined populations according to their region/superpopulation (i.e. AFR/African, AMR/American) in our haplotype file, traits file (haplotype counts per region), and haplotype network because of the large number of subpopulations we included in our analysis. In total, AHRR has 1379 haplotypes, so to make our haplotype network legible, we included only the top 50 most frequent haplotypes, and their ties, along with the haplotypes associated with the archaic samples, and any Oceanic haplotypes with frequencies of at least two. This left us with 55 haplotypes to plot. We used Relate 44 to generate ancestral recombination graphs. To do this, we first generated a new VCF from the gnomAD 1KGP + HGDP phased callset 115 with all of the modern human populations present and converted and filtered the VCF to the proper format using the functions native in Relate, including the RelateFileFormats --mode ConvertFromVcf and PrepareInputFiles.sh commands. Once filtering was completed, we were left with 857 variants for further analysis within the AHRR core haplotype region. To run Relate, we used the recommended settings of a mutation rate of 1.25e –8 , a generation time of 28 years, and an effective population size of 30,000. The trees we selected for further review included derived branches originating at least 1,000,000 years ago, with subsequent long, non-recombining branches, followed by rapid expansion very recently. We used Relate to extract subpopulation trees (RelateExtract --mode SubTreesForSubpopulation) to investigate the genealogies in more detail for variants with an archaic-looking tree that also were labelled as matches in map_arch 41 . Following this, we plotted mutational trees (TreeViewMutation.sh) using the built-in Relate functions with the YRI as our outgroup. We were interested in determining if any of the archaic variants in our results show genome-wide significant associations ( p  = 5 × 10 –8 or less) with complex traits and diseases. We compared the variants in our results against the IEU OpenGWAS project database 89 , 90 . For any matches that were found at genome-wide significance, we downloaded the associated GWAS summary statistics and compared these to our archaic alleles. If the effect allele did not match the archaic allele, we converted the direction of effect (beta) by flipping the sign and then updated the non-effect allele. We explored whether archaic variants showing genome-wide significant associations in our results were involved in regulating gene expression. Using FUMA GWAS’s SNP2GENE function 47 , 48 , we input the variants compiled in our OpenGWAS 89 , 90 results ( n  = 327) and then downloaded the eQTL table from the output. This process was repeated again just for significant variants that overlapped our core haplotypes ( n  = 114). SNP2GENE identifies eQTLs by matching the input SNP genomic position, alternative, and reference alleles against SNPs showing significant expression 47 in 17 curated cis -eQTL databases 139 . Variant annotation was done using the GRCh37 search in SNPnexus v4 45 , 46 to understand the functional consequences of archaic variants identified in our data ( n  = 880). Further, we used BioMart 49 to access GO Consortium information 56 , 57 for the genes with high-frequency archaic segments identified in our results ( n  = 118). We ran several analyses to determine if any of the genes found within our results were enriched in any pathways. We used ShinyGO 50 to test for significantly enriched KEGG pathways 55 for the genes that contain genome-wide significant markers ( n  = 47). We supplemented this with results from Enrichr 51 – 53 , which scans the Reactome 2022 54 and KEGG 2021 human 55 pathways, and the GO biological processes and GO molecular functions 56 , 57 databases for any other significant pathways not outlined in our first analysis. We repeated these analyses just for the genes overlapping core haplotype segments with genome-wide significant markers ( n  = 10). Lastly, we also performed enrichment analysis on the set of genes regulated by archaic eQTLs identified by FUMA GWAS 47 , 48 for all genes ( n  = 176) and just those regulated within the core haplotypes ( n  = 44). Further information on research design is available in the  Nature Portfolio Reporting Summary linked to this article.

Supplementary Information Description of Additional Supplementary Files Supplementary Data 1 Supplementary Data 2 Supplementary Data 3 Supplementary Data 4 Supplementary Data 5 Supplementary Data 6 Supplementary Data 7 Supplementary Data 8 Supplementary Data 9 Supplementary Data 10 Supplementary Data 11 Supplementary Data 12 Supplementary Data 13 Supplementary Data 14 Supplementary Data 15 Supplementary Data 16 Supplementary Data 17 Reporting Summary Transparent Peer Review file Supplementary Information Description of Additional Supplementary Files Supplementary Data 1 Supplementary Data 2 Supplementary Data 3 Supplementary Data 4 Supplementary Data 5 Supplementary Data 6 Supplementary Data 7 Supplementary Data 8 Supplementary Data 9 Supplementary Data 10 Supplementary Data 11 Supplementary Data 12 Supplementary Data 13 Supplementary Data 14 Supplementary Data 15 Supplementary Data 16 Supplementary Data 17 Reporting Summary Transparent Peer Review file