Multi-locus Analysis of Genomic Time Series Data from Experimental Evolution

Summary

A growing number of experimental biologists are generating “evolve-and-resequence” (E&R) data in which the genomes of an experimental population are repeatedly sequenced over time. The resulting time series data provide important new insights into the dynamics of evolution. This type of analysis has only recently been made possible by next-generation sequencing, and new statistical procedures are required to analyze this novel data source. We present such a procedure here, and apply it to both simulated and real E&R data.

Introduction

A common study design in population genetics consists of collecting genomic variation data from living organisms to make inferences about unobserved evolutionary and biological phenomena. The many areas where this design has been applied include demographic inference (see [1] for a recent review), recombination rate estimation [2–6], and detection of natural selection [7–13]. Recently, there has been much interest in utilizing time series genetic data—e.g., from ancient DNA [14–21], experimental evolution of a population under controlled laboratory environments [22–26], or direct measurements in fast evolving populations [27]—to enhance our ability to probe into evolution. In particular, understanding the genetic basis of adaptation to changes in the environment can be significantly facilitated by such temporal data. Specifically, the dynamics of allele frequencies in an evolving population potentially convey added information about how the genome functions [28], information which is inaccessible to methods which operate only on a static snapshot of that genome.

An experimental methodology which serially interrogates the genomes of an controlled population over time could potentially yield new insights. In fact, this methodology can now be realized thanks to the advent of next-generation sequencing. By sequencing successive generations of model organisms raised in a controlled environment, genetic time series data can be generated which describe evolution at nucleotide resolution [24, 25, 28, 29]. This so-called evolve-and-resequence (henceforth, E&R) methodology is fundamentally different than the observational approach described above, and new inference procedures are needed to analyze this type of data.

In this paper, we present such a procedure and study its ability to perform a number of testing and estimation tasks relevant to population genetics. Our method is based on an approximation to the multi-locus Wright-Fisher process, and is well-suited to the small population, discrete generation, and random mating setting in which many E&R experiments are conducted. Furthermore, because it is based on a canonical population genetic model of genome evolution, our method can directly estimate population genetic quantities such as fitness, dominance, recombination rate, and effective population size. It can also be used to design future experiments with sufficient power to reliably infer these quantities.

We first use simulated data to demonstrate the utility of our method. Then, we apply our method to analyze genome-wide data from a real E&R experiment of D. melanogaster, designed to study the adaptation to a novel laboratory environment over tens of generations.

Related work

There is a small but growing literature on the analysis of evolve-and-resequence data. Feder et al. [30] present a statistical test for detecting selection at a single biallelic locus in time series data. (Although it is not a major focus, their method can also be used to estimate the selection parameter.) Similar to our method, they model the sample paths of the Wright-Fisher process as Gaussian perturbations around a deterministic trajectory in order to obtain a computable test statistic. However, their aim is slightly different from ours in that they analyze yeast and bacteria data sets where the population size is both large and must be estimated from data. Here we focus on population sizes which are smaller and more typical of experiments performed on higher organisms, for example mice or Drosophila. We generally assume that the effective population size is known but also test our ability to estimate it from data. Also, because of the increased amount of drift present in the small population regime, we necessarily restrict our attention to selection coefficients which are somewhat larger than those considered by Feder et al. Finally, although Feder et al. do study the performance of their method when time series data are corrupted by noise due to finite sampling (as in e.g. a next-generation sequencing experiment), they do not model this effect. Here we properly account for the effect of sampling by integrating over the latent space of population-level frequencies when computing the likelihood.

Another related work is Baldwin-Brown et al. [31], which presents a thorough study of the effects of sequencing effort, replicate count, strength of selection, and other parameters on the power to detect and localize a single selected locus segregating in a 1 Mb region. Results are obtained by simulating data under different experimental conditions and comparing the resulting distributions of allele trajectories under selection and neutrality using a modified form of t-test. Because it is not model-based, this method is incapable of performing parameter estimation. As a result of their study, Baldwin-Brown et al. present a number of design recommendations to experimenters seeking to attain a given level of power to detect selection. In a related work, Kofler and Schlötterer [32] carried out forward simulations of whole genomes to provide guidelines for designing E&R experiments to maximize the power to detect selected variants.

