key: cord-0820651-3bsgx94z
authors: Souilmi, Yassine; Lauterbur, M. Elise; Tobler, Ray; Huber, Christian D.; Johar, Angad S.; Enard, David
title: An ancient coronavirus-like epidemic drove adaptation in East Asians from 25,000 to 5,000 years ago
date: 2020-11-16
journal: bioRxiv
DOI: 10.1101/2020.11.16.385401
sha: 1424e87ec5e234036ada88900dd51753182595d7
doc_id: 820651
cord_uid: 3bsgx94z

The current SARS-CoV-2 pandemic has emphasized the vulnerability of human populations to novel viral pressures, despite the vast array of epidemiological and biomedical tools now available. Notably, modern human genomes contain evolutionary information tracing back tens of thousands of years, which may help identify the viruses that have impacted our ancestors – pointing to which viruses have future pandemic potential. Here, we apply evolutionary analyses to human genomic datasets to recover selection events involving tens of human genes that interact with coronaviruses, including SARS-CoV-2, that started 25,000 years ago. These adaptive events were limited to ancestral East Asian populations, the geographical origin of several modern coronavirus epidemics. An arms race with an ancient corona-like virus may thus have taken place in ancestral East Asian populations. By learning more about our ancient viral foes, our study highlights the promise of evolutionary information to combat the pandemics of the future.

The current SARS-CoV-2 pandemic has emphasized the vulnerability of human 14 populations to novel viral pressures, despite the vast array of epidemiological and 15 biomedical tools now available. Notably, modern human genomes contain evolutionary 16 information tracing back tens of thousands of years, which may help identify the viruses 17 that have impacted our ancestors -pointing to which viruses have future pandemic 18 potential. Here, we apply evolutionary analyses to human genomic datasets to recover 19 selection events involving tens of human genes that interact with coronaviruses, 20 including SARS-CoV-2, that started 25,000 years ago. These adaptive events were Organization, 2019). The most recent outbreak began in late 2019 when SARS-CoV-2 -a less 36 virulent but far more contagious strain than those behind the two previous epidemics -emerged 37 in mainland China before spreading rapidly across the rest of the world, triggering an ongoing 38 pandemic (COVID-19) that so far has infected 45 million people and resulted in over one million 39 deaths worldwide (Dong et al., 2020) . 40 Enrichment is then calculated for each set of top ranked genes using a block-randomization 135 approach that accounts for the genomic clustering of neighboring CoV-VIPs (Enard and Petrov, 136 2020; STAR Methods). For our control gene set, we use protein-coding genes situated at least 137 500kb from CoV-VIPs to avoid overlapping the same sweep signals. Additionally, genes in the 138 control sets are chosen to have similar characteristics as the CoV-VIPs (e.g. similar 139 recombination rates, density of coding and regulatory sequences, percentage of immune genes, 140 percentage of genes that interact with bacteria; see STAR Methods for the complete list of 141 factors) to ensure that any detected enrichment is virus-specific rather than due to a 142 confounding factor (Enard and Petrov, 2020) . Finally, we also exclude the possibility that 143 functions other than viral interactions might explain our results by running a Gene Ontology 144 analysis (Gene Ontology Consortium, 2015; SI; Tables S2, S3 and Figure S1 ). 145

Applying this approach to each of the 26 human populations from the 1,000 genomes dataset, (Table S1 ). The enrichment is most substantial for the top-ranked gene sets ranging 157 between the top 10 and top 1,000 loci ( Figure 1 ; whole enrichment curve P=3.10 -6 for nSL, 158 P=4.10 -3 for iHS, P=6.10 -5 for iHS and nSL combined), and is particularly strong for the top 200 159 loci in large windows (1 Mb) where a four-fold enrichment is observed for both nSL and iHS 160 statistics (pertaining to between 10 to 13 selected CoV-VIPs amongst the top 200 ranked 161 genes; Table S4 ). This suggests that strong selection targeted multiple CoV-VIPs in the 162 common ancestors of modern East Asian populations. That the selected haplotype structures 163 are detected by both the iHS and nSL methods suggests that they are unlikely to have occurred 164 prior to 30,000 years ago, as both nSL and iHS have little power to detect adaptive events 165 arising before this time point in human evolution (Sabeti et al., 2006) . GTEx eQTLs proximal to CoV-VIP genes than expected by chance (iSAFE peak proximity test, 200 P<10-9; STAR Methods). Therefore, for each CoV-VIP gene, we chose a variant with a low 201

Relate p-value (<10 -3 ; STAR Methods) that is situated at or close to a GTEx eQTL associated 202 with the focal gene to estimate the likely starting time of selection for that gene (STAR Methods; 203 Figure S5 ). 204

Using this approach, we observe 42 CoV-VIPs (Table S5 and Figure S5 ) with selection starting 205 times clustered around a peak 870 generations ago (~200 generations wide, potentially due to 206 noise in our estimates; Figure 2 ). While this amounts to about four times more selected CoV-207 VIP genes than were detected using either nSL or iHS (both detected around ten CoV-VIPs 208 amongst the top 200 ranked genes; Table S4 ) this is not unexpected as Relate has more power 209 to detect selection events than nSL and iHS when the beneficial allele is at intermediate starting times forms a highly significant peak (peak significance test P=2.3.10 -4 ; Figure 2 ) when 213 comparing the observed clustering of CoV-VIPs start times with the distribution of inferred start 214 times for randomly sampled sets of genes (STAR Methods). Note that this peak significance test 215 is gene clustering-aware (STAR Methods). Further, this significance test is not biased by the 216 fact that CoV-VIPs are enriched for sweep signals, as the test remains highly significant 217 (P=1.10 -4 ) when using random control sets with comparable high-scoring nSL statistics (STAR 218 Methods). This suggests that the tight temporal clustering of selection events is a specific 219 feature of the CoV-VIPs, rather than a confounding aspect of any gene set similarly enriched for 220 sweeps. The figure shows the distribution of selection start times at CoV-VIPs (pink distribution) 227

compared to the distribution of selection start times at all loci in the genome (blue distribution). 228

Details on how the two distributions are compared by the peak significance test, and how the 229 selection start times are estimated with Relate, are provided in STAR Methods. 230

The genes with clustered selection starting times around 900 generations ago are enriched in 231 strong nSL signals, as shown by running the peak significance test using only CoV-VIPs and 232 controls with strong nSL signals ( Figure S6 ). Conversely, the peak disappears when restricting 233 this test to weaker nSL signals (P=0.53 when using the lowest 50% of nSL statistics; STAR 234 Methods). Importantly, our estimates of the timing of selection are not biased by our use of 235 methods that rely on selected variants not being fixed in the population at the time of genome 236 sampling (i.e. nSL and iHS). When rerunning our analytical pipeline focusing only on strong 237 candidate loci according to Tajima's D (Tajima, 1989 Using this two-interval configuration, CLUES infers that beneficial mutations experienced an 274 exponential increase in frequency between 800 and 500 generations ago for most of the CoV-275

VIPs, which was preceded by a more moderate establishment phase ( Figure 3A Beijing Chinese s = 0.004; Figure 4A , B). These patterns are consistent with the appearance of 282 a strong selective pressure that triggered a coordinated adaptive response across multiple 283 independent loci, which waned through time as the host population adapted to the viral pressure 284 and/or as the virus became less virulent. Although selective pressures other than a coronavirus 285 or related virus might also contribute to these patterns, we note that the signal is restricted 286 specifically at CoV-VIPs and none of 17 other viruses that we tested exhibit the same temporal 287 clustering ~900 generations ago in East Asia (peak significance test P>0.05 in all cases; STAR 288 Methods). Further, this test remained highly significant when retesting the temporal clustering of 289

CoV-VIPs using only other RNA VIPs as the control set (P=4.10-4; Table S1 Beijing CHB 1,000 Genomes population. Left: average selection coefficients between 0 and 500 306 generations ago. Right: average selection coefficients between 500 and 1,000 generations ago. 307

To further clarify that an ancient viral epidemic caused the strong burst of selection we observe 309 ~900 generations ago in the ancestors of East Asians, and not another ecological pressure 310 acting on the same set of genes, we test if the 42 selected CoV-VIPs are enriched for genes 311 with anti-or proviral effects relative to other CoV-VIPs (i.e. loci that are known to have a 312 detrimental or beneficial effect on the virus, respectively). Because the relevant literature for 313 coronaviruses is currently limited -which also applies to the relatively recent SARS-CoV-2 virus 314 -we extend our set of anti-and proviral loci beyond those associated with coronaviruses to 315 include loci reported for diverse viruses with high confidence from the general virology literature 316 (see SI: Host adaptation is expected at VIPs; Table S1 ). We find that 21 (50%) of the 42 CoV-317

VIPs that came under selection ~900 generations ago have high-confidence anti-or proviral 318 effects (vs. 29% for all 420 CoV-VIPs), a significant inflation in anti-and proviral effects 319 (hypergeometric test P=6.10 -4 ) that further supports our claim that the underlying selective 320 pressure was most likely a viral epidemic. 321

Coronavirus infections in humans are known to have pathological consequences for specific 323 bodily tissues, whereby we investigate if the genes targeted by selection in the ancestors of 324

East Asians are also enriched for regulatory functions in similar tissues. In light of our finding 325 that many putative causal mutations in CoV-VIPs were proximal to eQTLs, we investigate 326 whether selected mutations are situated closer to eQTLs for a specific tissue than expected by 327 chance, as this would indicate that the tissue was negatively impacted by the virus (prompting 328 the adaptive response). Briefly, we estimate a proximity-based metric that quantifies the 329 distance between the location of the causal mutation estimated by iSAFE and the tissue-specific 330 eQTLs for the 42 loci that likely started coming under selection ~900 generations ago, and 331 compare this to the same distances observed amongst randomly sampled sets of CoV-VIPs 332 ( Figure 5 ; STAR Methods). 333

Using this approach, we find that GTEx lung eQTLs lie closer to predicted causal mutations 334 amongst the 42 putative selected loci than for any other tissue (P=3.10 -5 ; Figure 5 ). Several 335 additional tissues known to be negatively affected by coronavirus -blood and arteries (Bao et tissue-specific eQTLs than expected by chance ( Figure 5) . Interestingly, the spleen shows no 339 tendency for eQTLs to lie closer to selected loci than expected around 900 generations ago 340 compared to other evolutionary times, perhaps because the spleen is replete with multiple types 341 of immune cells that might be more prone to more regular adaptation in response to diverse 342 The histogram shows how close selection signals localized by iSAFE peaks are to the GTEx 355 eQTLs from 25 different tissues, at peak-VIPs compared to randomly chosen CoV-VIPs (STAR 356

Methods). How close iSAFE peaks are to GTEx eQTLs compared to random CoV-VIPs is 357 estimated through a proximity ratio. The proximity ratio is described in the STAR Methods. It 358 quantifies how much closer iSAFE peaks are to eQTLs of a specific GTEx tissue, compared to 359 random expectations that take the number and structure of iSAFE peaks, as well as the number 360 and structure of GTEx eQTLs into account (STAR Methods). are unable to precisely identify the causal variants for the selected CoV-VIP genes observed in 377 the ancestors of East Asians -nor would these variants necessarily occur as outliers in a 378 GWAS conducted on the British population -we note that it is possible that other variants in the 379 same CoV-VIP genes may also produce variation in SARS-CoV-2 susceptibility and severity 380 amongst modern British individuals. 381

By contrasting variants in CoV-VIPs against those in random sets of genes, we find that variants 382

in CoV-VIPs have significantly lower p-values for both the susceptibility GWAS and severity 383 GWAS than expected (simple permutation test P<10-9 for both GWAS tests; STAR Methods). (Table S6) . This raises the possibility that such drugs could be 402 repurposed for therapeutic use in the current SARS-CoV-2 pandemic. Indeed, an additional six 403 of the 42 selected CoV-VIPs have been identified by (Finan et al., 2017) as part of the 404 "druggable genome" (Table S6) information on a specific gene, the evolutionary approach adopted here is able to leverage 476 evolutionary information embedded in modern genomes to identify candidate genomic regions 477 of interest. This is similar to the information provided by GWAS -i.e. lists of variants or genes 478 that are potentially associated with a particular trait or disease -though we note that the 479 information provided by evolutionary analyses comes with an added understanding about the 480 historical processes that created the underlying population genetic patterns. 481

The current limitation shared by population genomic approaches such as GWAS and the 482 evolutionary analyses presented here, is that they identify statistical associations, rather than 483 causal links, between genomic regions and traits, thereby necessitating additional research to 484 confirm causality. In addition to the various forms of empirical information that we provide here, 485 further evidence of causal relationships between the CoV-VIPs and COVID-19 etiology could be 486 obtained by examining which viral proteins the selected CoV-VIPs interact with, thereby 487 establishing the specific viral functions that are affected. As a preliminary observation, we find 488 that the 35 of the 42 selected SARS-CoV-2 VIPs tend to interact with more viral proteins than 489 expected by chance (SI). Such information will help establish genetic causality and will also 490 improve our understanding how hosts adapt in response to viruses. 491

The ultimate confirmation of causality requires functional validation that the genes interact with 492 the virus, or that drugs targeting these genes have a knock-on impact for the virus. Notably, 493

several CoV-VIP genes are existing drug targets showing the functional importance of these 494 particular loci (Table S6) , several of which are currently being investigated or used to treat 495 severe cases in the current COVID-19 pandemic. It remains to be established if the other genes 496

we have identified in this study might also help guide drug repurposing efforts and provide a 497 basis for future drug and therapeutic development to combat COVID-19 and related 498 pathologies. 499 Conclusion 500 By leveraging the evolutionary information contained in publicly available human genomic 501 datasets, we were able to infer ancient viral epidemics impacting the ancestors of contemporary 502

East Asian populations, which initially arose around 25,000 years ago, resulting in coordinated 503 adaptive changes across at least 42 genes. Importantly, our evolutionary genomic analyses 504 have identified several new candidate genes that might benefit current efforts to combat COVID-505 19, either by providing novel drug targets or by repurposing currently available drugs that target 506 these candidate genes (Tables S4 & S6 ). More broadly, our findings highlight the utility of 507 incorporating evolutionary genomic approaches into standard medical research protocols. 508

Indeed, by revealing the identity of our ancient pathogenic foes, evolutionary genomic methods 509 may ultimately improve our ability to predict -and thus prevent -the epidemics of the future. 510 

We used a dataset of 5,291 VIPs (Table S1 ). Of these, 1,920 of these VIPs are high confidence 518 VIPs identified by low-throughput molecular methods, while the remaining VIPs were identified 519 by diverse high-throughput mass-spectrometry studies. For a more detailed description of the 520 VIPs dataset, please refer to SI: Host adaptation is expected at VIPs. 521

To detect signatures of adaptation in various human populations, we used the 1,000 Genome values within a window of fixed size, centered at the genomic center of genes, halfway between 535 the most upstream transcription start site and the most downstream transcription end site. We 536 then rank the genes according to the average iHS or nSL (more precisely their absolute values) 537 in these windows. We get six rankings for six different fixed window sizes: 50kb, 100kb, 200kb, 538 500kb, 1,000kb and 2,000kb. We do this to account for the variable size of sweeps of different 539 strengths. We then estimate the sweep enrichment at CoV-VIPs compared to controls over all 540 these different window sizes considered together, or at specific sizes, as described below and in 541

Enard & Petrov (Enard and Petrov, 2020) . 542

Estimating the whole ranking curve enrichment at CoV-VIPs and its statistical 543 significance 544

To estimate a sweep enrichment in a set of genes, a typical approach is to use the outlier 545 approach to select, for example, the top 1% of genes with the most extreme signals. Here we 546 use a previously described approach to estimate a sweep enrichment while relaxing the 547 requirement to identify a single top set of genes. Instead of, for example, only estimating an 548 enrichment in the top 100 genes with the strongest sweep signals, we estimate the enrichment 549 over a wide range of top X genes, where X is allowed to vary from the top 10,000 to the top 10 550 with many intermediate values. This creates an enrichment curve as in Figure 1 . Figure 1 shows 551 the estimated relative fold enrichments at CoV-VIPs compared to controls, from the top 1,000 to 552 the top 10 nSL. The statistical significance of the whole enrichment curve can then be estimated 553 by using block-randomized genomes, as described in Enard & Petrov (Enard and Petrov, 2020) . 554

In brief, block-randomized genomes make it possible to generate a large number of random 555 whole enrichment curves while maintaining the same level of clustering of genes in the same 556 candidate sweeps as in the real genome, which effectively controls for gene clustering. 557

Comparing the real whole enrichment curve to the random ones then makes it possible to 558 estimate an unbiased false-positive risk (also known as False Discovery Rate in the context of 559 multiple testing) for the observed whole enrichment curve at CoV-VIPs. A single false positive 560 risk can be estimated for not just one curve but by summing over multiple curves combined, 561 thus making it possible to estimate a single false positive risk over any arbitrary numbers of rank 562 thresholds, window sizes, summary statistics, and populations. For instance, we estimate the 563 false-positive enrichment risk of P=2.10-4 at CoV-VIPs for rank threshold from the top 10,000 to 564 top 10, over six window sizes, for the five East Asian populations in the 1,000 Genomes data, 565 and for both nSL and iHS, all considered together at once. This makes our approach more 566 versatile and sensitive to selection signals ranging from a few very strong sweeps, to many, 567 more moderately polygenic hitchhiking signals. The entire pipeline to estimate false-positive 568 risks with block-randomized genomes is available at 569 https://github.com/DavidPierreEnard/Gene_Set_Enrichment_Pipeline (Enard and Petrov, 2020) . 570

To estimate a sweep enrichment at CoV-VIPs, we compare CoV-VIPs with random control sets 572 of genes selected far enough (>500kb) from CoV-VIPs that they are unlikely to overlap the 573 same large sweeps. We do not compare CoV-VIPs with completely random sets of control 574 genes. Instead, we use a previously described bootstrap test to build random control sets of 575 genes that match CoV-VIPs for a number of potential confounding factors that might explain a 576 sweep enrichment, rather than interactions with viruses. The bootstrap test has been described 577 in detail (Enard and Petrov, 2020) , and is available at 578 https://github.com/DavidPierreEnard/Gene_Set_Enrichment_Pipeline. estimates that the probability is not too small that the frequency of the mutation is different from 608 zero. For our purpose of estimating when selection started, the low time estimate is the best 609 suited, because it provides an estimate of when the frequency of a selected mutation was 610 already high enough to distinguish from zero, for those mutations where selection started from a 611 very low frequency. For cases where selection started with standing genetic variants that were 612 already distinguishable from zero, the Relate low time estimates for the emergence of mutations 613 do not provide a good proxy for when selection actually started. Thus, if we were able to 614 estimate when selection started for standing genetic variants, we might be able to observe an 615 even stronger peak than the one we see when just relying on those variants where selection 616 started from low frequencies. 617

Using the low Relate time estimates is also justified due to the fact that the sweep establishment 618 phase can take very variable amounts of time before the start of the sweep exponential phase. 619

During the establishment phase, selected alleles are still mostly governed by drift which makes 620 pinpointing the actual starting time of selection difficult. In this context, the low Relate time 621 estimates provide an estimate of the time when the selected alleles were no longer at very low 622 frequencies not statistically different from zero, and closer to entering the exponential phase, 623 which provides a more certain time estimate for when selection started for certain. 624

An important step is then to choose at each CoV-VIP locus, and all the other control loci, which 625

Relate mutation to use to get a single time estimate for each locus. Note that here we make an 626 assumption that each locus has experienced only one single adaptive event. Given our finding 627 that iSAFE peaks at CoV-VIPs are much closer to GTEx V8 eQTLs than expected by chance, it 628 is likely that the selected adaptive mutations are regulatory mutations at, or close to annotated 629 eQTLs for a specific gene. They are not necessarily exactly located at eQTLs, because current 630 eQTLs annotations may still be incomplete, and in our case we use eQTLs identified in GTEx 631 V8 using mostly European individuals, even though we analyse selection signals in East Asian 632 populations. Because of these limitations, we use the Relate estimated time at the mutation 633

where Relate estimates the lowest positive selection p-value within 50kb windows centered on 634 eQTLs. We also only consider variants with a minor allele frequency greater than 20%, given CoV-VIPs compared to the null expectation. 676

We then repeat the sliding of a 200 generations window to identify the maximum peak and 677 measure the same difference, but this time for random sets of Relate selected variants of the 678 same size (152 selected variants out of the 1,771 selected variants). To estimate p-values, we 679 then compare the actual observed difference with the distribution of differences generated with 680 one million random samples. 681

As mentioned in the Results, one potential issue is that we run the peak significance test after 682 we already know that CoV-VIPs are enriched for iHS and nSL top sweeps, and especially 683 enriched for nSL top sweeps. This enrichment may skew the null expectation for the distribution 684

of Relate times at CoV-VIPs. In other words, there is a risk that any set of genes with the same 685 sweep enrichment might exhibit the same peak as CoV-VIP. As a result, comparing CoV-VIPs 686 with randomly chosen non-CoV-VIPs may not be appropriate. To test this, we repeat the peak 687 significance test, but this time comparing the peak at CoV-VIPs with the peaks at random sets 688 of non-CoV-VIPs that we build to have the same distribution of nSL ranks as CoV-VIPs. To do 689 this, we define nSL bins between ranks 1 and the highest rank with a rank step of 100 between 690 each bin, and we count how many Relate selected variants fall in each bin (each gene has one 691 nSL rank and one Relate selected variant). To build the random set, we then fill each of the 100 692 bins with the same number of random non-CoV-VIPs, as long as their nSL rank falls within that 693 bin. We use the average nSL rank over the five East Asian populations, and the lower 694 population-averaged rank of either 1 Mb or 2Mb window sizes (where we observe the strongest 695 enrichment at CoV-VIPs, see Results). The results of the peak significance test are unchanged 696 when using the matching nSL distribution (peak significance test P=1.10-4 vs. P=2.3.10-4 697 without matching nSL distribution). 698

In further agreement with the fact that the sweep enrichment does not confound the peak 699 significance test, the peak at CoV-VIPs stands out more when repeating the peak significance 700 test using a smaller nSL top rank limit ( Figure S6) Indeed, the genomic regions at or close to CoV-VIP GTEx eQTLs are likely enriched for CoV-714 VIP regulatory elements, and therefore the most likely place to find CoV-VIP-related adaptations 715 in the genome. To localize where adaptation occurred, we use the iSAFE method that was 716 specifically designed for this purpose (Akbari et al., 2018) . iSAFE scans the genome and 717 estimates a score that increases together with proximity to the actual selected mutation. The 718 higher the score, the higher the odds that the scored variant is itself the selected one, or close 719 to the selected one. An important caveat is that iSAFE is designed to localize where selection 720 happened right after it happened, or as selection is still ongoing. In our case, we have evidence 721 that selection was strong at CoV-VIPs only more than 500 generations (~14,000 years) ago, 722 and then much weaker more recently (Figure 4 ). This could be an issue, because we expect 723 that recombination events that occurred after the strong selection might have deteriorated the 724 iSAFE signal that relies on haplotype structure. This is because recombination mixes together 725 the haplotypes that hitchhiked with the selected mutation, with those that did not. In line with 726 this, we often do not observe simple, clean iSAFE score peaks, but instead, iSAFE score 727 plateaus and more rugged peaks ( Figure S5 ). For this reason, we designed an approach to not 728 only identify the top of simple iSAFE peaks, but also more rugged peaks or plateaus. First, to 729 measure iSAFE scores, we combine all the haplotypes from the five East Asian populations 730 together as input, since we found that the selection signal at CoV-VIPs is common to all these 731 populations (iSAFE parameters: --IgnoreGaps --MaxRegionSize 250000000 --window 300 --732 step 100 --MaxFreq .95 --MaxRank 15). We then use a 500kb window sliding every 10kb to 733 identify the highest local iSAFE value in the 500kb window ( Figure S8 ). Once we have the 734 highest local iSAFE value and coordinate, we define a broader iSAFE peak as the region both 735 upstream and downstream where the iSAFE values are still within 80% of the maximum value 736 ( Figure S8 ). This way, we can better annotate iSAFE plateaus and rugged peaks, and take into 737 account the fact that they can span more than just a narrow local maximum ( Figure S5) . 738

Once the local iSAFE peaks are identified, we can ask how close GTEx eQTLs are to these 739 peaks compared to random expectations. We first measure the distance of each CoV-VIP GTEx 740 eQTL to the closest iSAFE peak. To avoid redundancy, we merge eQTLs closer than 1kb to 741 each other into one test eQTL at the closest, lower multiple of 1,000 genomic coordinates (for 742 example 3,230 and 3,950 would both become 3,000). We then measure the average of the log 743 of the distance between all CoV-VIPs and their closest iSAFE peak. We use the log (base 10) of 744 the distance, because it matters if the eQTL/iSAFE peak distance is 100 bases instead of 745

200kb, but it does not really matter if the distance is 200kb or 600kb, because the iSAFE peak 746 at 300kb is likely not related to the eQTL more than the peak at 600kb. Once we have the 747 average of log-distances, we compare it to its random expected distribution. To get this random 748 distribution, we measure the log-distance between each CoV-VIP eQTL and the iSAFE peaks, 749 but after shifting the iSAFE scores left or right by a random value between 1Mb and 2Mb ( Figure  750 S8; less, or no shift at all if this falls within telomeres or centromeres). We shift by at least 1Mb 751 to make sure that we do not rebuild the original overlap of iSAFE peaks with eQTLs again and 752 again (some iSAFE peaks, or more precisely rugged peaks and plateaus can be wide and 753 include several hundred kilobases; see Figure S5 ). The random shifting effectively breaks the 754 relationship between eQTLs and iSAFE peaks, while maintaining the same overall eQTL and 755 peak structure (and thus variance for the test). The random log-distance distribution then 756 provides an overall random average log-distance to compare the observed average long-757 distance with, as well as estimate a p-value. 758

Then, to more specifically ask if lung eQTLs at CoV-VIPs or the eqTLs of other specific tissues 759 are closer to iSAFE peaks than expected by chance, we can do the same but only using the 760 eQTLs of that specific tissue. The analysis represented in figure 5 is however more complicated 761 than just testing if CoV-VIP eQTLs for a specific tissue are closer to iSAFE peaks than expected 762 by chance by randomly sliding iSAFE values. Instead, what we ask is whether the 42 peak-VIPs 763 have eQTLs for a given tissue that are even closer to iSAFE peaks than the eQTLs of all CoV-764 VIPs in general. To test this, for example with lung eQTLs, we first estimate how close lung 765 eQTLs are to iSAFE peaks at peak-VIPs, compared to random expectations, by measuring the 766 difference between the observed and the average random log-distance, just as described 767 before. We then count the number of peak-VIPs with lung eQTLs (19 out of 25 peak-VIPs with 768

GTEx eQTLs), and we randomly select the same number of any CoV-VIP (which may randomly 769 include peak-VIPs) as long as the random set of CoV-VIPs has the same number of lung eQTLs 770 (plus or minus 10%) as the set of peak VIPs with lung eQTLs (the same gene can have multiple 771 eQTLs for one tissue). We make sure that the tested and the random sets have similar numbers 772 of genes and eQTLs so that the test has the appropriate null variance. We then measure the 773 difference between the observed log-distance, and the randomly expected average log-distance 774 for the random set of CoV-VIPs, exactly the same way we did before for the actual set of peak-775

VIPs. We then measure the ratio of the observed difference in log-distance between peak-VIPs 776 and the random expectation after many random shiftings (1,000), divided by the average of the 777 same difference measured over many random sets of CoV-VIPs. The final ratio tells us how 778 much closer lung eQTLs are to iSAFE peaks at peak-VIPs compared to CoV-VIPs in general, 779 and still takes the specific eQTLs and iSAFE peak structures at each locus into account, since 780 we compare differences in log-distances expected while preserving the same eQTL and iSAFE 781 peak structure (see above the description of the random coordinate shifting). One important last 782 detail about the test is that because we already found that the 50% of loci with the lowest nSL 783 signals do not show a peak of selection at CoV-VIPs around 900 generations ago (see Results), 784 we do not use these loci in this test since any iSAFE peak there is much more likely to represent 785 random noise, not actual selection locations, and thus likely to dilute genuine signals. Using this 786 test, we find that lung and other tissues' eQTLs at peak-VIPs are much closer to iSAFE peaks 787 than they are at CoV-VIPs in general. This test thus specifically tells that adaptation happened 788 closer to lung eQTLs, specifically around 900 generations ago compared to other evolutionary 789 times. By estimating the same ratio for 24 other tissues with at least 10 peak-VIPs with the 790 specific tested tissue eQTLs, we can finally rank each tissue for its more pronounced 791 involvement in adaptation ~900 generations ago, as done in figure 5. It is particularly interesting 792 in this respect that the tissue with least evidence for being more involved in adaptation at that 793 time more than other evolutionary times is spleen. Spleen indeed likely represents a good 794 negative control as a tissue strongly enriched in immune cell types and likely to have evolved 795 adaptively for most of evolution. 796

To compare the UK Biobank GWAS p-values at different loci, we assigned one p-value for each 798 gene, either CoV-VIPs, peak-VIPs or other genes, even though each gene locus can have many 799 variants with associated GWAS p-values. To assign just one single GWAS p-value to each 800 gene, we selected the variant with the lowest p-value at or very close (<1kb) to GTEx eQTLs for 801 a specific gene, in line with the fact that GWAS hits tend to overlap eQTLs (Hormozdiari et al., 802 2016) , and to remain consistent with the rest of our manuscript. We then compared the average 803 p-value between different sets of genes using classic permutations (one billion iterations). 804

We queried the databases DGIdb (Cotto et al., 2017) , and PanDrugs (Piñeiro-Yáñez et al., 806 2018) for drugs targeting CoV-VIPs and peak-VIPs. For hits from PanDrugs we limited the 807 results to only genes that are in direct interaction with the designated drug. Drugs targeting 808 peak-VIPs are presented in Table S6 . In addition, we present a list of peak-VIPs that are not 809 currently drug targets, but have been previously identified in (Finan et al., 2017) 

A global reference for human genetic variation

Identifying the favored mutation in a positive selective sweep

Disparities in Cancer Outcomes Due to 826 COVID-19-A Tale of 2 Cities

SARS-CoV-2 828 induced thrombocytopenia as an important biomarker significantly correlated with abnormal 829 coagulation function

Evolutionary dynamics of human Toll-like receptors and 833 their different contributions to host defense

Coevolution of host and virus: the pathogenesis of virulent 835 and attenuated strains of myxoma virus in resistant and susceptible European rabbits

DGIdb 3.0: a redesign and expansion of the drug-839 gene interaction database

An interactive web-based dashboard to track COVID-843 19 in real time

Enhancing GTEx by bridging the gaps between genotype, gene 845 expression, and disease

Genomewide Association Study of Severe Covid-848 19 with Respiratory Failure

Digestive Manifestations in Patients 851 Hospitalized with COVID-19

Evidence that RNA Viruses Drove Adaptive Introgression 853 between Neanderthals and Modern Humans

Ancient RNA virus epidemics through the lens of recent 855 adaptation in human genomes

Genome-wide signals of positive selection in 857 human evolution

Viruses are a dominant driver of 859 protein adaptation in mammals

On detecting 861 incomplete soft or hard selective sweeps using haplotype structure

The druggable genome and support for target 865 identification and validation in drug development

Gene Ontology Consortium: going forward. Nucleic Acids 867 Res

A SARS-CoV-2 protein interaction map 870 reveals targets for drug repurposing

Analysis 872 of cardiopulmonary findings in COVID-19 fatalities: High incidence of pulmonary artery thrombi 873 and acute suppurative bronchopneumonia

Polygenic adaptation after a sudden change in 875 environment

The landscape of recombination in African 878 Americans

SARS Reference

Colocalization of GWAS and eQTL 882 Signals Detects Target Genes

Ethnicity-dependent allele frequencies are correlated with COVID-19 885 case fatality rate

Gene expression levels are a target of recent natural selection in the human genome

Recent positive selection has acted on genes encoding proteins with more interactions within 891 the whole human interactome

hapbin: An Efficient Program 893 for Performing Haplotype-Based Scans for Positive Selection in Large Genomic Datasets

SARS-CoV-2 infection and obesity: Common inflammatory 896 and metabolic aspects

A 898 genetic method for dating ancient genomes provides a direct estimate of human generation 899 interval in the last 45,000 years

Genetic Ancestry and Natural Selection Drive 902 Population Differences in Immune Responses to Pathogens

Characterization of spike glycoprotein of SARS-CoV-2 on virus entry and its immune cross-905 reactivity with SARS-CoV

PanDrugs: a novel method to prioritize anticancer drug treatments according to individual 909 genomic data

Genetic Adaptation and Neandertal Admixture 912 Shaped the Immune System of Human Populations

Human Immunology through the Lens of Evolutionary Genetics

Clinical Virology

AncestryDNA COVID-19 Host Genetic 918 Study Identifies Three Novel Loci

Positive natural selection in the human 921 lineage

Obesity Is a Risk Factor for 923 Severe COVID-19 Infection

Positive selection of primate 925

TRIM5α identifies a critical species-specific retroviral restriction domain

A multimethod approach for county-scale geospatial analysis of emerging infectious diseases: a 929 cross-sectional case study of COVID-19 incidence in Germany

Background Selection Does Not Mimic the Patterns of Genetic Diversity 931 Produced by Selective Sweeps

Evolutionarily conserved elements in 934 vertebrate, insect, worm, and yeast genomes

A method for genome-wide genealogy 936 estimation for thousands of samples

An approximate full-likelihood method for 938 inferring selection and allele frequency trajectories from DNA sequence data

Disentangling selection on 941 genetically correlated polygenic traits using whole-genome genealogies

UK biobank: an open access resource for identifying the 944 causes of a wide range of complex diseases of middle and old age

selscan: an efficient multithreaded program to 946 perform EHH-based scans for positive selection

Statistical method for testing the neutral mutation hypothesis by DNA 948 polymorphism

The COVID-19 Host Genetics Initiative, a global 950 initiative to elucidate the role of host genetic factors in susceptibility and severity of the SARS-951

CoV-2 virus pandemic

Exploiting selection at linked sites to infer the 953 rate and strength of adaptation

A map of recent positive 955 selection in the human genome

Global Epidemiology of Bat 957 Coronaviruses

Middle East respiratory syndrome coronavirus (MERS-CoV)

The major genetic risk factor for severe COVID-19 is 960 inherited from Neandertals