Modeling the spatiotemporal spread of beneficial alleles using ancient genomes
Abstract
Ancient genome sequencing technologies now provide the opportunity to study natural selection in unprecedented detail. Rather than making inferences from indirect footprints left by selection in presentday genomes, we can directly observe whether a given allele was present or absent in a particular region of the world at almost any period of human history within the last 10,000 years. Methods for studying selection using ancient genomes often rely on partitioning individuals into discrete time periods or regions of the world. However, a complete understanding of natural selection requires more nuanced statistical methods which can explicitly model allele frequency changes in a continuum across space and time. Here we introduce a method for inferring the spread of a beneficial allele across a landscape using twodimensional partial differential equations. Unlike previous approaches, our framework can handle timestamped ancient samples, as well as genotype likelihoods and pseudohaploid sequences from lowcoverage genomes. We apply the method to a panel of published ancient West Eurasian genomes to produce dynamic maps showcasing the inferred spread of candidate beneficial alleles over time and space. We also provide estimates for the strength of selection and diffusion rate for each of these alleles. Finally, we highlight possible avenues of improvement for accurately tracing the spread of beneficial alleles in more complex scenarios.
Editor's evaluation
This is an important manuscript that presents an elegant framework to infer the dynamics of beneficial alleles over time and space. The authors present a new method and show, convincingly, its utility and great potential to reconstruct the evolutionary history of beneficial alleles. The method is also applied to loci that likely mediated human genetic adaptations, contributing to our understanding of human recent evolution. The work will be of broad interest to evolutionary biologists who seek to understand the dynamics of beneficial mutations in populations.
https://doi.org/10.7554/eLife.73767.sa0eLife digest
Analyzing the genomes of our ancient ancestors can reveal how certain traits spread through the human population over the course of evolution. Mutations that make individuals better equipped to survive their environment are more likely to be passed on to the next generation and become more common. For example, a genetic variant that enables adult people to digest sugars in dairy products has become more common in humans over time. Yet evolution does not only happen across time: it transverses space as well.
Modeling the geographic spread of such genetic mutations is challenging using existing methods. To overcome this, Muktupavela et al. developed a new computational method that uses modern and ancient human genomes to study the evolution of specific genetic variants across space and time. The tool can determine where certain variants first emerged, how quickly they spread across geographic areas, and how rapidly they became prevalent in human populations.
Muktupavela et al. applied their new method, which was based on a previously published framework, to track the spread of two common genetic variations that have previously been reported to be subject to natural selection: one that allows adult humans to digest dairy products, and another associated with skin pigmentation. They found that the mutation that enabled dairy consumption originated around what is now southwestern Russia or eastern Ukraine. The variation then spread westward, becoming increasingly more common over the course of the Holocene.
The mutation related to skin pigmentation emerged further south than the dairyrelated variation, and then also spread westward. Massive human migrations during the Neolithic and Bronze Age eras may have helped disperse both variants.
The model developed by Muktupavela et al. could help scientists track the geographic spread of other genetic variants in human populations, as well as provide new insights into how humans adapt to changing environmental conditions. Incorporating major events into the model, like mass migrations or glacial retreats, may lead to even more insights.
Introduction
Understanding the dynamics of the spread of a beneficial allele through a population is one of the fundamental problems in population genetics (Ewens, 2012). We are often interested in knowing the location where an allele first arose and the way in which it spread through a population, but this is often unknown, particularly in natural, nonexperimental settings where genetic sampling is scarce and uneven.
Patterns of genetic variation can be used to estimate how strongly natural selection has affected the trajectory of an allele and to fit the parameters of the selection process. The problem of estimating the age of a beneficial allele, for example, has yielded a rich methodological literature (Slatkin and Rannala, 2000), and recent methods have exploited finescale haplotype information to produce highly accurate age estimates (Mathieson and McVean, 2014; Platt et al., 2019; Albers and McVean, 2020). In contrast, efforts to infer the geographic origins of beneficial mutations are scarcer. These include Novembre et al., 2005, who developed a maximum likelihood method to model the origin and spread of a beneficial mutation and applied it to the CCR5Δ32 allele, which was, at the time, considered to have been under positive selection (Stephens et al., 1998; Sabeti et al., 2005; Novembre and Han, 2012). Similarly, Itan et al., 2009 developed an approximate Bayesian computation (ABC) approach using demic simulations, in order to find the geographic and temporal origins of a beneficial allele, based on presentday allele frequency patterns.
As ancient genome sequences become more readily available, they are increasingly being used to understand the process of natural selection (see reviews in Malaspinas et al., 2012; Dehasque et al., 2020). However, few studies have used ancient genomes to fit spatial dynamic models of the spread of an allele over a landscape. Most spatiotemporal analyses that included ancient genomes have used descriptive modeling in order to learn the spatiotemporal covariance structure of allele frequencies (Segurel et al., 2020) or hidden ancestry clusters (Racimo et al., 2020b), and then used that structure to hindcast these patterns onto a continuous temporally evolving landscape. In contrast to descriptive approaches, dynamic models have the power to infer interpretable parameters from genomic data and perhaps reveal the ultimate causes for these patterns (Wikle et al., 2019).
Dynamic models can also contribute to ongoing debates about the past trajectories of phenotypically important loci. For example, the geographic origin of the rs4988235(T) allele—upstream of the LCT gene and associated with adult lactase persistence in most of Western Eurasia (Enattah et al., 2002)—remains elusive, as is the way in which it spread (an extensive review can be found in Ségurel and Bon, 2017). The allele has been found in different populations, with frequencies ranging from 5% up to almost 100%, and its selection coefficient has been estimated to be among the highest in human populations (Bersaglieri et al., 2004; Enattah et al., 2008; Tishkoff et al., 2007). However, the exact causes for its adaptive advantage are contested (Szpak et al., 2019), and it has been suggested that the selection pressures acting on the allele may have been different in different parts of the continent (Gerbault et al., 2009). Ancient DNA evidence shows that the allele was rare in Europe during the Neolithic (Burger et al., 2007; Gamba et al., 2014; Allentoft et al., 2015; Mathieson et al., 2015) and only became common in Northern Europe after the Iron Age, suggesting a rise in frequency during this period, perhaps mediated by gene flow from regions east of the Baltic where this allele was more common during the onset of the Bronze Age (Krüttli et al., 2014; Margaryan et al., 2020). Itan et al., 2009 deployed their ABC approach to model the spatial spread of the rs4988235(T) allele and estimated that it was first under selection among farmers around 7500 years ago possibly between the Central Balkans and Central Europe. Others have postulated a steppe origin for the allele (Allentoft et al., 2015), given that the rise in frequency appears to have occurred during and after the Bronze Age migration of steppe peoples into Western Eurasia (Haak et al., 2015; Allentoft et al., 2015). However, the allele is at low frequency in genomes of Bronze Age individuals associated with Corded Ware and Bell Beaker assemblages in Central Europe who have high steppe ancestry (Mathieson et al., 2015; Margaryan et al., 2020), complicating the story further (Ségurel and Bon, 2017).
The origins and spread dynamics of largeeffect pigmentationassociated SNPs in ancient Eurasians have also been intensely studied (Ju and Mathieson, 2020). Major loci of large effect on skin, eye, and hair pigmentation have been documented as having been under recent positive selection in Western Eurasian history (Voight et al., 2006; Sabeti et al., 2007; Pickrell et al., 2009; Lao et al., 2007; Mathieson et al., 2015; Alonso et al., 2008; Hudjashov et al., 2013). These include genes SLC45A2, OCA2, HERC2, SLC24A5, and TYR. While there is extensive evidence supporting the adaptive significance of these alleles, debates around their exact origins and spread are largely driven by comparisons of allele frequency estimates in population groups, which are almost always discretized in time and/or space. Among these, selection at the TYR locus is thought to have occurred particularly recently, over the last 5000 years (Stern et al., 2019), driven by a recent mutation (Albers and McVean, 2020) that may have spread rapidly in Western Eurasia.
Here, we develop a method to model the spread of a recently selected allele across both space and time, avoiding artificial discretization schemes to more rigorously assess the evidence for or against a particular dispersal process. We begin with the model proposed by Novembre et al., 2005 and adapt it in order to handle ancient lowcoverage genomic data and explore more complex models that allow for both diffusion and advection (i.e., directional transport) in the distribution of allele frequencies over space, as well as for a change in these parameters at different periods of time. We apply the method to alleles in two of the aforementioned loci in the human genome, which have been reported to have strong evidence for recent positive selection: LCT/MCM6 and TYR. We focus on Western Eurasia during the Holocene, where ancient genomes are most densely sampled, and infer parameters relevant to the spread of these alleles, including selection, diffusion and advection coefficients.
Results
Summary of model
We based our statistical inference framework on a model proposed by Novembre et al., 2005 to fit allele frequencies in two dimensions to presentday genotype data spread over a densely sampled map. We extend this model in several ways:
We incorporate temporally sampled data (ancient genomes) to better resolve changes in frequency distributions over time.
We make use of genotype likelihoods and pseudohaploid genotypes to incorporate lowcoverage data into the inference framework.
We permit more general dynamics by including advection parameters.
We allow the selection, advection, and diffusion parameters to be different in different periods of time. Specifically, to reflect changes in population dynamics and mobility before and after the Bronze Age (Loog et al., 2017; Racimo et al., 2020a), we partitioned the model fit into two time periods: before and after 5000 years BP.
We explored the performance of two different spread models, which are extensions of the original model by Novembre et al., 2005, hereby called model A. This is a diffusion model containing a selection coefficient $s$ (determining the rate of local allele frequency growth) and a single diffusion term ($\sigma $). A more general diffusion model—hereby model B—allows for two distinct diffusion parameters for latitudinal (${\sigma}_{y}$) and longitudinal (${\sigma}_{x}$) spread. Finally, model C is even more general and includes two advection terms (v_{x} and v_{y}), allowing the center of mass of the allele’s frequency to diverge from its origin over time. The incorporation of advection is meant to account for the fact that population displacements and expansions could have led to allele frequency dynamics that are poorly explained by diffusion alone.
In order to establish a starting time point for our diffusion process, we used previously published allele age estimates obtained from a nonparametric approach leveraging the patterns of haplotype concordance and discordance around the mutation of interest (Albers and McVean, 2020). In the case of the allele in the LCT/MCM6 region, we also used age estimates based on an approximate Bayesian computation approach (Itan et al., 2009).
Performance on deterministic simulations
To characterize the accuracy of our inference method under different parameter choices, we first generated deterministic simulations from several types of diffusion models. First, we produced an allele frequency surface map with a specified set of parameters from which we drew 1040 samples matching the ages, locations, and genotype calling format (diploid vs. pseudohaploid) of the 1040 genomes that we analyze below when studying the rs1042602(A) allele.
We generated six different simulations with different diffusion coefficients and afterward ran our method assuming model B. The results (simulations B1−B6) are summarized in Figure 1, Figure 1—figure supplements 1–5, and Appendix 2—table 1. Overall, the model is more accurate at correctly inferring the parameters for the time period before 5000 years BP (Figure 1b), with decreased performance when longitudinal diffusion is high (Figure 1—figure supplement 5).
Next, we investigated the performance of model C, which includes advection coefficients. We generated four different simulations including advection (simulations C1−C4: Figure 2, Figure 1—figure supplements 1–3, and Appendix 2—table 2). We found that our method is generally able to estimate the selection coefficient accurately. However, in some of the simulations, we found discrepancies between the estimated and true diffusion and advection coefficients, often occurring because of a misestimated origin forcing the other parameters to adjust in order to better fit the allele frequency distribution in later stages of the allele’s spread (Figure 2). Despite the disparities between the true and inferred parameter values, the resulting surface plots become very similar as we approach the present, suggesting that different combinations of parameters can produce similar presentday allele frequency distributions.
Advection model applied to nonadvection simulations
We assessed model performance when we apply model C, which includes advection coefficient estimates, to simulations generated without advection (see Figure 1—figure supplements 6 and 7). We can observe that the advection coefficients are inferred to be nonzero (Figure 1—figure supplements 6b and 7b); however, the inferred allele frequency dynamic plots closely resemble the ones obtained with true parameter values (Figure 1—figure supplements 6a and 7a). This shows that complex interactions between the diffusion and advection coefficients can result in similar outcomes even when only diffusion is considered in the model.
The inference of the origin of the allele also differs when we compare results when using models B and C. In order to understand better how the model estimates the allele origin, we highlighted the first individual in simulations B1 and B4 that contains the derived allele. We can see that in the case of simulation B1 the inferred origin of the allele is close to the first observance of the derived allele in the model that includes advection. In contrast when the advection is not included, the origin of the allele is inferred to be closer to where it is initially rising in frequency (Figure 1—figure supplements 1a and 4a). However, this is not always the case. For instance, if we look at the results from the advection model on simulation B4, we can see that the origin of the allele is inferred relatively far from the sample known to have carried the first instance of the derived allele. Therefore, if there is a relatively large interval between the time when the allele originated and when the first ancient genomes are available, the beneficial allele can spread widely, but as this spread is not captured by any of the data points, inference of the precise origin of the selected allele is nearly impossible.
Impact of sample clustering on parameter estimates
We evaluated the impact of different sampling and clustering schemes on our inferences that could potentially arise by aggregating aDNA data from studies with different sampling schemes. We used a deterministic simulation to create three different degrees of clustering, which we will refer to as ‘homogeneous,’ ‘intermediate,’ or ‘extreme’ by varying the area from which we sample individuals to be used in our inferences (Figure 3—figure supplement 1). Additionally, we also tested the impact of biased temporal sampling in the periods before and after 5000 years BP by oversampling in the ancient period (75%/25%), equally sampling in the two periods (50%/50%), and oversampling in the recent period (25%/75%). Because we evaluated this temporal bias for each of the three spatial clustering sampling scenarios, this resulted in a total of nine different sampling scenarios. We note that the third ‘extreme’ spatial clustering scenario is completely unrealistic and one would not expect inferences of any degree of accuracy from it, but we believe it gives a good idea of the behavior of our method in the limiting case of extremely restricted spatial sampling.
A comparison of allele frequency maps generated using true parameter values and using parameter estimates from the different sampling schemes is shown in Figure 3—figure supplements 2–9. In Figure 3 we show the allele frequency map generated using the ‘intermediate 75%/25%’ clustering scheme. Parameter estimates used to generate all these figures are summarized in Appendix 2—table 3. Overall we can see that the allele frequency maps inferred from these scenarios closely resemble the maps generated using the true parameter values, despite the challenges in finding accurate values for the individual point estimates of some of the parameters, highlighting that various combinations of diffusion and advection coefficients can produce similar underlying frequency maps (as discussed in the section ‘Performance on deterministic simulations’). This suggests that the joint spatiotemporal information encoded in the inferred maps (not just the individual parameters estimates) should be used in interpreting model outputs, particularly when it comes to the advection and diffusion parameters. The selection coefficient estimates are inferred highly accurately, regardless of the sampling scheme chosen, and lie close to the true value, with only a slight underestimation in the time period after 5000 years BP (with the exception of the ‘extreme 25%/75%’).
Spatially explicit forward simulations
In addition to drawing simulated samples from a diffusion model, we used SLiM (Haller and Messer, 2019) to perform spatially explicit individualbased forwardintime simulations of selection acting on a beneficial allele by leveraging an R interface for spatial population genetics now implemented in the R package slendr (Petr, 2021).
We introduced a single beneficial additive mutation in a single individual and let it evolve across the European landscape. Before applying our method on the simulated data, we sampled 1040 individuals whose ages were loguniformly distributed to ensure that there were more samples closer to the present, as in the real data. We transformed the diploid genotypes to pseudohaploid genotypes by assigning a heterozygous individual an equal probability of carrying the ancestral or the derived genotype. The parameter values estimated by our model to the simulations described in this section are summarized in Appendix 2—table 4.
We can see that the origin of the allele inferred by the model closely corresponds to the first observation of the derived allele in the simulation (Figure 4). The inferred selection coefficient is only slightly higher than the true value from the simulation (0.0366 vs. 0.030). In general, the model accurately captures the spread of the allele centered in Central Europe, though we observe some discrepancies due to differences between the model assumed in the simulation (which, e.g., accounts for local clustering of individuals, Figure 4—figure supplement 1), and that assumed by our diffusionbased inference.
Dynamics of the rs4988235(T) allele
Having tested the performance of our method on simulated data, we set out to infer the allele frequency dynamics of the rs4988235(T) allele (associated with adult lactase persistence) in ancient Western Eurasia. For our analysis, we used a genotype dataset compiled by Segurel et al., 2020, which amounts to 1434 genotypes from ancient Eurasian genomes individuals, and a set of 36,659 genotypes from presentday Western and Central Eurasian genomes (Ségurel and Bon, 2017; Heyer et al., 2011; Marchi et al., 2018; Liebert et al., 2017; Gallego Romero et al., 2012; Itan et al., 2010; Charati et al., 2019). After filtering out individuals falling outside of the range of the geographic boundaries considered in this study, we retained 1332 ancient individuals. The locations of ancient and presentday individuals used in the analysis to trace the spread of rs4988235(T) are shown in Figure 5.
We used a twoperiod scheme by allowing the model to have two sets of estimates for the selection coefficient and the diffusion and advection coefficients in two different periods of time: before and after 5000 years ago, reflecting the change in population dynamics and mobility before and after the Bronze Age transition (Loog et al., 2017; Racimo et al., 2020a). We used two allele age estimates as input: a relatively young one (7441 years ago) obtained by using the estimated start of selection onset from Itan et al., 2009 (though we note this is necessarily a lower bound of the age of mutation origin), and a relatively old one (20,106 years ago) obtained from the age estimate from Albers and McVean, 2020. The results obtained for fitting the model on rs4988235(T) are summarized in Appendix 2—table 5 and Appendix 2—table 6, and in Figure 6b (younger age) and Figure 6—figure supplement 1 (older age).
Assuming the mutation age estimate is equivalent to the start of selection onset from Itan et al., 2009, the origin of the allele is estimated to be north of the Caucasus, around what is now southwestern Russia and eastern Ukraine (Figure 6b). Given that this age is relatively young, our method fits a very strong selection coefficient ($\approx 0.1$) during the first period in order to accommodate the early presence of the allele in various points throughout Eastern Europe, and a weaker (but still strong) selection coefficient ($\approx 0.03$) in the second period. We also estimate stronger diffusion in the second period than in the first, to accommodate the rapid expansion of the allele throughout Western Europe, and a net westward advection parameter, indicating movement of the allele frequency’s center of mass to the west as we approach the present.
Assuming the older age estimate from Albers and McVean, 2020, the origin of the allele is estimated to be in the northeast of Europe (Figure 6—figure supplement 1), which is at a much higher latitude than the first occurrence of the allele, in Ukraine. Due to the deterministic nature of the model, the frequency is implicitly imposed to expand in a region where there are no actual observed instances of the allele. The model compensates for this by placing the origin in an area with a lower density of available aDNA data and thus avoiding an overlap of the increasing allele frequencies with individuals who do not carry the derived rs4988235(T) allele (see Figure 6a). As the model expands rapidly in the southern direction (Appendix 2—table 6) it eventually reaches the sample carrying the derived variant in Ukraine.
Dynamics of the rs1042602(A) allele
Next, we investigated the spatiotemporal dynamics of the spread of an allele at a pigmentationassociated SNP in the TYR locus (rs1042602(A)), which has been reported to be under recent selection in Western Eurasian history (Stern et al., 2019). For this purpose, we applied our method to the Allen Ancient DNA Resource data (Reich and Mallick, 2019), which contains randomly sampled pseudohaploid genotypes from 1513 published ancient Eurasian genomes (listed in Supplementary file 1), from which we extracted those that had genotype information at this locus in Western Eurasia. We merged this dataset with diploid genotype information from highcoverage presentday West Eurasian genomes from the Human Genome Diversity Panel (HGDP) (Bergström et al., 2020), which resulted in a total of 1040 individuals with genotype information at rs1042602, which were used as input to our analysis. Geographic locations of individuals in the final dataset are shown in Figure 7.
Similarly to our analysis of the spread of the allele in rs4988235(T), we inferred the dynamics of the rs1042602(A) allele separately for the time periods before and after 5000 years BP and assuming the age of the allele to be 26,361 years (Albers and McVean, 2020). The inferred parameters for both time periods are summarized in Appendix 2—table 7, and the allele frequency surface maps generated using these parameters are shown in Figure 8b. The origin of the rs1042602(A) corresponds closely to the region where the allele initially starts to segregate in the time period between 7500 and 10,000 years BP as seen in Figure 8a. Estimates of the selection coefficient for both time periods (0.0221 and 0.0102 for the period before and after 5000 years BP, respectively) suggest that selection acting on the allele has decreased after 5000 years BP.
Robustness of parameters to the inferred geographic origin of allele
We carried out an analysis to characterize how sensitive the selection, diffusion, and advection parameters are to changes in the assumed geographic origin of the allele. For the rs4988235(T) allele, we forced the origin of the allele to be 10° away from our inferred origin in each cardinal direction, while assuming the allele age is equal to the inferred start of selection onset from Itan et al., 2009 (Appendix 2—table 8). In Figure 6—figure supplements 2–5, we can see the allele frequency dynamics of these four scenarios, respectively. We also forced the allele origin to be at the geographic origin estimated in Itan et al., 2009 (Figure 6—figure supplement 6, Appendix 2—table 9), which is westward of our estimate. In all five cases during the period prior 5000 years BP, the allele is inferred to expand in the direction of the first sample that is observed to carry the rs4988235(T) allele and is located in Ukraine. During the time period after 5000 years BP, the patterns produced by the model are rather similar, although the parameters associated with diffusion and advection differ, in order to account for the different starting conditions.
We also investigated how the results are affected when the estimated geographic origin of the rs1042602(A) allele is moved with respect to the initial estimate. We set the allele to be 10° east, 10° north, and 10° south of the original estimate as shown in Figure 8—figure supplements 1–3, respectively (for parameter estimates see Appendix 2—table 10). We did not look at a scenario in which the origin of the allele is moved to the west since it would either end up in the Black Sea or more westward than 10°. The selection coefficient remains similar to the original estimate throughout all three scenarios. The way the allele spreads across the landscape is also similar in all cases and, as in the case of rs4988235(T), the model accounts for the different origins of the allele by adjusting the diffusion and advection coefficients in the time period after 5000 years BP.
Robustness of parameters to the assumed age of the allele
In order to investigate how sensitive our inferences are to the point estimates of allele ages we obtained from the literature (Albers and McVean, 2020; Itan et al., 2009), we also fitted our model using the upper and lower ends of the 95% confidence intervals or credible intervals for each age estimate (depending on whether the inference procedure in the literature was via a maximum likelihood or a Bayesian approach). For the rs4988235(T) allele, the reported credible intervals for the (Itan et al., 2009) age are 8683 and 6256 years BP. For the rs1042602(A) allele, the reported confidence intervals for the age are 27,315 and 25,424 years BP (Albers and McVean, 2020).
When refitting the model for the rs4988235(T) allele, we found that the inferred selection coefficient is slightly lower when the allele age is assumed to be at the lower bound of the 95% credible interval (0.0867 vs. 0.0993 before 5000 years BP and 0.0321 vs. 0.0328 after 5000 years BP) and slightly higher when assumed to be at the upper bound (0.0994 vs. 0.0993 before 5000 years BP and 0.0572 vs. 0.0328 after 5000 years BP) (Appendix 2—table 5 and Figure 6—figure supplements 7 and 8). This occurs because the selection intensity must be higher or lower when there is more or less time, respectively, for the allele to reach the allele frequencies observed in the data. In the case of the rs1042602(A) allele, this only affects the earlier time period (Appendix 2—table 7). The rs4988235(T) allele’s geographic distribution in the more recent time periods is also less extended geographically when the age is assumed to be young. The inferred geographic origin of both alleles slightly differs under different assumed ages (Figure 8—figure supplements 4 and 5).
In addition, we assessed the likelihood of the bestfitted models with varying the ages of the rs4988235(T) and rs1042602(A) alleles (Figure 6—figure supplement 9 and Figure 8—figure supplement 6, respectively). We can see that in the case of rs4988235(T) allele the allele age used in this study (7441 years) gives the most likely solution among the explored ages. In case of the rs1042602(A) allele, we found that there are multiple nearly equally likely ages when looking at ages at least as old as 15,000 years.
Discussion
A spatially explicit framework for allele frequency diffusion can provide new insights into the dynamics of selected variants across a landscape. We have shown that under the conditions of strong, recent selection, our method can infer selection and dispersal parameters using a combination of ancient and presentday human genomic data. However, when allowing for advection, the inferred location tends to become less accurate. This suggests that migration events early in the dispersal of the selected allele could create difficulties in finding the true allele origin if net directional movement (i.e., via major migratory processes) had a large effect in this dispersal. This issue could be alleviated with the inclusion of more ancient genomes around the time of the origin of the mutation, perhaps in combination with a more finescaled division into periods where advection may have occurred in different directions.
The inferred geographic origin of the rs4988235(T) allele reflects the best guess of our framework given the constraints provided by its input, namely, the previously inferred age of the allele and the observed instances of this allele throughout Western Eurasia. We are also assuming that the allele must have arisen somewhere within the bounding box of our studied map. When assuming a relatively young allele age (7441 years ago, equal to the start of selection onset in Itan et al., 2009), the origin of the allele is placed north of the Caucasus, perhaps among steppe populations that inhabited the area at this time (Haak et al., 2015; Allentoft et al., 2015). This origin is further east than the geographic origin estimate from Itan et al., 2009, likely reflecting additional ancient DNA information that is available to us, and indicates an early presence of the allele in Eastern Europe. When assuming a relatively old allele age (20,106 years ago, Albers and McVean, 2020), the age is placed in Northeast Europe, perhaps among Eastern huntergatherer groups that inhabited the region in the early Holocene. We note that the number of available genomes for Eastern and Northeastern Europe during the early Holocene is scarce, so the uncertainty of the exact location of this origin is relatively high. Regardless of the assumed age, we estimate a net westward displacement of the allele frequency’s center of mass, and a rapid diffusion, particularly in the period after 5000 years ago.
Various studies have estimated the selection coefficient for the rs4988235(T) allele, and these range from as low as 0.014 to as high as 0.19 (Enattah et al., 2008; Mathieson and Mathieson, 2018; Mathieson, 2020; Stern et al., 2019; Burger et al., 2020; Peter et al., 2012; Gerbault et al., 2009; Itan et al., 2009; Bersaglieri et al., 2004). Recent papers incorporating ancient DNA estimate the selection coefficient to be as low as 0 (in certain regions of Southern Europe) and as high as 0.06 (Mathieson and Mathieson, 2018; Mathieson, 2020; Burger et al., 2020). It is also likely that the selection coefficient was different for different regions of Europe, perhaps due to varying cultural practices (Mathieson, 2020). In our case, the estimated selection coefficient during the first period—before 5000 years ago—depends strongly on the assumed allele age (s = 0.0993 vs. s = 0.0285). As in the case of the geographic origin, these estimates should be taken with caution as the number of available allele observations in the early Holocene is fairly low. The estimates for the second period—after 5000 years ago—are more robust to the assumed age: s = 0.0328 (95% CI: 0.0327–0.0329) if we assume the younger allele age (7441 years ago) and s = 0.0255 (95% CI: 0.0252–0.0258) if we assume the older allele age (20,106 years ago). These estimates are also within the range of previous estimates.
In the case of the rs1042602(A) allele, our estimated selection coefficients of 0.0221 (95% CI: 0.0216–0.0227) and 0.0102 (95% CI: 0.0083–0.0120) for the time periods before and after 5000 years BP, respectively, are generally in agreement with previous results. Wilde et al., 2014 used a forward simulation approach to infer a point estimate of 0.026. Another study using an approximate Bayesian computation framework (Nakagome et al., 2019) estimated the strength of selection acting on rs1042602 to be 0.013 (0.002–0.029). Although both studies utilized ancient DNA data, the estimates were obtained without explicitly modeling the spatial dimension of the selection process.
Our estimates of the longitudinal advection parameter are negative for both the SNPs in the TYR and LCT loci: the mutation origins are always to the east of the center of mass of the allele frequency distribution seen in presentday data. This perhaps reflects common migratory processes, like the largescale Neolithic and Bronze Age population movements from east to west, affecting the allele frequencies at these loci across the Eurasian landscape (Allentoft et al., 2015; Haak et al., 2015). As a form of regularization, we kept the range of explored values for the advection parameters to be small (−2.5–2.5 km per generation), while allowing the diffusion parameters to be explored over a much wider range of values. In certain cases, like the second period of the rs4988235(T) spread when the allele age is assumed to be young (Appendix 2—table 5), we find that the advection parameters are fitted at the boundary of the explored range, because the allele needs to spread very fast across the landscape to fit the data.
A future improvement to our method could include other forms of regularization that better account for the joint behavior of the advection and diffusion processes, or the use of priors for these parameters under a Bayesian setting, which could be informed by realistic assumptions about the movement of individuals on a landscape. Bayesian parameter fitting would likely provide a more robust understanding of the uncertainty of the estimates as well as an opportunity to formally compare different models using Bayes factors, although at the cost of an increase of computational intensity.
When investigating the robustness of the geographic origin of both rs4988235(T) and rs1042602(A), we found that parameters related to the beneficial allele’s expansion change in response to different assumed origins of the allele. The resulting allele frequency surface plots, however, appear very similar throughout the later stages of the process, showing that the model tends to adjust the diffusion and advection coefficients in a way such that the allele will end up expanding into the same areas regardless of the origin.
As we apply these methods to longer time scales and broader geographic areas, the assumptions of spatiotemporal homogeneity of the parameters seem less plausible. There may be cases where the allele may have been distributed over a wide geographic area but remained at low frequencies for an extended period of time, complicating the attempts to pinpoint the allele’s origin. In our study, we estimated diffusion and selection coefficients separately for two time periods before and after 5000 years ago to account for changes in mobility during the Bronze Age, but this approach may still be hindered by uneven sampling, especially when the allele in question exists at very low frequencies. Notably, our results for the spread of the rs4988235(T) allele during the older time period should be interpreted with caution, since they may be affected by sparse sampling in the early Holocene.
Potential future extensions of our method could incorporate geographic features and historical migration events that create spatially or temporally varying moderators of gene flow. An example of this type of processes is the retreat of glaciers after the last Glacial maximum, which allowed migration of humans into Scandinavia (Günther et al., 2018). These changing geographic features could lead to changes in the rate of advection or diffusion across time or space. They could also serve to put more environmentally aware constraints on the geographic origin of the allele, given that it cannot have existed in regions uninhabitable by humans, and to extend our analyses beyond the narrow confines of the Western Eurasian map chosen for this study. One could also envision incorporating variation in population densities over time or known migration processes in the time frames and regions of interest. These might have facilitated rapid, longrange dispersal of beneficial alleles (Bradburd et al., 2016; Hallatschek and Fisher, 2014) or caused allelic surfing on the wave of range expansions (Klopfstein et al., 2006). Additional information like this could come, for example, from previously inferred spatiotemporal demographic processes (e.g. Racimo et al., 2020b).
As described above, our model only accounts for diffusion in two directions. Further extension of our model could therefore incorporate anisotropic diffusion (Othmer et al., 1988; Painter and Hillen, 2018). Another possibility could be the introduction of stochastic process components in order to convert the partial differential equations into stochastic differential equations (Brown et al., 2000). Stochastic components could serve to induce spatial autocorrelation and capture local patterns of allele frequency covariance in space that might not be well modeled by the deterministic partial differential equations (PDEs) (Cressie and Wikle, 2015). They could also serve to induce stochasticity in allele frequency changes over time as a consequence of genetic drift (Crow and Kimura, 1970), allowing one to model the dynamics of more weakly selected variants, where drift plays an important role. Eventually, one could perhaps combine information across loci to jointly model the spatiotemporal frequency surfaces at multiple loci associated with the same trait. This could help clarify the dynamics of polygenic adaptation and negative selection on complex traits (IrvingPease et al., 2021), and perhaps hindcast the genetic value of traits across a landscape.
The availability of hundreds of ancient genomes (Marciniak and Perry, 2017) and the increasing interest in spatiotemporal method development (Bradburd and Ralph, 2019), such as the one described in this article, will likely lead researchers to posit new questions and hypotheses about the behavior of natural selection. In the case of a beneficial allele spreading on a landscape, new ontologies and vocabulary for describing positive selection in time and space will be needed. Abundant terms exist to classify the initial conditions and dynamics of a selective sweep in a single population (hard sweep, multipleorigin soft sweep, singleorigin soft sweep, partial sweep) (Hermisson and Pennings, 2005; Pritchard and Di Rienzo, 2010; Hermisson et al., 2017). In contrast, there is a lack of vocabulary for distinguishing between a scenario of strong selection that is locally constrained in space from a scenario of widespread selection extended over a landscape, or a model of neutral diffusion in space followed by parallel nonneutral increases in frequency at multiple locations. For example, Ralph and Coop, 2010 showed how multiple localized hard sweeps may be seen as a soft sweep at a larger populationwide scale. Existing vocabulary for spatiotemporal genetic processes is clearly not enough, limiting the types of questions or hypotheses we can pose about them.
Population genetic models that explicitly account for space and time are an important area of future methodological development (Bradburd and Ralph, 2019). We believe that methods such as the one described in this study show great promise at broadening the horizon of our understanding of natural selection across space and time in humans and other species. As in the case of demographic reconstruction (Ray and Excoffier, 2009), spatiotemporal information can greatly help improve our knowledge of how natural selection operated in the past.
Methods
The model
To describe the allele frequency dynamics in time and space, we first begin by using a deterministic model based on a twodimensional PDE (Fisher, 1937; Kolmogorov et al., 1937; Novembre et al., 2005). This PDE represents the distribution $p(x,y,t)$ of the allele frequency across a twodimensional ($x,y$) landscape at time $t$:
where
Here, $\sigma $ is the diffusion coefficient, $s$ is the selection coefficient, and $d$ is the dominance coefficient (Novembre et al., 2005). We assumed an additive model and fixed $d=2s$ in all analyses below. We call this ‘model A,’ but we also evaluated the fit of our data under more complex models that are more flexible, and are described below.
Model B is a more general diffusionreaction model, which incorporates distinct diffusion terms in the longitudinal and latitudinal directions (${\sigma}_{x}$ and ${\sigma}_{y}$, respectively):
Model C is a generalization of model B that incorporates advection terms in the longitudinal and latitudinal directions (see, e.g. Cantrell and Cosner, 2004 for a motivation of this type of model in the context of spatial ecology):
Here, v_{x} and v_{y} represent the coefficients for advective velocity along the longitude and latitude respectively.
In Appendix 1, we motivate the construction of these equations using model C as an example and show that Equation 4 can be obtained by taking an infinitesimal limit of a random walk on a twodimensional lattice, after including a reaction term due to selection. Models A and B are then shown to be special cases of model C.
For evaluating the likelihood of the observed data, we use a binomial genotype sampling model. Let ${g}_{i}\in 0,1,2$ be the genotype of individual $i$ at the locus of interest, let a_{i} be the number of reads carrying ancestral alleles, and let d_{i} be the number of reads carry derived reads. Let $({x}_{i},{y}_{i})$ be the coordinates of the location from which individual $i$ was sampled, and t_{i} its estimated age (e.g., from radiocarbon dating). Then, the likelihood for individual $i$ can be computed as follows:
Here, $p({x}_{i},{y}_{i},{t}_{i})$ is the solution to one of the partial differential equations described above (Equation 1, Equation 2, or Equation 4, depending on the process model chosen), evaluated at location $({x}_{i},{y}_{i})$ and time t_{i}. In turn, $P[{d}_{i},{a}_{i}{g}_{i}=h]$ is the likelihood for genotype $i$. Furthermore, $P[{g}_{i}=hp({x}_{i},{y}_{i},{t}_{i})]$ is a binomial distribution, where $n$ represents the ploidy level, which in this case is 2:
Then, the likelihood of the entire data can be computed as
where M is the total number of individuals for which we have data, d is the vector containing the derived read count for each individual, and a is the vector containing the ancestral read count for each individual. We computed genotype likelihoods directly on the BAM file read data using the SAMtools genotype model (Li, 2011) implemented in the software ANGSD (Korneliussen et al., 2014).
When only randomly sampled pseudohaploid allele counts are available, we used a Bernoulli sampling likelihood (conditional on the genotype g_{i}) on the righthand side of Equation 6 instead. Briefly, assuming that the probability of an individual having genotype $g$ at a particular locus given the underlying allele frequency $p$ follows a binomial distribution and that the probability of sampling a read given the genotype of an individual follows a Bernoulli distribution with probability of success $\frac{1}{2}g$, then the probability of sampling a read given the genotype follows a Bernoulli distribution with probability of success $p$.
Map
We restricted the geographic area explored by our model fit to be between 30°N to 75°N, and between 10°W and 80°E. For numerical calculations, we used a grid constructed using a resolution of approximately one grid cell per latitude and longitude. We used Harvesine functions in order to transform the distance from degrees to kilometers between two geographic points. The diffusion of the allele frequency was disallowed in the map regions where the topology is negative (i.e., regions under water), based on ETOPO5 data (NOAA, National Geophysical Data Center, B. C, 1988). For this reason, we added land bridges between the European mainland and Sardinia, and between the mainland and Great Britain, in order to allow the allele to diffuse in these regions (see Appendix 2—figure 1).
Parameter search
Parameter optimization was done via maximum likelihood estimation with a twolayer optimization setup. The first layer consists of a simulated annealing approach (Bélisle, 1992) starting from 50 random points in the parameter space. The initial 50 points are sampled using Latin hypercube sampling to ensure an even spread across the parameter space. The output of this fit was then fed to the LBFGSB algorithm to refine the parameter estimates around the obtained maximum and obtain confidence intervals for the selection, diffusion and advection parameters (Byrd et al., 1995).
The parameters optimized were:
The selection coefficient ($s$), restricted to the range 0.001–0.1.
Two dispersal parameters ${\sigma}_{x}$ and ${\sigma}_{y}$ in the longitudinal and latitudinal directions, respectively, restricted to the range of 1–100 square kilometers per generation.
The longitudinal and latitudinal advection coefficients v_{x} and v_{y}, respectively. As a form of regularization, we set the range of explored values to be narrowly centered around zero: –2.5 to 2.5 kilometers per generation.
The geographic origin of the allele, which is randomly initialized to be any of the 28 spatial points shown in Appendix 2—figure 2 at the start of the optimization process.
We chose to construct our method in a way that uses the age of the allele as an input parameter rather than estimating it. We do this since there are multiple equally possible solutions with various combinations of allele age and selection coefficient values as shown in Appendix 2—figure 3. The latitude and longitude are discretized in our model in order to solve the differential equations numerically, thus the origin of a mutation is measured in terms of discrete units. For this reason, when using the LBFGSB algorithm, we fixed the previously estimated origin of the allele and did not explore it during this second optimization layer. For numerical calculations, we used the Livermore Solver for Ordinary Differential Equations (Hindmarsh, 1983) implemented in R package ‘deSolve’ (Soetaert et al., 2010), which is a generalpurpose solver that can handle both stiff and nonstiff systems. In case of stiff problems, the solver uses a Jacobian matrix. Absorbing boundary conditions were used at the boundaries of the map. For visualization purposes, we masked the allele frequencies from areas with negative topology (i.e., areas covered by large bodies of water). Time was measured in generations, assuming 29 years per generation. During the optimization, we scaled the time and the parameters by a factor of 10, which allowed us to decrease the execution time of the model.
We initialized the grid by setting the initial allele frequency to be p_{0} in a grid cell where the allele originates and 0 elsewhere. p_{0} was calculated as $1/(2*D*A)$, where $D$ is the population density and is equal to 2.5 inhabitants per square kilometer, which is the estimated population density in Europe in 1000 BC (Colin McEvedy, 1978; Novembre et al., 2005). In the equation, $D$ is multiplied by 2 because we assume that the allele originated in a single chromosome in a diploid individual. $A$ is the area in square kilometers of the grid cell where the allele emerged.
Asymptotic 95% confidence intervals for a given parameter ${\theta}_{j}$ were calculated using equation
where $F(\mathit{\theta})$ is an estimate of the observed Fisher information matrix (Fisher, 1922; Efron and Hastie, 2016; Casella and Berger, 2001).
Implementation
The abovedescribed model was implemented in R version 3.6. To numerically solve the differential equations and obtain maximum likelihood estimates, we used the libraries deSolve (Soetaert et al., 2010), ReacTran (Soetaert and Meysman, 2012), and bbmle (Bolker, 2020). Scripts containing the code used in this article are available on GitHub: https://github.com/RasaMukti/stepadna, (copy archived at swh:1:rev:d024767648d873f329a8e17fcaf6034c99157120; Muktupavela, 2021).
Individualbased simulations
For the individualbased spatiotemporal forward simulations, we first defined a spatial boundary for a population spread across a broad geographic region of Europe. In order to ensure a reasonably uniform distribution of individuals across this spatial range throughout the course of the simulation, we set the maximum distance for spatial competition and mating choice between individuals to 250 km (translated, on a SLiM level, to the interaction parameter maxDistance), and the standard deviation of the normal distribution governing the spread of offspring from their parents at 25 km (leveraged in SLiM’s modifyChild() callback function) (Haller and Messer, 2019). We note that we have chosen the values of these parameters merely to ensure a uniform spread of individuals across a simulated landscape. They are not intended to represent realistic estimates for these parameters at any time in human history.
After defining the spatial context of the simulations and ensuring the uniform spread of individuals across their population boundary, we introduced a single beneficial additive mutation in a single individual. In order to test how accurately our model can infer the parameters of interest, we simulated a scenario in which the allele appeared in Central Europe 15,000 years ago with the selection coefficient of the beneficial mutation set to 0.03. Over the course of the simulation, we tracked the position of each individual that ever lived together with its location on a twodimensional map, as well as its genotype (i.e., zero, one, or two copies of the beneficial allele). We then used this complete information about the spatial distribution of the beneficial allele in each time point to study the accuracy of our model in inferring the parameters of interest.
Appendix 1
Here, we motivate the construction of model C as a largescale limit of a random walk model on a lattice (Karlin and Taylor, 1975; Cantrell and Cosner, 2004). We think of the allele frequency as a variable $p$ that can increase in magnitude due to its inherent advantage (selection), spread across a landscape (diffusion) or move directionally as a consequence of migration (advection). We imagine a lattice composed of small square cells of size $\mathrm{\Delta}x$ x $\mathrm{\Delta}y$, where a certain amount of allele frequency $p$ can occur at a given time point $t$. At each small time step (of duration $\mathrm{\Delta}t$), inflow and outflow of p can occur in the xdirection with probability h or in the ydirection with probability 1h, and the magnitude of these flows depend on the amount of $p$ present in neighboring cells. If flow of p is along the xaxis, it does so in the positive direction with probability $\alpha $ and in the negative direction with probability $1\alpha $. If flow of p is along the yaxis, it does so in the positive direction with probability $\beta $ and in the negative direction with probability $1\beta $. The allele frequency can also increase in magnitude locally via a function $\gamma ()$ that depends on its dominance (d), selection coefficient (s) and current magnitude ($p(x,y,t)$). Then, we obtain:
We can also write this as:
If we divide both sides by $\mathrm{\Delta}t$ and take the limit of infinitesimally small $\mathrm{\Delta}x$, $\mathrm{\Delta}y$, and $\mathrm{\Delta}t$, while assuming that, in this limit, $\frac{\mathrm{\Delta}{x}^{2}}{\mathrm{\Delta}t}$ and $\frac{\mathrm{\Delta}{y}^{2}}{\mathrm{\Delta}t}$ are finite (Okubo, 1980), we obtain:
where ${u}_{x}=\frac{\mathrm{\Delta}x}{\mathrm{\Delta}t}$, ${u}_{y}=\frac{\mathrm{\Delta}y}{\mathrm{\Delta}t}$, ${\lambda}_{x}=\frac{\mathrm{\Delta}{x}^{2}}{\mathrm{\Delta}t}$, ${\lambda}_{y}=\frac{\mathrm{\Delta}{y}^{2}}{\mathrm{\Delta}t}$.
If we let ${\sigma}_{x}^{2}=h{\lambda}_{x}$, ${\sigma}_{y}^{2}=(1h){\lambda}_{y}$, ${v}_{x}=h(12\alpha ){u}_{x}$, ${v}_{y}=(1h)(12\beta ){u}_{y}$, then we obtain Equation 4. Thus, we can see that the squared diffusion coefficient ${\sigma}_{x}^{2}$ depends on the square of the length of the cells in the xaxis relative to the duration of a time step (${\lambda}_{x}$), and on the probability that flows occurs in the xaxis at a given time step ($h$). Similarly, the squared diffusion coefficient ${\sigma}_{y}^{2}$ depends on the square of the length of the cells in the yaxis relative to the duration of a time step (${\lambda}_{y}$), and on the probability that flows occurs in the yaxis at a given time step ($1h$). The advection coefficient v_{x} depends on the advective velocity along the xaxis (u_{x}) as well as on the probability of flow occurring along the xaxis ($h$) and the directional bias $12\alpha $, which depends on the probability that flow occurs in the positive xdirection ($\alpha $). Finally, the advection coefficient v_{y} depends on the advective velocity along the yaxis (u_{y}), as well as on the probability of flow occurring along the yaxis ($1h$) and the directional bias $12\beta $, which depends on the probability that flow occurs in the positive ydirection ($\beta $).
We can recover model B as a special case of model C if we fix $\alpha =\beta =\frac{1}{2}$, assuming isotropy in the two directions, so that $\mathrm{\Delta}x=\mathrm{\Delta}y$. We can also recover model A if we additionally fix $h=\frac{1}{2}$.
Appendix 2
Data availability
The current manuscript is a computational study, therefore no data have been generated for this manuscript. Software code along with publicly available data used for this study are deposited to GitHub: https://github.com/RasaMukti/stepadna/tree/main/reproducibles (copy archived at swh:1:rev:d024767648d873f329a8e17fcaf6034c99157120).

Allen Ancient DNAID AADR. Allen Ancient DNA Resource (AADR).

EMBLEBIID HGDP. Human Genome Diversity Panel.
References

Complex signatures of selection for the melanogenic loci Tyr, TYRP1 and DCT in humansBMC Evolutionary Biology 8:1–14.https://doi.org/10.1186/14712148874

Convergence theorems for a class of simulated annealing algorithms on rdJournal of Applied Probability 29:885–895.https://doi.org/10.2307/3214721

Genetic signatures of strong recent positive selection at the lactase geneAmerican Journal of Human Genetics 74:1111–1120.https://doi.org/10.1086/421051

Spatial population genetics: it’s about timeAnnual Review of Ecology, Evolution, and Systematics 50:427–449.https://doi.org/10.1146/annurevecolsys110316022659

Blurgenerated nonseparable spacetime modelsJournal of the Royal Statistical Society 62:847–860.https://doi.org/10.1111/14679868.00269

A limited memory algorithm for bound constrained optimizationSIAM Journal on Scientific Computing 16:1190–1208.https://doi.org/10.1137/0916069

BookSpatial Ecology via Reaction‐Diffusion EquationsJohn Wiley & Sons.https://doi.org/10.1002/0470871296

BookAtlas of World Population HistoryGreat Britain: Penguin Books Lyd. and Allen Lane.

BookStatistics for SpatioTemporal DataHoboken, New Jersey, USA: John Wiley & Sons.https://doi.org/10.1111/j.15384632.2012.00859.x

Inference of natural selection from ancient DNAEvolution Letters 4:94–108.https://doi.org/10.1002/evl3.165

BookComputer Age Statistical InferenceCambridge University Press.https://doi.org/10.1017/CBO9781316576533

Identification of a variant associated with adulttype hypolactasiaNature Genetics 30:233–237.https://doi.org/10.1038/ng826

BookMathematical Population Genetics 1: Theoretical IntroductionSpringer Science & Business Media.https://doi.org/10.1007/9780387218229

On the mathematical foundations of theoretical statisticsPhilosophical Transactions of the Royal Society of London. Series A, Containing Papers of A Mathematical or Physical Character 222:309–368.https://doi.org/10.1098/rsta.1922.0009

The wave of advance of advantageous genesAnnals of Eugenics 7:355–369.https://doi.org/10.1111/j.14691809.1937.tb02153.x

Herders of Indian and European cattle share their predominant allele for lactase persistenceMolecular Biology and Evolution 29:249–260.https://doi.org/10.1093/molbev/msr190

Slim 3: forward genetic simulations beyond the WrightFisher modelMolecular Biology and Evolution 36:632–637.https://doi.org/10.1093/molbev/msy228

Soft sweeps and beyond: understanding the patterns and probabilities of selection footprints under rapid adaptationMethods in Ecology and Evolution 8:700–716.https://doi.org/10.1111/2041210X.12808

ODEPACK, A systematized collection of ODE solversScientific Computing 1:55–64.https://doi.org/10.12691/ajmo411

The origins of lactase persistence in europePLOS Computational Biology 5:e1000491.https://doi.org/10.1371/journal.pcbi.1000491

A worldwide correlation of lactase persistence phenotype and genotypesBMC Evolutionary Biology 10:36.https://doi.org/10.1186/147121481036

BookA First Course in Stochastic Processes (Second)New York: Academic Press.https://doi.org/10.1016/C20091285698

The fate of mutations surfing on the wave of a range expansionMolecular Biology and Evolution 23:482–490.https://doi.org/10.1093/molbev/msj057

ANGSD: analysis of next generation sequencing dataBMC Bioinformatics 15:356.https://doi.org/10.1186/s1285901403564

Harnessing ancient genomes to study the history of human adaptationNature Reviews. Genetics 18:659–674.https://doi.org/10.1038/nrg.2017.65

Demography and the age of rare variantsPLOS Genetics 10:e1004528.https://doi.org/10.1371/journal.pgen.1004528

Fads1 and the timing of human adaptation to agricultureMolecular Biology and Evolution 35:2957–2970.https://doi.org/10.1093/molbev/msy180

Inferring the model and onset of natural selection under varying population size from the site frequency spectrum and haplotype structureProceedings. Biological Sciences 286:20182541.https://doi.org/10.1098/rspb.2018.2541

BookData Announcement 88MGG02,, Digital Relief of the Surface of the EarthPublisher: U.S. Department of Commerce.

Human population structure and the adaptive response to pathogeninduced selection pressuresPhilosophical Transactions of the Royal Society of London. Series B, Biological Sciences 367:878–886.https://doi.org/10.1098/rstb.2011.0305

BookDiffusion and ecological problems: mathematical modelsIn: Andersson S, editors. Biomathematics. Berlin, Heidelberg & New York: SpringerVerlag. pp. 254–272.https://doi.org/10.1016/B9780444502735.X50253

Models of dispersal in biological systemsJournal of Mathematical Biology 26:263–298.https://doi.org/10.1007/BF00277392

Beyond broad strokes: sociocultural insights from the study of ancient genomesNature Reviews. Genetics 21:355–366.https://doi.org/10.1038/s415760200218z

On the evolution of lactase persistence in humansAnnual Review of Genomics and Human Genetics 18:297–319.https://doi.org/10.1146/annurevgenom091416035340

Estimating allele ageAnnual Review of Genomics and Human Genetics 1:225–249.https://doi.org/10.1146/annurev.genom.1.1.225

Solving differential equations in R: package desolveJournal of Statistical Software 33:1–25.https://doi.org/10.18637/jss.v033.i09

Reactive transport in aquatic ecosystems: rapid model prototyping in the open source software REnvironmental Modelling & Software 32:49–60.https://doi.org/10.1016/j.envsoft.2011.08.011

Dating the origin of the CCR5delta32 AIDSresistance allele by the coalescence of haplotypesAmerican Journal of Human Genetics 62:1507–1515.https://doi.org/10.1086/301867

BookSpatioTemporal Statistics with RBoca Raton, Florida : CRC Press: CRC Press.https://doi.org/10.1201/9781351769723
Decision letter

Aida M AndrésReviewing Editor; University College London, United Kingdom

George H PerrySenior Editor; Pennsylvania State University, United States

Isabel AlvesReviewer; l'institut du thorax, INSERM/CNRS, University of Nantes, France

Anders ErikssonReviewer; Institute of Genomics, University of Tartu, Estonia
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Modelling the spatiotemporal spread of beneficial alleles using ancient genomes" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and George Perry as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Isabel Alves (Reviewer #1); Anders Eriksson (Reviewer #2).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
The reviewers and editor found the manuscript elegant and interesting and appreciate the novelty and promise of your new method. They also value your new inferences on two important targets of positive selection in humans. You will find the individual reviews below, to be addressed individually. Here we provide a summary of the desired revisions, to assist you in the revision process.
1) We believe that the manuscript needs a more systematic test for the sensitivity of the method to 1) misspecification of the age of the aDNA samples, and 2) the geographic location/clustering of the samples. This could be done for example by exploring the effect of these factors on the accuracy of the parameter estimates.
2) Related, we believe that further discussion is needed on the extent to which differences observed across time periods (e.g. before and after 5000 BP) reflect particularities of the alleles (e.g. selection), populations (e.g mobility) or the density of the samples across time and space.
3) Where models are applied to specific loci, results are often reported in the text without enough quantitative information and we think that clearer and more detailed information about the inferences must be presented in the text. As an example, where the model is run for different allele age estimates the (approximate maximum) likelihoods of parameter estimates for the different scenarios should be clearly shown in the text.
4) We believe that the authors should clarify some of their choices for the method. We would like to see a discussion of whether the method may be suitable for treating allele age as a free parameter, or whether it is possible (and superior) to report Bayesian posterior distributions for the parameters and Bayes factors to compare the fit of different models to the data.
5) The loci studied represent unusual cases of very strong positive selection in particular geographic areas. We would welcome more discussion regarding the extent to which the current version of the method can be applied to other alleles, to other geographic regions and even to other species.
6) Finally, the reviewers raised issues regarding the presentation of the numerical implementation that we hope that the authors can easily address.
Reviewer #1 (Recommendations for the authors):
I very much appreciated reading this manuscript and the effort of the authors in developing such an exciting tool. However, I still have some concerns about the accuracy and the meaning of the estimates.
1) We know that ancient DNA samples are clustered in time and space because they are an aggregate of samples coming from different studies often with different scopes in terms of time and geographical area covered. What's the impact of having clustered samples in the parameter estimates? The authors comment on that a little bit when comparing the accuracy tests performed with the deterministic vs forward in time simulations. However, I think a more systematic analysis of different sampling schemes would help clarifying the robustness of the parameter inference. Moreover, I think this would greatly help understanding the differences in behavior of the statistical method before and after 5000 BP. Indeed the authors allow the model parameters to differ across different time periods (before and after 5000 BP) but it is not clear whether differences in the estimates reflect any particularity of the populations (e.g mobility) or they just reflect differences in the density of the samples across time and space.
2) How would the method perform if simulations would be performed without advection (numerical solution to equation 3) but parameter estimates would be inferred from a model C, which includes advection ?
3) Finally and more of a general question: what's the point of estimating the geographical origin of an allele if there are no samples dated from around the same time as the assumed allele age? I have the feeling this can be easily misinterpreted but maybe I am missing something. The authors discuss this (lines 199210) when discussing the results obtained for the dynamics of the rs4988235(T) for the two assumed allele ages but it is not clear to me if the estimate for the local of origin is actually meaningful under similar situations. I wonder if with simulation one could try to understand what the method is retrieving…
With respect to the presentation of the work I found that the manuscript is globally well explained but there are some parts a bit difficult to follow:
1) it should be clearly stated in the main that the age of the allele is provided as input and not inferred contrarily to the geographical origin of the allele and the value used to perform the inference should also be provided in the corresponding part of the main text.
2) The comparison of the results for the two ages assumed for the rs4988235(T) mutation (Itan et al. 2009 and Albers and McVean 2020), mainly between lines 199210, is very confusing. From Figure 5a it seems there is no information on the allele for the time period between the allele age inferred by Itan et al. 2009 and Albers and McVean 2020. However, from the text it is not clearly stated whether there is or not information. If there is no information at all how does the algorithm infer an origin in Northeastern Europe?
3) It would also help to follow the results if the authors would clearly state in the main text (lines: 212230) the age of the allele assumed for the rs1042602(A).
Reviewer #2 (Recommendations for the authors):
Given the ability to of the model to sample at any point in space and time, it would have been very interesting to test the model under actual geographical and temporal distributions of data, such as the Allen Ancient DNA Resource used for one of the empirical investigations, perhaps even taking the observed pattern of missingness in the data into account.
I would recommend using latitude and longitude systematically throughout the manuscript instead of mixing x, y, latitude, longitude.
In the methods: Please fix the typo in Eq (9). Also please clarify which boundary conditions where used, and how the equations were solved numerically (using Eq. (8) or (9)?), including how boundary conditions were implemented.
L437: " the lefthand side of equation 5". This doesn't look right, do you mean the righthand side of Eq. (6)?
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Reviewer #1 (Recommendations for the authors):
It is now the second time I review the manuscript "Modelling the spatiotemporal spread of beneficial alleles using ancient genomes" by Muktupavela et al. and I am glad to see that the authors addressed most of the reviewers' concerns. There is now a more systemic assessment of the advantages/limitations of the method which I truly believe to be crucial for future implementations. That being said, I would like to recommend this manuscript for publication.
Reviewer #2 (Recommendations for the authors):
The authors have made substantial improvements to the manuscript, and is ready for acceptance except for one remaining point:
Contrary to what the authors claim in their rebuttal, Eqs (8) and (9) in Appendix 1 still contain errors.
First, p(x+\δ x,y,t) and p(x+\δ x,y,t) appears twice on line 827.
Second, the coefficients of p(x,y,t) on both sides of Eq (9) sum to 1, which is correct. The term should therefore be absent (have coefficient zero) in Eq (8).
https://doi.org/10.7554/eLife.73767.sa1Author response
Reviewer #1 (Recommendations for the authors):
I very much appreciated reading this manuscript and the effort of the authors in developing such an exciting tool. However, I still have some concerns about the accuracy and the meaning of the estimates.
1) We know that ancient DNA samples are clustered in time and space because they are an aggregate of samples coming from different studies often with different scopes in terms of time and geographical area covered. What's the impact of having clustered samples in the parameter estimates? The authors comment on that a little bit when comparing the accuracy tests performed with the deterministic vs forward in time simulations. However, I think a more systematic analysis of different sampling schemes would help clarifying the robustness of the parameter inference. Moreover, I think this would greatly help understanding the differences in behavior of the statistical method before and after 5000 BP. Indeed the authors allow the model parameters to differ across different time periods (before and after 5000 BP) but it is not clear whether differences in the estimates reflect any particularity of the populations (e.g mobility) or they just reflect differences in the density of the samples across time and space.
We thank the reviewer for highlighting the importance of taking into account clustering effects that could potentially arise by aggregating aDNA data from studies with different sampling schemes. To elaborate on this issue in more detail, we have included a new analysis (chapter “Impact of sample clustering on parameter estimates” in the manuscript) of the expected impact of different sampling and clustering schemes on our inferences. We used a deterministic simulation to create three different degrees of clustering which we will refer to as “homogeneous”, “intermediate” or “extreme” by varying the area from which we sample individuals to be used in our inferences, now included as Figure 3—figure supplement 1 in the paper. Additionally, we also tested the impact of biased temporal sampling in the periods before and after 5000 year BP by oversampling in the ancient period (75%/25%), equal sampling in the two periods (50%/50%), and oversampling in the recent period (25%/75%). Because we evaluated this temporal bias for each of the three spatial clustering sampling scenarios, this resulted in a total of 9 different sampling scenarios. We note that the third “extreme” spatial clustering scenario is completely unrealistic and one would not expect inferences of any degree of accuracy from it, but we believe it gives a good idea of the behaviour of our method in the limit case of extremely restricted spatial sampling.
A comparison of allele frequency maps generated using true parameter values and using parameter estimates from the different sampling schemes are shown in the manuscript Figure 3—figure supplements 29. In Figure 3—figure supplement 6 and Figure 3 in the manuscript we show the allele frequency map generated using the “intermediate 75%/25%” clustering scheme. Parameter estimates used to generate all these figures are summarised in Appendix 2–Table 3 in the manuscript. Overall we can see that the allele frequency maps inferred from these scenarios closely resemble the maps generated using the true parameter values, despite the challenges in finding accurate values for the individual point estimates of some of the parameters, highlighting that various combinations of diffusion and advection coefficients can produce similar underlying frequency maps (as discussed in the manuscript section “Performance on deterministic simulations”). This suggests that the joint spatiotemporal information encoded in the inferred maps (not just the individual parameters estimates) should be used in interpreting model outputs, particularly when it comes to the advection and diffusion parameters. The selection coefficient estimates are inferred highly accurately, regardless of the sampling scheme chosen, and lie close to the true value, with only a slight underestimation in the time period after 5000 years BP (with the exception of “extreme 25%/75%”).
2) How would the method perform if simulations would be performed without advection (numerical solution to equation 3) but parameter estimates would be inferred from a model C, which includes advection ?
The following text was added to the main text as section “Advection model application to nonadvection simulations”.
“We assessed the performance of the model C, which includes advection coefficient estimates on simulations generated without advection (results shown in Figure 1—figure supplement 6 and Figure 1—figure supplement 7). We can observe that the advection coefficients are inferred to be nonzero (Figure 1—figure supplement 6b and Figure 1—figure supplement 7b), however the inferred maps of spatial allele frequency dynamics closely resemble the ones obtained with true parameter values (Figure 1—figure supplement 6a and Figure 1—figure supplement 7a). This shows that complex interactions between the diffusion and advection coefficients can result in similar outcomes even when only diffusion is considered in the model.
The inference of the origin of the allele also differs when we compare the results for using model B and model C. In order to understand better how the model estimates the allele origin, we highlighted the first individual in simulations B1 and B4 that carries the derived allele. We can see that in the case of simulation B1 the inferred origin of the allele is close to the first observance of the derived allele in the model which includes advection. In contrast, when the advection is not included, the origin of the allele is inferred to be closer to where it is initially rising in frequency (Figure 1—figure supplement 1a and Figure 1—figure supplement 4a). However, this is not always the case. For instance, if we look at the results from the advection model on simulation B4, we can see that the origin of the allele is inferred relatively far from the sample known to have carried the first instance of the derived allele. Therefore, if there is a relatively large interval between the time when the allele originated and when the first ancient genomes are available, the beneficial allele can spread widely, but as this spread is not captured by any of the data points, inference of the precise origin of the selected allele is nearly impossible.”
3) Finally and more of a general question: what's the point of estimating the geographical origin of an allele if there are no samples dated from around the same time as the assumed allele age? I have the feeling this can be easily misinterpreted but maybe I am missing something. The authors discuss this (lines 199210) when discussing the results obtained for the dynamics of the rs4988235(T) for the two assumed allele ages but it is not clear to me if the estimate for the local of origin is actually meaningful under similar situations. I wonder if with simulation one could try to understand what the method is retrieving…
Indeed, inference about spatiotemporal dynamics of an allele in a time period which is not covered by any ancient samples whatsoever makes little sense beyond making extremely broad conclusions (such as, perhaps, to get an idea where it is very unlikely to have originated). Specifically, our new results described above show that in situations when the advection is disallowed in the model, the allele origin tends to be inferred to lie within the region of highest allele frequencies. When the advection is considered, the model infers the origin of the allele to be in close proximity to the very first carrier of the derived in the data set of available sequenced genomes. In other words, in the absence of any useful information, difficulties may arise when inferring this parameter when the age of the allele is much larger than the age of the oldest sample in the data set.
With respect to the presentation of the work I found that the manuscript is globally well explained but there are some parts a bit difficult to follow:
1) it should be clearly stated in the main that the age of the allele is provided as input and not inferred contrarily to the geographical origin of the allele and the value used to perform the inference should also be provided in the corresponding part of the main text.
On lines 124128 we mention that we use previously published allele age estimates as starting points for the diffusion process.
We mention the exact age estimates of the rs4988235(T) allele that we use as input on lines 239240 under the section “Dynamics of the rs4988235(T) allele”.
The allele age estimate for the rs1042602(A) allele has now been added on line 275. Thank you for bringing up this issue.
2) The comparison of the results for the two ages assumed for the rs4988235(T) mutation (Itan et al. 2009 and Albers and McVean 2020), mainly between lines 199210, is very confusing. From Figure 5a it seems there is no information on the allele for the time period between the allele age inferred by Itan et al. 2009 and Albers and McVean 2020. However, from the text it is not clearly stated whether there is or not information. If there is no information at all how does the algorithm infer an origin in Northeastern Europe?
We rephrased lines 252260. We hope that the information regarding the inference of the Northeastern European origin is now conveyed more clearly.
“Assuming the older age estimate from Albers and McVean (2020), the origin of the allele is estimated to be in the Northeast of Europe (Figure 6—figure supplement 1), which is at a much higher latitude than the first occurrence of the allele, in Ukraine. Due to the deterministic nature of the model, the frequency is implicitly imposed to expand in a region where there are no actual observed instances of the allele. The model compensates for this by placing the origin in an area with a lower density of available aDNA data and thus avoiding an overlap of the increasing allele frequencies with individuals who do not carry the derived rs4988235(T) allele (see Figure 5a). As the model expands rapidly in the southern direction (Appendix 2–Table 6) it eventually reaches the sample carrying the derived variant in Ukraine.”
3) It would also help to follow the results if the authors would clearly state in the main text (lines: 212230) the age of the allele assumed for the rs1042602(A).
We have now added in the main text the age of the rs1042602(A) allele in the line 275.
Reviewer #2 (Recommendations for the authors):
Given the ability to of the model to sample at any point in space and time, it would have been very interesting to test the model under actual geographical and temporal distributions of data, such as the Allen Ancient DNA Resource used for one of the empirical investigations, perhaps even taking the observed pattern of missingness in the data into account.
In the paragraph on lines 130134 we describe a simulation scheme in which we sampled individuals from the diffusion allele frequency maps to match the locations and ages of the individuals from the Reich lab ancient DNA data and the HGDP panel.
I would recommend using latitude and longitude systematically throughout the manuscript instead of mixing x, y, latitude, longitude.
Thank you for pointing this out, in the new version of the text, we now always refer to geographic coordinates exclusively as longitude and latitude.
In the methods: Please fix the typo in Eq (9). Also please clarify which boundary conditions where used, and how the equations were solved numerically (using Eq. (8) or (9)?), including how boundary conditions were implemented.
The following explanation was added in the chapter “Parameter search” in lines 539544.
“For numerical calculations we used the Livermore Solver for Ordinary Differential Equations (Hindmarsh, 1983) implemented in R package “deSolve” (Soetaert et al., 2010), which is a general purpose solver that can handle both stiff and nonstiff systems. In case of stiff problems the solver uses a Jacobian matrix.
Absorbing boundary conditions were used at the boundaries of the map. For visualisation purposes we masked the allele frequencies estimated in areas with negative topology (i.e. areas covered by large bodies of water).”
L437: " the lefthand side of equation 5". This doesn't look right, do you mean the righthand side of Eq. (6)?
Thank you for pointing out the typo. The line has now been corrected.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Reviewer #2 (Recommendations for the authors):
The authors have made substantial improvements to the manuscript, and is ready for acceptance except for one remaining point:
Contrary to what the authors claim in their rebuttal, Eqs (8) and (9) in Appendix 1 still contain errors.
First, p(x+\δ x,y,t) and p(x+\δ x,y,t) appears twice on line 827.
Second, the coefficients of p(x,y,t) on both sides of Eq (9) sum to 1, which is correct. The term should therefore be absent (have coefficient zero) in Eq (8).
Thank you for the comments and my apologies for the inconvenience caused by our previous response regarding this. We now recognize the error and have corrected it by removing the term p(x,y,t) from the equation (8).
https://doi.org/10.7554/eLife.73767.sa2Article and author information
Author details
Funding
Villum Fonden (00025300)
 Rasa A Muktupavela
 Fernando Racimo
Lundbeckfonden (R30220182155)
 Martin Petr
 Fernando Racimo
Novo Nordisk Fonden (NNF18SA0035006)
 Martin Petr
 Fernando Racimo
Carlsbergfondet (CF190712)
 Thorfinn Korneliussen
National Institutes of Health (R01 GM132383)
 John Novembre
Novo Nordisk Fonden (NNF22OC0076816)
 Fernando Racimo
European Research Council (951385)
 Fernando Racimo
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Graham Gower, Evan IrvingPease, Montgomery Slatkin, the members of the Racimo group, and two anonymous reviewers for helpful comments and advice. FR and RM were funded by a Villum Fonden Young Investigator award to FR (project no. 00025300). FR was also supported by a Novo Nordisk Fonden Ascending Investigator Award (NF22OC0076816) and a ERC Synergy grant (ID 951385). Additionally, MP and FR were supported by a Lundbeckfonden grant (R30220182155) and a Novo Nordisk Fonden grant (NNF18SA0035006) to the GeoGenetics Centre. TSK was funded by a Carlsberg grant (CF190712). JN was funded by NIH grant R01 GM132383.
Senior Editor
 George H Perry, Pennsylvania State University, United States
Reviewing Editor
 Aida M Andrés, University College London, United Kingdom
Reviewers
 Isabel Alves, l'institut du thorax, INSERM/CNRS, University of Nantes, France
 Anders Eriksson, Institute of Genomics, University of Tartu, Estonia
Version history
 Preprint posted: July 22, 2021 (view preprint)
 Received: September 9, 2021
 Accepted: November 21, 2022
 Version of Record published: December 20, 2022 (version 1)
 Version of Record updated: December 21, 2022 (version 2)
Copyright
© 2022, Muktupavela et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 674
 Page views

 83
 Downloads

 1
 Citations
Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Evolutionary Biology
 Genetics and Genomics
In many species, meiotic recombination events tend to occur in narrow intervals of the genome, known as hotspots. In humans and mice, double strand break (DSB) hotspot locations are determined by the DNAbinding specificity of the zinc finger array of the PRDM9 protein, which is rapidly evolving at residues in contact with DNA. Previous models explained this rapid evolution in terms of the need to restore PRDM9 binding sites lost to gene conversion over time, under the assumption that more PRDM9 binding always leads to more DSBs. This assumption, however, does not align with current evidence. Recent experimental work indicates that PRDM9 binding on both homologs facilitates DSB repair, and that the absence of sufficient symmetric binding disrupts meiosis. We therefore consider an alternative hypothesis: that rapid PRDM9 evolution is driven by the need to restore symmetric binding because of its role in coupling DSB formation and efficient repair. To this end, we model the evolution of PRDM9 from first principles: from its binding dynamics to the population genetic processes that govern the evolution of the zinc finger array and its binding sites. We show that the loss of a small number of strong binding sites leads to the use of a greater number of weaker ones, resulting in a sharp reduction in symmetric binding and favoring new PRDM9 alleles that restore the use of a smaller set of strong binding sites. This decrease, in turn, drives rapid PRDM9 evolutionary turnover. Our results therefore suggest that the advantage of new PRDM9 alleles is in limiting the number of binding sites used effectively, rather than in increasing net PRDM9 binding. By extension, our model suggests that the evolutionary advantage of hotspots may have been to increase the efficiency of DSB repair and/or homolog pairing.

 Evolutionary Biology
 Microbiology and Infectious Disease
Drug resistance remains a major obstacle to malaria control and eradication efforts, necessitating the development of novel therapeutic strategies to treat this disease. Drug combinations based on collateral sensitivity, wherein resistance to one drug causes increased sensitivity to the partner drug, have been proposed as an evolutionary strategy to suppress the emergence of resistance in pathogen populations. In this study, we explore collateral sensitivity between compounds targeting the Plasmodium dihydroorotate dehydrogenase (DHODH). We profiled the crossresistance and collateral sensitivity phenotypes of several DHODH mutant lines to a diverse panel of DHODH inhibitors. We focus on one compound, TCMDC125334, which was active against all mutant lines tested, including the DHODH C276Y line, which arose in selections with the clinical candidate DSM265. In six selections with TCMDC125334, the most common mechanism of resistance to this compound was copy number variation of the dhodh locus, although we did identify one mutation, DHODH I263S, which conferred resistance to TCMDC125334 but not DSM265. We found that selection of the DHODH C276Y mutant with TCMDC125334 yielded additional genetic changes in the dhodh locus. These double mutant parasites exhibited decreased sensitivity to TCMDC125334 and were highly resistant to DSM265. Finally, we tested whether collateral sensitivity could be exploited to suppress the emergence of resistance in the context of combination treatment by exposing wildtype parasites to both DSM265 and TCMDC125334 simultaneously. This selected for parasites with a DHODH V532A mutation which were crossresistant to both compounds and were as fit as the wildtype parent in vitro. The emergence of these crossresistant, evolutionarily fit parasites highlights the mutational flexibility of the DHODH enzyme.