Illingworth et al. [33] derive a probabilistic model for time series data generated from large, asexually reproducing populations. The population size is sufficiently large (on the order of ∼ 10⁸) that population allele frequencies evolve quasi-deterministically. The deterministic trajectories are governed by a system of differential equations describing the effect of a selected (“driver”) mutation on nearby linked neutral (“passenger”) mutations. Randomness arises due to the finite sampling of alleles by sequencing. The main difference between the setting of Illingworth et al.’s and our own concerns genetic drift. While drift may be ignored when studying a large population of microorganisms, we show that it confounds our ability to detect and estimate selection in populations of order ∼ 10³. Thus, for E&R studies on (smaller) populations of macroscopic organisms, methods which assume that allele frequencies evolve deterministically may not perform as well as those which explicitly take drift into account.

Topa et al. [34] present a Bayesian model for single-locus time series data obtained by next-generation sequencing. In each time period, the allele count is modeled as a draw from a binomial distribution with number of trials equal to the depth of sequencer coverage, and success probability equaling the population-level allele frequency. The posterior allele frequency distribution is used to test for selection by comparing a neutral model to one in which unobserved allele frequencies to depend on time. In the non-neutral case, a Gaussian process is used to allow for directional selection acting on the posterior allele frequency distributions.

Finally, Lynch et al. [35] derive a likelihood-based method for estimating population allele frequency at a single locus in pooled sequencing data. The method allows for the possibility of sequencing errors as well as subsampling the population prior to sequencing. Using theoretical results as well as simulations, the authors give guidelines on the (subsampled) population size and coverage depth needed to reliably detect a difference in allele frequency between two populations. Unlike the other methods surveyed here, the approach of Lynch et al. is not designed to analyze time series data. Hence the data requirements needed to reliably detect allele frequency changes using their method—for example, sequencing coverage depth of at least 100 reads—are potentially greater than for methods are informed by a population-genetic model of genome evolution over time.

Novelty of our method

Our method differs from the above-mentioned approaches in several regards. To the best of our knowledge, ours is the first method capable of analyzing time series data from multiple linked sites jointly. We find that this is advantageous when studying selection in E&R data. Furthermore, it enables us to analyze features of these data which cannot be studied using single-locus models, such as local levels of linkage disequilibrium and the effect of a recombination hotspot. Additionally, because our model is based on a principled approximation to the Wright-Fisher process, it can numerically estimate the selection coefficient, dominance parameter, recombination rates, and other population genetic quantities of interest. In this way it is distinct from the aforementioned simulation-based methods [31, 32], methods which only focus on testing for selection [30, 31, 34], or methods based on general statistical procedures which are not specific to population genetics [34, 35].

Software and data availability

Source code implementing the method described in this paper is included in S1 Code. The experimental data analyzed in Analysis of a real E&R experiment data are from Franssen et al. [36] and are available on the Dryad digital repository http://dx.doi.org/10.5061/dryad.403b2.

Results

As described above, the primary methodological advance of this paper is to derive a tractable approximation to the discrete, multi-locus Wright-Fisher model with selection. This approximation enables us to perform statistical inference on time-series data generated in E&R experiments. Before studying how our approximation performs on both simulated and real data, we give a brief overview of its motivation and derivation.

A brief overview of the method

We consider the following model of an E&R experiment. A sexually reproducing population of N diploid individuals is evolved in discrete, non-overlapping generations. Pooled DNA sequencing [37, 38] is performed T times at generations t₁ < t₂ < ⋯ < t_T. At each segregating site in the resulting data set, we assume that there are two alleles, denoted A₀ and A₁. (As will be seen below, up to a change in the sign of the selection coefficient associated with each site, the model is agnostic to which allele is called A₀ or A₁.) Let L and R denote the number of loci and the number of experimental replicates, respectively. The array D ∊ [0, 1]^T×L×R counts relative frequency with which the A₁ allele was observed for each combination of generation, locus and replicate.

Given D and a vector of underlying population-genetic parameters θ, let ℙ(D∣θ) denote the model likelihood. In an idealized E&R experiment, generations are discrete and non-overlapping, mating is random, and the population size is fixed, so the likelihood is well approximated by the classical Wright-Fisher model of genome evolution [39]: