## Abstract

The majority of empirical population genetic studies have tried to understand the evolutionary processes that have shaped genetic variation in a single sample taken from a present-day population. However, genomic data collected over tens of generations in both natural and laboratory populations are increasingly used to find selected loci underpinning adaptation over these short timescales. Although these studies have been quite successful in detecting selection on large-effect loci, the fitness differences between individuals are often polygenic, such that selection leads to allele frequency changes that are difficult to distinguish from genetic drift. However, one promising signal comes from polygenic selection’s effect on neutral sites that become stochastically associated with the genetic backgrounds that lead to fitness differences between individuals. Previous theoretical work has established that the random associations between a neutral allele and heritable fitness backgrounds act to reduce the effective population size experienced by this neutral allele. These associations perturb neutral allele frequency trajectories, creating autocovariance in the allele frequency changes across generations. Here, we show how temporal genomic data allow us to measure the temporal autocovariance in allele frequency changes and characterize the genome-wide impact of polygenic selection. We develop expressions for these temporal autocovariances, showing that their magnitude is determined by the level of additive genetic variation, recombination, and linkage disequilibria in a region. Furthermore, by using analytic expressions for the temporal variances and autocovariances in allele frequency, we demonstrate that one can estimate the additive genetic variation for fitness and the drift-effective population size from temporal genomic data. We also show how the proportion of total variation in allele frequency change due to linked selection can be estimated from temporal data. Overall, we demonstrate that temporal genomic data offer opportunities to identify the role of linked selection on genome-wide diversity over short timescales, and can help bridge population genetic and quantitative genetic studies of adaptation.

ADAPTATION can occur over remarkably short ecological timescales, with dramatic changes in phenotypes occurring over just a few generations in natural populations. This rapid pace of adaptive change has been known to be mirrored at the genetic level since the early work of Fisher and Ford (1947) testing whether the rapid decline in a coloration polymorphism was consistent with natural selection or genetic drift. Since then, researchers have continued to use temporal data to detect selection on polymorphisms over short timescales in natural populations (Dobzhansky 1943, 1971; Fisher and Ford 1947; Kettlewell 1958, 1961; Mueller *et al.* 1985b), as well as quantify the rate of genetic drift (Prout 1954; Wallace 1956; Nei and Tajima 1981; Pollak 1983; Mueller *et al.* 1985a; Waples 1989; Wang and Whitlock 2003). However, this line of work in sexual populations has been partially eclipsed by a vast body of work examining large-scale population genetic and genomic data sets from a single contemporary timepoint. More recently, studies have applied similar temporal approaches to whole-genome data to discover selected loci in contemporaneous natural populations (Bergland *et al.* 2014; Rajpurohit *et al.* 2018), evolve and resequence studies (Teotónio *et al.* 2009; Burke *et al.* 2010; Johansson *et al.* 2010; Turner *et al.* 2011; Orozco-terWengel *et al.* 2012; Turner and Miller 2012; Franssen *et al.* 2017; Barghi *et al.* 2019), and ancient DNA (Mathieson *et al.* 2015; Fu *et al.* 2016). Furthermore, numerous methods have been developed to estimate effective population size (Nei and Tajima 1981; Pollak 1983; Waples 1989) and detect selected loci (Malaspinas *et al.* 2012; Mathieson and McVean 2013; Feder *et al.* 2014; Terhorst *et al.* 2015) from time-series data.

Overall, these approaches have identified compelling examples where selection has driven extreme allele frequency change at particular loci that is inconsistent with drift alone. However, most adaptation on ecological timescales likely involves selection on phenotypes with polygenic architecture and abundant standing variation (Endler 1986; Hendry and Kinnison 1999; Kinnison and Hendry 2001; Kopp and Hermisson 2009b). We know from theory that adaptation on such traits can result from very subtle allele frequency changes across the many loci that underlie the trait (Bulmer 1980), at least for the short-term evolutionary response (Hermisson and Pennings 2005; Chevin and Hospital 2008; Jain and Stephan 2015, 2017; Thornton 2018; Höllinger *et al.* 2019). These changes may be individually indistinguishable from genetic drift in temporal data. This poses a challenge for population genetic approaches to quantify selection: rapid phenotypic adaptations occurring on ecological timescales may leave a signal on genome-wide patterns of diversity that is undetectable by methods focused on individual loci.

Here, we explore an alternative: rather than aiming to find selected loci, we can use temporal data to quantify the genome-wide effects of linked selection during polygenic adaptation. Linked selection introduces a new source of stochasticity into evolution, as a neutral allele’s frequency change depends on the fitness of the set of random genetic backgrounds it finds itself on (Gillespie 2000). The impact linked selection has on neutral loci is mediated by associations [linkage disequilibria/disequilibrium (LD)] with selected loci and hence their recombination environment; neutral loci tightly linked to selected sites experience greater average reductions in diversity than more loosely coupled sites. Studies using a single timepoint have long exploited this idea, with some of the first evidence of pervasive natural selection being the correlation between diversity and recombination in *Drosophila* (Aguade *et al.* 1989; Begun and Aquadro 1992). Various forms of linked selection give rise to such patterns with much attention focusing on the hitchhiking (positive selection; Smith and Haigh 1974) or background selection models (negative or purifying selection; Charlesworth *et al.* 1993; Hudson and Kaplan 1995). Recent genomic studies have modeled patterns of genome-wide diversity considering substitutions, functional constraints, and recombination environments to estimate parameters of hitchhiking and background selection models, and have begun to differentiate between these models (McVicker *et al.* 2009; Hernandez *et al.* 2011; Elyashiv *et al.* 2016). Across-taxa comparisons have shown that signals of linked selection are present in many sexual organisms, and that in some species a proportion of the stochastic change in allele frequencies is due to the randomness of linked selection instead of genetic drift (Cutter and Payseur 2013; Corbett-Detig *et al.* 2015; Coop 2016). Likewise, in asexual and facultatively sexual organisms, both theory and empirical work show that linked selection and interference are primary determinants of the levels of genetic diversity (Neher and Shraiman 2011; Neher 2013; Good *et al.* 2014, 2017).

In this paper, we extend this well-established approach of quantifying genome-wide selection through its indirect impact on linked neutral sites to temporal genomic data. We show that during rapid polygenic selection, linked selection leaves a signal in temporal genomic data that can readily be differentiated from neutral processes. Specifically, selected alleles perturb the allele frequency trajectories of neighboring neutral loci, increasing the variance of neutral allele frequency change and creating covariance between the neutral allele frequency changes across generations. Earlier work has modeled this effect on neutral alleles as a long-term reduction in the effective population size (Wray and Thompson 1990; Santiago and Caballero 1995, 1998; Woolliams *et al.* 1993; Robertson 1961), but the increasing availability of genome-wide frequency data across multiple timepoints allows us to directly quantify the extent of linked selection over short ecological timescales (tens of generations). We develop theory for the variances and covariances of neutral allele frequency change under selection, and show that analogous to diversity in a single timepoint study, their magnitudes depend on the local fitness variation, recombination, and LD. Furthermore, we show that our theory can (1) directly partition the variation in genome-wide frequency change into the components caused by drift and selection, (2) estimate the additive genetic variance for fitness and how it changes over time, and (3) detect patterns of fluctuating selection from temporal data. Overall, we believe that our approach to modeling temporal genomic data will provide a more complete picture of how selection shapes allele frequency changes over ecological timescales in natural populations, potentially allowing us to understand short-term effects of linked selection that would otherwise not be perceptible from studies using a single timepoint.

## Outline of Temporal Autocovariance Theory

Our goal is to understand how linked selection affects the frequency trajectories of neutral sites by modeling the variances and covariances of neutral allele frequency changes (, where is the population frequency at time *t*; see Table 1 for a list of all notation used). We assume a closed population with discrete, nonoverlapping generations. When there are no heritable fitness differences between individuals, genetic drift is the only source of stochasticity of allele frequency change due to two sources of variation: random nonheritable or environmental differences in offspring number, and Mendelian segregation of heterozygotes. Both are directionless such that when averaged over evolutionary replicates, the expected change in allele frequency due to drift alone is and quantified by the variance in allele frequency change (as quantified by the variance effective population size; Wright 1938; Crow and Kimura 1970; Charlesworth 2009).

When there are heritable fitness differences between individuals in the population, a third source of stochasticity affects a neutral allele’s frequency change: neutral alleles can become randomly associated with the genetic backgrounds that determine the fitness differences between individuals (Santiago and Caballero 1995, 1998; Robertson 1961). Even though the neutral alleles do not impact fitness, their frequency trajectories are perturbed by their fitness background, as those on advantageous backgrounds leave more descendants, while those on disadvantageous backgrounds leave fewer. We can partition a neutral allele frequency’s change into these three uncorrelated stochastic components [following Santiago and Caballero (1995), see Appendix section *Decomposition of Allele Frequency Change* for proof],(1)where, , and are the neutral allele’s frequency changes due to nonheritable variation in fitness between diploid individuals, Mendelian segregation of heterozygotes into offspring, and heritable variation in fitness (we refer to this as the *heritable change* in neutral allele frequency), respectively. Note that while throughout the paper we consider the allele frequency change between adjacent generations, the same approach can be extended to situations where the study system cannot be observed every generation. Like the stochastic components of drift, the allele frequency change due to heritable fitness differences is directionless . Additionally, since each component is uncorrelated with the others, the variance in allele frequency change is(2)The terms and capture the variance due to the random reproduction process, and the former can accommodate extra nonheritable variance in offspring number [as long as individuals are exchangeable with respect to their genotype (Cannings 1974)], while the term captures heritable fitness variation due to systematic differences in the fitnesses of individuals caused by their genotypes.

In addition to inflating the within-generation variance in allele frequency change, heritable fitness variation has another profound effect on neutral alleles: while the stochastic components of drift have independent effects on frequency change each generation, heritable variation in fitness creates *temporal autocovariance* in neutral allele frequency changes across generations. The contribution of temporal autocovariance is evident by writing the total cumulative allele frequency change as the sum of allele frequency changes each generation,(3)These covariance terms are expected to be nonzero only when there is heritable variation in fitness [assuming there is neither non-Mendelian segregation, nor covariance between the parental and offspring environment, *e.g.*, as in Heyer *et al.* (2005)].

Temporal autocovariance is caused by the persistence over generations of the statistical associations (LD) between a neutral allele and the fitnesses of the random genetic backgrounds it finds itself on; as long as some fraction of association persists, the heritable variation for fitness in one generation predicts the change in later generations, as illustrated by the fact that (see Figure 1A). Ultimately, segregation and recombination break down haplotypes and shuffle alleles among chromosomes, leading to the decay of autocovariance with time.

