Impacts of Recurrent Hitchhiking on Divergence and Demographic Inference in Drosophila

Abstract In species with large population sizes such as Drosophila, natural selection may have substantial effects on genetic diversity and divergence. However, the implications of this widespread nonneutrality for standard population genetic assumptions and practices remain poorly resolved. Here, we assess the consequences of recurrent hitchhiking (RHH), in which selective sweeps occur at a given rate randomly across the genome. We use forward simulations to examine two published RHH models for D. melanogaster, reflecting relatively common/weak and rare/strong selection. We find that unlike the rare/strong RHH model, the common/weak model entails a slight degree of Hill–Robertson interference in high recombination regions. We also find that the common/weak RHH model is more consistent with our genome-wide estimate of the proportion of substitutions fixed by natural selection between D. melanogaster and D. simulans (19%). Finally, we examine how these models of RHH might bias demographic inference. We find that these RHH scenarios can bias demographic parameter estimation, but such biases are weaker for parameters relating recently diverged populations, and for the common/weak RHH model in general. Thus, even for species with important genome-wide impacts of selective sweeps, neutralist demographic inference can have some utility in understanding the histories of recently diverged populations.


Introduction
The advancement of DNA sequencing technology, along with computational capacity and methodology, continues to revolutionize the field of population genetics. Harnessing the power of whole genome data sets, researchers have begun to explore a wider variety of evolutionary models. One such model that has received considerable attention recently is a model of recurrent hitchhiking, where genetic diversity at neutral regions is reduced due to repeated selective sweeps at nearby loci. This reduction in diversity has been explored theoretically (Kaplan et al. 1989;Stephan et al. 1992;Wiehe and Stephan 1993) showing that the expected reduction in diversity can be approximated as a function of RHH model parameters: c ¼ 2N e s and k, where N e is the effective population size, s is the selection coefficient, and k is the rate of positively selected substitutions. Subsequent studies have examined such RHH models using forward simulation, focusing attention on how Hill-Robertson interference (HRI; Hill and Robertson 1966) between linked beneficial mutations on different haplotypes reduces the probability of fixation (Gerrish and Lenski 1998;Chevin and Hospital 2008).
The impact of natural selection on genomic diversity may be particularly significant for species with very large population sizes, such as Drosophila melanogaster (e.g., Sella et al. 2009;Langley et al. 2012). In abundant taxa, the population adaptive mutation rate is elevated and the weak influence of genetic drift may allow natural selection to favor alleles with modest selection coefficients. By estimating RHH parameters, Jensen et al. (2008) suggested that selective sweeps may reduce genomic diversity in D. melanogaster to half of neutral levels. While this study implicated a model of relatively strong and infrequent sweeps, the study of Andolfatto (2007) instead favored a model of substantially weaker but more frequent adaptive substitutions. Though both of these studies utilized the same genomic data (synonymous polymorphism data at 137 X-linked loci in D. melanogaster) to infer selection strength and adaptive mutation rate, the methods of the two studies led to distinct conclusions. The Jensen study utilized an ß The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com Approximate Bayesian Computation method to jointly infer adaptive mutation rate and selection strength. The Andolfatto study, however, used a maximum likelihood approach to estimate the product of the adaptive mutation rate and selection strength, followed by a McDonald-Kreitman approach to separate the two (McDonald and Kreitman 1991). While these models should imply strongly different proportions of substitutions driven by positive selection, their alignment with estimates of this quantity from genome-wide Drosophila data is unclear. And likewise, the predictions of each model for the role of HRI and for the fixation of neutral variants have not been investigated. We therefore sought to clarify the relationship between published Drosophila RHH models and adaptive divergence.
If linkage to natural selection substantially impacts Drosophila genetic diversity at neutral sites, the accuracy of demographic inference methods that assume neutrality is not assured. Most sites in the fly genome experience direct functional constraint (Halligan and Keightley 2006), which may lead to an excess of rare alleles from deleterious polymorphisms. Many sites that are not under as much direct selection pressure, such as synonymous sites and middles of short introns, are by definition very close to nonsynonymous sites and other functional sites that may experience natural selection. Selective sweeps could skew the genome-wide allele frequency spectrum, in particular by generating a skew toward rare alleles (Braverman et al. 1995) that may resemble the predictions of recent population growth. In line with these concerns, Schrider et al. (2016) found that the presence of positive selection can bias demographic parameter estimates for a single population's history, and can lead to misidentification of demographic models. However, much interest centers on the inference of demographic parameters between recently diverged populations, and it remains unclear whether Drosophila-like RHH on shorter time-scales is sufficient to bias parameters concerning population divergence times, populationspecific size changes, and migration rates. We therefore use RHH simulations to investigate the impact of RHH on estimation of these parameters.

