## Abstract

Genome-wide association studies (GWAS) have revealed that many traits are highly polygenic, in that their within-population variance is governed, in part, by small-effect variants at many genetic loci. Standard population-genetic methods for inferring evolutionary history are ill-suited for polygenic traits: when there are many variants of small effect, signatures of natural selection are spread across the genome and are subtle at any one locus. In the last several years, various methods have emerged for detecting the action of natural selection on polygenic scores, sums of genotypes weighted by GWAS effect sizes. However, most existing methods do not reveal the timing or strength of selection. Here, we present a set of methods for estimating the historical time course of a population-mean polygenic score using local coalescent trees at GWAS loci. These time courses are estimated by using coalescent theory to relate the branch lengths of trees to allele-frequency change. The resulting time course can be tested for evidence of natural selection. We present theory and simulations supporting our procedures, as well as estimated time courses of polygenic scores for human height. Because of its grounding in coalescent theory, the framework presented here can be extended to a variety of demographic scenarios, and its usefulness will increase as both GWAS and ancestral-recombination-graph inference continue to progress.

SOME of the most compelling examples of phenotypic evolution come from time courses that reveal the pace of evolution, either through observations across generations (Cook *et al.* 1986; Grant and Grant 2002) or through changes in the fossil record (Gingerich 1983; MacFadden 2005; Bell *et al.* 2006). For many traits and species, it can be difficult to ascertain whether these changes reflect genetic change. For example, we have fairly detailed knowledge of human height through time, but some changes in height are likely driven by environmental and dietary changes (Stulp and Barrett 2016). Thanks to ancient DNA, we can now sometimes obtain a partial picture of long-term genetic changes involving relatively simple traits like pigmentation (Ludwig *et al.* 2009) or more complex traits (Mathieson *et al.* 2015). However, we are usually not fortunate enough to have access to genotype data from across time, and even when ancient DNA are available, the resulting time courses will necessarily be incomplete.

One alternative to direct measurement of phenotypes through time is to reconstruct the history of a phenotype using contemporary genetic data. Positive selection on simple genetic traits drives large allele-frequency changes at the causal loci and linked neutral alleles (Smith and Haigh 1974). There are many procedures for detecting this kind of selection on individual alleles and for dating and modeling their spread through populations (Tajima 1989; Fay and Wu 2000; Sabeti *et al.* 2002; Voight *et al.* 2006; Ronen *et al.* 2013; Garud *et al.* 2015; Crawford *et al.* 2017).

One obstacle to understanding the evolutionary basis of phenotypes is the polygenic architecture of many traits. Complex traits (traits affected by many genetic loci and by environmental variation) are ill-suited to study by single-locus methods. In recent years, genome-wide association studies (GWAS) have made it possible to aggregate subtle evolutionary signals that are distributed across the many genetic loci that are associated with a trait of interest (Turchin *et al.* 2012; Berg and Coop 2014; Robinson *et al.* 2015; Field *et al.* 2016; Racimo *et al.* 2018; Uricchio *et al.* 2018). For example, Field *et al.* (2016) developed the singleton density score (SDS) to infer recent selection on a variety of traits among the ancestors of people in the United Kingdom. [As discussed below, some empirical findings of these studies have not replicated using effect-size estimates from less-structured GWAS samples, raising the possibility that the selection tests are sensitive to residual population stratification (Berg *et al.* 2018; Sohail *et al.* 2018). Nonetheless, given correct effect-size estimates, these methods are useful.]

Field *et al.* relied on the fact that selection distorts the gene genealogies, or coalescent trees, at genetic loci under selection. In particular, loci under positive selection will have increased in frequency in the recent past, leading to relatively faster coalescence of lineages than if the allele frequency had been constant [the method of Palamara *et al.* (2018) also capitalizes on this idea]. The principle on which the SDS relies is quite general: selection, even when its effect is spread over many loci, leaves systematic signals in coalescent trees at loci underlying trait variation.

The ancestral recombination graph (ARG) (Griffiths and Marjoram 1997) collects coalescent trees at loci along a recombining sequence, encoding information about allele-frequency changes at each site, as well as recombination events between sites. Recently, computational methods for inferring ARGs have advanced considerably (Rasmussen *et al.* 2014; Mirzaei and Wu 2016), allowing a range of applications (Palacios *et al.* 2015).

In this work, we consider ways in which ARGs—and in particular, the coalescent trees of sites associated with a phenotype—might be used to reconstruct the history of the phenotype with which they are associated. The ARG-based approaches we consider are motivated by polygenic traits, and the population-mean level of a polygenic score (a prediction of phenotype from an individual’s genotype) is the target of estimation. We present methods for estimating the time course of a population-mean polygenic score through the past, as well as a test for assessing whether an estimated time course is consistent with neutral evolution alone.

We begin by describing estimators and a hypothesis test for phenotypic time courses based on previous theory. Next, we apply these procedures to simulated data, using both true and reconstructed ARGs. Finally, we apply our methods to some human heights in the GBR (Great Britain) subset of the 1000 Genomes panel (1000 Genomes Project Consortium *et al.* 2012), using ARGs inferred by RENT+ (Mirzaei and Wu 2016).

## Theory

### Background and motivation

The ARG expresses the shared genealogical history of a sample of individuals at a set of genetic loci, accounting for correlations among neighboring loci that arise because of linkage. The ARG contains a coalescent tree for every site in the genome—these trees for specific sites are *marginal trees*. Our methods make use of the marginal trees at a set of sites that are associated with a phenotype. In particular, we concentrate on the information revealed by the number of coalescent lineages that remain (*i.e.*, that have not yet coalesced) at a time *t* in the past. The number of lineages at a given time in the past is described by a stochastic process called the ancestral process (*e.g.*, Tavaré 1984).

The intuition behind the methods we present here is shown in Figure 1. If an allele has been selected upward in frequency in the recent past, then the number of chromosomes carrying the allele will likely have increased. Looking backward in time, the number of carriers in the recent past is less than in the present, which forces an excess of recent coalescence events. Similarly, if an allele has been selected downward recently, then there will tend to be fewer recent coalescent events compared with the neutral expectation. If the trait that has been under selection is polygenic, then the signal at each locus associated with the trait will be smaller, and its strength and direction will depend on the effect size at the locus. In Appendix A, we derive the relationship between the rate of coalescence and selection on the phenotype. We show that phenotypic selection acting to increase the population-mean trait value increases the rate of coalescence for alleles that increase the trait and lowers the coalescence rate for alleles that decrease the trait. In contrast, stabilizing selection acting on a trait for which the population mean is at the fitness optimum does not have a systematic directional effect on the coalescent rates.

Our target of estimation is the population-average polygenic score for a trait going backward through time. By polygenic score, we mean a weighted sum of an individual’s genotypes, where the weights are the additive effect sizes of each allele. In our case, we are interested in the population-average polygenic score, so we take a weighted sum of the allele frequencies:(1)Here, _{} is the population-average polygenic score at time *t* in the past; is the population frequency of one of the two alleles at locus *i* at time *t* in the past, and is the additive effect size of the allele whose frequency is , where the effect sizes have been scaled so that the other allele has an effect size of zero. (In practice, effect sizes will be estimated with error, but in this article we treat the effect sizes as known.) The 2 arises because of diploidy.

If the *k* loci included in the calculation of include all the causal loci, then changes in (the population-average polygenic score) reflect changes in the population-average phenotypic value in the absence of changes in the distribution of environmental effects on the trait, changes in the effect size, epistasis, and gene-by-environment interaction. Even if these strong assumptions are not met, rapid changes in could provide evidence that natural selection has acted on the trait.

Our strategy for estimating is to estimate the historical allele-frequency time courses, for the loci associated with a trait. Given an estimator of the allele-frequency time courses, we estimate polygenic scores as(2)If the loci contributing to the polygenic score are independent, then the variance of the polygenic-score estimator is(3)We present three methods for estimating historical allele-frequency time courses. A number of authors have investigated estimating allele-frequency time courses from coalescent genealogies (Slatkin 2001; Coop and Griffiths 2004; Chen and Slatkin 2013) or by applying Wright–Fisher diffusion theory to time-series data (Bollback *et al.* 2008; Schraiber *et al.* 2016). Our approaches are cruder than some of these but have the advantage of being fast enough to be applicable to thousands of GWAS loci.

### Estimating the allele-frequency time course at a single locus

We present several methods for estimating the historical allele-frequency time course at a specific biallelic locus given a coalescent tree at the locus (Figure 2). (Our procedures could be generalized to loci with multiple alleles.) In each case, the goal is to estimate the frequency of an allele of interest (*e.g.*, the effect allele) at locus *i* at time *t* in the past, or

#### Proportion-of-lineages estimator:

The simplest way to estimate a historical allele frequency is to treat the lineages ancestral to the sample at time *t* as representative of the population at time *t*. If the locus has evolved neutrally between the present and time *t* in the past, then the lineages ancestral to the sample at time *t* are a random sample—with respect to allelic type—from the population at time *t*. If the lineages ancestral to the sample are a random sample from the population at time *t*, then a reasonable estimator of is the proportion of lineages at time *t* that carry the allele of interest,(4)where is the total number of lineages ancestral to the present-day sample at locus *i* at time *t* and is the number of lineages at time *t* that carry the allele of interest. Assuming that the mutation distinguishing the alleles has only appeared once in the history of the sample, the lineages that carry the derived allele are those ancestral to contemporary copies of the derived allele, which must coalesce to one lineage before coalescing with the rest of the tree. If the tree for the locus is known, then the branch on which the mutation must have appeared is known, but the exact time of the mutation is not. In practice, we assume that the mutation occurred in the middle of the branch connecting the derived subtree to the rest of the tree. (We make this assumption when implementing all estimators.)

If the population size at time *t* is large compared with the number of ancestral lineages at time *t*, then conditional on the number of ancestral lineages carrying the allele of interest, is distributed approximately as a random variable. Thus, Equation 4 is the maximum-likelihood estimator of and, conditional on its sampling variance can be estimated as (dropping the subscript *i*’s and parenthetical *t*’s for compactness)(5)In Appendix B, we give a Bayesian interpretation of the proportion-of-lineages estimator that relies on connections between neutral diffusion and the ancestral process (Tavaré 1984).

If chromosomes carrying the two alleles have differed in fitness between the present and time *t*, then Equation 4 will in general be a biased and inconsistent estimator of _{}. Chromosomes carrying alleles that have been favored by selection will be more likely to leave offspring in the present-day population than will chromosomes carrying unfavored alleles. Favored alleles will thus be overrepresented among the lineages ancestral to the sample compared with their actual frequency at time *t*. However, even if Equation 4 is an inconsistent estimator of the population’s allele frequency, it retains the interpretation of being the allele frequency among lineages ancestral to the sample. Thus, when Equation 2 is applied to allele frequencies estimated by Equation 4, the result is the mean polygenic score among chromosomes ancestral to the present-day sample at some time in the past.

#### Waiting-time estimator:

The proportion-of-lineages estimator proposed in Equation 4 works well under neutrality, but under selection it tends to underestimate the degree of allele-frequency change experienced by selected alleles. One potential solution is to consider the chromosomes in the population as two separate subpopulations (Hudson and Kaplan 1988)—one for the carriers of each of the two alleles at the locus—and to estimate the sizes of those two populations over time. At locus *i*, denote the sizes of these two subpopulations at time *t* as (for the allele of interest) and (for the other allele). The frequency of the allele of interest at time *t* in the past, or is(6)Here, we propose to estimate and on the basis of properties of the coalescent trees for the two alleles. We then estimate _{} by plugging these estimates, _{} and into Equation 6, giving(7)The estimator in Equation 7 is not unbiased in general, even if the estimates of and are unbiased. But separating the problems of estimating the two subpopulations confers an advantage: this estimator does not assume, as the proportion-of-lineages estimator does, that the two allelic types have had equal fitness between *t* and the present.

Assuming that the mutation that distinguished the two alleles occurred only once in the history of the sample, the chromosomes carrying the two alleles can be treated as two distinct subsamples between the time of the mutation and the present. Among the ancestors of the allele subsample carrying the allele of interest, coalescent time *τ* accrues according towhere *t* is measured in generations. It follows thatsuggesting that can be estimated by assessing the rate at which coalescent time accrues. can then be estimated analogously and can be estimated by Equation 7. We present two related approaches to estimating the rate of accrual of coalescent time in each subsample: one approach in which estimates are made with respect to waiting times between coalescent events, and another in which estimates are made with respect to the number of lineages ancestral to the subsample at a specified time in the past. In both approaches, we assume that and are piecewise constant, but it is possible to modify these estimators under other assumptions about how and change between time points.

The number of coalescence events in a time interval depends on the amount of coalescent time passed. Our approaches amount to inverting this relationship to estimate the amount of coalescent time passed on each lineage (and thus their relative population sizes) in a method-of-moments framework. In this subsection, we assess and according to the time passed between fixed numbers of coalescent events.

Suppose that assumed to be constant from a time point of interest (defined here to be ) until ℓ coalescences have occurred within the subsample, which at consists of lineages. Define as the waiting time to the *k*th coalescence, starting from the coalescence event (or from if ). Define the total waiting time from to the ℓth coalescence as Then, in units of generations, each is an independent, exponentially distributed random variable with rate Thus,(8)and(9)For large and intermediate values of *Y* is approximately normally distributed with expectation and variance given by Equation 8 and Equation 9 (Chen and Chen 2013). One estimator of which is both a method-of-moments estimator and the maximum-likelihood estimator under the asymptotically normal distribution of *Y*, is(10)With this is the “skyline” estimator of Pybus *et al.* (2000), and with it is the generalized skyline estimator of Strimmer and Pybus (2001). It is unbiased under the assumption of constant between coalescence events. In our implementation, after the tree has coalesced down to one lineage, we assume that remains constant into the past if its allele is ancestral and that remains constant before dropping to zero in the middle of the branch on which the mutation arose if its allele is derived. (In principle, derived *vs.* ancestral status may be determined by the tree topology or “forced” on the basis of prior knowledge if the tree topology may be in error. We use the tree topology.)

One may estimate by(11)where is estimated analogously to In principle, the estimates of and may be based on waiting times for different numbers of coalescences and they will likely be constant for different intervals of time. By a first-order Taylor approximation argument described in Appendix C,but in practice, can be substantially biased. Its approximate variance is given in Appendix C.

One important decision in using this estimator is the choice of ℓ—*i.e.*, the number of coalescence events to wait for—for each allele. Small values of ℓ lead to variable estimates. On the other hand, the estimator assumes that the size of the subpopulation is constant until ℓ coalescent events have occurred, which may be increasingly unrealistic for large In this article, we use and defer investigation of different choices of ℓ for future work.

#### Lineages-remaining estimator:

The next estimator, which we term the “lineages-remaining” estimator, operates on a principle similar to the waiting-time estimator. The coalescent time passed on each allele’s background is evaluated by looking at the local rate of coalescence, and the relative numbers of carriers of each allele are estimated to form an allele-frequency estimate. The difference is that whereas the waiting-time estimator estimates population sizes for each allele between coalescent events, the lineages-remaining estimator estimates the population sizes between prespecified times by comparing the number of lineages of each type that remain (*i.e.*, have not coalesced) at the more ancient end of a time interval with the number present at the more recent end of the interval.

Suppose that during an interval of generations. At the end of the interval closer to the present there are lineages, and at the end of the interval further into the past there are lineages. The expected number of lineages remaining at the end of the interval can be approximated as(12)where *τ* is the amount of coalescent-scaled time elapsed during the interval (Griffiths 1984; Maruvka *et al.* 2011; Chen and Chen 2013; Jewett and Rosenberg 2014). Here, is a haploid population size, so For large and intermediate lengths of time, is approximately normally distributed (Griffiths 1984; Chen and Chen 2013). An estimator for which is both a method-of-moments estimator—one based on the approximate value of —and the maximum-likelihood estimator under the limiting normal distribution, is(13)As the tree coalesces to few lineages, there are several edge cases in which the estimator is undefined. Our methods for handling these cases are discussed in Appendix C.

It may be impractical to estimate because in generations may be unknown. However, cancels in the estimator of the allele frequency (14)where is the number of carriers of the reference allele in the population at locus *i* at time *t*, and and are the numbers of lineages carrying the reference allele at the ends of the time interval closer to and more distant from the present, respectively. Approximate expressions for the variance of the estimator in Equation 14 are in Appendix C.

In practice, the estimator in Equation 14 will be evaluated at a set of times. Here, we evaluate the estimator every 0.001 coalescent units, starting at the present and extending back four coalescent units. The grid of times at which changes in the number of lineages of each type are considered will influence the estimate. Finer grids will lead to estimates that are more variable but less biased by the assumption of constant and between time points.

### Testing time courses for selection

Once a polygenic-score time course has been constructed, it can be tested for selection. Whereas the relevance for trait evolution of the estimators proposed in the previous sections depends on the proportion of trait variance accounted for by the polygenic score, polygenic-score time courses can be tested for selection even if they account for a small proportion of the trait variance. Further, although we apply the test in this section to estimated polygenic-score time courses, it is also applicable to measured time-series data on allele frequencies when available.

We propose a test for selection that amounts to a modification of the framework of Berg and Coop (2014), an analog of – tests for phenotypic selection (Whitlock 2008). Berg and Coop proposed to test for overdispersion of polygenic scores among population samples, relative to neutral expectations. Here, we check for overdispersion among a set of time points along one population branch. We denote our time-based modification of as

Suppose that for each time in a sequence of times, we observe a population-level polygenic score, as in Equation 1. Assume that changes across time at distinct loci are independent. Following the framework, we model polygenic scores using a multivariate normal distribution. Specifically, we posit that at each locus, over short timescales, allele-frequency changes between time points follow a normal distribution (Cavalli-Sforza *et al.* 1964; Nicholson *et al.* 2002). (Approximate normality fails over longer timescales, in part because allele frequencies are bounded by 0 and 1.) Here, *f* is the coalescent time that has passed between time points and is the allele frequency at locus *i* at one end of the interval. (In practice, we choose to be the allele frequency at the end of the interval closer to the present.) The parameter *f* is constant across loci.

If allele-frequency changes at each locus are independent and normally distributed, then changes between time points in the polygenic scores are normal with expectation 0 and variance The variance of polygenic-score changes can also be written as where is the additive genetic variance of the trait.

Imagine we have a time course of polygenic scores with Under neutrality, for each time point the statistic(15)has a normal distribution. ( can be recomputed for each value of *j* using allele frequencies at In practice, we use the same —the one computed from allele-frequency estimates closest to the present—for each time interval.) Moreover, under neutrality, allele-frequency changes in distinct time intervals are independent, so values of are independent for distinct *j*. Thus, the sum across time points(16)has a distribution. [For details on the equivalence of the statistic in Equation 16 to the statistic of Berg and Coop (2014) in this scenario, see Appendix D.] In contrast, under directional selection, changes in allele frequency across time depend on the effect size of the locus, leading to large changes in the polygenic score and values larger than predicted by the distribution. This test is an analog of time-course tests for phenotypic selection based on neutral Brownian motion (Lande 1976; Turelli *et al.* 1988), but with the advantage that we know

can be compared with a distribution to test for significance, or it can be compared with a distribution obtained by permuting either effect sizes or their signs. If the distribution is used, then coalescent times elapsed between polygenic scores can be estimated by assessing estimated allele-frequency changes between time points, either at putatively neutral loci or trait-associated loci. Using estimated times produces type-I error rates closer to the nominal value because the estimators we use tend to change at a systematically different rate than the actual allele frequencies: the proportion-of-lineages estimator changes more slowly than do population allele frequencies, and the other estimators change more quickly than population allele frequencies. In practice, we use the sample variance of as an estimate of coalescent time passed between time points, where is an estimated allele-frequency change at a variable locus and is the estimated allele frequency at the more recent end of the time interval. Assessing elapsed time on the basis of trait-associated loci may lead to a power decrement, but it will not be large—most of the signal comes from coordination of small shifts across loci and not from larger-than-expected allele-frequency changes (Berg and Coop 2014).

### Data availability

Supplemental material, including code for the article, available at Figshare: https://doi.org/10.25386/genetics.6955367. Additionally, code used for running the simulations and implementing all data analyses presented here is available at http://github.com/mdedge/rhps_coalescent. The version of the code used here is permanently archived at doi: 10.5281/zenodo.1461077.

## Simulation Results

We examined the performance of our methods in coalescent simulations. In particular, we simulated coalescent trees for unlinked loci associated with a phenotype. The simulated loci evolve neutrally or under directional selection. In particular, if an allele at locus *i* has effect size on a trait, and the trait experiences a selection gradient at time *t*, then the selection coefficient on the allele—representing the fitness of the heterozygote minus the fitness of the ancestral homozygote—is (Charlesworth and Charlesworth 2010, equation B3.7.7). The coalescent simulations were run in mssel (Berg and Coop 2015), a version of ms (Hudson 2002) that takes allele-frequency time courses that may be produced by selection, and assumed a constant population size. Effect sizes for the derived allele are drawn from a normal distribution centered at zero. For details on the simulations, see Appendix E.

We consider the performance of the methods when (1) the true trees are provided as input, and (2) when the trees must be reconstructed from sequence data. We use the software RENT+ (Mirzaei and Wu 2016) for tree reconstruction.

Before considering systematic results over many simulations, we present estimated time courses for one representative simulation (Figure 3). In the simulation shown, the trait was under selection upward in the past but has evolved neutrally recently.

The proportion-of-lineages estimator (Equation 4) estimates allele frequencies as the proportion of lineages ancestral to the sample carrying each allele. It is expected to perform well under neutrality, and it does here: in the neutral period between the offset of selection and the present, it tracks the true polygenic score closely. During the period of selection, looking backward in time, the proportion-of-lineages estimator strays off target, slowly recovering in the period before the onset of selection. Looking forward in time, it is as if the proportion-of-lineages estimator “anticipates” shifts due to selection. As mentioned in the *Proportion-of-lineages estimator* subsection, the apparent anticipation occurs because, if there has been selection between the present and time *t* in the past, then the ancestors of the present-day sample at time *t* are a biased sample from the population at time *t*. For example, if the trait has been selected upward, then the ancestors of a present-day sample will have had high trait values compared with their peers.

The waiting-time estimator (Equation 11) and the lineages-remaining estimator (Equation 14) do not rely on an explicit neutrality assumption. Instead, they track the relative passage of coalescent time—measured, roughly, in terms of coalescence events—for each allelic type. These estimators track the rapid change in the polygenic score during the period of selection much more closely than the proportion-of-lineages estimator. (Only the lineages-remaining estimator is shown in Figure 3.) At the same time, these estimators rely on a highly stochastic signal—in the case of the waiting-time estimator (with ), single coalescence events—and they are noisier than the proportion-of-lineages estimator as a result.

The patterns seen in Figure 3 reflect the performance of the methods over many simulations, as detailed in the next subsections.

### Estimator performance: bias and mean squared error

Figure 4 shows bias and mean squared error (MSE) of our estimators of the historical polygenic score across three scenarios: one in which the trait has evolved neutrally, one in which there has been recent directional selection on the trait, and a third in which there has been directional selection on the trait in the past but the trait has recently evolved neutrally. These estimators are also compared with a “straight-line” estimator: a straight line that goes from the present value to the ancestral state (*i.e.*, all derived-allele frequencies zero) in two coalescent units. In the neutral case, none of the estimators show marked bias and the proportion-of-lineages estimator has the lowest variance (and thus lowest MSE). Estimators formed from trees reconstructed by RENT+ (dashed lines) rather than the true trees (solid lines) are noisier and they do not outperform the straight-line estimator under neutrality.

In the presence of selection, the proportion-of-lineages estimator is badly biased, and the severity of the bias increases during the interval of selection (looking backward in time). The waiting-time and lineages-remaining estimators are less strongly biased in the presence of selection and they achieve similar MSEs under selection and neutrality. Again, estimators formed from RENT+ trees perform worse than estimators formed from the true trees, but in the presence of selection, they outperform the straight-line estimator.

### Interval estimation: coverage probabilities

Figure 5 shows the coverage probabilities of nominal 95% confidence intervals formed on the basis of the proportion-of-lineages, waiting-time, and lineages-remaining estimators and their (approximate) variances.

Under neutrality, and when using the true trees, all confidence intervals have approximately the correct coverage, although coverage for confidence intervals for both the lineages-remaining and waiting-time estimators decays further into the past. The decay of the coverage probability far in the past makes sense for the lineages-remaining and waiting-time estimators: both of these estimators implicitly assume that the number of carriers of each allele in the population remains constant between coalescent events. This assumption may be a reasonable approximation in the recent past, when coalescent events are frequent, but become untenable in the distant past, when coalescence times are longer.

Confidence intervals computed on the basis of RENT+ trees only achieve the nominal coverage in the very recent past and become anticonservative further back in time. This behavior is expected; the variances we use incorporate stochasticity in the coalescent process but do not account for randomness arising from errors in tree estimation. [Stochasticity in the coalescent process has been called “coalescent error” and contrasted with randomness from errors in tree estimation, or “phylogenetic error” (Ho and Shapiro 2011)].

Under selection, confidence intervals from the proportion-of-lineages estimator have very low coverage—this arises from the bias documented in Figure 4. The coverage probabilities of the waiting-time and lineages-remaining estimators are less changed by selection.

### Power of

We assessed the performance of the statistic (Equation 16) as a test statistic for detecting selection in the simulations shown in Figure 4 and Figure 5. We constructed the test statistic from the allele-frequency estimates produced by each of the proposed estimators and compared it against both the theoretical distribution and against a permutation distribution. Table 1 shows the results. Under neutrality (first two rows), comparing against a distribution formed by randomly permuting the effect sizes produces acceptable type-I error rates. (There are 100 simulations using RENT+ and the values in Table 1 do not differ significantly from 0.05.) When the theoretical distribution is used, RENT+ type-I error rates are unacceptably high, but the type-I error rates produced from the true trees are acceptable.

Of the methods with acceptable type-I error, tests using the allele frequencies estimated by the proportion-of-lineages estimator have by far the highest power. It may seem paradoxical that the proportion-of-lineages estimator is the best of our estimators at detecting selection, given that the estimated time courses produced by the proportion-of-lineages estimator are biased in the presence of selection. However, in our simulations, the proportion-of-lineages estimator generally moves in the correct direction in the presence of selection, albeit more slowly than it should. In contrast, the other two allele-frequency estimators are highly variable, leading to wide null distributions and decreased power. The proportion-of-lineages estimator can also be thought of as the mean polygenic score among lineages ancestral to the sample, and the test for selection responds to changes in the mean polygenic score of the ancestors that are faster than would be expected under the null hypothesis of neutral evolution.

With the proportion-of-lineages estimator, using true trees unsurprisingly gives better power than RENT+ trees, but RENT+ trees still have substantial power.

In Table 1, power is higher when selection occurs closer to the present. To explore the relationship between the timing of selection and present-day sample size, we conducted additional simulations. In these simulations, we assessed the power of the test (using the proportion-of-lineages allele-frequency estimates from true trees, and comparing with a permutation distribution) to detect an approximate one-standard-deviation shift in the population-mean polygenic score. We varied the timing of the shift and the present-day sample size. Figure 6 shows the results. For detecting selection close to the present, power increases with sample size. However, for selection further in the past, power reduces to the type-I error rate, regardless of the present-day sample size. This is because power to detect selection depends on unusual coalescent times during the period of selection, and by 0.1 coalescent units in the past, most coalescent events have already occurred, even in large samples. For example, even extremely large present-day samples have, in expectation, ancestors tracing back 0.01 coalescent units in the past and ancestors 0.1 coalescent units in the past (Maruvka *et al.* 2011; Jewett and Rosenberg 2014). Thus, it will likely be impossible to detect all but the strongest selective events by their signatures in coalescent trees if they are over 0.1 coalescent units in the past. ’s power to detect selection up to coalescent units into the past represents an extension of the SDS statistic (figure S6 in Field *et al.* 2016), which has excellent power in the very recent past but very little power beyond the expected length of a terminal branch [ in coalescent units, where *n* is the present-day sample size (Fu and Li 1993)]. In Appendix F, we show empirical power for a test statistic analogous to SDS computed from the lengths of the terminal branches (Figure F1). This SDS analog has similar power to near the present, but its power decays more rapidly for selection further in the past.

### Simulations with larger numbers of loci

All of the above simulations are of polygenic scores that incorporate 100 loci each. Reassuringly, the performance is extremely similar for all estimators and tests if the number of loci is increased, using the estimators based on the true trees. In Supplemental Material, Figure S1, Figure S2, and Table S1, we show results analogous to Figure 4, Figure 5, and Table 1 with 1000 loci per polygenic score.

## Empirical Application: Human Height

We applied our proposed estimators to human polygenic scores for height. Genetic variation within Europe related to human height has been studied by many investigators interested in polygenic selection (Turchin *et al.* 2012; Berg and Coop 2014; Robinson *et al.* 2015; Field *et al.* 2016; Berg *et al.* 2017; Racimo *et al.* 2018; Uricchio *et al.* 2018). A recent pair of articles compared the results produced by existing tests for polygenic selection when applied to human height using GWAS effect sizes from two different studies (Berg *et al.* 2018; Sohail *et al.* 2018). Most previous work has used GWAS effect sizes from the Genetic Investigation of Anthropometric Traits (GIANT) consortium (Wood *et al.* 2014), whereas the new work uses GWAS effect sizes from the larger and presumably less structured UK Biobank sample (Sudlow *et al.* 2015). Tests for polygenic selection on height provide much less evidence for selection when UK Biobank effect sizes are used than when effect sizes from GIANT are used. One possible explanation is that GIANT effect sizes are contaminated by some degree of population stratification.

In Figure 7, we show estimated population-mean polygenic-score time courses among populations ancestral to the GBR (British in England and Scotland) subsample of the 1000 Genomes Project (1000 Genomes Project Consortium *et al.* 2012). In the top panel GIANT effect sizes are used, and in the bottom panel UK Biobank effect sizes are used. Polygenic scores were constructed by taking the top locus in each of approximately independent genetic regions. [These polygenic scores are identical to those used in Berg *et al.* (2018); our “UK Biobank” is their “UKB-GB.”] Coalescent trees for these loci were estimated in RENT+. (Details in Appendix G.)

When the time courses are constructed using GIANT, all three estimators suggest that the population-mean polygenic score for height has increased in the recent past. Using the proportion-of-lineages approach, an increase of approximately three present-day polygenic-score standard deviations is estimated. In contrast, time courses estimated using the UK Biobank effect sizes show little apparent change in the recent past: the proportion-of-lineages estimator suggests a recent decrease of ∼0.34 standard deviations.