The effect that heritable variation has on neutral alleles has traditionally been modeled in a quantitative genetics framework where a large number of loosely linked polymorphisms contribute to heritable fitness differences between individuals, and the impact of heritable fitness variation on a neutral allele is quantified as a reduction in its long-run effective population size (Santiago and Caballero 1995, 1998; Robertson 1961). This form of linked selection can be contrasted with classic population genetic hitchhiking theory (Smith and Haigh 1974), which considers how neutral alleles closely linked to a new beneficial mutation are affected as it sweeps to fixation. While classic population genetic linked selection models consider how neutral variation is affected by strong associations caused by tight linkage to an advantageous site, quantitative genetic models of linked selection consider the weakest forms of associations: those between unlinked loci within an individual (Morley 1954; Santiago and Caballero 1995; Robertson 1961) [see Barton (2000) for more on the two models of linked selection]. These quantitative genetic models of linked selection match the expressions for the loss of diversity found in classic hitchhiking models under a steady flux of loosely linked advantageous alleles entering the population [see page 2110 in Santiago and Caballero (1998)] and genome-wide background selection [see equation 12 of Santiago and Caballero (1998) and Nordborg *et al.* (1996)]. The distant associations considered by these models are quickly established but are rapidly broken down by segregation and independent assortment, yet still can have a marked effect on diversity (Santiago and Caballero 1995; Robertson 1961). However, since the impact of heritable fitness variation has traditionally been modeled as causing a reduction in the effective population size, there has been no direct way to separately estimate its effects from those of drift. We show that with temporal genomic data, one can directly measure the levels of temporal variances and autocovariances of allele frequency change in the population. Additionally, we that show temporal autocovariance is created under both tight and loosely linked selection, and below, develop expressions for its magnitude that are applicable to both, bridging the two regimes of linked selection.

## A Model for Multilocus Temporal Autocovariance

Here, we develop theory for the temporal autocovariance in a neutral allele’s frequency changes through time, generated by the presence of heritable fitness in the population. We measure the temporal autocovariance at a single diallelic neutral locus. Since only allele frequency changes due to heritable variation in *fitness* contribute to temporal autocovariance, we can focus exclusively on the behavior of in deriving our expressions for the autocovariance across timepoints. We imagine that an individual *i* has fitness , *i.e.*, that their expected number of children is . We assume a constant population size, and so the population average fitness . Additionally, we assume that all fitness variation has an additive, polygenic architecture. Then, with *L* loci contributing to fitness, we can write individual *i*’s fitness as , where is the effect size in generation *t* and is individual *i*’s gene content at locus *l*. Here, each is analogous to a selection coefficient acting at locus *l*, since the fitnesses for genotypes , , and are and , respectively. This formulation is approximately equivalent to exponential directional selection on some additively determined trait, implying that selection does not create LD between unlinked loci; see *Temporal Variance and Autocovariance Under Multilocus Selection* in the Appendix for more detail.

When fitness variation exists in the population [(that is, ], the frequency of the neutral allele changes stochastically, as fitter individuals leave more descendants that inherit the neutral allele they carry and less-fit individuals leave fewer. Across the population, stochastic associations can form between the genetic components of an individual’s fitness and the neutral allele they carry, leading the neutral allele frequency to change due to fitness differences across individuals. The total heritable change in neutral allele frequency can then be partitioned into each individual’s contribution to this change based on their fitness and the number of the tracked neutral alleles they carry, , giving us(4)(Santiago and Caballero 1995). Substituting each fitness with its genetic basis and simplifying [see *Temporal Variance and Autocovariance Under Multilocus Selection* in the Appendix for derivation and equation 10 of Kirkpatrick *et al.* (2002)] gives(5)where is the *gametic* LD between the neutral allele and the allele at the selected site *l* on the same gamete, whereas is the *nongametic* LD, or the covariance *across* the neutral and selected allele on the two different gametes forming an individual [see p. 121 of Weir (1996) for details and Appendix Figure A2A for an illustration]. Intuitively, this expression tells us that the heritable change in neutral allele frequency is determined by the gametic and nongametic LD between the neutral site and all sites that affect an individual’s fitness, scaled by the magnitude of each selected locus’s effect. Alternatively, we can see this as the multivariate breeder’s equation, where the neutral allele is a correlated trait responding to selection on other traits/loci (Lande 1979). This expression is the multilocus analog of the change in a neutral site’s frequency due to hitchhiking at a single linked site [*e.g.*, see equations 2 and 3 in Stephan *et al.* (2006)].

Because the effects of the nongametic LD are relatively weak compared to the gametic LD for tightly linked loci (see *The Strength of Unlinked and Nongametic Associations* in the Appendix for an expression of their strength), we ignore these, and hereafter omit the primes in our notation so that refers to . Since , we can write the covariance as . Hereafter, we also omit the subscript *H* since . Expanding these terms, the covariance between the allele frequency changes at generations *t* and *s* can be written as This statement for temporal autocovariance is fairly general, as it can handle fluctuating selection (*e.g.*, when varies with time *t*) and any additive multilocus evolution (as long as the LD dynamics can be specified). Looking at the first term in this sum, we see that the temporal autocovariance is determined in part by the terms . These expected LD products reflect the degree to which the association between the neutral locus and a selected site persists from generation *t* to *s* (here, ). Intuitively, the higher the initial association between the neutral and selected loci, and the slower the rate of decay of LD between sites, the greater temporal autocovariance will be.

### The multilocus temporal autocovariance model with directional selection

Thus far, in reaching Equation 6 we assume only that fitness is additive across loci. In this section, we develop a model of how temporal autocovariance behaves specifically under directional selection beginning at a specific time. We make three assumptions to simplify our expressions. First, we assume that the effect size remains constant through time, such that for all *t* (we relax this assumption in *Fluctuating selection*). Second, we ignore the contribution of the second term of Equation 6, (for ), to temporal autocovariance. Under the case where the population is initially at mutation–drift–recombination equilibrium, we expect this product to be zero as there is no directional association between the two selected sites and the neutral site. However, we note that interaction between selected sites [Hill–Robertson interference (HRi)] will cause this term to become negative (Barton and Otto 2005), a point we return to later. Third, we assume that the selected sites increase in frequency independently, such that the dynamics of the LD between the neutral and selected site pairs can be modeled using two-locus dynamics. Using a deterministic continuous-time model for the dynamics of the LD between the selected and neutral site (Smith and Haigh 1974; Barton 2000), we rewrite the terms in the expression for temporal autocovariance as(6)(7a)(7b)where is the square of the correlation between the neutral site and selected site *l* at time *t* (a common measure of LD; Hill and Robertson 1968), and is the recombination fraction between the neutral site and selected site *l*.

We can further simplify this expression by assuming that there is no covariation between the additive genic variation at a selected site, and the LD between that selected site and the neutral site (see *Temporal Variance and Autocovariance Under Multilocus Selection* in the Appendix for more detail). This allows us to factor out the average additive genic variation for fitness at time *s* and write the covariance as(8)where is the additive genic variance for fitness, which is the additive genetic variance for fitness without the contribution of LD between selected sites, . In part, our expression not relying on the LD between selected sites is a result of ignoring the second term in Equation 6; we revisit the consequences of this assumption further on in *Comparing theory to simulation results*.

This expression allows us to calculate the temporal autocovariance in cases where we know the vector of recombination fractions between the neutral and each of the selected sites, . Often we do not know the exact positions of these sites, but we can treat these positions as randomly placed on a chromosome and further simplify our model to understand the factors that determine temporal autocovariance.

### Temporal autocovariance for an average neutral polymorphism

In the second part of our derivation, we develop a simple intuitive model of how temporal autocovariance is determined by a few key parameters when we make two additional assumptions. First, we assume that selected sites are randomly and uniformly distributed along the chromosome, such that a site’s position on the genetic map is a random variable (where *R* is the region’s length in morgans), and the focal neutral site with which we calculate temporal autocovariance lies in the middle of this idealized chromosome at the origin (as depicted in Figure 1B). Then, the recombination fraction between the focal neutral site and a selected site at random position *g* is given by the mapping function , which maps the position *g* to a recombination fraction. A simple choice for is Haldane’s mapping function, (Haldane 1919) (note we take the absolute value of *g* to translate the position *g* to a distance to the focal neutral site), and we use that here. Second, we assume that the LD between each selected site and the focal neutral site depends only on the recombination fraction between the two loci, and not their absolute positions or effect sizes; then, we rewrite as the function . For example, if the population was initially at drift–recombination balance, this would be where (Hill and Robertson 1968; Ohta and Kimura 1969). These assumptions allow us to conceptually understand the factors that determine temporal autocovariance; in practice, in temporal studies with LD data and recombination maps, one can directly calculate the sum in Equation 8 (see *Empirically Calculating the Average LD Persisting Across Generations* in the Appendix). We then write the temporal autocovariance experienced by a neutral allele in a region *R*-M long containing fitness variation at time *s* as(9)(see *Temporal Variance and Autocovariance Under Multilocus Selection* in the Appendix for details).

This integral is the sum of the initial LD between a typical neutral locus and a selected site, weighted by the decay of LD due to recombination over generations. Selection enters here through the total additive genic variance for fitness for the region divided by the genetic map length of the region . Thus, a key compound parameter in describing the temporal covariance is the additive genic variance per morgan, a quantity somewhat similar to the ratio of new adaptive mutations per base pair to recombination per base pair, , that occurs in models of recurrent sweeps (Stephan *et al.* 1992) and models of the limits of selection with linked loci (Robertson 1970, 1976). Note that this does not include the effects of genome-wide fitness variation, *e.g.*, the impacts that unlinked selected sites have on the neutral site due to the associations created when the sites sort within the same individuals. We quantify the magnitude of these in *The Contribution of the Rest of the Genome to Temporal Autocovariance at a Locus* in the Appendix.

To validate our theory, we simulate a fixed region of *R* M and calculate the covariance in allele frequency changes by averaging over many uniformly distributed neutral sites within this region. Then, the random distance between a neutral site’s position *n* and a selected site’s position *g* is , where ; this random variable *c* has a triangle distribution, . Averaging over the positions of both randomly placed neutral and selected sites, the temporal autocovariance is(10)where indicates an expectation taken over the position of the randomly placed neutral sites, and we define as the average LD between selected and neutral sites that persists from generations *t* to *s* . As is common with estimating the expected values of other ratios like (Bhatia *et al.* 2013), we use a ratio of expectations rather than the expectation of the ratio.

We can also use this expression to calculate the variance of allele frequency change. The standardized variance has two components: the drift term and the heritable variance in offspring number. Adding these independent contributions, the standardized variance is(11)where is the nonheritable variance in offspring number. Under a Wright–Fisher model of reproduction, , this simplifies to(12)When combined, this expression for the variance in allele frequency change and our expression for temporal autocovariance are in agreement with Robertson (2009) and Santiago and Caballero (1995; 1998) when predicting the total variance in allele frequency change (see *Connecting our Model with the Models of Robertson and Santiago and Caballero* in the Appendix). With the above expressions for the variances and covariances, we have a complete set of theoretical expressions for the variance–covariance matrix of allele frequency change, which we call **Σ**, with the diagonal variance elements given by Equation 12, and the upper- and lower-triangle covariance elements given by Equation 10.

### Modeling the dynamics of additive genic and genetic variation

Our expressions for temporal autocovariance (Equation 10) require an expression for , the additive genic variation through time. However, we lack general expressions for the dynamics of the additive genic variation during selection, as these dynamics are quite complex for a few reasons. First, since our theory considers polygenic selection at a finite number of loci in a region, additive genic variation is not as constant as it would be under an infinitesimal model (Bulmer 1980). Second, we allow for arbitrary levels of recombination from very tight linkage to loose linkage. Previous work has shown that predicting the dynamics of additive genic variation in a system with an arbitrary level of recombination is difficult, as both the additive genic and genetic variances depend on the higher-order moments of LD [see Barton and Turelli (1987), p. 607 of Turelli and Barton (1990), and Barton (1991)].

A primary determinant of the additive genic variation is the heterozygosity of the selected sites. Assuming effect sizes are constant through time and across loci, we can rewrite the additive genic variation, , as(13)(14)where is the sum of site heterozygosity at time *t*. Ideally, we would directly use in a region; however, this would require knowing *a priori* which sites are being selected. Instead, we assume that the trait is sufficiently polygenic that frequency changes due to selection are weak, and that the change in heterozygosity at neighboring neutral polymorphic sites approximately mirrors that at selected polymorphisms (this is the case under the infinitesimal model, where the change in frequencies due to selection is no different from the change due to drift) (Robertson 1960; Bulmer 1980, Kimura 1984). Then, using the sum of site heterozygosity at neutral sites, , as a proxy for the sum of site heterozygosity at selected sites,(15)where we define as the factor by which decreases at time *s*, approximated by neutral sites’ allele frequency changes. Under this approximation, the dynamics of genic variation are determined by one free parameter, , and the directly measurable sum of site heterozygosity at neutral sites through time.

Our focus here is on the short-term response of a population, and so we look at the decay of genetic backgrounds present at the onset of directional selection. In reality, new mutations consistently create additive genetic variation for fitness; thus, an equilibrium level of additive genetic variance in the population can be maintained. The long-run effect of linked selection under this equilibrium model is handled by Santiago and Caballero (1995, 1998); see *Connecting our Model with the Models of Robertson and Santiago and Caballero* in the Appendix.

### Multilocus simulation details

To test our theoretical expressions, we have conducted extensive forward simulations of directional selection on a polygenic trait. We vary four critical parameters in these simulations: (1) the level of additive genetic variance at the onset of selection , (2) the level of recombination (*R* in morgans), (3) the number of selected sites in the region (*L*), and (4) the population size (*N*). We choose our grid of the selection and recombination parameters based on the levels we would expect across a wide variety of organisms; see the *Multilocus Simulation Details* section in the Appendix for details. We used three different population sizes , but note that we use and a subset of the other parameters in our figures.

Before the onset of selection, we create the initial diploid population from a pool of gametes created by msprime (Kelleher *et al.* 2016), such that the initial allele frequency distribution and LD between sites is at mutation–drift–recombination balance. Details of how msprime was called are available this paper’s code repository (https://github.com/vsbuffalo/tempautocov) in R/simpop.r. Then, we pass this pool of gametes into a forward Wright–Fisher-with-recombination simulation routine and let it evolve for four generations neutrally before initiating selection on the fifth generation. These first four generations of neutral evolution (without mutation) serve as a control to validate that the variance in neutral allele frequency change is as expected under a Wright–Fisher model, and that temporal autocovariance between a generation before selection and during selection is zero.

We generate genetic variation for fitness by choosing *L* random loci from the neutrally evolved sites, and randomly assign an effect size of or , such that the expected total amount of additive genic variation is (note that the initial additive genic and genetic variance are equal, , as the LD contribution is zero for randomly chosen sites). The details of this are given in the *Multilocus Simulation Details* section in the Appendix. This approach creates some additional variance around the target level of additive genic variation, as the sum of site heterozygosities will vary stochastically across simulation replicates. At the onset of selection, an individual *i*’s trait value is calculated as where is their number of alleles with effect size *α* at locus *l*. Then, their absolute fitness is calculated using an exponential fitness function [see p. 17 in Turelli and Barton (1990)]. Under our Wright–Fisher model, we sample the parents of the next generation according to a multinomial distribution, where the probability of individual *i* being a parent is .

We record 50 generations of simulated evolution, after which we compute the standardized sample temporal variance–covariance matrix **Q** (this is the sample analog of our theoretical variance–covariance matrix **Σ**), for each replicate as follows. First, we mark frequencies reaching fixation or loss as missing values. This allows the frequency changes before fixation/loss to contribute to the measured covariance, rather than removing the entire locus’s trajectory, which would act to condition the covariance on more intermediate frequencies. Note that one cannot ignore fixations or losses, as these have and thus an autocovariance of zero, which would act to underestimate the true level of autocovariance at segregating sites. Having marked fixations/losses as missing, we take the frequency matrix and calculate a vector of allele frequency changes using each neutral locus *n*’s observed generations. Finally, we calculate the sample standardized variance–covariance matrix **Q**, averaging over *M* neutral loci such that element is calculated as(16)though see *Accounting for Allele Frequency Sampling Noise* in the Appendix for a bias-corrected version when samples rather than population allele frequencies are used. Sums over missing values only use pairwise-complete observations, implemented by R’s cov() function’s use = ’pairwise.complete’ argument.

We have extensively validated our simulation procedure in a neutrally evolving population, ensuring that the decay of LD and the allele frequency change match expectations (see Supplemental Material, Figures S1.1, S1.2, and S1.3).

### Comparing theory to simulation results

To validate our expressions for temporal autocovariance, we compare the levels of autocovariance and variance predicted by Equation 10 and Equation 12 to the average levels observed across simulation replicates. To calculate the theoretical values of temporal autocovariance and variance, our expression requires the additive genic variation at *s*, ; however, as described in *Modeling the dynamics of additive genic and genetic variation*, we lack an analytic expression for the dynamics of genic variation to plug into . Following the approach of others in evolutionary quantitative genetics [see p. 930 in Turelli and Barton (1994)], we substitute the numerical values calculated directly from the simulation data for . Additionally we consider two other numerical values related to the additive genic variance: the observed additive genetic variance from our simulations [ at time *s*, which includes the contribution of LD between selected sites], and the additive genic variation at time *s* as approximated by the observed decay in the sum of site heterozygosity at neutral sites [, as described in *Modeling the dynamics of additive genic and genetic variation*].

Figure 2 compares the fit of our theory with differing additive genetic variances with the empirical covariances from our multilocus simulations. In each panel, we plot the level of temporal autocovariance between the allele frequency change across the first two generations of selection and some later allele frequency change where *s* varies along the *x*-axis. Each point represents the temporal autocovariance (calculated across all sites in a region according to Equation 16) averaged across 100 replicate simulations, with the color of the point indicating the number of selected sites in the region. Within each panel, the temporal autocovariance predicted by Equation 10 is plotted as a set of three lines, one for each of the three different types of variance we have substituted in for . Overall, the fit is close but varies depending on the type of variance used for ; we discuss each in turn below.

Using empirical additive genic variation (solid lines), our theory provides a good fit to the simulation results for a short period after selection is initiated (around five generations) in regions with tighter linkage across a range of additive genetic variation parameters (; see Figure S3.2 for varying over orders of magnitude). With looser linkage ( M), our theory using the empirical additive genic variation fits much more closely over a longer duration (∼10–15 generations). Note that some variability is caused by the noise of each simulation replicate around the target initial additive genetic variation (*Multilocus simulation details*), as each replicate samples sites from a neutral coalescent. Our theory also accurately predicts the temporal autocovariance for different choices of reference generation, *i.e.*, varying *t* (see Figure S3.1).

When we use the sum of site heterozygosity at neutral sites (, shown as a short-dashed line) as a proxy for additive genic variation, the theory fits simulations over the same time span as using the empirical additive genic variation. This is because: (1) the sum of site heterozygosity at neutral sites closely matches the sum of site heterozygosity at selected sites and (2) both closely follow the dynamics of additive genic variation through time (see Figure S2.1). Using has the advantage that we can directly measure the neutral sum of site heterozygosity, which proves useful later in *Estimating Linked-Selection Parameters from Temporal Autocovariance*, as we use this approach to help infer the initial additive genic variation at the onset of selection.

Finally, we find that using the additive genetic variance accurately predicts the dynamics of temporal autocovariance over tens of generations (see the long-dashed lines in Figure 2). Furthermore, calculating the temporal autocovariance using the empirical additive genetic variation better fits simulation data in regimes with tight recombination, where using genic variation performs poorly after the first few generations (*e.g.*, the column of panels where ). Thus, using the additive genetic variance in our framework provides a good fit to the temporal dynamics over relatively long time spans.

What differentiates from that could explain this better fit? The additive genic variation ignores the contribution of LD between selected sites. We can write the additive genetic variance as(17)where is the LD between selected sites *i* and *j* at time *s*. At the onset of selection, there is no expected LD between selected sites since the sites and effect sizes were randomly sampled; in other words, . We see this in Figure 2, as the temporal autocovariances predicted with match those of when (see also Figure S2.1, which plots the empirical additive genic and genetic variances over time). Over time, these two quantities diverge as negative LD build up. While negative LD between selected sites build up due to epistasis under some forms of selection (known as the Bulmer effect) (Bulmer 1971, 1980), this is known not to happen under multiplicative selection [see p. 50 and p. 177 in Bürger (2000)] that is equivalent to the exponential directional selection fitness surface we have used in our simulations. Instead, the buildup of negative LD between selected sites is likely due to HRi between selected sites (Hill and Robertson 1966), which affects the total additive genetic variation that selection is acting on. HRi refers to the creation of negative LD among beneficial alleles in a finite population resulting from the fact that beneficial alleles that are on the same haplotype move more quickly through the population than beneficial alleles on deleterious backgrounds, resulting in negative LD. This negative LD among beneficial alleles lowers compared to the genic (Hill and Robertson 1966; Barton and Otto 2005; Good *et al.* 2014; Crouch 2017). In the derivation of our expression for temporal autocovariance, we greatly simplified the multilocus dynamics by ignoring the second term in Equation 6. This term includes the expected product of two LD terms; each is the LD between the neutral site and a selected site. Using full multilocus theory, one may find that, by including these LD products, rather than factors out the expression in Equation 7, but we leave this for future work. Importantly, our simulation results suggest that the negative LD created by selective interference only affects the temporal autocovariances through the variance term , and that the actual variance determining temporal autocovariance is the additive genetic variance, .

In addition to modeling autocovariance through time, our theory can predict the total temporal variance in allele frequency, , when there is heritable variation for fitness. Furthermore, from Equation 3 recall that we can decompose into variance and covariance components. The variance components are determined by both the magnitude of drift and selection according to Equation 11, and the covariance components are determined solely by selection according to Equation 10 (assuming no inheritance of environmental factors). Using our theory, we have predictions for each of these components given the amount of additive genic/genetic variation for fitness, the population size (*N*), and the amount of recombination (*R*). In Figure 3, we compare the magnitudes of these components (averaged over the replicates of our simulations) to our theoretical predictions. We depict the predictions for the variance and covariance components using both the empirical additive genetic variance and the neutral sum of site heterozygosity proxy as adjacent bars, each around a point range with the point representing the average value over simulation replicates, and the bars indicating the lower and upper quartiles over simulations.

Finally, we have found that across a wide range of recombination and additive genetic variation parameters, the temporal autocovariance is largely determined by the compound parameter , and the number of generations between *t* and *s*, which is a factor in Equation 9. We show in Figure 4 that the temporal autocovariance from simulations across a wide range of and *R* parameters fall roughly on the same curve for each number of elapsed generations.

### Data availability

All code to reproduce these results is available on GitHub at https://github.com/vsbuffalo/tempautocov. Larger simulation data sets used to create figures are in the Supplemental material available at FigShare: https://doi.org/10.6084/m9.figshare.7709930.

## Estimating Linked-Selection Parameters from Temporal Autocovariance

Our multilocus theory provides analytic expressions for the expected variances and covariances of a neutral allele’s frequency; thus, a natural approach to parameter estimation is to equate these expectations to averages from the data and apply the method of moments. We describe a method-of-moments procedure below to estimate the initial additive genetic variance at the onset of selection in the first generation and the drift-effective population size (*N*) from temporal data within a single region *R*-M long, and then show a simple extension that allows this to be applied to genome-wide data. Our basic approach is to first calculate the sample variances and covariances of the *τ* observed generation-to-generation allele frequency changes, averaging over all of the putatively neutral sites in a region. We then equate these sample variances and covariances to our analytic expressions for the variances and covariances, leaving us with an overdetermined system of equations, which we solve using least squares. We demonstrate that this simple estimation procedure provides accurate estimates of initial additive genetic variance and the drift-effective population size. We focus on this procedure, as it is simple and handles incomplete trajectories due to missing data or fixation/loss well. Calculating pairwise-complete covariances can leave sample covariance matrices nonpositive definite, which makes maximum likelihood estimation perhaps much more difficult. Throughout, we use population allele frequencies (*i.e.*, there is no sampling noise). which simplifies the description of the method; in Appendix section *Accounting for Allele Frequency Sampling Noise* we describe how the method is changed by finite sampling of chromosomes from a population.

From our multilocus theory, we have analytic expressions for each element of the covariance matrix of allele frequency changes in a region. To model the additive genic variance through time, we use the empirical neutral sum of site heterozygosity approximation as described in *Modeling the dynamics of additive genic and genetic variation*. This approximates the rate that the additive genic variation decreases through time from some initial level, , which we wish to estimate. In total, we have unique moment equations, which for the variance and covariance are defined as(18)(19)Here, the first line gives the form of *τ* equations for the variance of allele frequency changes between subsequent generations, which includes the effect of genetic drift, . The second line gives the form of the covariances of allele frequency changes among different generations. The term is the average level of LD after the generations that have elapsed, given there are *R* M of recombination. In our multilocus theory section and simulations, this is equal to the integral in Equation 10. However, we can also directly calculate a sample from observed LD in a region (for details, see Equation 55 in the Appendix).

Following the method of moments, we equate each of these independent equations for to the observed sampling moments, the elements of the upper triangle of the observed heterozygosity-normalized covariance matrix described in Equation 16. This yields equations with two unknown parameters: and . We solve this overdetermined system of equations using least squares, an approach similar to the generalized method of moments in econometrics (Hansen 1982). This approach finds parameter estimates that minimize the squared error between the moment-based parameter estimate and the true parameter value, with respect to the true parameter value. We write the elements of the upper triangle of the observed covariance matrix in the vector , and write the method-of-moments equations as,(20)where the elements of and , in the same order as , are given by(21)where is an indicator variable that is one when and zero otherwise.

Then, we can readily estimate the parameters and using least squares. We then obtain an estimate of by taking . Since these equations are not statistically independent, we cannot assume . However, this does not affect our estimates and , as the least squares procedure is unbiased regardless of the covariance structure between the error terms [see p. 26 in Christensen (2011)].

Using this method-of-moments approach, we sought to infer the parameters of 20 of the replicates across the 254 parameter combinations (the same as used in Figure 2). We use the first five generations after the onset of selection to infer and , as for this short time span the additive genetic variance is well approximated using the sum of neutral site heterozygosity approach (see *Comparing theory to simulation results*). Each simulation replicate includes ∼500 neutral sites (the exact number is random, see *Multilocus Simulation Details* in the Appendix for details).

Applying our approach to these simulations, we find that we can infer both the initial level of additive genetic variation and the effective population size *N* from multilocus temporal data. In Figure 5A, we show that our method of moments gives reasonable estimates for the initial level of additive genetic variance over orders of magnitude of additive genetic variation, and different recombination regimes. As additive genetic variation for fitness becomes weaker (the left side of the figure), our estimates become more noisy. In Figure 5B we show the simultaneously estimated population size *N* against the true population size value. We also plot the estimated (not accounting for selection) from a simple temporal estimator, (Krimbas and Tsakas 1971; Waples 1989), where *F* is Wright’s standardized variance (Wright 1931). While with high and low *R* the method-of-moments approach still underestimates *N*, it performs far better than a standard temporal estimator that does not account for selection. In Appendix Figure A3, we include a version of this figure calculated using the method of moments on sample allele frequencies for a sample of size chromosomes.

We can extend this approach to whole-genome data by imagining partitioning the genome into *B* nonoverlapping windows of length in base pairs (*e.g.*, megabase windows). We first assume that windows contribute uniformly to the genome-wide level of additive genetic variance for fitness, and show how our method-of-moments approach can be used to estimate a global . Assuming a uniform distribution of genetic variance across base pairs, the total additive variance is , across our *B* windows, where is the additive genetic variance per base pair. As each window *i* contributes to , our least squares approach given by Equation 20 becomes(22)However, we expect *a priori* that windows containing more coding bases might disproportionately contribute to the total additive genetic variance. This suggests an alternative model to fit where partitions of the additive genetic variance across windows are proportional to the number of coding bases, similar to background selection and other linked selection models (McVicker *et al.* 2009; Rockman *et al.* 2010; Corbett-Detig *et al.* 2015). Thus, we could write total where is the number of coding or exonic base pairs in window *i* (this could be any quantifiable annotation feature in the window), and is the total number of coding bases in the genome. With window *i* contributing to the additive genetic variance and having map length , we now define , and as having elements given by the equations(23)Again, the parameters of this model, and , can be estimated with least squares. When analyzing genome-wide data, these various models could potentially be compared to an out-of-sample procedure, using inferred parameters to estimate the mean-squared predictive error between the two models for the remaining windows (Elyashiv *et al.* 2016). The C.I.s for our method-of-moments estimates could be obtained through bootstrapping genomic windows since the errors are not identically and independently distributed.

### Estimating the proportion of variance in frequency change due to linked selection

We can also estimate what fraction of allele frequency change over *t* generations is due to linked selection acting to perturb the frequency trajectories of neutral alleles. We have developed two approaches: first, a more conservative approach that considers only the contribution of selection to the temporal autocovariance, and second, a more exact approach that uses the estimated effective population size to include the contribution of selection to both variances and covariances of allele frequency change.

First, a simple estimate of the total fraction of the variance in allele frequency change (*G*) caused by linked selection is(24)However, this estimator is conservative because it ignores the contribution that linked selection has on the variance in allele frequency change across a single generation [the term in Equation 2]. If we include these variance terms, we have a less-conservative estimator that we call :(25)We can think of the numerator of the second term in Equation 25 as the variance in allele frequency change in a Wright–Fisher population without selection. Recall that under a Wright–Fisher model, the standardized variance across *t* generations is approximately , where this second approximation works for short time spans . This suggests that we can use our method-of-moments estimate of the effective population size without the effects of selection, , and compare the fraction of standardized variance we expect under this rate of drift to the empirical standardized variance,(26)Figure 6 shows the estimated from the method-of-moment estimates using 20 replicates of simulated data. We learn three important points about our estimator. First, for low or , the estimator is quite noisy. Second, although the signal can be noisy for low , the relationship between and the level of recombination is consistent with selection affecting the total variance in allele frequency changes across the genome. Finally, this suggests that a negative relationship between and recombination rate, calculated in windows across the genome, is a robust signal of linked selection impacting the total variance in allele frequency change. Furthermore, we would expect a positive relationship between the number of coding base pairs per window (when such information is available) and , which could serve as another robust signal of linked selection impacting the total variance in allele frequency change.

### Fluctuating selection

Thus far, we have assumed that fitness effect sizes are constant through time, that is , for all . In natural populations, changes in the environment or composition of the population may cause these effect sizes to change through time, due to changing selection pressures and changes in the epistatic environment experienced by alleles. If these changes occur within the timeframe of recorded allele frequency changes, the levels of temporal autocovariance will differ from the levels predicted from our directional selection theory. However, from Equation 6 we can see that the magnitude of temporal autocovariance is determined in large part by .

Here, we discuss how temporal autocovariance behaves under an example of strong fluctuating selection: when selection on a trait changes direction at some point. Specifically, we change the fitness function to after some timepoint ; this is equivalent to changing if , and otherwise for all other .

When such a strong change in the direction of selection occurs, the temporal autocovariance between timepoints before and after the change becomes negative, since temporal autocovariance is determined by the product for (here we are holding effects constant across loci). We have validated this using the same simulation procedure as described in *Multilocus simulation details*, except that on generation 15 we reverse the direction of selection on the trait by changing the fitness function to . In Figure 7A, we show the temporal autocovariance for varying *s* along the *x*-axis (in this case, , and ). During the first five generations, the temporal autocovariance behaves as it does under directional selection, decaying due to the decrease in additive genetic variance and the breakdown of LD. Then, on the 15th generation, the direction of selection on the trait with breeding value reverses, and temporal autocovariance becomes negative since for all and . Under this simple flip in the direction of selection pressure, the genic variance, , in the expressions for the temporal autocovariance can be replaced with , akin to a genetic/genic covariance. In Figure 7A, the gray line is our predicted level of temporal autocovariance proportional to given by Equation 8 before generation 15, and after that generation it is proportional to (using the empirical additive genetic variance). However, note that the dynamics of additive genetic variance under fluctuating selection are more complex than under directional selection. Whereas under directional selection the genetic variance decays as selection proceeds, under fluctuating selection there can be a transient increase in the additive genetic variance (seen in Figure 7A between generations 15–24). This transient inflation of the additive genetic variance is caused by the increase in the heterozygosity of haplotypes that experience reduced heterozygosity due to directional selection. With the direction of selection reversed, previously selected haplotypes move to more intermediate frequencies, which increases the additive genetic variance until selection proceeds and this variance decays (generations 25 and onwards). Overall, the dynamics of additive genetic variance under fluctuating selection are more complicated than under directional selection, which makes inference using our sum of site heterozygosity approximation infeasible.

Since fluctuating selection can create *negative* temporal autocovariance, the total amount of autocovariance over time [*e.g.*, ] can misrepresent the actual amount that linked selection is affecting allele frequencies on shorter timescales. In turn, this leads our estimator *G*, for the fraction of variance in allele frequency due to linked selection, to underestimate the total contribution of linked selection to variation in allele frequencies over the time period. We show an example of this in Figure 7B, which depicts the total variance in allele frequency change through time, partitioned into variance and covariance components, and colored according to whether the generation was before or after the reverse in directional selection. The covariances increase as they would under directional selection, but after generation 15, the contribution of covariance begins to decrease as negative autocovariances accumulate. By generation 20, the total covariance terms have a net negative effect, and actually act to decrease the total variance for a few generations below the constant level expected under drift and heritable variation.

To more fully capture the contribution of selection allele frequency change, we modify *G*, using the absolute value of the covariances(27)which prevents negative temporal autocovariances from canceling out the effects of positive temporal autocovariances, since both are a reflection of linked selection acting on neutral allele frequency changes. We show in Figure 7C, where even after the change in the directional selection we see a steady accumulation of covariance contributing to total variance . This also suggests that one plausible way to check for genomically widespread fluctuating selection would be to test if . However, we note that in contrast to our simulations, it is likely that in natural populations only a subset of alleles change their relationship to fitness. This may act to dampen the magnitude of, but not completely reverse, the direction of genome-wide temporal autocovariance, and so different approaches may be needed to identify fluctuations.

## Discussion

Currently, the prevailing empirical approach to studying linked selection relies on using samples from a single timepoint, and modeling the patterns of diversity subject to different functional constraints and in different recombination environments. The early theoretical work underpinning empirical analyses of linked selection’s effects on diversity were primarily full sweep, recurrent hitchhiking models, where beneficial mutations arise in the population and then sweep to fixation (Smith and Haigh 1974; Kaplan *et al.* 1989; Stephan *et al.* 1992). Furthermore, by looking at the patterns of diversity around amino acid substitutions or the site frequency spectrum in low-recombination regions, researchers have teased apart the effects of background selection and hitchhiking in *Drosophila* (Begun *et al.* 2007; Elyashiv *et al.* 2016), humans (McVicker *et al.* 2009; Hernandez *et al.* 2011), and some plant species (Nordborg *et al.* 2005; Schmid *et al.* 2005; Williamson *et al.* 2014; Beissinger *et al.* 2016). Yet, as other theoretical models of hitchhiking incorporating changes in the environment (Kopp and Hermisson 2007, 2009a,b), sweeps originating from standing variation (Hermisson and Pennings 2005), and multiple competing beneficial haplotypes (Pennings and Hermisson 2006) have been developed, it has become rather more difficult to detect the signals of these hitchhiking phenomena in empirical data. We have proposed here that temporal autocovariance offers a unique and measurable signal of linked selection over shorter timescales that provides a fuller picture of the ways in which genome-wide diversity has been affected by these other hitchhiking phenomena.

### Empirical applications and future directions

Here, we have developed expressions for temporal variances and autocovariances, and applied these to model temporal variances and covariances during directional selection on a trait. We have demonstrated how one can: (1) estimate the additive genetic variance for fitness and the drift-effective population size, (2) estimate the fraction of variance in allele frequency change due to linked selection, and (3) evaluate whether fluctuating selection is operating from these temporal variances and autocovariances. However, we recognize a series of limitations and difficulties when applying these methods to empirical data and natural populations.

First, one difficulty with temporal sampling of natural populations is the risk that the genetic composition may change drastically due to migration or biased sampling across timepoints. Since our theory assumes a constant-sized and isolated population, migration into the sampled population presents a serious potential confounder. For example, seasonal migration could create an influx of new alleles that could at best dampen signals of directional selection across seasons, or at worst create a signal of artificial covariance between like seasons. Similar biases could occur if a sampling method incidentally preferred certain subgroups in a stratified population. This could occur, for example, if individuals differed in some behavior affecting their likelihood of being sampled across different temporal environments, which could cause spurious covariances. While such sampling issues and migration might be able to be detected *post hoc* from exploratory data analysis approaches like PCA, studying isolated natural populations and carefully designing sampling schemes would lead to the best inference. Importantly, the effects of gene flow and other temporal inhomogeneities could differ across recombination environments, as among-population differentiation will be more pronounced in regions of low recombination and high functional density (Keinan and Reich 2010; Nachman and Payseur 2012; Burri 2017). In situations where migration is a factor, one way forward might be to study the contribution of linked selection after partitioning out the effects of gene flow across recombinational and functional environments, by extending admixture inference approaches that estimate the genome-wide effect of drift and admixture.

Second, our method-of-moments approach relies on assuming that we can approximate the dynamics of the decay of the additive genetic variance that was present at a reference timepoint by using the observed changes in the sum of site heterozygosity in a region as a proxy for its decay. While we have shown that this model works under directional selection in a relatively idealized setting, inference in natural populations may be complicated by changes in the environment that induce the effect sizes across loci to vary across timepoints. While our fluctuating selection results show that our directional selection theory extends to changes in the direction of selection with minor adjustments (*e.g.*, Figure 7), having the effect sizes vary across each generation would complicate the dynamics of additive genetic variance through time and make inference of difficult. Similarly, we assume that effect sizes are constant across sites. Variance in effect sizes across a region will not bias our results unless there is covariance between effect size and local recombination rate. Further work is required to develop statistical methods to test for violations of our assumptions about effect sizes remaining constant across time, and to potentially incorporate these complications into inference.

Along similar lines, we have focused on directional selection under a multiplicative model. However, selection experiments, a natural place to apply our method, often use truncation selection, which generates systematic epistasis for fitness and thus LD among loci (Bürger 2000; Walsh and Lynch 2018). Similarly, in natural populations stabilizing selection will act on many traits, which can also generate LD among loci. These selective processes will act to rapidly reduce the additive genetic variance for fitness across time, especially in low-recombination regions, and reduce the initial additive genetic variance in low-recombination regions. Simulating the effect of truncation selection and stabilizing selection on temporal covariances in allele frequencies seems a useful direction. Our model may prove to be a useful null model of selection that the complications of epistasis and different dynamics of the additive genetic variance could be tested against.

Finally, while we have demonstrated that our method-of-moments estimation approach is simple and leads to unbiased estimators, we see opportunities for simple extensions and other inference procedures. Differences in recombination rates and coding density are relatively easily accommodated (Equation 23). In fitting our model, we assumed a parametric form to the initial LD between neutral and selected loci; however, in practice the initial LD between neutral polymorphisms and putatively functional sites could be estimated empirically. Using the empirical LDs in Equation 23 could make the inference somewhat analogous to LD score regression (Bulik-Sullivan *et al.* 2015). The LD score of a SNP is simply the sum of the s of the focal SNP to all SNPs within some large physical genomic window, which could be used in place of the integral in Equation 10. Using this equation, could be estimated by regressing the temporal covariance of a SNP on its LD score. An alternate approach would be to use likelihood methods to model each neutral site’s frequency changes using the set of pairwise LD and recombination distances between the neutral site and all neighboring polymorphic sites. Then, genome- or region-wide estimates could be found via composite likelihood methods, in a similar manner to McVicker *et al.* (2009) and Elyashiv *et al.* (2016). Furthermore, one could include different parameters for neighboring polymorphic sites with specific functional annotations—such as those in genic regions, introns, and exons—to see how different classes of sites contribute to the additive genetic variance for fitness. Our hope is that statistical methods to quantify the effects of linked selection over short timescales will improve and be combined with measures of phenotypic change, leading to a more synthetic view of how selection on ecological timescales occurs at the genetic and phenotypic levels.

As the number of empirical temporal genomic studies continues to increase, it is worth mentioning how our study of temporal autocovariance suggests a few ways to optimize experimental design to increase the power to differentiate the effects of selection from drift. First, one should ideally sample frequencies from consecutive generations for at least some timepoints during the duration of the experiment. This is because the variance in allele frequency change between adjacent generations is only impacted by the heritable variance in offspring number, (see Equation 2), and not by the accumulation of temporal autocovariance terms (*e.g.*, Equation 3); this allows for more accurate estimates of *G*. In cases where a long study duration is needed but sequencing is limited to only a subset of generations, a mixed-duration sampling design, such as sampling generations 1 through 4, then 10 through 14, and so forth could serve as a compromise. Second, as described in Accounting for Allele Frequency Sampling Noise in the Appendix, the shared sampling noise between adjacent timepoints creates a negative bias in autocovariance that must be corrected for. As described in Equation 98 and Equation 99, we can estimate this bias from data, but this introduces additional uncertainty into our parameter estimates. In cases where the experimenter suspects *a priori* that fluctuating selection is occurring, *e.g.*, between two seasons, we recommend at least two temporal samples per season. This allows one to differentiate negative covariance occurring from the bias correction procedure underestimating bias from negative covariance caused by fluctuating selection through comparing nonadjacent timepoints that differ in season. Finally, one can directly remove the effects of the technical sampling noise created by variation in sequencing by dividing up temporal samples and barcoding them into two groups (*e.g.*, A and B). Then, the sample covariance estimate does not share the technical sampling noise, reducing the bias (but note some bias remains due to the sampling process where individuals are sampled from the population.)

### Connecting temporal linked selection with single-timepoint studies

Our goal in this paper is to suggest that quantifying variance and autocovariance using temporal data sets can help us understand the impact that linked selection has across the genome on short timescales, which supplements our current view informed mainly by single-timepoint studies. A range of approaches to estimate the parameters and impacts of models of linked selection from a single contemporary timepoint have been developed (Begun and Aquadro 1992; Wiehe and Stephan 1993; Hudson 1994; McVicker *et al.* 2009; Sella *et al.* 2009; Elyashiv *et al.* 2016). These estimates necessarily reflect linked selection over tens to hundreds of thousands of generations. One question is whether these estimates of the proportion of allele frequency change due to linked selection should line up with those over shorter time periods? Some forms of linked selection may be fairly uniform over time, whereas rare, strong sweeps will have a huge impact on long-term patterns of variation, but may be hard to catch in temporal data. Conversely, as we discuss below, fluctuating selection may lead to stronger signals of linked selection on short timescales than seen in long-term snapshots.

Studies of contemporary data have revealed multiple lines of evidence for the effect of linked selection in a variety of taxa. If linked selection is pervasive across the genome, diversity could be severely dampened as most sites would be in the vicinity of selected sites, thus reducing the genome-wide level of diversity without leaving strong local signals differentiated from the background. This is one proposed resolution of Lewontin’s paradox, the observation that diversity levels occupy a narrow range across taxa with population sizes that vary by orders of magnitude (Lewontin 1974; Smith and Haigh 1974; Gillespie 2001; Leffler *et al.* 2012). Elyashiv *et al.* (2016) estimated a 77– reduction in neutral diversity due to selection on linked sites in *Drosophila melanogaster*, and concluded that no genomic window was entirely free of the effect of selection. Similarly, Corbett-Detig *et al.* (2015) found evidence of a stronger relative reduction in polymorphisms due to linked selection in taxa with larger population sizes. However, these reductions fall short of the many orders of magnitude required for linked selection to explain Lewontin’s paradox (Coop 2016).

One limitation of these approaches is that they require estimating , the level of diversity in the absence of linked selection, usually from the diversity in high-recombination regions with low gene content. The average genome-wide reduction of diversity can then be judged relative to . Ideally, would be a measure of the average diversity due entirely to drift and demographic history, *i.e.*, unaffected by heritable fitness variation. However, there are two complications with this. First, as Robertson (2009) first showed, even a site completely unlinked from sites creating heritable fitness variation experiences a reduced effective population size due to the total additive genetic variance for fitness at these unlinked sites, and thus lower diversity [see also Santiago and Caballero (1995)]. The second complication is that if linked selection is sufficiently strong, the bases used to measure may not be sufficiently unlinked from fitness-determining sites to plateau to the Robertson (2009) level of diversity, a known potential limitation (Coop 2016; Elyashiv *et al.* 2016). Overall, the empirical studies relying on present-day samples from a single timepoint could be underestimating the effects that pervasive linked selection has on diversity. If linked selection can be observed over suitable timescales in temporal data, we might be able to disentangle some of these effects. For example, if high-recombination regions still show temporal autocovariance in allele frequency change, we would have evidence that even these regions are not free of the effect of linked selection and we might be able to estimate its long-term impact on levels of diversity.

Temporally or spatially fluctuating selection has long been discussed as an explanation for abundant, rapid phenotypic adaptation over short timescales, yet over longer timescales both phenotypic changes and molecular evolution between taxa are slow (Gingerich 1983; Hendry and Kinnison 1999; Messer *et al.* 2016). However, most of our approaches to population genomic data are built on simple models with constant selection pressures, as typically we have not had the data to move beyond these models (Messer *et al.* 2016). Currently, many approaches to quantify the impact of linked selection due to hitchhiking assume classic sweeps, where a consistent selection pressure ends in the fixation of a beneficial allele (Wiehe and Stephan 1993; Sella *et al.* 2009; Hernandez *et al.* 2011). However, fluctuating selection can have a larger effect on reducing diversity than classic sweeps (Barton 2000) depending on the timescales over which such fluctuations occur. In fact, as Barton (2000) points out, the total effect of classic Maynard–Smith and Haigh-type sweeps on diversity is limited by the relatively slow rate of substitutions. We show that when the direction of selection on a trait abruptly reverses, this creates negative autocovariance between the allele frequency changes before and after the reverse in direction. We can observe the shift by plotting autocovariances over time and noting when they become negative, indicating a negative additive genetic covariance between fitness at two timepoints. Here, we assume a simple form of fluctuating selection, where selection pressures on all of our sites flip at some timepoint. In reality, selection pressures will change on only some traits, and some of the genetic response will be constrained by pleiotropy, thus only some proportion of the additive genetic variance will change. Still, we expect some level of negative covariance after a reversal in the direction of selection, and there is an additional signal of fluctuating selection by comparing how the strength of temporal autocovariance varies with recombination and the initial level of LD in the genome.

### Connecting estimates of from temporal genomic data and quantitative genetic studies

The temporal covariance of allele frequencies potentially offers a way to estimate the additive genetic variance for fitness, as illustrated by our method-of-moments approach across genomic windows. The additive genetic variance for fitness can, like any other trait, be estimated through quantitative genetics methods, which exploit the phenotypic resemblance between relatives and their known kinship coefficients (Kruuk 2004; Shaw and Shaw 2013) [see Hendry *et al.* (2018) for a review], and these methods have been applied to estimate the additive genetic variance for fitness from natural populations (Mousseau and Roff 1987; Burt 1995). Ideally, one could reconcile quantitative genetic measures of fitness variance with estimates from allele frequency covariance. For example, Charlesworth (2015) undertook a similar analysis in *D. melanogaster*, comparing population genetic estimates of fitness variance to quantitative genetics estimates, highlighting a discordance potentially consistent with undetected large-effect alleles that are likely maintained by some form of balancing selection. By allowing us to directly measure fitness variation from population genetic data over very short timescales, temporal data could help untangle the causes of this discordance. A natural extension of this would be to see which regions contain the greatest inferred levels of additive genetic variance for fitness and test for functional covariates such as the number of coding bases, etc. Whereas previous temporal studies have focused on finding loci under selection, inferring the level of additive genetic variance could provide a more complete view of how much selection operates over short timescales.

### Conclusions

With temporal data, we can directly partition the total variance in allele frequency change across generations, , into components according to the underlying process governing their dynamics: drift and linked selection. Since the trajectory of a neutrally drifting polymorphism does not autocovary, evidence of temporal autocovariance across neutral sites in a closed population is consistent with linked selection perturbing these sites’ trajectories. If we consider drift to be the process by which nonheritable variation in reproductive success and Mendelian segregation cause allele frequencies to change, then this is estimable from and separable from the effects of linked selection using temporal data. This helps frame the long-running debate about the roles that neutral drift and linked selection have in allele frequency dynamics into a problem that can potentially be directly quantified by the contribution of each distinct process with temporal data.

## Acknowledgments

We thank Nick Barton, Doc Edge, Matt Osmond, Enrique Santiago, Michael Turelli, and two anonymous reviewers for feedback on previous versions of the manuscript, and Aneil Agrawal, Dave Begun, Sarah Friedman, Bill Hill, John Kelly, Tyler Kent, Chuck Langley, Sally Otto, Jonathan Pritchard, Kevin Thornton, Anita To, and members of the Coop laboratory for helpful conversations. This research was supported by an National Science Foundation (NSF) Graduate Research Fellowship grant awarded to V.B. (1650042), and National Institutes of Health (R01 GM-108779) and NSF (1353380) grants awarded to G.C.

## Appendix

### Decomposition of Allele Frequency Change

This decomposition of *neutral* allele frequency change between two consecutive generations is based on that of Santiago and Caballero (1995). We imagine a closed Wright–Fisher population of *N* diploids, where each diploid *i* contributes gametes to the next generation. We assume that the population size is constant, such that one diploid begets one diploid and thus . The neutral allele frequency in the next generation can be thought of as each of the *N* parents passing their average genotype (where is the number of tracked neutral alleles individual *i* carries) to their gametes, plus a random Mendelian deviation to each offspring *j*. Then, the frequency in the next generation can be written as(28)where , and is an indicator function that is one when the individual *i* is a heterozygote (*i.e.*, ), and zero otherwise.

If we further decompose the number of offspring of individual *i* into the genetic and nongenetic contributions, , then(29)and the change of the neutral allele’s frequency is the difference where . Then,(30)These *d*’s broadly capture nonheritable variation in an individual’s offspring number, with . In a quantitative genetics framework, can include nongenetic variation in the lifetime reproductive success of individuals , while in population genetic models can accommodate the sampling of parents to form the next generation, *e.g.*, multinomial sampling of individuals from fitnesses (Santiago and Caballero 1995). From now forwards, we assume that variation in these *d*’s is nonheritable. We note that in nonpanmictic populations, chance covariances could be created between a neutral polymorphism and environmental component of their phenotype, especially as a population expands its range into new environments that affect the phenotype and variants “surf” to higher frequencies (Edmonds *et al.* 2004; Excoffier and Ray 2008; Hallatschek and Nelson 2008).

Note that by construction, the allele frequency change components and are orthogonal within each individual, given the neutral allele frequency *x*. Partitioning the allele frequency change for an individual *i* into its components,(31)(32)where .

Two random variables are uncorrelated if , and are orthogonal if either has an expected value of zero such that . We show briefly that taking expectations *over conceptual evolutionary replicates*, the terms are orthogonal. First, the terms and are orthogonal (dropping *i* subscripts),(33)since across evolutionary replicates due to the assumption that population size is constant, and , as across all evolutionary replicates there is no dependence between a particular neutral allele an individual carries and their fitness (though in particular replicates, such associations occur). Similarly, for the case (other cases are all zero and can be ignored), it can be shown using the law of total expectation that [and likewise with and ]. Note that across individuals within a population, there are weak covariances in their number of offspring as the total number of offspring must sum to *N*; under a multinomial offspring distribution, these are of order .

### Temporal Variance and Autocovariance Under Multilocus Selection

We assume the phenotype of an individual *i* has an additive polygenic basis, such that their breeding value is , which deviates around a mean of zero, and is the additive effect size at locus *l* in generation *t* and is individual *i*’s allele count at this locus (note that the effect of nonheritable environmental noise affecting the trait is accounted for in the terms above). We impose directional selection on this trait using an exponential fitness function, such that individual *i*’s fitness is (assuming ). If we assume individuals’ phenotypic values do not deviate too far from their mean value of zero, we can approximate as: . Then, we can write the change in neutral allele frequency due to only heritable variation in fitness as a covariance between fitness and the neutral allele frequency across individuals in generation *t*,(34)(35)(36)which is the is the Robertson–Price covariance (Robertson 1966; Price 1970; Lynch *et al.* 1998; Walsh and Lynch 2018).

Now, we break up the genotypic value into the contributions of each of the two gametes that formed individual *i*, , and likewise with the trait locus , where and are all indicator variables. Expanding out the covariances, we have(37)Each of these covariances is between the neutral allele and a selected allele, either on the same gamete (either maternal or paternal) or across gametes. These covariances can be written as LD terms,(38)where is the LD between alleles on the same gamete (the gametic LD), and the is the across-gamete LD (nongametic LD) [see p. 121 in Weir (1996)]. This equation also appears in Kirkpatrick *et al.* (2002) (equation 10).

We ignore nongametic LD as these are weak under random mating, and write the multilocus temporal covariance between the allele frequency changes and as(39)since .

#### Modeling the dynamics of LD between selected and neutral sites

Here, we outline a model of the changes in LD between the focal neutral site and selected sites, which allows us to derive an expression for the first term in Equation 39. Typically, models of multilocus selection track the genetic changes in a population by transforming between a representation of haplotype frequencies to a representation of allele frequencies, LD, and higher-order LD (Barton and Turelli 1987, 1991; Turelli 1988; Turelli and Barton 1990). We for the moment avoid the difficulty of a full multilocus treatment by assuming that the linkage between selected sites is loose enough that one selected site’s frequency change is independent of the change at other selected sites, *e.g.*, there is no selective interference; this was an assumption in past treatments [Santiago and Caballero (1995, 1998); see Barton (2000) for discussion of this]. Specifically, we model the dynamics of the LD between the neutral site and each selected site as they would behave under a single-sweep model.

We adapt Barton’s (2000) model for LD dynamics during a sweep. We imagine a polymorphic neutral locus has alleles and with frequencies *p* and . We partition the allele frequency of by conditioning on which allele at the selected site (either or ) is carried on the same background, *e.g.*, and , such that . Then, the LD between neutral and selected sites can be expressed as . To simplify notation, we denote the frequency of the neutral allele on the different fitness backgrounds at time *t* by and , and the selected allele at locus *l* frequency as , then the LD between the focal neutral site and selected site *l* at time *t* is(40)Selection changes the frequency through time and recombination acts to disassociate the neutral allele with its backgrounds. The expected difference is maintained with probability each generation (Barton 2000). If the initial generation is *t* and the future generation is *s* , and is the recombination rate between the focal neutral locus and selected site *l*, this leads to(41)Then, we can use this to describe the dynamics of to generation *t*,(42)[compare with equations 30 and 31 in Stephan *et al.* (2006)]. Now, we can find the expected product by multiplying both sides by and taking expectations. We treat the allele frequency trajectory as deterministic, giving us,(43)Then, we simplify this by replacing with [where is the square of the correlation between the neutral site and selected site *l* at time *t* Hill and Robertson (1968)],(44)Returning to Equation 39 and replacing the terms with our expression above,(45)again ignoring the cross-associations between selected sites. Dividing our temporal covariance by the neutral site’s , we can write the multilocus temporal covariance in a standardized form (analogous to Wright’s F),

#### Using average additive genetic variation

We can approximate Equation 46 by noticing that the terms are similar to an additive genic variation if effect sizes remain constant through time. We make that assumption here, writing , leading to(47)We can further simplify this by assuming that there is no covariance between the additive genic variation at a selected site, and the LD between that selected site and the neutral site. We write the additive genic variation at site *l* at time *s* as , and the average additive genic variation across loci as . Then, each locus’s additive genic variation can be expressed as: . Substituting this, the autocovariance is(48)We assume that this last term, which is nonzero in expectation only if there is covariance between the additive genic variation at a selected site, and the expected LD between the selected and the neutral sites is zero. Rewriting the total genic variation as ,

#### Continuous approximation to chromosomes

Currently, we have treated the positions of the selected sites as fixed, conditional on knowing that a selected site *l* has a recombination fraction away from the focal neutral site. Here, we make a few further assumptions. First, we assume the *L* selected loci are each independently and identically uniformly distributed along a *continuous* region of *R* M in length , and we now calculate the covariance at a focal neutral site at the origin by taking expectations over the random positions of these sites. Second, we assume that the LD between the neutral and selected site *l* only depends on the recombination fraction between the sites. Since *g* is the genetic distance, the recombination fraction is now provided by a mapping function that maps genetic distances to recombination fractions. Throughout the paper and in the simulations, we use Haldane’s mapping function, (note the absolute value translates positions on to distances to the focal neutral site). Next, we assume that the LD between two sites can be completely determined by the distance between the focal neutral site and a random selected site, allowing us to rewrite as the function . Now, letting represent the expectation taken over the random positions of selected sites on the genetic map,(50a)(50b)(50c)since the cancels with the *L* from the sum of expectations.

As our trait is neutrally evolving before directional selection starts, we use the expected neutral LD, , where (Hill and Robertson 1968; Ohta and Kimura 1969) when *t* is the first generation when selection begins.

### The Contribution of the Rest of the Genome to Temporal Autocovariance at a Locus

Additionally, we can consider the impact that other unlinked selected sites have on the neutral allele’s frequency trajectory. In Equation (50c), we are modeling the temporal autocovariance at a focal neutral site in a chromosome with total length *R*. Here, we assume the whole genome can be modeled in a similar fashion, as a single very large chromosome with genetic length *C*, the entire map length of the genome, and the focal neutral falls in the center of this. Then, we can take a piecewise integral over all linked sites (those less than a recombination fraction of one-half away) and all unlinked sites with respect to the genetic distance (those at a recombination fraction of one-half away). Our piecewise integral gives us(51)(52)where is the map distance at which a selected site becomes approximately unlinked to the neutral site, *e.g.*, the recombination fraction is . As we move away from the neutral site, the first term in the brackets accounts for the accumulation of linked selected sites. Eventually, the genetic distance away from the neutral site becomes unlinked and the second term accounts for the contribution of these unlinked sites on both the same chromosome, as well as sites on other chromosomes. Note that here we assume that the density of additive genetic variance per morgan is constant. As in the main text, this ignores the contribution of nongametic LD, formed as two gametes unite in an individual, and can be converted into gametic LD via recombination [see Figure A2; this process is elaborated in *The Strength of Unlinked and Nongametic Associations* in the Appendix, and see also p. 521 in Tenesa *et al.* (2007)]. Selected sites that are unlinked or loosely linked to the neutral site (*e.g.*, ) quickly become associated with the neutral site, but also quickly decay; their contribution acts to further decrease the population size by a factor of two [see the discussion on p. 2115 in Santiago and Caballero (1998)].

As a simple thought experiment, we might ask for what value of *M* does the contribution from unlinked sites dominate the contribution from linked sites? For and assuming that the level of LD is that under mutation–drift–recombination balance (*e.g.*, using the equation of Tomoko Ohta), we plot the relative contributions of linked selected sites (on the focal chromosome) and unlinked selected sites (on other chromosomes) for various spans of the covariance (*e.g.*, ), and the size of the remaining portion of the genome in morgans (*M*) in Figure A1.

### Averaging Covariance Across Multiple Loci

Thus far, our covariance assumes that a single neutral site is positioned in the center of a region *R*-M long, with selected sites uniformly distributed along this region. However, in our simulations we simulate a region that contains many neutral sites, which we average over in calculating the temporal autocovariance. In this case, we average over the random distance between a neutral site’s position *n* and a selected site’s position *g*, which is , where . This random variable *c* is distributed according to the triangle distribution, ; we replace the uniform probability density function (PDF) in Equation (50b) with the triangle density PDF and average over the distance between sites,(53)(54)where indicates we take the expectation also over neutral sites, and we use to denote the average LD between selected and neutral sites that persists from generations *t* to *s* . Note that in calculating the standardized covariance above, we use a ratio of expectations rather than the expectation of the ratio (Bhatia *et al.* 2013).

### Empirically Calculating the Average LD Persisting Across Generations

In the previous expressions for temporal autocovariance, we stepped through a conceptual model for the average levels of LD between neutral and selected sites that persists across generations , where the positions of selected and neutral sites are randomly distributed along a chromosomal region. In systems with a known recombination map and studies where LD can be calculated, we have the recombination fraction , and the pairwise LD between two loci *i* and *j* (where is the matrix of pairwise LD calculated at time *t*). Since we do not *a priori* know whether a site is selected or not, we sum over all polymorphic *M* loci, thus characterizing the average LD in a region as(55)This sum is the empirical analog to the integral in Equation 10.

### The Strength of Unlinked and Nongametic Associations

Here, we characterize the contribution of completely genetically unlinked loci segregating for fitness variation to the change in frequency of our neutral allele. Across evolutionary replicates, there is no expected covariance between the neutral allele an individual carries and their fitness ; rather, for unlinked loci, chance associations are created from the variance around this sampling process of neutral alleles into individuals with varying fitness . As the neutral allele and fitness variation independently assort themselves into individuals, the chance associations that form have a variance given by . This has the form of the sampling variance of a covariance, which for random variables *X* and *Y* is given on p. 472 Kendall *et al.* (1994)),(56)(57)where and are the means of *X* and *Y*, respectively, the variance is taken over conceptual replicate populations, and the covariance is calculated over the individuals in a population. Then, applying this to our covariance ,(58)Thus the chance covariances that form between the neutral alleles individuals carry and their fitness have a variance proportional to .

### Nongametic LD’s Contribution to Temporal Autocovariance

Throughout the paper, we ignore the effects of nongametic LD, , the disequilibria that occurs between the two gametes (maternal and paternal) at two loci (see Figure A2A for an illustration of gametic LD and nongametic LD ). Following the equation for the sampling variance of nongametic LD in Weir (1996) (see p. 124), the chance nongametic disequilibrium that builds up sampling gametes into *N* individuals is,(59)Note the similar form to Equation 58 as both the unlinked and nongametic LD arise from the random sampling of alleles at different loci into individuals.

There is no expected covariance between our gametic and nongametic LD within a generation , assuming random mating. However, for , as a fraction of the nongametic LD may be converted into gametic LD in the next generation. Specifically, following Santiago and Caballero (1995), we can write the product of the nongametic LD in generation *t* with the gametic LD in generation *s* as(60)where a proportion *r* the nongametic LD in generation *t* is converted into gametic LD, and a proportion of this is carried forward unbroken by recombination over the remaining generations (see Figure A2B for an illustration of this process).

In our analysis in the main text, we ignore these terms as is expected to be small due to its inverse dependence on *N*. However, these terms are necessary for the analysis of looser linkage (Santiago and Caballero 1995, 1998).

### Connecting our model with the models of Robertson and Santiago and Caballero

Here, we describe the models of Santiago and Caballero (1995, 1998), relating their work on the long-run effective population size experienced by a neutral allele where there is (1) unlinked heritable fitness variation (Santiago and Caballero 1995) or (2) linkage, where fitness-determining sites are randomly scattered along a chromosome (Santiago and Caballero 1998). Overall, their models are formulated in a quantitative genetics tradition, where the population genetic dynamics at the selected loci are not explicitly modeled (although these links are made more explicitly in their 1998 paper). In contrast, in deriving our expressions for temporal autocovariance and variance, we use a population genetic approach, modeling the dynamics at selected sites (though we simplify from the full multilocus treatment, *e.g.*, we assume selected loci experience independent sweeps and we ignore the LD between selected sites). We show that we can reconcile the two approaches, and demonstrate that the temporal autocovariance expressions we develop in our models are implicit in their model. We also work through their expressions for with heritable variance, because it represents a quite useful result but their original presentation was spread across two papers (and a change in notation).

#### Santiago and Caballero’s 1995 and 1998 models for

While our goal in the main text of our paper was to develop expressions for the temporal variances and autocovariances in allele frequency change when there is heritable fitness variation in the population, the goal of both the 1995 and 1998 Santiago and Caballero papers was to derive an expression for the long-run when there is heritable variation for fitness in the population. In their 1995 paper, they found the effective population size for large *t* (see p. 1018, equation 16) to be(61)(62)where is the heritable variation for fitness ( in our notation), and is the nonheritable variation in offspring number (*i.e.*, under a Wright–Fisher model, ). For a neutral locus completely unlinked from fitness variation [the situation first considered by Robertson (1961)], [see equation 17 in Santiago and Caballero (1995)]. Here, *G* represents the decay rate of the additive genetic variance associated with a particular haplotype. Note that we have simplified their expressions by assuming no assortative mating, and that we try to follow their notation as closely as possible (consequently, the *Q* here is unrelated to the of the main text). Santiago and Caballero (1995) assume that continual artificial selection maintains a constant level of fitness variation in the population each generation, yet the *particular* fitness backgrounds the neutral allele is stochastically associated with only contribute a fraction of *G* in the next generation, in the generation after, and so on, as selection reduces genetic variation for fitness (note: in their 1998 paper, they use *Z* instead of *G*). Similarly, the associations between the neutral and fitness backgrounds decay at a rate due to independent assortment. Note that Robertson (1961) assumed that the fitness backgrounds that become stochastically associated with the neutral allele do not experience any decay in their fitness variation ; in this case, as Robertson’s work found. With an arbitrary amount of linkage between the focal neutral and fitness backgrounds, Santiago and Caballero (1998) show that only *Q* is affected, and derive an expression for *Q* that only depends on *G* and the size of the genome in morgans, *L* [see equation 6 in Santiago and Caballero (1998)].

In our main text, we model the temporal autocovariance created by heritable fitness variation, which also impacts the cumulative variance in allele frequency change . To illustrate how our model connects with theirs, below is the cumulative variance in allele frequency change for three generations in their 1995 notation, with the corresponding changes in allele frequency below:(63)Grouping terms by the generation that the initial association was formed, we see how Santiago and Caballero (1995) define the terms in their notation,(64)since the associations created in generation *i* persist with probability , with proportion of its original fitness variation in generation *t*. In general, the cumulative impact of the associations formed *i* generations ago has coefficient . Using these terms simplifies this equation to(65)or, in general,(66)Then, Santiago and Caballero (1995) note that assuming , , and population size *N* are constant across generations, the magnitudes of all of the effects , and are constant across all generations (for all *i*, so we omit the *i* subscript for these terms), except for a geometric decay due to drift at a rate per generation that effects all terms. Such that, when we include the decay in the variance due to drift,(67)[compare with p. 1018 of Santiago and Caballero (1995)]. In the long run, the variance in the neutral allele’s frequency change hits a balance. Many copies of the neutral allele segregating in the population are on fitness backgrounds that it has recently become stochastically associated with, as segregation and recombination have not broken these associations apart. A few copies of the neutral allele are on fitness backgrounds they became associated with many generations ago, that have by chance survived to remain associated. In all cases, the effect that these associations have on present-day allele frequency change is weakened by the fact that natural selection has reduced the genetic variance of these fitness backgrounds. Since the long-run variance in allele frequency under drift in a Wright–Fisher population is(68)one can estimate the effective population size using the observed difference in variances and . Note that this is a different long-run effective population size to that used by others, (Crow and Kimura 1970). Santiago and Caballero use Equation 68, taking the difference and rearranging to end up with the large *t* estimator of ,(69)[compare with p. 1018 in Santiago and Caballero (1995)]. Rearranging,(70)This very conveniently simplifies the sum in Equation 67, as we can show with the case of ,or, generally,(71)Inserting this into Equation 70, the long-run effective population size can be written as(72)where in Robertson’s (1961) model, and in Santiago and Caballero’s (1995) model. Note that *r* here represents the recombination fraction between fitness variation and neutral sites, which differs from equation 17 in Santiago and Caballero (1995), where *r* represents the correlation between parental fitness).

Then, Santiago and Caballero (1995) show,(73)(74)(75)[compare with equation 11 in Santiago and Caballero (1995)]. Inserting these into Equation 72, we have(76)which is a simplified version of equation 16 in Santiago and Caballero (1995), and which further simplifies to when and , the effective population size of a Wright–Fisher population of *N* hermaphroditic individuals.

#### The covariances caused by fitness associations

With an understanding of the basics of Santiago and Caballero’s (1995, 1998) models, how they connect to our notation, and how they reach their expression for the long-run effective population size, we turn now to finding the temporal autocovariances implicit in their model. We start by looking at the variance in allele frequency between generations 0 and 4 (Equation 63), including an additional generation so the pattern is clearer later,The cross terms like , and are all expected products of independent random variables, where the expectation of each random variable is zero, and consequently are all zero. The only nonzero cross terms are products of . When we look at the covariances with the allele frequency change in the initial generation and a later generation *s*, ,(77)(78)(79)(80)where the coefficient comes from the fact that, in Santiago and Caballero’s work, the products represent *both* the and terms, so a single temporal autocovariance in our notation is one-half their joint covariance term.

However, if our reference generation is different, say two, the associations from earlier generations that have persisted to that generation can also lead to covariances to later generations. Looking at the covariances(81)(82)(83)Likewise, the covariances include the associations that persist from earlier generations. In general,(84)This expression is more complex than our expression for temporal autocovariance because it is modeling the LD in generation *t* as it builds up from generations 1 to *t*. In contrast, our expressions for covariance incorporate all of this buildup of LD as the initial LD term (for a single-locus case). This expression for autocovariance implied by Santiago and Caballero’s (1995) work matches our expression for autocovariance (for a single locus) when the arbitrary first generation is *t* rather than one(85)Using the expression for [equivalent to in our notation] derived in *The Strength of Unlinked and Nongametic Association* in the Appendix,(86)This is analogous to Equation 8 for a single locus, where is the additive variation in generation *s* [equivalent to our ], and the factor represents the chance buildup of LD between the neutral site and an unlinked fitness background. In our expression, we condition on existing LD between the neutral site and its fitness background, whereas they assume a buildup of LD to a drift–recombination equilibrium. We can see this by returning to the term of Equation 63,(87)(88)where following Santiago and Caballero’s (1995) approach (see p.1018), we can replace each with and let as we focus on the buildup of LD. This gives us the general equation,(89)and taking this geometric series to infinity converges [since ) and replacing ] with the chance associations that build up gametes sampled into individuals (Equation 58),(90)When we assume , , and is a constant, this gives us(91)which is analogous to the measure of LD, standardized to rescale the fitness variation. The right-hand side is identical to the identity-by-descent equilibrium under drift–recombination balance Sved (1971). Note that our expression can be recovered from Equation 84 when the reference generation such that the LD hits its equilibrium level,(92)which is identical to our expression for temporal autocovariance when initial LD is due to a neutral drift–recombination balance, and the change in during selection is modeled as a geometric decay at rate *G*. Note that the terms in underbraces indicate the corresponding terms in Equation 8 for a single locus.

