1,000 ancient genomes uncover 10,000 years of natural selection in Europe

Ancient DNA has revolutionized our understanding of human population history. However, its potential to examine how rapid cultural evolution to new lifestyles may have driven biological adaptation has not been met, largely due to limited sample sizes. We assembled genome-wide data from 1,291 individuals from Europe over 10,000 years, providing a dataset that is large enough to resolve the timing of selection into the Neolithic, Bronze Age, and Historical periods. We identified 25 genetic loci with rapid changes in frequency during these periods, a majority of which were previously undetected. Signals specific to the Neolithic transition are associated with body weight, diet, and lipid metabolism-related phenotypes. They also include immune phenotypes, most notably a locus that confers immunity to Salmonella infection at a time when ancient Salmonella genomes have been shown to adapt to human hosts, thus providing a possible example of human-pathogen co-evolution. In the Bronze Age, selection signals are enriched near genes involved in pigmentation and immune-related traits, including at a key human protein interactor of SARS-CoV-2. Only in the Historical period do the selection candidates we detect largely mirror previously-reported signals, highlighting how the statistical power of previous studies was limited to the last few millennia. The Historical period also has multiple signals associated with vitamin D binding, providing evidence that lactase persistence may have been part of an oligogenic adaptation for efficient calcium uptake and challenging the theory that its adaptive value lies only in facilitating caloric supplementation during times of scarcity. Finally, we detect selection on complex traits in all three periods, including selection favoring variants that reduce body weight in the Neolithic. In the Historical period, we detect selection favoring variants that increase risk for cardiovascular disease plausibly reflecting selection for a more active inflammatory response that would have been adaptive in the face of increased infectious disease exposure. Our results provide an evolutionary rationale for the high prevalence of these deadly diseases in modern societies today and highlight the unique power of ancient DNA in elucidating biological change that accompanied the profound cultural transformations of recent human history.

also include immune phenotypes, most notably a locus that confers immunity to Salmonella infection at a 23 time when ancient Salmonella genomes have been shown to adapt to human hosts, thus providing a 24 possible example of human-pathogen co-evolution. In the Bronze Age, selection signals are enriched near 25 genes involved in pigmentation and immune-related traits, including at a key human protein interactor of 26 SARS-CoV-2. Only in the Historical period do the selection candidates we detect largely mirror 27 previously-reported signals, highlighting how the statistical power of previous studies was limited to the 28 last few millennia. The Historical period also has multiple signals associated with vitamin D binding, 29 providing evidence that lactase persistence may have been part of an oligogenic adaptation for efficient 30 calcium uptake and challenging the theory that its adaptive value lies only in facilitating caloric 31 supplementation during times of scarcity. Finally, we detect selection on complex traits in all three 32 periods, including selection favoring variants that reduce body weight in the Neolithic. In the Historical 33 period, we detect selection favoring variants that increase risk for cardiovascular disease plausibly 34 reflecting selection for a more active inflammatory response that would have been adaptive in the face of 35 increased infectious disease exposure. Our results provide an evolutionary rationale for the high 36 prevalence of these deadly diseases in modern societies today and highlight the unique power of ancient 37 DNA in elucidating biological change that accompanied the profound cultural transformations of recent 38 human history. 39 Main 40 Gene-culture co-evolution-whereby cultural adaptations including technological developments 41 lead to new lifestyles that change selection pressures-have been widely discussed as a potential major 42 driver of genetic adaptation 1 . To date, however, there have been few empirical examples, possibility due 43 to the lack of ancient DNA data in sufficient sample sizes to reveal changes in allele frequencies before 44 and after cultural change. This deficiency can be addressed with large ancient DNA datasets. Several 45 central hypotheses have been put forward regarding how human cultural evolution may have driven 46 human biological evolution 2 . 47 48 The first hypothesis relates to metabolic traits. The advent of agriculture induced a shift toward 49 starch-rich and less diverse diets, which would be expected to lead to selection for loci that more 50 effectively metabolize such diets and address their deficiencies of key nutrients 3 . Farming may have 51 paradoxically also contributed to food scarcity. In times of plenty and food stability, population growth 52 occurred at much faster rates than in the hunting and gathering period. However, these larger populations 53 could also have been subject to periods of famine due to drought, agricultural disease outbreaks, or poor 54 food distribution which might lead to additional selection for reduced caloric demand or more efficient 55 energy metabolism. 56 57 The second hypothesis relates to gene-culture co-evolution associated with immunity. As humans 58 began living in closer proximity to domesticated animals in the Neolithic, they would have been exposed 59 to disease affecting those animals. In the Bronze Age and Historical periods, larger increases in 60 population size as well as population movement occurred due to improved technology and mobility. 61 However, this would also have radically increased the opportunity for transmission of infectious disease 62 and pressures on the immune system to more effectively combat them. The immune system has innate 63 aspects associated with inflammatory processes and adaptive aspects associated with recognition of 64 specific antigens. Making both these arms of the immune system more active can have deleterious coverage ~0.9x) and could have different haplotype structure 21 . To avoid additional biases associated with 125 misestimating allele frequencies with heterogeneous data, we did not include ancient shotgun data or 126 modern data in our analysis. 127 128 We leveraged a model of the demographic history of Europe over the past 10,000 years that has 129 been inferred from ancient DNA studies 4,5,22-29 . Broadly, these studies conclude that most Europeans in 130 this period derive the great majority of their ancestry from three primary ancestral sources which came 131 together in the course of two major demographic transitions corresponding to significant shifts in the 132 archeological record. The first is the transition from hunting and gathering to farming, which was 133 accompanied by a major population transition in central and western Europe. In this period, ancestry from 134 the Mesolithic inhabitants of central and western Europe was largely displaced by ancestry from farmers 135 whose ancestors originated in Anatolia, and who were amongst the first peoples in the world to use 136 agriculture a few thousand years before. This economic and demographic shift began in southeastern 137 Europe after around 6,500 BCE but had spread to the far reaches of Europe as well as Britain by 4,000 138 BCE. The second major demographic transition occurred during the shift from the Neolithic to the Bronze 139 Age with the arrival of Steppe Pastoralists from the Eurasian Steppe. In the subsequent millennia leading 140 to the Historical period, there were subtle shifts in the proportion of Steppe ancestry that largely arose 141 from the homogenizing of populations with different Steppe ancestry proportions. 142 143 We assigned individuals to different groupings based on f4-statistics, time period (based on direct 144 radiocarbon dates or well understood archaeological contexts), and geographic location. We removed 145 individuals that were outliers from each time period and were found to have atypical ancestry of that 146 period based on f4-statistics. . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made wide expectation can be identified as evidence for selection. In this case, the frequency of the allele in 210 population C has risen to a frequency that is ~50% above the expected frequency based on mixture 211 proportions and suggests that natural selection has elevated the frequency of this allele in the time period 212 since admixture. 213 A scan for selection at individual loci 214 To identify candidate selected loci in our dataset, we used our three epoch model and applied a 215 method that utilizes the admixture events that occurred in each epoch. Under neutrality, the allele 216 frequency of an admixed population is expected to be the weighted average of the allele frequencies in the 217 source populations that contributed to the admixture. Significant deviations from this expectation suggest 218 that natural selection has acted at a particular locus (Fig. 2). After correcting for inflation of the test 219 statistic independently in each of the three epochs, we used a cutoff of 5ｘ10 -8 as a genome-wide 220 significance threshold. This is a common significance threshold in genome-wide association studies 221 (GWAS) 30 , and also roughly corresponds to a P value of <0.05 after Bonferroni correction for a 1.2 222 million SNP target set (Methods). Previous work has examined the impact of sample size, the strength of 223 selection, the time that selection has acted, mis-specification of the mixture proportions, and additional 224 unmodeled mixtures in detecting selection using this method and has shown that, after applying a 225 correction for genomic inflation, these issues result in reduced power but not an increased rate of false-226 positives 4,31 . Additional work on the same method with slightly different statistical formulation has 227 confirmed this robustness to deviations from the model 32 . To further study the effect of model 228 misspecification as well as the effect of sample size on our power to detect signals, we carried out two 229 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint additional analyses. First, we examined our model's robustness to mis-estimates of the admixture 230 proportions and found that deviations on the order of 15% resulted in little reduction in power (Extended 231 Data Fig. 2). Second, we found that reduced sample size below 80% of the dataset size used for analysis 232 has a major effect on power to detect selection signals (Extended Data Fig. 3). 233 234 Following a previous strategy used to mitigate false-positives in ancient DNA scans of selection 235 due to biases affecting the sequences aligning to a particular variant, we considered loci to be candidates 236 for selection if at least 3 alleles within 1 Mb of each other and the causal gene significantly deviated from 237 their expected frequency 4 . This distance is also in agreement with a recent study examining the optimal 238 window size for linking GWAS-associated SNPs to causal genes 33 . To determine if functional categories 239 of genes were significantly associated with selection signals, we carried out enrichment analysis using 240 FUMA 34 , which maps SNPs to genes and performs gene set enrichment analysis for GWAS and GO 241 annotations incorporating LD information as well as gene matching by length and conservation scores 242 (Methods). Across all epochs, we discovered a total of 25 regions containing alleles with frequency changes 247 that significantly deviated from genome-wide expectation (Fig. 3 In the Bronze Age, we do not detect evidence for continued selection on candidate variants that 326 are directly associated with a change in diet. Instead, we found evidence for selection at or near genes that 327 affect skin and eye pigmentation. 328 329 The strongest signal is at the allele rs16891982, in the gene SLC45A2, which is known to play a 330 major role in light skin pigmentation, and for which there has been previous evidence for selection 51 . The 331 second strongest signal based on our analysis is in the allele rs11636232 in OCA2/HERC2, which is a 332 primary determinant of light eye color in Europeans 4,52 . 333 334 As in the Neolithic period, we also found several candidate genes involved in immunity beyond 335 those seen in the HLA region. We observed selection at rs10797666 in the major histocompatibility 336 complex class I-related gene MR1, which is an immune sensor of microbial ligands, including 337 Mycobacterium tuberculosis, Streptococcus pyogenes, Salmonella enterica and Escherichica coli 53 . We 338 also find evidence of selection at a number of alleles in genes in the killer-cell immunoglobulin-like 339 receptor locus (KIR gene family), which are expressed on the cell membrane of natural killer cells. KIR 340 receptors interact with major histocompatibility molecules to detect pathogen-infected cells and have a 341 crucial role in host defense. This locus is highly polymorphic across human populations worldwide, and 342 diversity at this locus has been correlated with pathogen load across populations 54 . We also find evidence 343 for selection at the MRGPRX3-4 locus, which includes genes that are key physiological and pathological 344 mediators of itch and related mast cell-mediated hypersensitivity reactions, as well as potential targets to 345 reduce drug-induced adverse reactions 55,56 . A final immune-related candidate is the gene MARK3, which 346 is a host protein that is one of the key interactors with SARS-CoV-2 and is important in mediating the 347 maladaptive host response to COVID-19. The allele under selection appears to be linked to a lead signal 348 for monocyte count, which has now been shown to be important in the pathology of COVID-19 57-59 . 349 There is direct ancient DNA evidence for pathogen infection being a major source of mortality in the 350 Bronze Age. The earliest evidence for Yersinia pestis infections in Europe ascertained from ancient DNA 351 comes from the Bronze Age at times particularly close to or after the arrival of pastoralists from the 352 Eurasian Steppe, from where both of these pathogens have been recovered from humans several millennia 353 prior to their first evidence in Europe 60,61 . Thus, pathogens entering Europe along with Steppe Pastoralists 354 in the Bronze Age could have been a driving force behind changes in these immune related genes. 355 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint 356 We also observed significant frequency deviation from the expectation due to genetic drift in 357 alleles lying on several genes related to cardiovascular disease. One candidate is in the angiotensin gene, 358 AGT, which causes vasoconstriction and increases blood pressure 62 . The protein encoded by AGT is a 359 frequent antagonist in drugs that treat heart disease. Additionally, we also observed another locus that 360 reached significance, rs915843, which is a missense allele in ABCG1, a gene that controls tissue lipid 361 levels and the efflux of cellular cholesterol to HDL 63 . 362 363 Finally, we observed candidates in genes where mutations have been linked to reproduction. One 364 of our significant variants is at rs7188473, a splice mutation in the gene HYDIN. Homozygous carriers of 365 this allele suffer from primary ciliary dyskinesia-5, which affects sperm motility and leads to male 366 infertility 64 . 367 368 More broadly, across this epoch, we find a statistically significant enrichment of signals near 369 genes related to skin, hair, and eye pigmentation (Supplementary Table 3 The variant with the strongest significant deviation from expectation is in the LCT gene, which is 375 responsible for conferring the ability to digest lactose in adulthood in Europeans. This is also consistently 376 the strongest signal of natural selection detected in scans in modern Europeans, and in line with findings 377 in previous publications 65 , this allele appears to have experienced its major change in frequency primarily 378 in the past few thousand years, and not during the Bronze Age when the allele was first introduced in 379 central and western Europe. 380 381 We found a selective signal in DHCR7 (the focal SNP that deviates most from expectation is in 382 an enhancer region several kb upstream of the gene), which governs availability of 7-dehydrocholesterol 383 for conversion to vitamin D3 by the action of sunlight on the skin. Milk is rich in 7-dehydrocholesterol, 384 suggesting that selection on this locus as well as LCT might have been related to the need for increased 385 production of vitamin D 66 . This locus has also been linked to several auto-immune diseases. We also 386 detect evidence for deviation in allele frequency from expectation in the missense variant rs653178 in the 387 gene SH2B3. This allele doubled in frequency from the Bronze Age to the Historical period and is a major 388 risk locus for Celiac disease. The allele we identify as a signal of selection has recently also been shown 389 to be fine-mapped in a GWAS for vitamin D binding 67 . Functional investigation of the effect of the 390 SH2B3 genotype in response to lipopolysaccharide and muramyl dipeptide revealed that carriers of the 391 SH2B3 allele showed stronger activation of the NOD2 recognition pathway. This suggests that SH2B3 392 also plays a role in protection against bacterial infection 68 . 393 394 395 The second strongest signal was in the gene encoding the toll-like receptor locus TLR, which is 396 expressed on the membranes of leukocytes. Variants at this locus have been associated with host immune 397 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made By separating our analysis into different epochs, we were able to examine overlaps in candidate 420 alleles across epochs as well as a previous scan examining modern Europeans from the 1000 Genomes 421 Project. Outside of the HLA region, we found no overlaps of any of the loci discovered in the Neolithic 422 period with any of the other epochs. All of the candidates we discovered in the Historical period had also 423 been discovered in a scan comparing modern Europeans to ancient Europeans 4 . As expected, this scan 424 was largely blind to selection during the European Neolithic, showing the value of direct comparison of 425 groups of ancient DNA samples to study the selection that occurred in this time (Fig. 3). 426 427 428 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made In summary, we tested for polygenic selection in a way that reduces statistical power but is more 501 robust to confounding. We identify a set of variants significantly associated with a trait in a Japanese 502 sample, whose population structure and thus potential for population stratification is uncorrelated to that 503 in our ancient European sample. We then carried out a test for selection by comparing the chi-squared 504 statistic of trait-associated alleles, considered alongside with the direction of change in frequency, to those 505 of random SNPs (Methods). 506 507 Signals of polygenic selection that differ across epochs of European history 508 509 In each of the three epochs, we tested for selection by comparing our selection statistic in trait-510 associated SNPs to a distribution of matched controls resampled 10,000 times. The control was composed 511 of an equal number of SNPs matched for deciles of derived allele frequency, recombination rate 87, and 512 intensity of background selection 88 (Methods). We restricted the 220 total traits in the Biobank of Japan to 513 only those that had at least 20 SNPs significantly (P<1×10 -6 ) associated with the trait. To assess the 514 directionality of genetic change, we polarized our selection statistic to the direction of the effect allele in 515 the GWAS (polarized chi-squared statistic) and asked whether the mean observed polarized statistic for a 516 trait was below the 2.5 or above the 97.5 percentiles of all the matched null samples. In total, we 517 identified 39 traits that reach this level of significance across the three epochs (Fig. 4, Supplementary 518 Table 8). In carrying out this analysis, we checked that null distributions for all traits were approximately 519 normally distributed (Extended Data Fig. 10), and that we had enough variants to prevent the same variant 520 being sampled multiple times (Supplementary Table 9). 521 522 In the hunter-gatherer to farming transition, we detected evidence for selection favoring body 523 mass-decreasing and cholesterol-increasing alleles. We also found evidence for selection on traits related 524 to blood cell biomarkers such as platelet and hemoglobin concentration. 525 526 In the Bronze Age, we detected evidence for selection on alleles associated with some disease 527 endpoints such as hepatitis and ulcerative colitis, as well as several blood and blood-pressure-related 528 traits. We also observed evidence for selection favoring alleles increasing triglyceride levels. 529 530 The vast majority of polygenic adaptation signals we observe were in the Historical period, 531 though several of the traits we identify could be genetically correlated. Importantly, we detect selection 532 favoring alleles that increase rates of heart disease due to myocardial infarction (heart attacks). In 533 addition, we detected evidence for selection on alleles associated with related phenotypes like angina, as 534 well as biomarkers for cardiovascular disease and cardiovascular prescriptions such as beta-blockers. 535 Finally, we observed signals of selection favoring alleles that increase risk for several common auto-536 immune diseases. 537 538 To investigate alternative ways of carrying out these analyses, we repeated these studies by 539 including effect sizes in the polarized chi-squared statistic; we found that the majority of our signals were 540 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint replicated using the magnitude of effect as well as direction (Fig. 4, Supplementary Table 10). The cases 541 of non-replication were largely sub-significant (for example, ~2% vs ~5% for body weight). We also 542 carried out an additional sensitivity analysis by carrying out the same scan, but this time removing SNPs 543 that were in the lowest 27.5% of the chi-squared statistic distribution (Extended Data Fig. 11). These are 544 positions where the direction of frequency change may have been mis-estimated. Therefore, by restricting 545 to sites that are deviating more significantly from expectation, we might increase our confidence in our 546 estimates of the effect direction but reduce power as the number of SNP positions that we could use in the 547 analysis would decrease. This analysis replicated 90% of the original signals across all epochs 548 (Supplementary Table 11). Importantly, we did not detect evidence for natural selection on height, 549 perhaps because of lack of power or because previous analyses may have been confounded by population 550 stratification 4,89 . A final sensitivity analysis we carried out was to repeat the polygenic selection analysis 551 using summary statistics that were identified from a large consortia study of 25 phenotypes carried out 552 using a within-sibling GWAS design 90 . In theory, family-based GWAS designs can control for 553 demographic and indirect genetic effects, but even with the relatively large number of samples, only 6 out 554 the 25 traits had more than 20 SNPs that met our inclusion P value threshold. Out of these 6, only 3 traits 555 overlapped with traits that were also seen in the Biobank of Japan dataset. Largely, the within-sib GWAS 556 data agreed with our previous analysis, with the exception of selection favoring BMI-decreasing alleles in 557 the Neolithic. However, the within-sib analysis was considerably underpowered compared with using the 558 BBJ dataset, as the total number of SNPs were different by an order of magnitude (Supplementary Table  559 13). 560 561 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint composite approach integrating multiple classical selection tests on modern European genomes 91 . We 578 found only one candidate shared between our study and with modern genomes from the ancient DNA 579 based scan in the EN (the HLA region) and two in the BA period (the HLA region and SLC45A2), but 580 9/12 of Mathieson et al.'s 4 candidates were shared with our Historical epoch candidates (Fig. 3). 581 Similarly, Pybus et al. 91 , found none of the candidates we found in the EN epoch, one in the BA epoch, 582 but 7 that match our candidates from the Historical period. A possible explanation for this is that the 583 admixture in Europe over the past 10,000 years has obscured signals of selection that occurred before the 584 immediate past 14 . 585 586 Our approach looking at time trajectories of alleles over a 10,000-year old period also made it 587 possible to assess the impact of alleles in the germline over long time scales. As an example of this, we 588 studied the frequency trajectory of the CCR5-∆32 variant, which in homozygous form provides protection 589 against HIV in European individuals. We studied the frequency of this allele using a proxy SNP 590 (rs73833033) that is in high LD with it. Across the 3 epochs, we find no evidence for selection of this 591 allele (p=0.55, 0.05, 0.34 for the EN, BA, and H epochs respectively) in line with the evidence from 592 modern samples, despite previous reports of selection at this locus 92,93 (Extended Data Fig. 4). 593 594 It is important to recognize that the number of candidates we find in each epoch does not reflect 595 the intensity of natural selection in that time. Rather, many factors feed in to epoch-specific statistical 596 power. Consider the example of LCT: it is possible that 6,000-3,000 years ago, selective pressures on 597 lactase persistence have been stronger than in the Historical period. Here, selection overcame genetic drift 598 and drove the very rare allele to a frequency of several percentage points of the population. Yet, the 599 largest change in allele frequency, from a few percent to the majority allele in northern Europe, only 600 occurred in the Historical period. These are the changes that we are most powered to detect. Another 601 important caveat is that the genomic control and null model we rely on may not be equally adequate in all 602 parts of the genome, particularly in the HLA region where mutation rates, recombination rates, natural 603 selection, and genetic drift are highly atypical 94 . Nevertheless, our estimation of genome-wide admixture 604 proportions using qpAdm suggests that our expected frequencies broadly capture the allele frequency 605 shifts associated with admixture. 606 607 Our results also allow us to interpret our signals of selection in light of archaeological, 608 evolutionary, and biological evidence. In particular, they allow us to test theories about gene-culture co-609 evolution, specifically with regard to hypotheses about how major changes in lifestyle in the last ten 610 millennia in Europe have or have not resulted in signals of genetic adaptations. 611 612 One important set of insights relates to the genomic impact of the major transition from hunting 613 and gathering to farming. A set of alleles that were targets of selection in this period have to do with 614 decreased body weight/size. Complementarily, the archeological record also shows an overall decrease in 615 body size during the Neolithic transition 95 . One hypothesis is that a reduction in overall calorie intake, a 616 trait that would be genetically correlated with reduced body weight, was advantageous in the Neolithic 617 when famines and resource instability became more frequent 96 . Similarly, selection for more efficient 618 storage and use of glucose in tissues during periods of famine or pathogen outbreaks might also underlie 619 several of our selection signals associated with insulin secretion, regulating glucose in the blood stream. 620 621 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made Our results also allow us to re-examine the hypothesis for selection on the lactase persistence 622 allele. A recent study suggested that milk exploitation and consumption started well before the lactase 623 persistence allele began to be selected, and that the selected allele did not show consistent associations 624 with improved fitness or health in modern individuals, perhaps suggesting that the ability to digest lactose 625 into adulthood was only selected for in conditions of food scarcity 97 . While this study was exclusively 626 restricted to just this allele and phenotype, our genome-wide scan adds additional perspective on the 627 rationale for selection at this locus by connecting it to selection at other candidates in the same time 628 period. In the Historical period, along with the LCT locus, we detect selection candidates in SH2B3 and 629 DHCR7-two genes that are directly related to vitamin D binding as well as a candidate in SLC45A2, a 630 major locus of light skin pigmentation in Europeans, a phenotype which also promotes vitamin D 631 synthesis from sunlight. Taking these results together, our results may suggest an alternative to the caloric 632 supplementation hypothesis; namely, that selection acted to increase calcium uptake-via improved 633 vitamin D absorption as well as increased dietary uptake through the consumption of milk-a finding that 634 has also been discussed in another recent study 13 . Vitamin D is almost entirely absent in a plant-based 635 diets, and these results might also help explain the differences in both lactase persistence and skin 636 pigmentation between Northern and Southern Europe, with increased sun exposure in Southern Europe 637 allowing for sufficient vitamin D synthesis despite a similar dietary transition. increasing cardiovascular disease and auto-immune disease risk can also be interpreted in this light. All 653 else being equal, our findings suggest that today, people with hunter-gatherer genomes would have been 654 at lower risk for cardiovascular and autoimmune disease. The high prevalence of cardiovascular disease 655 in modern societies could be in part due to past selection for heightened inflammatory response-the 656 immune system's primary response to harmful stimuli including pathogens 105 . That is, beginning in the 657 Neolithic and intensifying in the subsequent periods, humans were subject to a greater infectious load, 658 and selection for proinflammatory genes and a strong inflammatory function due to the secretion of 659 adipokines, which underlie cardiometabolic diseases, may have resulted in increased risk for 660 cardiovascular disease. 661 662 While we find evidence for two hypotheses concerning gene-culture co-evolution in the last ten 663 thousand years-with selection for traits related to metabolism as well as immune response-we did not 664 have power to detect selection for cognitive or neuro-psychiatric disease traits, due to the limited data and 665 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint relatively small sample size for these traits in the Biobank of Japan data. There is no evidence in the 666 genetic data for selection on such traits, but future larger studies might provide power to detect such 667 signals. 668 669 While our work offers some methodological improvements compared to previous efforts, the 670 greatest improvement in resolution comes from the quality and quantity of ancient DNA data. More 671 ancient genomes are becoming available from different geographic regions and time periods. Extending 672 the type of analysis we report here to these datasets has the potential to enrich our understanding of the 673 history of natural selection on humans beyond what could be learned through analyses of contemporary 674 sample alone, where ancient selective events are obscured by admixture and drift, and where their timing 675 cannot be directly determined. 676 677