Further, when the change from the present to the most recent time point (0.001 units in RENT+) is assessed by the test, both sets of effect sizes yield some evidence of selection, but the evidence is stronger with GIANT effect sizes than with UK Biobank effect sizes. Specifically, with GIANT, from 10,000 permutations, and with UK Biobank, . [We do not claim that the difference between the values across the data sets is itself significant (Gelman and Stern 2006), merely that the pattern of weakened evidence in the UK Biobank matches that observed by recent work (Berg *et al.* 2018; Sohail *et al.* 2018).] However, using both data sets, the evidence for selection is limited to periods very close to the present. If the same set of times evaluated in the simulations is evaluated for the height polygenic scores (approximate coalescent times before the present), neither polygenic-score time course provides evidence for selection ( in both cases). GIANT effect sizes produced much lower *P*-values for recent selection on height in the UK in a recent article (Field *et al.* 2016), but that work used a sample of 3195 genomes, whereas the GBR subset of the 1000 Genomes sample contains only 91 genomes.

Thus, our estimators broadly recapitulate the pattern of other methods for detecting polygenic selection, finding evidence suggestive of selection when GIANT effect sizes are used but much weaker evidence when UK Biobank effect sizes are used.

## Discussion

We have proposed a set of estimators and tests for population-mean polygenic scores over time, given (additive) effect sizes for a trait at independent trait-associated loci and coalescent trees for the trait-associated loci. Estimation of the population-mean polygenic-score time course is most effective when the trait (and its associated loci) evolves neutrally and the ancestors of the sample are representative of the ancestral population. When the trait has been under selection, estimation is still possible, but the estimates obtained are noisier. Tests for polygenic selection that are based on coalescent trees have the potential to be powerful in the recent past.

In terms of practical applications, we have produced one estimator that produces good estimates of population-mean polygenic-score time courses under neutrality and that is also well powered to detect departures from neutrality (the proportion-of-lineages estimator). The other two estimators are less biased by selection, but they are variable and less useful for detecting selection. At this writing, one sensible procedure for fitting these methods to data would be to form initial estimates using the proportion-of-lineages estimator and test them for selection using the statistic. If the test suggests selection, then the ancestors of the sample may not be representative of the ancient population and polygenic-score time courses from the proportion-of-lineages estimator may be biased. In that case, the waiting-time or lineages-remaining estimators might be applied.

These methods add to a set of methods that use GWAS information to study the history of complex traits (Berg and Coop 2014; Field *et al.* 2016; Berg *et al.* 2017; Racimo *et al.* 2018; Uricchio *et al.* 2018). Many of these methods have been applied to human height, and our methods produce similar conclusions when applied to the same data (Berg *et al.* 2018; Racimo *et al.* 2018; Sohail *et al.* 2018; Uricchio *et al.* 2018). Our methods add to previous work by estimating the historical time courses of mean polygenic scores and by leveraging ARGs.

As with other population-genetic methods for studying polygenic traits, results from our methods are accompanied by many qualifiers to interpretation (Novembre and Barton 2018).

In general, estimates arising from the methods presented here should not be viewed as necessarily reflecting the historical time course of trait values within a population. In our simulations, the association between genotype and phenotype was assumed to remain constant over time. In contrast, polygenic scores estimated in practice are better thought of as functions that encode present-day associations between genotype and phenotype. The genotype–phenotype association captured by a polygenic score may be due to causal effects of the included genotypes, but it might also be due to linkage disequilibrium between tag SNPs and ungenotyped causal SNPs, to indirect genetic effects (Kong *et al.* 2018; Walsh and Lynch 2018, Chap. 22), or to environmental effects that covary with genotype for other reasons. Any of these sources of association between genotype and phenotype might change as the environment and genetic background of the population change over time, causing time courses estimated by our method to deviate from the history of average trait values in the population. For example, the genetic architecture may change across the time period over which estimates are made, for example because of changes in linkage disequilibrium between tag loci and causal loci (Martin *et al.* 2017), or because loci that explained trait variation in the past have since fixed or been lost. Further, real-world estimates of effect size will be subject to noise in estimation and possibly bias due to stratification. [Even small amounts of stratification can seriously mislead tests for selection (Berg *et al.* 2018; Sohail *et al.* 2018).] Particularly in case-control studies, ascertainment biases may also lead to confounding between the evolutionary status of an allele (*i.e.*, derived or ancestral status) and power to detect trait associations (Chan *et al.* 2014). Also, importantly, changes in the environment may drive changes in mean levels of the trait that either amplify or oppose changes in population-mean polygenic scores, either via their direct effects or via gene–environment interactions. Finally, for most human traits, current polygenic scores explain relatively small proportions of the trait variance. In addition to all the issues above that might lead to bias in an inferred time course, using a weakly predictive polygenic score adds measurement error to the inferred time course—much of the variance in the trait will not be reflected by the polygenic score. This added variance would not be expected to increase the type-I error rate of tests of neutrality, but it will make the inferred time courses less likely to reflect the history of the trait closely.

Beyond these general caveats, the methods we propose here have several limitations that suggest directions for future work. The first three limitations concern outstanding statistical issues. First, the polygenic scores we estimate here are weighted sums of effect sizes estimated under additive models. Our variance estimates also assume that the loci incorporated in the polygenic score are in linkage equilibrium. Because our estimators work by estimating historical allele frequencies at the loci contributing to the polygenic score, they can in principle be adapted to estimate any function of allele frequencies, including trait predictions that account for dominance, epistasis, and linkage among loci, and also to single-locus trajectories. However, the strategies we use here for variance estimation and hypothesis testing may need to be modified for more general functions of allele frequencies. Second, any applications of these methods to real data will entail noise that is not accounted for by the variance estimates we propose. In particular, effect sizes will be estimated with error and so will the coalescent trees for sites included in the polygenic score. It will be important to incorporate these sources of variance in future estimates of sampling variation. Third, as suggested in the *Theory* section, the waiting-time and lineages-remaining estimators implicitly contain smoothing parameters. Fully characterizing the effects of these smoothing parameters and of alternative smoothing strategies—such as those used to smooth coalescent-based estimates of population-size history (Drummond *et al.* 2005; Minin *et al.* 2008)—will reveal the potential of these estimators, which have high variance in the forms in which they are used here.

The next three possible extensions are suggested by biological applications and by the coalescent framework in which we work. First, the theory we develop here is for a single population, but our setting within a coalescent framework suggests the possibility for extension to multiple populations, perhaps by developing multivariate analogs of our statistics within a coalescent-with-migration framework (Kaplan *et al.* 1991). Similarly, whereas we work with polygenic scores for a single trait, our methods can be extended to consider polygenic scores for multiple, correlated traits. In a similar vein, Berg *et al.* (2017) have recently extended the statistic to multiple correlated traits, drawing inspiration from the framework of Lande and Arnold (1983). Working with multiple traits will allow us to distinguish hypotheses in which a trait is directly subject to selection from hypotheses in which a correlated trait is the target of selection. Finally, because the coalescent framework explicitly represents the evolution of the sample backward in time, it will be productive to incorporate ancient samples.

The methods we propose are promising in part because they capitalize on illuminating descriptions of genetic variation that will become increasingly widely available. ARGs encode all the coalescence and recombination events reflected in the present-day sample, and thus are richly informative about the history of the sample’s ancestors. Statistics computed from these ARGs have the potential to capture all the information about an allele’s frequency time course that is available in a present-day sample. Approaches to traits that are based on sample ARGs will improve with the development of our understanding of the architecture of complex traits and of our ability to reconstruct ARGs.

## Acknowledgments

We thank members of the Coop laboratory and Arbel Harpak for useful discussions, Jeremy Berg for providing the height loci and effect sizes used in Berg *et al.* (2018), and Sajad Mirzaei for support with RENT+ software. Joshua Schraiber and Jay Taylor provided helpful comments on the manuscript. The authors acknowledge support from National Institutes of Health (R01 GM-108779) and National Science Foundation (1262327 and 1353380).

## Appendix A: The Relationship Between Coalescent Rates and Phenotypic Selection

In this appendix, we describe the consequences of directional and stabilizing selection on a trait for coalescence rates at loci associated with the trait.

We begin by considering the general effects of selection on coalescence before considering the specific cases of directional and stabilizing selection. If the frequency of the allele of interest at locus *i* at time *t* is in a population of chromosomes, then the probability of any one of pairs coalescing in the preceding generation is

(Hudson and Kaplan 1988). Writing the frequency in the preceding generation as then(A2)assuming that is small so that higher-order terms in the Taylor expansion can be ignored. We can decompose the change in frequency from the preceding generation into , the contributions of random genetic drift (*N*) and selection (*S*), respectively. Then, taking the expectation over the frequency in the preceding generation,(A3)because and is the deterministic change due to selection.

If the allele of interest affects a trait that is under selection, then the form of depends on the allele’s role in the trait’s architecture and on the specifics of selection on the trait. The simplest case is one in which directional selection acts on a phenotype whose genetic architecture is purely additive. Denote the selection gradient—that is, the slope of fitness regressed on phenotype—as *α*. If the effect size of the *i*th locus on the phenotype is and the loci that affect the phenotype are unlinked, then(A4)(equation 3.17 in Charlesworth and Charlesworth 2010). Therefore, the expected coalescent rate over the possible allele-frequency changes in the preceding generation is(A5)Therefore, recent directional selection increases the coalescence rate among alleles that move the phenotype in the direction in which selection acts. For example, if larger trait values have been selected recently, the coalescence rate will be higher among alleles associated with larger values of the trait and lower among alleles associated with smaller values of the trait. This is true regardless of the frequency of the allele.

If the phenotype is subject to both directional selection and stabilizing selection, the allele frequency matters. For example, one common and reasonably general model incorporating directional and stabilizing selection is the quadratic optimum model, where the fitness of an individual with phenotype *G* is There is directional selection on the phenotype when the population mean deviates from the optimum (When the population mean is below the optimum, and selection will favor larger trait values.) Under this model,(A6)see Barton (1986), Bürger (2000), and Simons *et al.* (2018) for a recent discussion. Note the similarity of the directional selection term—with the selection gradient now controlled by —to Equation 20, resulting in similar coalescent rates to Equation 21. The Δ term for stabilizing selection matches that of a disruptive-selection model with an unstable equilibrium at

Incorporating both directional and stabilizing selection, the expected coalescent rate is approximately(A7)The stabilizing selection term decreases the coalescent rate of whichever allele has frequency < The average coalescent rate among “up” alleles depends on the value of averaged over the frequency distribution of Thus, when directional selection acts on the phenotype, we again expect an increased coalescence rate among alleles whose effects have been favored. If the population is at the optimum, *i.e.*, then averaging across loci there is no net effect on the coalescent rates if there is a symmetric distribution of ∼ (across all effect sizes). This symmetric distribution will result when the population is at equilibrium and there is no mutational bias (*i.e.*, there is no difference in the mutational input of up and down alleles).

## Appendix B: A Bayesian View of the Proportion-of-Lineages Estimator

To arrive at a Bayesian view of the estimator proposed in Equation 4, define the frequency of an allele of interest at a biallelic locus at time *t* in the past as If the allele of interest is equally likely to be derived or ancestral, a reasonable prior distribution would be proportional to the limit (as the mutation rates approach 0) of the stationary distribution of the neutral diffusion with mutation, or (equation 5.70 in Ewens 2004). Assume that the two alleles at the locus have experienced equal reproductive success between the present and time *t*—that is, assume that the locus has evolved neutrally in the time interval. In that case, the lineages ancestral to a random sample of chromosomes in the present are themselves a random sample from the population at time *t*, and if there are *r* ancestral lineages at time *t*, then the number carrying the allele of interest is a random variable. If we observe that *j* out of *r* ancestral lineages at time *t* carry the allele of interest, then the posterior distribution for *x* is proportional toNoticing that is the kernel of a distribution gives the posterior,

Equation B1 suggests the posterior mean of as an estimate of the frequency (when and ), which is equal to Equation 4. The posterior variance is(B2)which differs from the binomial sampling variance in Equation 5 by a factor of

It may be interesting to note that there is a connection between this Bayesian view of the estimator and the neutral Wright–Fisher diffusion process, reviewed by Tavaré (1984). Under a neutral Wright–Fisher diffusion model without mutation, the probability that an allele frequency changes forward in time from to in time *t* can be written as (equation 7.18 in Tavaré 1984)(B3)here is the probability that *r* lineages of an initially infinite number of lineages survive to time *t*. Because the neutral diffusion is time reversible (Ewens 2004),and the expression on the right side of Equation B3 is thus also an expression for the allele frequency in the past given the allele frequency in the present, One way to interpret Equation B3 is that if *r* lineages survive to time *t*, then the number *j* of lineages carrying a focal allele at time *t* is a random variable. Then, conditional on *j*, the frequency has a distribution.

## Appendix C: Mathematical Details for the Waiting-Time and Lineages-Remaining Estimators

### Taylor-Series Expansions for Expectations and Variances of Ratios

Both the waiting-time and lineages-remaining estimators take the form(C1)where and and are population-size estimators that are independent of each other. We assume that Equation C1 is equivalent to Equation 7 from the main text, with the subscripts and time notation removed for compactness.

In this subsection we present general Taylor-series approximations for the expectation and variance of the estimator in Equation C1 in terms of the expectation and variance of and Similar presentations are available in other references (sections 10.5–10.6 in Stuart and Ord 1987), but we present the argument here for completeness.

The first-order Taylor expansion of Equation C1 evaluated at is(C2)where and are the partial derivatives of *f*, and *Z* represents the contribution of higher-order terms. The expectation of is thereforewith the second step coming from the linearity of expectation and the fact that is not random. Choosing and to be equal to the expectations of and gives(C3)where the approximation depends on the higher-order terms represented by *Z* being small. Equation C3 underlies the claim in the main text that the waiting-time estimator is approximately, to first order, unbiased (under the unrealistic assumption that the number of carriers of each allele is constant between coalescent events).

To obtain the approximate variance of the estimator in Equation C1, we use Equation C3 and note that(C4)Substituting the first-order Taylor expansion in Equation C2 for ignoring the higher-order terms represented by *Z*, gives(C5)The last step comes from the definitions of variance and covariance and the fact that we choose and to be equal to the expectations of and . In our case, , implying that and Plugging these values into Equation C5 gives(C6)with the useful alternative form(C7)Finally, in our setting, and and are independent conditional on the true but unknown allele-frequency time course (Kaplan *et al.* 1988; Coop and Griffiths 2004). Conditional independence gives and Using these identities in Equation C7 gives

### Approximate Variance of the Waiting-Time Estimator

As noted in the main text, the waiting-time estimator has the formwhere*Y* is the waiting time associated with a set of ℓ coalescent events on the allele of interest’s background that envelop time *t*, and is the number of allele-of-interest lineages that exist before any of the ℓ coalescent events have occurred. is analogous, substituting in a waiting time and number of lineages from the other allele’s background, and possibly a different value of

It is also noted in the main text (Equation 9) that the waiting time *Y* has varianceThus, the variance of is(C9)If Equation C9 reduces to If thenand By Equation C8 and the fact that and are unbiased, the approximate variance of the waiting-time estimator is(C10)where and are the variances of and respectively. If is used to estimate both and then and Equation C10 reduces to(C11)If and are estimated using the same value of and ℓ is much smaller than the starting number of lineages for both loci, so that and then Equation C9 can be rewritten as

### Approximate Variance of the Lineages-Remaining Estimator

In the main text, we propose to estimate as(C12)where is the number of generations elapsed between two ends of a predefined time interval, is the number of lineages of the allele of interest at the end of the time interval closer to the present, and is the number of lineages of the allele of interest at the end of the interval further into the past. The only random component of is which depends on the number of coalescent events that occur during the interval. This estimator of is leads to an estimator for the allele frequency of interest,(C13)(The expression in Equation C13 is the same as in Equation 14.)

To obtain an approximate sampling variance for we use Equation C8, making the substitutionsandTo use Equation C8, we need (at least approximate) expectations and variances for both these terms. To compute an approximate expectation, we replace and with their expectations. Taking the case of notice that by Equation 12,(C14)In Equation C14, *τ* represents the coalescent time passed on the allele of interest’s background, or and so(C15)and analogously for (C16)Next we compute the approximate variances of these terms. Because is a constant,(C17)By a first-order Taylor approximation, the variance is approximately(C18)where is the derivative of the right side of Equation C17 with respect to evaluated at (The argument for this approximation is entirely parallel to the one in Equation C5, but for a function of a single random variable.) The derivative of the right side of Equation C17 with respect to evaluated at is(C19)To write the approximate variance of we use the asymptotic variance derived by equation 5 in Griffiths (1984), using the version for and (see also Chen and Chen 2013). To change the expression into our notation, we make the replacements (because we assume only a single mutation distinguishing the two allelic types, so conditional on that mutation), and (by Griffiths’ equation 4). Doing so and simplifying yields(C20)Plugging the expressions for and into Equation C18 gives an approximate variance,(C21)Analogously, for the term,(C22)Now, to obtain an approximate variance for we plug the expressions in Equation C21, Equation C16, Equation C21, and Equation C22 into Equation C8. Doing so gives the expressionNoticing that expanding the products in the denominators and factoring out common terms gives:Finally, noticing that and and defining the coalescent time passed with respect to the whole sample, gives(C23)To estimate the approximate variance in Equation C23, we replace with its estimated value from the lineages-remaining estimator, with the (known or estimated) coalescent time elapsed between time points, and the expectations of and with their realized values.

The approximate variance in Equation C23 relies on an asymptotic variance of the number of ancestral lineages from Griffiths (1984). In practice, Equation C23 gives variance estimates that are much too large when the number of lineages are small and the time passed is short. For example, with and Equation C23 produces Because its maximum variance is 0.25. Our strategy is to use the minimum of Equation C23 and the estimated variance of the waiting-time estimator with or (Equation C11). The rationale is that represents the variance obtained when one estimates from only one coalescent event on each background—the lineages-remaining estimator is always based on at least one coalescent event on each background.

### Computing the Lineages-Remaining Estimator for Edge Cases

If no coalescences have occurred in the interval and then the estimator in Equation 13 is undefined. In this case, we compute the estimator for a larger value of in which at least one coalescence event has occurred, extending to the next time on the predefined grid in which a coalescence has occurred. The estimator is also undefined if because the expected number of lineages remaining approaches but never reaches 1. In this case, we assume the coalescent time passed in the interval is which is the expected amount of time to coalesce from lineages to one lineage. Finally, once a subtree coalesces to one lineage, we assume that population size remains constant into the past for the ancestral allele and, for the derived allele, that it remains constant before dropping to zero at in the middle of the branch connecting the derived subtree to the rest of the tree.

The approximate variance in Equation C23 is undefined if the estimated allele frequency is equal to 0 or 1, in which case we define the estimated variance of to be 0. It is also undefined if only one lineage remains on either background at the more recent end of the time interval [ = 1 or =1]. If one lineage remains on one of the backgrounds, then we estimate the variance as the variance of the waiting-time estimator with If only one lineage remains on each of the two backgrounds, we define the estimated variance of to be which is the variance of a uniform random variable.

## Appendix D: The Correspondence of and

In this section, we explain the correspondence between the statistic proposed in Equation 16 and the statistic of Berg and Coop (2014). If the population history is tree-like, then can be thought of as checking for overdispersion (relative to neutral expectations) of polygenic scores at the tips of the population tree. In contrast, checks for overdispersion at a set of time points along one branch of a population tree (*i.e.*, a single population at different times).

Briefly, Berg and Coop (2014) start with a vector whose elements are population-mean polygenic scores for a set of populations minus an assumed value for the population-mean polygenic score in an ancestral population. Under neutrality, follows a multivariate normal (MVN) distribution, where is the additive genetic variance, and ** F** is a matrix reflecting the shared history of the populations. In particular, if the population history is tree-like, then entry

*F*in the

_{ij}**matrix reflects the branch length shared by populations**

*F**i*and

*j*since their descent from the ancestral population. If then the transformed vector (where

**comes from the Cholesky decomposition of**

*C***, such that ) obeys an MVN distribution. Thus, the sum of the squared entries of the transformed vector X (the quantity labeled ) obeys a distribution under neutrality, where**

*F**w*is the length of

In the setting of this article, we have population-mean polygenic scores for the population ancestral to one present-day population, assessed at various times. Each *t _{j}* represents some amount of time before the present and In practice, we set and treat the population-mean polygenic score in the present, as Berg and Coop treated the ancestral polygenic score, soIn this setting, with the present-day population serving the role of “ancestor,” the shared branch length for two time points along one branch of a population tree is their shared distance backward from the present day, which is to say the branch length from the present back to the more recent of the time points. Thus, for this instance of each entry

*F*in the

_{ij}**matrix is equal to , meaningThe Cholesky decomposition then givesThe inverse of this matrix isIf is a constant, then the entries of are equal to the values of Equation 15, and is equal to the quantity in Equation 16.**

*F*## Appendix E: Simulation Details

We ran simulations of polygenic scores, as well as coalescent trees and haplotypes associated with them, with and without selection.

For each trait, we simulated the frequency of the derived allele at each of *k* independent, trait-associated loci. We ensured that at all loci the minor allele frequency was at least 0.01, both at the onset of selection (if applicable) and at the present. Allele-frequency time courses were simulated using the normal approximation to the diffusion.

In particular, for each locus, an effect size for the derived allele was drawn from the normal distribution with expectation zero and variance equal to where is the desired heritability (because we are focused on polygenic scores rather than traits themselves, ), is the desired variance of the trait (*i.e.*, polygenic score; always 1 in our simulations), *w* is a modified version of Watterson’s constant for the population with *c* equal to the minimum minor allele frequency being drawn (here, ), and *k* is the number of independent loci affecting the trait. Drawing effect sizes from this distribution gives the property that, under neutrality, if the loci affecting the trait are independent and a trait is formed by adding the individual’s polygenic score to an independent “environmental” term with variance , then the resulting trait has variance in expectation.

Once the effect size was specified, we selected a derived-allele frequency from the population-level neutral site frequency spectrum, conditional on the requirement that the minor allele frequency be at least 0.01. That is, we set the probability of drawing a derived-allele frequency to be proportional to for all (and zero otherwise). Under neutrality, this randomly drawn allele frequency was treated as the derived-allele frequency at the present. In simulations including a period of selection before the present, this frequency was used as the allele frequency at the start of selection. In simulations with selection, the derived-allele frequency was simulated forward in time in steps of coalescent units (*i.e.*, one diploid generation). Conditional on frequency the frequency at the next time step, was drawn from a normal distribution with expectation and variance Here, *s* is the selection coefficient on the derived allele at time *t*, which is computed as (equation 3.17 in Charlesworth and Charlesworth 2010), where *α* is the selection gradient on the trait at time *t* and *β* is the effect size of the derived allele. Two clarifications about the selection coefficient: first, this is the selection coefficient for the heterozygote, meaning that if the fitness of the ancestral homozygote is proportional to 1 (marginalizing over the other loci that affect the trait), then the fitness of the heterozygote is and the fitness of the derived homozygote is Similarly, the effect size *α* is the expected trait value difference between the heterozygote and the ancestral homozygote, holding the genotype at all other loci constant. Second, the selection gradient is equal to where *S* is the selection differential on the trait—that is, the mean difference in phenotype between the parents of the next generation and the overall population (section 3.3.iii in Charlesworth and Charlesworth 2010)—and is the trait variance. Here, we used the expected trait variance at the onset of selection (which was 1) to set the value of *α*, which was retained over the course of the period of selection, regardless of the realized value of the trait variance.

Once the allele-frequency time course was simulated forward to the present, it was retained if the minor allele frequency in the present was at least 0.01. (Otherwise, the time course was discarded and the process was repeated.) To simulate the period from the origin of the mutation up until the onset of selection, we simulated steps backward in time using the normal approximation to the neutral diffusion, conditional on loss of the allele. In particular, given a derived-allele frequency the allele frequency coalescent time units further in the past was drawn from a normal distribution with expectation and variance (Przeworski *et al.* 2005; Berg and Coop 2015; Lee and Coop 2017).

In each set of simulations with selection, we retained polygenic-score time courses if their difference in value between the beginning and end of the period of selection was within 5% of the target value, which is given by where *N* is the effective population size, is the duration of selection in coalescent units, and *S* is the selection differential.

To simulate coalescent trees and haplotypes at each locus, these allele-frequency time courses were used as input to mssel (Berg and Coop 2015), a version of ms (Hudson 2002) that can incorporate by conditioning on user-specified allele-frequency time courses. In mssel, we set the sample size to 200 and the number of derived chromosomes to the value of a binomial random draw, where *p* is the present-day frequency of the derived allele. If the random draw was equal to 0 or 200, then it was set to 1 or 199, respectively, to force the locus to be segregating in the sample. We chose a population size of 10,000, a haplotype length of 200,000 (with the effect locus at position 100,000), a per-base-pair mutation rate of and a per-base-pair recombination rate of These values imply population-scaled ms inputs of -r 199.5 and -t 159.68. Here, and our slight differences from these values arise from the fact that we compute where *f* is the probability distribution function of either a binomial random variable or a binomial random variable, respectively, with That is, we define the locus-wide recombination (or mutation) rate as one minus the probability that no recombinations (or no mutations) occur anywhere in the locus. We used the flags -T and -L to include the true trees and times to the most recent common ancestor in the output.

The coalescent trees produced by mssel for the selected sites were used as the true trees, and the haplotypes produced were used as input to RENT+ to generate estimated trees (Mirzaei and Wu 2016). Both true and estimated trees were read into R (R Core Team 2013) and handled using the ape package (Paradis *et al.* 2004). All the estimators and tests used in this article were coded in R.

## Appendix F: Power Comparison with an Analog of trait SDS

Field *et al.* (2016) proposed the SDS as a test statistic for detecting very recent selection. Because recent selection distorts the terminal branches of the coalescent at a selected locus, Field *et al.* reasoned that recent selection would alter the distribution of singleton mutations around selected sites. Specifically, favored alleles are expected to have short terminal coalescent branches and thus relatively few singletons nearby compared with disfavored alleles.

Field *et al.* use distances to the nearest singletons on each allele’s background to estimate the difference in the mean length of the terminal branches for the two alleles. This estimated mean difference is then standardized by an empirical mean and standard deviation computed from neutral simulations at a similar derived-allele frequency. To test for selection on polygenic traits, they construct a trait SDS (tSDS) by summing each locus’s standardized SDS, signed so that positive values suggest selection for higher trait values.

Here, we compute power for an analog of tSDS in the simulations shown in Figure 6, with the differences that in our case the true terminal branch lengths are known, and we use the true effect sizes in computing tSDS. Both of these changes should enhance the power of tSDS. Specifically, at each locus the SDS value is

where is the mean terminal branch length for the ancestral allele, is the mean terminal branch length for the derived allele, is the sample mean of computed from 500 neutral simulations at the same sample size and derived-allele frequency (in bins of 0.005), and is the sample standard deviation of from the same 500 neutral simulations. [It is also possible to replace and with analytical values from Fu and Li (1993) under the assumption that the allele frequency does not change. Doing so gives similar results to those we report here.] To compute the tSDS statistic for a polygenic score, we use(F2)where is the derived-allele effect size at locus *i* and is the value of Equation F1 for locus *i*. [In Field and colleagues’ version, is replaced with which has lower power when the effect size is known.] A polygenic score will tend to have higher tSDS if higher polygenic-score values have been selected for recently. To test for significance, we form a permutation distribution by recomputing tSDS 10,000 times while randomly permuting the effect sizes.

The power of the tSDS analog as a function of the timing of a selective pulse leading to a change in the polygenic score of approximately one standard deviation is shown in Figure F1. We also show the power of with a sample size of 1000 chromosomes for comparison (from Figure 6). Compared with the tSDS has lower power, and its power drops off more rapidly as the timing of selection moves further into the past. This is not surprising because tSDS relies only on terminal branches, and uses all the intercoalescent times.

At this writing, tSDS has the advantage of scaling readily to large samples. In contrast, requires a reconstructed tree, which is currently error prone and computationally intensive. However, as tree reconstruction becomes faster and more accurate, the ability to use information from nonterminal branches will provide power benefits, especially beyond the very recent past.

## Appendix G: Height Analysis Details

We estimated the time courses of two of the polygenic scores for height studied in Berg *et al.* (2018). Information about the loci included in each polygenic score (including reference SNP ID, effect size, chromosome, and position) was provided by Jeremy Berg. Each polygenic score was initially constructed by choosing the locus with the lowest *P*-value for a test of association with height, conditional on a minor allele frequency of at least 5%, within each of ∼1700 approximately independent genomic regions defined by Berisa and Pickrell (2016). (For GIANT, the polygenic score includes 1697 loci, and for UK Biobank, 1700 loci are included.) For the GIANT polygenic score, effect sizes and *P*-values were taken from Wood *et al.* (2014), and for the UK Biobank polygenic score, effect sizes and *P*-values were taken from Neale Lab (2017). [Berg *et al.* (2018) includes several polygenic scores constructed from UK Biobank effect sizes; the one we use corresponds to their UKB-GB.]

We used sequence information from the 1000 Genomes GBR subsample (1000 Genomes Project Consortium *et al.* 2012, release 20130502), which includes 91 genomes. For each locus included in each polygenic score, we extracted phased sequence information in a window extending 100,000 bases from the locus on each side using tabix (Li 2011). The resulting .vcf file was processed into a form acceptable by RENT+ using vcftools (Danecek *et al.* 2011) and R (R Core Team 2013). We then used RENT+ (Mirzaei and Wu 2016) to estimate an ARG for the region, including the “-t” flag to estimate branch lengths. Our version of RENT+ was modified slightly to print branch lengths to greater precision and to print RENT+’s internal estimate of *θ* that it uses to scale the branch lengths.

For both polygenic scores, we rescaled the effect sizes so that their standard deviation in the present-day sample is one, assuming linkage equilibrium among loci. In particular, if the sample frequencies of the effect alleles at the *k* loci are and the effect sizes are then if the loci are independent, the variance of the polygenic-score *Z* isIf is computed from the original effect sizes, then we rescale the effect sizes as when computing the polygenic scores.

Finally, we rescaled the branch lengths estimated by RENT+. RENT+ estimates a per-nucleotide, population-scaled mutation rate using Watterson’s estimator (Watterson 1975). That is, for a sample of *n* haplotypes each covering a region of *w* base pairs, with *S* the number of segregating sites in the sample, *θ* is estimated asThe time in coalescent units separating a pair of haplotypes is then estimated as where *H* is the Hamming distance between a pair of haplotypes. For samples of size greater than two, RENT+ estimates branch lengths using a similar logic but using local distance matrices that take into account inferred recombinations.

We scaled RENT+’s estimated times because the Watterson estimator is biased downward if the population has been growing exponentially, and the human population has grown superexponentially (Keinan and Clark 2012). As a result, the inferred coalescent times from RENT+ were implausibly large. We first multiplied branch lengths at each locus by the *θ* estimate at each locus, giving branch lengths in units of twice the mutational distance per base pair at all sites. To convert these mutational distances into approximate coalescent units, we computed the time to the most recent common ancestor (tMRCA) at each locus, and then rescaled the branch lengths at all loci by a constant factor that set the mean tMRCA to two, the expectation in coalescent units under neutrality.

## Footnotes

Supplemental material available at Figshare: https://doi.org/10.25386/genetics.6955367.

*Communicating editor: R. Nielsen*

- Received August 10, 2018.
- Accepted October 23, 2018.

- Copyright © 2019 by the Genetics Society of America