### Multilocus Simulation Details

#### Targeting an initial level of additive genetic variation

We choose *θ* for the coalescent simulations to target a total number of segregating sites , where *L* is the number of selected sites (a parameter we vary in our multilocus simulations) and or more randomly placed neutral sites over which we can calculate the temporal autocovariance. Then, the total number of target sites is then inflated by a factor of 1.5 to ensure a sufficient number of sites given the random mutation process. The *L* selected sites are randomly chosen from the segregating sites, and all remaining mutations are neutral. Thus, using Watterson’s expression for the expected number of segregating sites under the coalescent (Watterson 1975), we have where is Euler’s γ. Each of the *L* selected sites is given a random effect size of with equal probability, where we choose by targeting a specific level additive genic variation where is the sum of site heterozygosity of the *L* neutrally evolving sites that will be selected once selection begins. Under neutrality, . Then, we set . We empirically validate that our target genic variation is close to the empirically observed level.

#### Choosing the simulation parameter range

We simulate over a grid of parameter ranges, varying the number of loci *L*, the target additive genic variation in the region (by varying *α*) and the recombination map length of the region *R* in morgans. We have varied these parameters over ranges that encompass a wide range likely to be encountered in natural populations. To do this, we found that additive genetic variation for lifetime reproductive success varied over orders of magnitude, from 0 [*e.g.*, in male red deer (Kruuk *et al.* 2000)] to 1.1 [in the male red-billed gull (Teplitsky *et al.* 2009)]. These values of additive genetic variation are for the entire genome (we write these as , where indicates genome-wide); our simulations model a region of varying map length. We expect most recombination map lengths to be roughly over the scale of M in length, and we chose to investigate how temporal autocovariance behaves across a spectrum of recombination, from a completed linked region , the scale over which a strong classic hard sweep could affect diversity (0.5 cM), to a region where the ends are approximately unlinked (4.5 M), overall giving us a parameter range of . Over this grid of our region recombination map lengths, and the total organism recombination map lengths, we get a rough estimate of the number of regions we would expect with this level of recombination in the organism’s total genome, imagining a homogeneous recombination rate, by taking . Using this estimate of the number of regions in the genome, we calculate the genetic variation per region over our grid of parameters as and target a level of additive genic variation per region equal to this regional additive genetic variation . From preliminary simulations, we found that we cannot detect much temporal autocovariance below with the initial level of LD from mutation-drift balance, so we ignore parameters less than this value (other than as a control). Additionally, to reduce the number of simulations, we exclude as this only excludes a small region of the parameter grid and preliminary simulations demonstrated that the behavior of temporal autocovariance with high is evident with the included values. Overall, this gives us a spectrum of target additive genic variation per region of . Note that to prevent our plots from being too dense, we include only a representative subset of this parameter grid in our figures.

