When viruses impose a selective pressure on a host, adaptation happens in response at 39 the host loci that interact with the virus (Enard et al., 2016 ). An epidemic can then be detected 40 through the enrichment in signals of adaptation at VIPs that interact with a virus. One obstacle 41 to utilizing this method is that we do not have a direct access to the interactions between a host 42 and an ancient virus. However, host-virus interactions are well known in present viruses. of knowledge of ancient viruses can be explained by the difficulty to recover ancient viral DNA 70 (Spyrou et al., 2019) , and by the fact that many viruses have their genomes coded by RNA that 71 is known to degrade much faster than DNA. In light of the scarcity of molecular remains, how 72 can we identify the viruses that drove ancient epidemics during human evolution?
In order to identify viruses that drove ancient epidemics, we can use present interactions 75 as proxies for the interactions with ancient viruses from the same phylogenetic family. Using this 76 approach, we previously found evidence that 50,000 years ago 
Among available haplotype-based statistics, the iHS or integrated Haplotype Score (Voight et 
We further ask whether specific types of viruses drove recent adaptation more frequently 122 than others. More specifically, our previous results on adaptive introgression from Neanderthals 123 to modern humans , and the known zoonoses from the recent past 124 (Geoghegan et al., 2017; Kreuder Johnson et al., 2015) , suggest that RNA viruses may have 125 driven more recent adaptation than DNA viruses because they jump more often from a species Table   139 S1), providing an evenly powered comparison.
First, using all VIPs compared to non-VIPs, we find a strong enrichment of iHS selective 142 sweeps at VIPs overall, suggesting that viruses were a major selective pressure during recent 143 human evolution that drove multiple strong adaptive events. Second, we further find that the 144 selective sweeps enrichment is much more pronounced at RNA VIPs compared to DNA VIPs.
These results validate previous results and are consistent with RNA viruses being a significant 146 selective pressure in the past 50,000 years of human evolution.
We use a dataset of ~4,500 VIPs that are human proteins known to interact physically with 152 viruses (Table S1 ). Of these, we previously manually curated 1,920 from the virology literature 153 (Methods). The articles that report these interactions also frequently report their effect on the 154 viral replication cycle, and 66% of the manually curated interactions also have reported proviral 155 (interaction is beneficial to viral proliferation), or antiviral effects (interaction is detrimental to the 156 virus) (Methods). While this percentage is high, it is still likely to be an underestimate due to the 157 fact that it is limited to the cases where the effect of the interaction has been investigated in the 158 first place. This makes VIPs prime candidates for host adaptation, with a strong likelihood of 159 functional mutations affecting the viral replication cycle. The remaining ~2,600 VIPs were 160 identified using high-throughput methods and were retrieved from the VirHostnet 2.0 database 161 as well as from a number of additional studies (Table S1) . For this analysis, we pooled the 162 manually curated and the high-throughput VIPs together into one group that we systematically 163 compared to the rest of the genome. In total, 20 viruses that infect humans have more than 10 164 VIPs and 14 viruses have more than 100 VIPs, with influenza virus (IAV) and HIV having the 165 highest numbers of known VIPs (1,505 and 1,209 respectively; Table S1 ). Note that VIPs and 
For instance, VIPs are more highly constrained, more highly expressed than non-VIPs, and 
Using this strategy, we find that VIPs are enriched for recent selected sweeps with Figure 1A,B) . Conversely, when we compare DNA VIPs with genes far from them, we do not observe any enrichment of strong selective sweeps (Fig 2B) .
We nevertheless observe some enrichment when using the rank thresholds from top 2,000 to 304 top 500, suggesting weaker sweeps signals at DNA VIPs. These results suggest that RNA with sweeps in control genes that belong to none of the VIP-enriched GO annotations (Fig. 3) . If 337 GO annotations rather than viruses explain our results, we expect the genes from the VIP-338 enriched GO annotations to be significantly enriched for sweeps compared to other genes.
Because the confounding effect should be independent from viruses, we should be able to 
To determine a window size that would more specifically detect strong adaptation in the human 439 genome, we ran population simulations of incomplete selective sweeps and measured the 440 power of iHS to detect weak or strong adaptation when using different window sizes. In 441 particular, we used two window sizes, with small 50kb windows and much larger 1,000kb 442 windows. We ran the population simulations using discoal . Each 
A. Incomplete selective sweeps from a de-novo mutation that reached a 50% frequency. B. Sweep from a 466 standing, 5% standing variant at the start of selection and that reached a 50% frequency. C. and D. same 467 as A. and B., respectively, but for sweeps that reached a 70% instead of 50% frequency. The y axis 468 represents the statistical power (true positive rate) at 0.1% false positive rate (FPR). The x axis 469 represents the range of simulated selection intensities, ranging from 2Ns=50 to 2Ns=1,000 where 470 N=10,000 in our simulations, and s is the selection coefficient. Blue curve: power with 50kb windows.
Orange curve: power with 1,000kb windows. Grey curve: ratio of the power with 1,000kb windows over 472 the power with 50kb windows.
The bootstrap test we used to match VIPs with control non-VIPs that match for multiple 476 confounding factors has already been described extensively in a previous manuscript that the 477 reader can refer to for more ample details (Castellano, 2019 ). An implementation of the 478 bootstrap test is available at https://github.com/DavidPierreEnard, as part of a larger pipeline 479 that also estimates the whole enrichment curve p-value (see below). In brief, the bootstrap test 
It is important to note that we match many (11) confounding factors, which inevitably restricts 527 the number of non-VIPs that can be used as controls. Smaller numbers of control non-VIPs 528 increase the false positive rate of the bootstrap test (see below). However, the False Discovery
Rate analysis strategy that we use to test the significance of the whole curve enrichment fully 530 takes this limitation into account (see below). Note also that in order to limit the number of false 531 positive sweeps due to low recombination regions included in the bootstrap test, we did not use 532 genes with recombination rates estimated to be lower than 0.2 cM/Mb (Hinch et al., 2011) .
Testing the significance of the whole enrichment curve
In addition to estimating a p-value for each iHS rank threshold from top 2,000 to top 20 with the 536 bootstrap test, we also estimated how significantly the whole enrichment curve stands above 537 null expectations. We do this for multiple reasons, some of which already mentioned above.
First, the different thresholds across the whole curve are not independent from each other. As a 539 consequence, an excess of significant p-values at multiple thresholds may reflect a correlation 540 of false positives due to the dependence between thresholds. Second and most importantly, it is 541 very likely that we used the bootstrap test in conditions that make it non-nominal (meaning a 542 non-uniform distribution of p-values under the null hypothesis of no enrichment). Indeed, as we 543 mentioned in the results, we had to exclude many non-VIPs too close to VIPs (<500kb) to avoid 544 counting VIP sweeps as non-VIP sweeps, and vice versa. We also used only non-VIPs that 545 could match VIPs for a large number of confounding factors. The issue then is that we only had 546 a limited number of non-VIPs that were far enough from VIPs and that could be used as 547 controls. For example, for the VIPs vs. non-VIPs test, we could use only 1,704 control non-VIPs.
A smaller set of controls means that in the bootstrap test, the same control non-VIP has to be 549 re-sampled more times (it is a bootstrap), thus resulting in a smaller set of control non-VIPs.
Smaller samples of control non-VIPs are inevitably related to (i) a higher variance of the control 551 null distribution, and (ii) also to a higher variance of the overall average of the null distribution.
Bias (i) can result in a higher rate of false negatives with a decreased power to detect an 553 significant sweep enrichment, but more alarmingly bias (ii) can on the contrary inflate the rate of 
This is a very serious limitation of our analysis, but fortunately there is a very simple, yet 558 computationally intensive remedy. Since we know that the bootstrap test is likely not nominal, 559 we can estimate the true False Discovery Rate associated with the whole approach, by re-560 running the entire analysis pipeline many times on randomized genomes where the iHS ranks 561 have been swapped randomly between genes. We specifically use 10,000 random genomes.
We can then estimate how significant the real enrichment curve from top 2,000 to top 20 is 563 compared to the same curved measured on the 10,000 random genomes. As a statistic to 564 estimate the significance of the whole curve, we use the difference of the observed number of 565 VIPs, minus the expected number according to controls, and we sum this difference over all 566 rank thresholds from top 2,000 to top 20 (Figure 1) . We then compare the real value of the 567 statistic for the real genome, with the distribution of 10,000 random values obtained from the 568 10,000 randomized genomes. For each randomized genome, we use exactly the same VIPs 569 and control non-VIPs as we did when testing the real genome, but these genes are now 570 associated with randomly swapped iHS ranks. The sample size of the control sets of non-VIPs 571 is exactly the same using random genomes as when using the real genome, meaning that the 572 results of the boostrap test reproduce exactly the same biases for the random genomes as 573 when testing the real genome. The p-value obtained for the whole rank thresholds enrichment 574 curve after 10,000 randomized genomes is thus an unbiased, nominal p-value that matches the 575 actual False Discovery Rate.
It is however very important to note that it is true that the randomized genomes provide an 578 unbiased, whole curve test, only because we do not randomize genomes in a completely 579 random fashion. If we just swapped randomly iHS ranks between genes, we would lose the very 580 important property of the real genome that genes with top iHS ranks are likely to be neighbors 581 and to form clusters. Indeed, selective sweeps, and especially strong large selective sweeps, 582 can overlap with multiple genes in the human genome, thus creating clusters of neighboring high-ranking iHS genes. In our case, it was therefore crucial to randomize iHS ranks between 584 genes in a way that conserved the same exact clustering structure between top iHS rank genes.
To achieve this, we cut the genes ordered as they are across human chromosomes in 100 586 blocks of contiguous genes, and then randomly shuffled these blocks. The size of the blocks is 587 much larger than the size of even very large sweeps of multiple megabases, thus ensuring that 588 the clustering structure of the top iHS rank genes is preserved within the blocks. Randomly 589 swapping the blocks still results in randomly swapping iHS ranks between genes within the 590 blocks. Because the randomized genomes obtained this way preserve the clustering of iHS 591 signals, we can still measure enrichments using the number of genes, instead of having to count 592 the number of sweeps. Note that the entire pipeline to get the whole enrichment curve p-value is 593 available at https://github.com/DavidPierreEnard with a user manual, together with the bootstrap 594 test.
To test whether host intrinsic functions can explain the sweep enrichment at VIPs and 598 particularly at RNA VIPs, we measured enrichment for sweeps among genes that belong to GO 599 functions that are over-represented in VIPs. Specifically, we tested whether genes that belong 600 to over-represented functions in VIPs are enriched for sweeps completely independently of 601 whether they interact with viruses or not, by completely excluding VIPs and genes at less than 602 500kb from VIPs (and more likely to overlap VIP sweeps) from the analysis. Thus, we used the 603 bootstrap test to compare iHS sweeps in 718 non-VIPs more than 500kb from VIPs that belong 604 to over-represented GO functions in VIPs (GO functions had to be found in 50 VIPs or more, 605 100% over-represented or more), with sweeps in 423 non-VIPs more than 500kb from VIPs that VIPs with no over-represented GO function. We ran the bootstrap test testing for a deficit of 611 sweeps in the 423 non-VIPs out of over-represented GO functions using the larger number of 612 non-VIPs with at least an over-represented function (the larger the pool of controls compared to 613 the tested set, the better). We also flipped the two groups when testing GO functions for RNA 614 VIPs for the same reason that it is always better to have the larger group being the control 615 group in the bootstrap test. We estimated the significance of the bootstrap test by measuring the 616 whole enrichment curve p-value as already described above. We used the same confounding factors to match during the bootstrap test as we previously did to match VIPs with control non- Table S1 . We found that 66% of the interactions with VIPs 