Modeling the effect of changing selective pressures on polymorphism and divergence

Abstract

The most common models of sequence evolution used to make inferences about adaptation rely on the assumption that selective pressures at a site remain constant through time. Instead, one might plausibly imagine that a change in the environment renders an allele beneficial and that when it fixes, the site is now constrained-until another change in the environment occurs that affects the selective pressures at that site. With this view in mind, we introduce a simple dynamic model for the evolution of coding regions, in which non-synonymous sites alternate between being fixed for the favored allele and being neutral with respect to other alleles. We use the pruning algorithm to derive closed forms for observable patterns of polymorphism and divergence in terms of the model parameters. Using our model, estimates of the fraction of beneficial substitutions α would remain similar to those obtained from existing approaches. In this framework, however, it becomes natural to ask how often adaptive substitutions originate from previously constrained or previously neutral sites, i.e., about the source of adaptive substitutions. We show that counts of coding sites that are both polymorphic in a sample from one species and divergent between two others carry information about this parameter. We also extend the basic model to include the effects of weakly deleterious mutations and discuss the importance of assumptions about the distribution of deleterious mutations among constrained non-synonymous sites. Finally, we derive a likelihood function for the parameters and apply it to a toy example, variation data for coding regions from chromosome 2 of the Drosophila melanogaster subgroup. This modeling work underscores how restrictive assumptions about adaptation have been to date, and how further work in this area will help to reveal unexplored and yet basic characteristics of adaptation.

Figure 1

Transition rate diagrams. (A) Synonymous site; (B) Non-synonymous site. Double circles indicate states that could be polymorphic.

Figure 2

Types of adaptations in the model. Each row depicts a process at a site, represented by a quadruplet; the fixed nucleotide is marked with a frame, light shaded nucleotides are neutral or equally beneficial, and dark shaded nucleotides are deleterious. (A) An unconstrained adaptation with a fixation event: A change in the environment affects an unconstrained site where A is fixed, and now favors only C. The intermediate state is short-lived, as eventually, through mutation and selection, C fixes at the site, which is now constrained. (B) An unconstrained adaptation without a fixation event: By chance, C is the fixed nucleotide before the environmental change. The change again favors C, turning the site into a constrained one, such that selection will keep C fixed. (C) A constrained adaptation: A change in the environment affects a constrained site with A fixed, favoring only C. Again, through mutation and selection, C eventually fixes at the site.

Figure 3

Estimates of the rate of adaptation. (A) Estimates of the total rate of adaptive events ν given the parameter g (obtained using the function FindRoot in Mathematica 8). Note that the rate is presented as its proportion with respect to the (constant) neutral mutation rate μ; assuming the value μ =0.058 mutations per Myr in Drosophia (Haag-Liautard et al. 2007), the estimate of ν varies between 1.83·10⁻³ and 3.14·10⁻³ adaptations per Myr. The shaded regions denote the approximate 95% confidence interval (see section 3.6). (B) Estimates of α for a given value of g (we note that, at this scale, the confidence region is barely perceptible). For comparison, we also estimate α according to the method of Smith and Eyre-Walker (2002) using the same parameters, shown here as a red dashed line.

Figure 4

Non-synonymous divergence. (A) Simulated curves of equal non-synonymous divergence between D. simulans and D. yakuba, given values of g and ν (darker shades represent lower divergence values). The measured value is shown as a dashed orange line. (B) The same graph with superimposed curves of equal intersection set size (that is, non-synonymous sites that are divergent between the two species and polymorphic in their ancestor).

Figure 5

Tree topologies. (A) Topology 0 shows a hypothetical data set, in which divergence is measured between two species and polymorphism is measured in their ancestral population. (B) Topologies 1 and Topology 2 (shown in C) represent more realistic data sets, in which divergence is measured between the species with red circles and polymorphism is measured in the one with orange circles.

Figure 6

The intersection sets in four cases. See text for details.

Figure 7

Probability that a site is polymorphic or divergent. The probability that a non-synonymous site is both divergent and polymorphic (b_n, left) and the probability that it is divergent (d_n, right), for a given value of g and for ν = 0.03μ, 0.04μ, 0.05μ, 0.06μ (darker to lighter lines, respectively); the rest of the parameters are taken from our estimates for the Drosophila melanogaster subgroup (see Supplementary Information). (A) For the hypothetical case (corresponding to topology 0 in Figure 5) in which divergence is measured between D. simulans and D. yakuba and polymorphism measured in their ancestral population (with parameters for the ancestral population taken from D. simulans). (B) Using divergence between D. yakuba and D. melanogaster and polymorphism in D. simulans (topology 1 in Figure 5) (C) Using divergence between D. yakuba and D. erecta, polymorphism in D. simulans (topology 2 in Figure 5).

Figure 8. The effects of weakly deleterious mutations on estimates of adaptive parameters

We plotted the estimates of model parameters as a function of δ, the ratio of the probability of observing a segregating site at unconstrained versus a constrained site (obtained using the function FindRoot in Mathematica 8). The graphs focus on the range of δ where the model would fit the data. The blue curves correspond to topology 1 and the red to topology 2. (A) f as a function of δ. (B) ν/μ as a function of δ. (C) g as a function of δ. (D) α as a function of δ.

Figure 9

Generalized definition of g. We define for each site the set F of (equally) beneficial nucleotides; mutations within F are neutral and mutations outside F are deleterious, such that the fixed nucleotide is always expected to be in F. Selective changes (adaptations and relaxations) are precisely those that alter this set. For example, in our model, we allow only #F = 1 or 4 (where #F denotes the number of nucleotides in the set): #F=1 for constrained sites, where there is only one preferred nucleotide, and #F=4 for unconstrained sites, where all nucleotides are equally favorable and thus all mutations are neutral. Here, unconstrained adaptations reduce the set F, constrained adaptations change its contents, and relaxations expand it. A general definition can be made in a similar fashion: let F₀ and F₁ denote the favorable set of a site before and after a selective event, respectively, so (A) if F₀ ⊂ F₁ (i.e., the set was expanded), the event is a relaxation; (B) if F₀ ⊃ F₁ (the set was reduced), the event is an unconstrained adaptation; and otherwise, (C) this event is a constrained adaptation. Now, g is defined like before as the fraction of unconstrained adaptations (B) out of all adaptive events (B and C). This definition accords with the notion of novelties and modifications (see Discussion).

Figure 10

Performance of inference of adaptive parameters with genome-wide data. We picked 24 combinations of adaptive parameters, ν/μ = 0.03, 0.04, 0.05, 0.06 and g = 0, 0.2, 0.4,…, 1, and calculated the region in which their inferred values would fall in 95% of cases (see details in Supplementary Information). The sample sizes (~3·10⁶ non-synonymous sites) and other parameters were chosen based on the Drosophila data from chromosome 2 (see Table 4). The resulting regions are shown in terms of ν/μ and g (A and C) and of α and g (B and D) for topologies 1 and 2, correspondingly.

Figure 11

Phylogenetic tree of the Drosophila melanogaster subgroup. See Methods for details.

Figure A1

Nodes and branches in a topology 1 tree.