### Accounting for Allele Frequency Sampling Noise

In practice, one will calculate the temporal variance–covariance matrix on allele frequency trajectories calculated from sampled chromosomes from the population. We assume a binomial sampling process, where *n* chromosomes are sampled from the population such that and . We can then write , and our covariances can be written as(93)Note this simplifies to(94)since in these cases, the sampling noise at a timepoint is not shared between the estimated allele frequency changes. However, if the sampling noise from timepoint is shared, biasing the sample estimate of covariance:(95)Similarly, the variance is biased, as it is impacted by the binomial sampling noise too,(96)In practice, these covariances are calculated over loci in a region or across the entire genome. We assume that the tracked allele has been randomly swapped (*e.g.*, the tracked allele frequency is not systematically the minor, major, or reference allele), such that for all *t* and *l*. Then, the unbiased covariance estimate as calculated over loci is(97)Then, we can use an unbiased plugin estimate of the frequency sampling variance [see p.191 in Nei (1987)] to estimate these bias terms, and add or subtract them from the estimator accordingly. Accounting for finite sampling, the unbiased sample variance–covariance matrix now has elements:(98)and variance(99)In *Comparing theory to simulation results*, we used population frequencies in introducing the method of moments estimators of and . Here, we discuss the performance of these estimators with sample allele frequencies. Our simulations are identical to those described in *Multilocus simulation details*, except we have increased the target number of neutral sites in each region so it is around . We mimic sampling of chromosomes from the population, and use the bias-corrected sample variance–covariance matrix in the method-of-moments approach described in *Estimating Linked-Selection Parameters from Temporal Autocovariance* to estimate and .