McDonald-Kreitman Analysis
To estimate the proportion of substitutions in the D. melanogaster genome fixed by natural selection, we applied a genome-wide asymptotic McDonald-Kreitman analysis (McDonald and Kreitman 1991;Messer and Petrov 2013) using the web tool from Haller and Messer (2017). Here, we surveyed 197 genomes from a Zambian population of D. melanogaster (Lack et al. 2015), which is believed to be within the ancestral range of the species (Pool et al. 2012). These genomes are masked for identity by descent, apparent heterozygosity, and recent cosmopolitan admixture. In our analysis, we require any given site to be called in at least 50% of the genomes. We also applied a more conventional McDonald-Kreitman analysis, in which we required the minor allele to be segregating above 10% frequency in our sample. Applying this filter to both putatively neutral and selected site classes should reduce bias from deleterious polymorphisms. To estimate the number of substitutions, we used a Drosophila simulans genome aligned to the D. melanogaster genome (Stanley and Kulathinal 2016).
In this analysis, we estimated the proportion of substitutions driven to fixation by natural selection as a ¼ 1 À DsPn D n P s (Smith and Eyre-Walker 2002). Here, D s is the number of synonymous substitutions, P n is the number of nonsynonymous polymorphic sites, D n is the number of nonsynonymous substitutions, and P s is the number of synonymous polymorphic sites. a was calculated for nine site classes (nonsynonymous, 2-fold synonymous, 3-fold synonymous, 5 0 untranslated regions, 3 0 untranslated regions, intron, intergenic, and RNA-coding) and individually for each major chromosome arm (2 L, 2 R, 3 L, 3 R, and X). Four-fold synonymous sites were evaluated as proxies for neutral evolution. Site classes were taken from flybase.org for release 5.56 of the D. melanogaster genome.

Simulations
In this study, we are interested in the effects of recurrent hitchhiking on demographic inference. To examine this, we ran forward simulations using SLIM version 2.5 (Haller and Messer 2016) to model recurrent hitchhiking. Because full-forward simulations are memory intensive and slow when simulating large populations, it is necessary to rescale simulation simulation parameters. We started by running test simulations to get an idea of the largest population size that we could simulate in a reasonable amount of time. We concluded that diploid populations of 50,000 individuals were a sensible target. This results in 50Â rescaling assuming an effective population size of roughly 2,500,000 (e.g., Duchen et al. 2013). Following the results of Uricchio and Hernandez (2014), we determined that under the RHH models of interest, a size reduction to 50,000 individuals should closely maintain the genetic variation of a nonrescaled population. Further, both algorithms provided in the cited paper yielded near identical scaled parameters when reducing the population size from 2,500,000 to 50,000 individuals. Because of this, we used the simpler Algorithm 1 of Uricchio and Hernandez (2014). The main idea behind this rescaling method is that patterns of genetic diversity are maintained when populationscaled parameters h ¼ 4N e l, q ¼ 4N e r, and c ¼ 2N e s are fixed while N e is varied. As such, if we decrease N e ¼2,500,000 by 50Â to N e ¼50,000, then h, q, and a must be increased by 50Â. The algorithm is laid out in step form below.
Here, the subscripts refer to before rescaling (subscript 0) and after rescaling (subscript 1). s 1 is the selection strength, N 1 is the population size, r 1 is the per base pair per generation per chromosome recombination rate, and L 1 is the simulated locus length.
We ran simulations under two different models of RHH, both of which were estimated from D. melanogaster data. Since the rate of adaptive substitutions and the average selective advantage are highly confounded in terms of their impact on diversity levels, we wanted to examine complementary models. The first model we chose to study is from Jensen et al. (2008). Here, the rate of incoming adaptive substitutions (k) is low (k ¼ 4.2E-11; 2Nk ¼ 2:1E-4) and the average strength of selection (s ¼ 0.002; 2Ns ¼ 10; 000) is high. The second model, from Andolfatto (2007), consisted of a high rate of adaptive substitutions (k ¼ 6.9E-10; 2Nk ¼ 3:45E-3) with a very low average selection strength (s ¼ 1.2E-5; 2Ns ¼ 60). For both models, we used a mutation rate l ¼ 3.27E-9 (4Nl ¼ 0:0327Þ (Schrider et al. 2013) and a recombination rate r ¼ 2.5E-8 (4Nr ¼ 0:25Þ. We also modeled gene conversion at a rate of 6.25E-8/bp/generation. The tract length of the gene conversion was drawn from a geometric distribution with a mean length of 518 base pairs (Comeron 2012). In forward simulations, it is not possible to directly specify the rate of adaptive substitutions. Instead, one must input the rate of beneficial mutations . In the absence of interference among selected mutations, this can be derived using k and the probability of fixation (Kimura 1962): In our simulations, each beneficial mutation had its selection coefficient drawn randomly from an exponential distribution with a mean equal to about half of the rescaled selection strength. The means of these distributions were chosen such that the average selection strength of a fixed mutation is the rescaled selection strength, s 1 (since the more strongly beneficial mutations drawn from this distribution are more likely to fix). The variation in selection coefficients helps to avoid the artificial scenario of interference between mutations with precisely identical fitness. We wanted to analyze ten kilobases (kb) from each simulation for demographic inference. Because a sweep can affect regions far from the target of a sweep, we simulated extra flanking regions for each side of the 10 kb that was used for the demographic inference, while analyzing only the middle region. For the common/weak sweep model, we simulated 480 base pairs on each side of the 10 kb for a total of 10,960 base pairs simulated. For the rare/strong model, we simulated 20 kb flanking loci for a total of 50 kb simulated. Using the formula (2N e sÞ À2r=s (Smith and Haigh 1974; Barton 2010), we expect a 1% reduction in neutral diversity from a sweep 20 kb away under the rare/strong model, while under the common/weak model a sweep should reduce diversity 480 bp away by only 0.03%. Hence, these simulations should incorporate a large majority of sweep effects on neutral diversity predicted by the associated RHH models, while maintaining computational tractability. RHH simulation parameters are provided in table 1.

Demographies Simulated
Forward simulations require a burn-in period to generate appropriate genetic variation. Thus, both recurrent hitchhiking models were run for 500,000 (10N e ) generations. These "trunk" simulations were then used for the demographic simulations. There are two relevant two-population demographic models that we were interested in. These demographies included a bottleneck model and an isolation with migration (IM) model. In our bottleneck model, the populations split and one population experiences a bottleneck. The parameters of the bottleneck model were taken from Thornton and Andolfatto (2006

Demographic Inference
To examine how recurrent hitchhiking affects demographic inference, we used dadi version 1.6.3 (Gutenkunst et al. 2009) to estimate demographic model parameters. We first attempted to fit a two-epoch and three-epoch size change model to the trunk simulations (where there were no size changes simulated) to examine whether recurrent hitchhiking can misidentify demographic models. For the bottleneck simulations, we fit two and three parameter bottleneck models.
The three parameter bottleneck model consisted of a population size reduction, a length of time as the reduced population size, and an instantaneous size change back to the original size that occurs some time in the past. In the two population bottleneck model, the length of the bottleneck is fixed and not optimized. We examined each bottleneck model both with and without fitting the ancestral size change models as well. In this way, we could better parse ancient parameters from more recent demographic parameters post population split. Finally, we tested the IM models with both shared and private selective sweeps. In these cases, the timing of the population split and the migration rate are estimated. As in the bottleneck cases, we tested both IM models with and without ancestral size changes. In total, 14 demographic models were investigated for both hitchhiking models and 11 demographic models were tested for the neutral simulations (since there is no shared/private sweep distinction in the neutral case).
We ran 1,000 simulations of both the common/weak RHH and the rare/strong RHH model. For any given demographic model tested, we randomly chose 50 simulations to generate a site frequency spectrum to run dadi on. dadi was then run ten times on each SFS. The parameters of the dadi run with the highest likelihood across the ten runs were chosen as the inferred parameters. We repeated this process 200 times for each demography tested.

Data Deposition
This study produced no empirical data. All scripts necessary to recapitulate the analyses presented can be found at http:// github.com/jeremy-lange/RHH_project; last accessed July 17, 2018.

Effects of Hill-Robertson Interference
In order to investigate the effects of recurrent hitchhiking on divergence and diversity in Drosophila, we performed forward simulations reflecting two published models of RHH representing relatively common/weak selection (Andolfatto 2007) and rare/strong selection (Jensen et al. 2008), respectively. Before proceeding with further analysis, we checked to see if our initial adaptive mutation rates based on these models were producing the prescribed rates of adaptive substitution, or if instead an important impact of interference must be accounted for. While the above studies assumed no interference between positively selected mutations, Andolfatto (2007) suggested that an interesting next step would be to examine how the presence of interference influences the observed versus expected adaptive substitution rate ðkÞ and selection coefficient of fixed beneficial substitutions (s) in the simulations.
Under a model of no interference, we can expect approximately kLgP fix adaptive substitutions as a product of the adaptive substitution rate at L sites across g generations and their probability of fixation. However, if multiple beneficial mutations are sweeping simultaneously, competition among sweeping haplotypes will lead fewer mutations to fix. Further, mutations that do fix will tend to have a higher selection strength than the input distribution of selective effects. The dynamics of the common/weak and rare/strong models that we tested are very different. In the rare/strong model, on average, we expect 554 beneficial mutations to occur per simulation (over the full 50 kb locus), with 52.2 fixing on average during the 500,000 generations. We expect one beneficial fixation approximately every 9,523 generations. It is unlikely, therefore, that any given beneficial mutation would experience interference from another. In the common/weak model, however, we expect 315,197 beneficial mutations per simulation. Under a model of no interference, we would expect on average 189 of these beneficial mutations to fix. This makes it very likely that more than one beneficial mutation would be sweeping at any given time. Thus, we expect interference to be more significant in the common/weak model relative to the rare/strong model.
To test this expectation, we ran both the common/weak and the rare/strong models as described in the methods section, using the published k and s to generate input parameters. We tracked every mutation that fixed across all simulations and recorded the average selection coefficient. To accurately reflect the RHH models that we were simulating, our goal for our simulations was to approximately match the expected number of fixed beneficial mutations as described above. Across the 1,000 simulations of the rare/strong model, an average of 51.8 beneficial mutations fixed per simulation compared to an expectation of 52.2. For the common/ weak RHH model, we found that, on average, 182 beneficial mutations fixed, corresponding to a modest 3.7% reduction in the expected number of adaptive fixations. We attribute this reduction to interference between positively selected mutations. In order to emulate the properties of the common/weak model, we increased the adaptive mutation rate by 3.7%. This recovered the desired rate of adaptive fixations, averaging 188 adaptive fixations per simulation.
We also conducted simulations to test whether our inclusion of gene conversion was crucial, and we found that it indeed had an important effect on the degree of Hill-Roberton interference. In simulations without gene conversion, the retuned common/weak model averaged 182 adaptive fixations. Thus, the absence of gene conversion reduced the rate of adaptive substitution by 3.3%.
We used the above retuned parameter and ran 1,000 simulations as described in the methods section to more accurately reflect the common weak sweep RHH model. Output of these simulations was used for the simulation analyses detailed below. Note that for a nonrescaled Drosophila population, our adjusted RHH mutational parameters for the Andolfatto (2007) model would correspond to a beneficial mutation rate of ¼ 5.5E-12.

Impacts of RHH on Adaptive Divergence
In these simulations, we can calculate the proportion of substitutions fixed by selection. In the rare/strong RHH model, 1.23% of all fixations are adaptive while in the retuned common/weak model, 17.5% of fixations are adaptive (table 2). These proportions can be compared with the same quantity (a) estimated from extended McDonald-Kreitman analyses of empirical data (Smith and Eyre-Walker 2002). Our genomewide analysis estimated that approximately 18.8% of substitutions in the D. melanogaster genome were driven to fixation by natural selection ( fig. 1; supplementary table S1, Supplementary Material online). Our estimates of a are largely concordant with other studies (Andolfatto 2005, Begun et al. 2007, but see Mackay et al. 2012). Perhaps surprisingly, a estimates for 2 and 3-fold synonymous site classes were very low and hence indicated evolutionary patterns similar to 4-fold synonymous sites. It is unclear why these sites may not have been frequent targets of adaptive protein evolution and, further, it is unclear why a estimates did not converge to a particular range when the minimum allele frequency threshold was increased in the asymptotic McDonald-Kreitman method (Messer and Petrov 2013). For full results, see supplementary table S1, Supplementary Material online for nonasymptotic a estimates with a simple 10% frequency threshold. Overall, this empirical analysis suggests that the common/weak RHH model seems more compatible with adaptive divergence estimates in D. melanogaster.

Hitchhiking Effects on Demographic Inference
The final goal of this study was to examine how models of recurrent hitchhiking affect inferences on demography. We tested whether selective sweep models involving substantial hitchhiking effects would violate assumptions of neutrality made by demographic inference tools and bias parameter estimation. We ran our simulations for 500,000 generations as a single population before we added a population split and distinct demographies. As previously demonstrated (Jensen et al. 2008), the rare/strong model entails a large impact on nucleotide diversity (p reduced by over 50% vs. neutrality; supplementary table S2, Supplementary Material online). This model also generates a notable excess of rare alleles (supplementary fig. S1, Supplementary Material online), in line with theoretical expectations (Braverman et al. 1995). By comparison, the common/weak model reduces p by just 15% and does not produce an excess of rare alleles.
For the "trunk" simulations, we first asked if dadi would fit a model with population size changes over the true model of constant size. Here, we tested both two and three epoch size change models. The two epoch model consisted of two parameters to maximize: a single size change (reduction or expansion) that occurs some time in the past. The three epoch model has three parameters: a size change (reduction or expansion), a length of time at this new size, and some time in the past when the population recovers. The parameters of these models were most affected by the RHH models. Under neutrality, dadi prefers a model of constant population size (the true model) for both the two and the three epoch models ( fig. 2). In the two epoch size change model, dadi infers a population expansion for the rare/strong sweep model, in line with the observed excess of rare alleles. For the three epoch model, a bottleneck ending 1.2N e generations ago is instead favored. Qualititatively similar single-population results have recently been reported (Schrider et al. 2016). However, for common/weak sweep model, estimated population size changes are subtle (<2-fold) and quite ancient (>6N e generations ago), in agreement with the lesser impact of RHH on genetic variation under this model.
Our primary interest was to assess whether the demographic parameters that relate recently diverged populations are similarly biased by RHH. We therefore attempted to infer two distinct bottleneck models occurring after the 500,000 generation burn-in: a two and a three parameter model. Both models consisted of a population size contraction, a bottleneck length, and a length in time since recovery back to original population size. All three parameters are inferred in the three parameter model while bottleneck length is fixed and not inferred in the two parameter model. Results showed that in the two parameter bottleneck model, bottleneck strength was accurately inferred under neutrality and both selection models (supplementary table S3, Supplementary Material online). In the three parameter bottleneck case, bottleneck strength and bottleneck duration were accurately inferred under neutrality and in the common/weak RHH model, with little upward bias in the rare/strong RHH case ( fig. 3; supplementary table S3, Supplementary Material online). Time since recovery was most affected by the presence of selection. Under neutrality, this parameter was accurately recapitulated for both the two and three parameter bottleneck models. Under selection, however, this parameter was overestimated by at least twofold. This result is in line with the fact that both postbottleneck growth and selective sweeps leave behind an excess of rare alleles. Thus, if both are occurring in our simulation but dadi assumes that only neutral events have occurred, it makes sense that this inference method is biased towards overestimating bottleneck recovery times in the presence of positive selection. Models were also run in which ancestral size change parameters and bottleneck parameters were estimated simultaneously, with qualitatively similar results (supplementary table S4, Supplementary Material online).
We simulated two distinct IM models: one with shared selective sweeps between populations, and one simulating In the two-epoch scenario, we allow dadi to infer a single population size change. In the three-epoch scenario, we allow dadi to infer a size change and require it to return to its original effective population size. local adaptation. For both cases, we estimated the time since the population split and the migration rate. It should be noted that the timescale simulated here, 0.5N e , is much longer than the estimated divergence of any D. melanogaster populations (Duchen et al. 2013;Kern and Hey 2017). This scenario, therefore, should be viewed as an extreme scenario in terms of divergence time. Under neutrality, dadi correctly estimated both the divergence time and the migration rate. Under both RHH models, however, these parameter estimates were biased in both the shared sweep and local adaptation cases. In the shared sweep model, dadi overestimated migration rates, with less bias for divergence time ( fig. 4; supplementary table S3, Supplementary Material online). In contrast, for the local adaptation case, the divergence time was greatly overestimated while the migration rate was less biased. As with the bottleneck models, we also made demographic inferences by combining the ancestral size change models with the recent divergence IM models, again with comparable results (supplementary table S4, Supplementary Material online).

Discussion
Simulation provides us with flexible tools for studying the impact of selective sweeps on genetic diversity and divergence. In this study, we have examined two published RHH models and their consequences for divergence and demographic inference. We found that a model of common and weak positive selection appears to fit the genomic estimates of adaptive divergence better than a rare/strong sweep model, and that this favored model appears to entail a slight degree of Hill-Robertson interference. After retuning the common/weak model to yield adaptive substitution rates similar to the published parameters, we investigated the impact of both RHH FIG. 4.-Demographic parameter estimates for simulations of the depicted isolation-migration model are shown. Some bias was observed for this relatively ancient population split time even under neutrality. However, differences from neutral estimates for RHH cases with shared sweeps or local sweeps were consistent with the effects of decreased or increased genetic differentiation, respectively. models on demographic inference. Here, we confirmed that RHH can bias demographic parameters, but we found that the magnitude of such bias was greater for long term population parameters and lesser for more recent parameters.
Although we find agreement between the common/weak model and McDonald-Kreitman estimates of adaptive divergence, we do not argue for a specific quantitative RHH model in light of some important caveats. First, the model that gives a estimates in line with our genome-wide estimates is based on an analysis of nonsynonynous sites specifically (Andolfatto 2007). While this model may share traits in common with the true genomic RHH model for this species, it is best viewed as a qualitative example of the type of model compatible with this aspect of data. Second, we are assessing RHH models based on their general agreement with McDonald-Kreitman-based estimates of a, but such estimates can be biased depending on the history of population size change (Eyre-Walker 2002) and recombination rate change (Comeron et al. 2014). Finally, estimates of a derived from RHH models depend not only on adaptive parameters, but also upon the neutral mutation rate that we use in simulations. The raw mutation rate may be somewhat higher than our simulated l (Schrider et al. 2013;Huang et al. 2016). However, more than half of these mutations should be prevented from fixing by selective constraint (Halligan and Keightley 2006). If the l that we use in our simulations is too high, the predictions for a may be too low, and vice versa. In spite of these quantitative uncertainties, we argue that our analyses provide general insight into the RHH models that are most plausible for this species.
Although this is not an inference study, our results suggest the importance of HRI in a common/weak RHH model like that of Andolfatto (2007). We found that this model led to a loss of approximately 3.7% of beneficial substitutions before retuning the adaptive mutation rate. These results raise the possibility that the effects of interference in D. melanogaster may be stronger than previously suggested (Castellano et al. 2016). In this 2016 study, the authors conclude that approximately 27% of beneficial substitutions are lost due to HRI. This estimate is based on the assumption that there is no interference in regions of recombination that exceed 2 cM/ Mb. However, our results raise the prospect of slight interference in regions of high recombination (2.5 cM/Mb). Therefore, the genomewide effects of interference on the adaptive substitution rate may be somewhat stronger than previously estimated.
Our results also shed light on the impact of RHH on linked neutral variation and its consequences for demographic estimation. It is clear that strong effects of natural selection across the genome can violate the neutralist assumptions of typical demographic inference methods, whether due to background selection (Ewing and Jensen 2016) or selective sweeps (Schrider et al. 2016). Our results support the biasing effect of RHH on ancient parameter estimation and demographic model choice. This is in line with a previous study using single population simulations (Schrider et al. 2016). However, we show that there can be relatively less impact on recent parameters in two population models. Thus, such methods may retain utility even for populous species like D. melanogaster, but the biases that do exist should be borne in mind for this species and investigated for other taxa.
Our simulations emulate a population of D. melanogaster, a populous species with a compact genome. It has been shown that selection may be prevalent in large fly populations (Sella et al. 2009;Karasov et al. 2010;Langley et al. 2012). Further, due to the compactness of the genome, a selective sweep can affect relatively large regions. It is less clear how much of an impact natural selection would have on demographic inference in smaller N e species, such as humans. It is presumed that much of the genetic variation in the human genome is most affected by genetic drift instead of natural selection. Recent studies, however, have argued that an appreciable rate of adaptive substitution has shaped genetic variation in humans (Enard et al. 2014) and that soft sweeps play the dominant role in adaption in human evolution (Schrider and Kern 2017). Thus, even though theory suggests that demographic inference should be more accurate in a less populous species where the effects of natural selection should be lessened, further study is needed to delimit the parameter space in which RHH biases demographic inference.
The simulations in this study entail specific, important caveats. First, we simulated data reflecting the highly recombining portion of the Drosophila genome. The inclusion of low recombination regions would presumably exacerbate the effects of selective sweeps, HRI, and background selection on genetic variation and hence demographic estimates. Our simulations also did not model selective constraint. The biases we observed, therefore, would presumably be worse if analyzed sites had an excess of rare alleles due to deleterious variants. The RHH models that we simulated invoke an arbitrary distribution around a point estimate of s. Importantly, we do not know the true distribution of fitness effects in nature. Theoretical work has suggested that advantageous mutations may be exponentially distributed (Gillespie 1984, Orr 2003, as we have done in this study, although processes such as migration-selection balance may lead to departures from this prediction (Yeaman and Whitlock 2011). There are important caveats to such a distribution in regards to our inference. For instance, the size of the flanking region was determined by the point estimate of the selection strength. Any selective sweep with a selective strength substantially greater than the mean of the distribution will not fully be captured by our analyzed region.
The ABC method used in the Jensen study attempts to capture the locus-to-locus variance in genetic variation to infer adaptive mutation rate and selection strength. Such a method is highly dependent on the size of the genomic regions used in the study. Longer loci entail higher discriminatory power to distinguish between a rare/strong RHH model and a common/ weak RHH model because common/weak selection reduces variance from one locus to the next even in small genomic windows. The loci used in the Jensen and Andolfatto studies were on average 680 base pairs, potentially biasing the ABC method of the Jensen study towards high values of s. The Andolfatto (2007) study utilized a McDonald-Kreitman approach to estimate s, which is not sensitive to the size of locus length. Our estimates of adaptive divergence are in accordance with the results of both the Andolfatto (2007) study and also a more recent study (Keightley et al. 2016). Using D. melanogaster polymorphism data, Keightley suggested a scaled selection strength of N e s ¼ 12. This is very close to our common/weak model's scaled selection strength of N e s ¼ 15 Further, the inferred probability that a new mutation is beneficial from the Keightley study (0.5%) is very close to the retuned estimate from this study (0.4%). Of course, while we draw these conclusions about weak/common RHH models, it is still possible that there could be a number of strong sweeps in the genome.
Our simulations only modeled complete sweeps from new mutations. It has been argued that this is not necessarily the dominant adaptive model in nature (Pritchard et al. 2010;Hernandez et al. 2011;Schrider and Kern 2017), so it is worth considering how other models of natural selection may alter conclusions drawn in this study. If there are two simultaneous soft sweeps, for instance, it is more likely that there exists a haplotype with both favorable variants prior to selection starting, reducing the impact of HRI. Likewise, because soft sweeps have more limited impacts on genetic variation (Pennings and Hermisson 2006), their impact on demographic estimation should be less severe as well.
It is also important to note that complete sweeps may only be a single component of Darwinian selection in nature. Other models of selection may impact genetic variation in fly populations with less effect on divergence, including fluctuating selection (Mustonen and L€ assig 2010;Bergland et al. 2014) and diminishing selection (Vy et al. 2017). Thus, depending on the modes of positive selection that are prevalent in nature, the total impact of hitchhiking on genetic variation and demographic inference may be greater or lesser than simulated here-underscoring the need for further investigation of this topic.

Supplementary Material
Supplementary data are available at Genome Biology and Evolution online.