678
Ancient DNA data curation 679 680 We obtained ancient DNA sequencing data from the Allen Ancient DNA Resource 681 (https://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-682 ancient-dna-data, version 51), and subsetted the data to only include samples that were enriched for ~1.2 683 million nuclear targets with an in-solution hybridization capture reagent. 684 685 To analyze the data, we began with the raw read data for all of these samples and sorted the read 686 pairs by searching for the expected identification indices and barcodes for each library, allowing up to one 687 mismatch from the expected sequence in each case. We removed adapters and merged together sequences 688 requiring a 15 base pair overlap (allowing up to one mismatch), taking the highest quality base in the 689 merged segment to represent the allele. We mapped the resulting sequences to the hg19 human reference 690 genome 106 using the same command of BWA 107 (version 0.6.1). We removed duplicate sequences 691 (mapping to the same position in the genome and having the same barcode pair), and merged libraries 692 corresponding to the same sample (merging across samples that the genetic data revealed were from the 693 same individual). For each individual, we restricted to sequences passing filters (not overlapping known 694 insertion/deletion polymorphisms, and having a minimum mapping quality 10), and trimmed two 695 nucleotides from the end of each sequence to reduce deamination artifacts. In addition, we also restricted 696 to sequence data with a minimum base quality of 20. 697 698 We assessed evidence for ancient DNA authenticity by measuring the rate of damage in the first 699 nucleotide, flagging individuals as potentially contaminated if they had less than a 3% cytosine-to-700 thymine substitution rate in the first nucleotide for a UDG-treated library and less than a 10% substitution 701 rate for a non-UDG-treated library. We used contamMix to test for contamination based on 702 polymorphism in mitochondrial DNA 108 and ANGSD to test for contamination based on polymorphism 703 on the X chromosome in males 109 , removing individuals with evidence for contamination. For population 704 genetic analysis to represent each individual at each SNP position, we randomly selected a single 705 sequence (if at least one was available). For the selection analysis, in order to obtain read count 706 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint information on a per sample basis, we used BCFtools 110 version 1.3.1 to obtain reference and alternate 707 counts at each genomic position. 708

Principal components analysis
We carried out PCA using the smartpca package of EIGENSOFT 7.2.1 111 . We used default 711 parameters and added two options (lsqproject:YES and numoutlieriter:0) to project the ancient individuals 712 onto the PCA space. We used 991 present-day West Eurasians 22,112,113 as a basis for projection of the 713 ancient individuals. We also computed FST between groups using the parameters inbreed:YES and 714 fstonly:YES. We restricted these analyses to the dataset obtained by merging our ancient DNA data with 715 the modern DNA data on the Human Origins Array and restricted it to 597,573 SNPs. We treated 716 positions where we did not have sequence data as missing genotypes. Fig. 1 shows the PCA of all ancient 717 samples. Extended Data Fig. 1 shows underlying modern samples used for the projection, along with the 718 ancient individuals. 719 720 Extended Data Fig. 1. PCA of ancient samples as well as the basis set of modern samples used in the 721 projection analysis (in grey). 722 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made To estimate the population allele frequencies at each site, we obtained the maximum likelihood 754 estimate from the likelihood of a given frequency p using an approach first described in Mathieson et al. 755 2015. Let p be a reference allele frequency, Ri be the number of reads with the reference allele, Ti be the 756 total number of sequences, N be the number of samples for the population, and ε be a probability of error. 757 Let the binomial probability mass function be denoted as ( , , ) = * ! " + # (1 − ) $%# . Then the 758 likelihood of a frequency p given the read data is 759 760 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made Minimizing the negative log-likelihood function produced the allele frequency estimate for each 763 population at every site. Samples with 0 reference and alternate reads at a site were excluded from the 764 calculation of the maximum likelihood estimate. We used the SLSQP solver from SciPy 115 to minimize 765 the negative log-likelihood function, setting the bounds for the allele frequency at 0.01 and 0.99. We also 766 removed all positions where all reads were missing in any of the populations used in the scan. 767 768 The expected frequency of the target population was also obtained given the mixing proportions frequency of the target population. We tested when the observed allele frequency deviated from 776 expectation using the likelihood ratio test. 777 778 This statistic was used to compute a P value from the χ , ' distribution. To address genomic 781 inflation, a control factor was applied to the statistics such that I median statistic 7.9:; J ' = 0.45494 after 782 removing 49,000 SNPs of functional importance 4,31 . We also removed genomic positions that were 783 covered by >15,000 reads (coverage >10x mean coverage) across our dataset, due to potential mis-784 mapping artifacts. 785 786 Previous work has examined the robustness of this particular approach to mis-estimation of allele 787 frequencies as well as sample sizes, but we added additional power calculations to our particular scenario. 788 First, we carried out an analysis where we modified the ancestry proportions in 5% increments from the 789 actual proportion and again examined the number of 1Mb regions that remained significant according to 790 our criterion after genomic inflation correction. Our results suggest that we are well-powered to detect the 791 majority of our signals even with mis-specification of the admixture proportion by over 30% (Extended 792 Data Fig. 2). 793 794 Second, we carried out a sub-sampling analysis where we down sampled the overall dataset in 5% 795 increments (that is, reducing the sample sizes of both the two source populations and the target population 796 across all 3 epochs in steps of 5%), and then examined the number of 1 Mb regions that remained 797 significant according to our criterion after genomic inflation correction. We see that with 90% of the data, 798 we are essentially recapturing most of our signals, though the lack of a clear plateau in our analysis 799 suggests that increasing sample sizes further is likely to continue to improve the power to discover new 800 loci (Extended Data Fig. 3). Small increases in power are seen even at slightly larger sample sizes, as our 801 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint sampling process is carried out at the level of individuals. This is as expected, as coverage varies greatly 802 by sample and across genomic position, but the overall upward trajectory of increased power with 803 increased sample size is clear. 804 805 806 Extended Data Fig. 2. The power to discover significant genomic regions as a function of admixture 807 proportion mis-specification shown in percent deviation from the actual mixture proportion. Grey line 808 shows the quadratic fitted estimate. 809 810 811 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made and enrichment analysis were performed on the human reference genome build GRCh37 116 . We include 847 significant loci in our selection scan and genomics annotations in a 1 Mb neighborhood using 848 LocusZoom 117 (Extended Data Fig. 6). 849 850 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made Biobank of Japan to test for selection in 220 different traits 120 . For each trait, we classified each allele as 864 trait-increasing or trait-decreasing using the effect direction. We then polarized our admixture scan 865 selection statistic such that a positive sign indicated directional selection of the trait-increasing allele. In 866 other words, for a given loci i, our polarized statistic was computed as | χ & | sgn(β < P ), where | χ & | is the 867 magnitude of the chi-squared statistic from the monogenic selection scan, and sgn(β < P ) is the sign of the 868 effect size for the allele that increased or decreased from expectation. For each trait, we compared the 869 mean polarized statistic of the GWAS significant SNPs (at significance level P<1×10 -6 ) to the distribution 870 of the mean polarized statistics of randomly sampled SNPs (Extended Data Fig. 7). The rationale for this 871 test is that trait-associated SNPs would be more likely than random to undergo short-term directional 872 selection. 873 874 875 Extended Data Fig. 7. Null distribution of random subsamples of matched controls (in blue bars) and the 876 observed statistic (the red line) for a single trait, height, in the Bronze Age. Observed statistics below the 877 2.5th percentile or above the 97.5th percentile of the null were considered significant. 878

879
To investigate the impact of population stratification on the GWAS effect sizes, we regressed 880 GWAS effect sizes on PC loadings on both the UK Biobank (UKB) and Biobank of Japan (BBJ) datasets. 881 In Extended Data Fig. 8, we show the results of this regression on the best studied and most heritable of 882 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made For each trait, we took 100 kb windows and chose only a single SNP with the lowest P value to 898 represent blocks of independent associations with the trait, and we computed the mean polarized statistic 899 of the set of SNPs that were also below a genome-wide significance threshold of 1 x 10 -6 . To match these 900 observed variants to controls, we binned the other variants in the genome based on derived allele 901 frequency, B statistic, and recombination rate. The derived allele frequency bins were separated into 8 902 equally sized bins that ranged from 0 to 1. The B statistic bins were divided by deciles 121 . The 903 recombination rate bin thresholds were computed such that the recombination rates of the SNPs used in 904 the admixture scans were evenly distributed across 8 bins. In carrying out this empirical sampling 905 procedure, we ensured that matched control distributions were normally distributed (Extended Data Fig.  906 9), and that this was the case across different bins (Extended Data Fig. 10). To ensure that we very rarely 907 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint resampled the same alleles each time in our random distributions, we ensured that there were a minimum 908 of 100 variants in each bin (Supplementary Table 9). 909 910 911 Extended Data Fig. 9. Null distributions across all the traits we tested, across different P value 912 thresholds of the GWAS, showed normality in the null distributions and were centered around 0. 913 914 915 Extended Data Fig. 10. Distribution of chi-squared statistics across different bins also show null 916 distributions centered around 0. 917 918 We then randomly sampled variants that were not associated with the trait 10,000 times to match 919 the profiles of the variants with the lowest association P values and computed the mean polarized 920 selection statistic for the sampled variants. We considered there to be directional selection in the trait-921 increasing direction if less than 2.5% of the sampling trials had a mean polarized selection statistic higher 922 than that of the variants with the lowest P values. Similarly, we considered there to be directional 923 selection in the trait-decreasing direction if less than 2.5% of the sampling trials had a mean polarized 924 selection statistic lower than that of the variants with the lowest P values. We only report significant 925 association on traits that had at least 20 SNPs that were significant at a GWAS P value threshold of less 926 than 1×10 -6 ). 927 928 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint We also repeated this analysis by multiplying each variant's chi-squared statistic by its effect size, 929 thereby computing a score that also includes the magnitude of its effect on the trait beyond just looking at 930 its direction. We report all results for this analysis in Supplementary Table 10. Finally, we carried out an 931 analysis reducing the SNPs used by removing sites where the chi-squared value was in the bottom 27.5% 932 of each epoch. We chose to use 27.5% as it was just past the modal value of the distribution at 25%. The 933 positions that were removed by this process are ones where we have a higher likelihood of mis-estimating 934 the direction of frequency change. This analysis replicated the majority (90%) of the original signals. 935 Traits that were not replicated were due to reduction in the overall number of SNPs being reduced to 936 below 20, a condition we required for the sampling process to be reasonable. Only 4 new traits that were 937 not significant previously were seen to be significant, but these were sub-significant (5%) in the original 938 analysis (Supplementary Table 11). 939 940 941 Extended Data Fig. 11. Distribution of the frequency (y-axis) of chi-squared statistic values (x-axis) 942 across epochs. Black lines show the 25th and 50th percentiles and the red line is the 27.5th percentile. 943 In addition to carrying out analysis with the Biobank of Japan dataset, we repeated the analysis 944 with the UK Biobank dataset Supplementary Table 12. We caution that the results of this analysis are not 945 readily interpretable or comparable with those from the Biobank of Japan in light of the issues we discuss 946 with the applications of GWAS with known population stratification artifacts to detect polygenic 947 selection, but we provide these results for completeness. Finally, we reran the polygenic selection analysis 948 using data from the within-sibling GWAS consortium. The GWAS estimates from within-sibling studies 949 are, in theory, better controlled for issues associated with stratification, but the total number of genome-950 wide significant loci that met our threshold limited this approach to only a handful of traits. We report the 951 results of this in Supplementary Table 13. 952 953

954
The code used to run the selection scans for individual alleles and polygenic traits is available at 955 https://github.com/Narasimhan-Lab/1000-genomes-natural-selection. 956 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made This file contains genome-wide selection scan results and allele frequencies for the Historical period. 1330 . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted August 26, 2022. ; https://doi.org/10.1101/2022.08.24.505188 doi: bioRxiv preprint