In Figure A3, we show the performance of our estimators in the case where chromosomes have been sampled from the population. Overall, there are two important differences compared with Figure 5 of the main text. First, while the estimator performs well for high levels around , the variance around the estimator increases significantly as grows weaker. As the covariances are proportional to , sampling noise grows larger than the theoretical temporal autocovariance as becomes weaker. Then, one cannot discriminate against the chance covariances formed by the sampling process from the temporal autocovariance created by linked selection without either a large sample size or more timepoints. Second, the approach underestimates for very loose linkage . This is another consequence of the first problem; as sampling noise grows equal to or larger than the magnitude of temporal autocovariance, the estimation procedure performs poorly. As the linkage becomes looser, the magnitude of temporal autocovariance grows weak relative to the sampling noise, and this noise can be partially absorbed by . This effect can be somewhat ameliorated by calculating the sample variance–covariance matrix over a shorter region of the genome such that *R* is smaller (as long as SNP density is sufficient) or by increasing the sample size. Finally, the estimation of effective population size shown in Figure A3 is also affected by , as discussed in *Estimating Linked-Selection Parameters from Temporal Autocovariance* (though the effect is obscured by the sampling noise): high levels of in regions with low recombination generally lead to underestimates of .

To understand how sample size affects the method-of-moments estimators, Figure A4 depicts median relative error of and for various sample sizes.

## Footnotes

Supplemental material available at FigShare: https://doi.org/10.6084/m9.figshare.7709930.

*Communicating editor: N. Barton*

- Received July 31, 2019.
- Accepted September 22, 2019.

- Copyright © 2019 by the Genetics Society of America