Faster than neutral evolution of constrained sequences: the complex interplay of mutational biases and weak selection

Abstract

Comparative genomics has become widely accepted as the major framework for the ascertainment of functionally important regions in genomes. The underlying paradigm of this approach is that most of the functional regions are assumed to be under selective constraint, which in turn reduces the rate of evolution relative to neutrality. This assumption allows detection of functional regions through sequence conservation. However, constraint does not always lead to sequence conservation. When purifying selection is weak and mutation is biased, constrained regions can even evolve faster than neutral sequences and thus can appear to be under positive selection. Moreover, conservation estimates depend also on the orientation of selection relative to mutational biases and can vary over time. In the light of recent data of the ubiquity of mutational biases and weak selective forces, these effects should reduce the power of conservation analyses to define functional regions using comparative genomics data. We argue that the estimation of true mutational biases and the use of explicit evolutionary models are essential to improve methods inferring the action of natural selection and functionality in genome sequences.

FIG. 1.—

Accelerated rate of evolution when selection is counteracting a mutational bias. In the shown scenario, purifying selection favors C/G bases over A/T bases. Mutations occur according to a standard HKY85 model specified by the mutational bias, π_A+T and a transition/transversion ratio of four. Mutations G ↔ C and A ↔ T are neutral and mutational biases are symmetric (π_A = π_T, π_C = π_G). (A) Ratio of equilibrium substitution rate, r, over the equilibrium neutral rate, r₀, calculated according to equation (2). Values of r/r₀ < 1 indicate sequence conservation; values r/r₀ > 1 indicate rate acceleration. Rate acceleration is observed if weak purifying selection counteracts mutationally preferred states. (B) Comparison of mutational bias and actual composition bias for the selection coefficients that yield maximal rate acceleration for the respective mutational bias. Weak purifying selection counteracting the mutational biases effectively lowers the resulting composition bias compared with that expected if no selection were acting. Gray areas denote the respective regions where |Δρ_A+T| < |Δπ_A+T|. (C) Increase of divergence, d, over neutral divergence, d₀, for the selection coefficients that yield maximal rate acceleration for the respective mutational biases. Time is measured in units of the expected number of neutral substitutions per site.

FIG. 2.—

Performance of ML branch-length estimation in the presence of mutational biases and weak selective constraint using different neutral inference models. In all four evolutionary scenarios, the true mutation model is an HKY85 model specified by an equilibrium A/T content π_A+T = 0.8 at neutrally evolving sites and a transition/transversion ratio of four. In (A) and (C), C/G bases are preferred over A/T bases at constrained sites. Purifying selection is thus acting in opposition to the mutationally preferred states. In (B) and (D), A/T bases are preferred over C/G bases; purifying selection is acting in unison with the mutational biases. ML branch-length estimates at constrained sites, t^*, and at neutral sites, t₀^*, were calculated using three different neutral inference models: the Jukes-Cantor (JC69) model, the HKY85 model with its mutational biases estimated from the data (HKY85:ρ) and the HKY85 model using the true neutral mutation parameters (HKY85:π). Each used the default value of four for the transition/transversion ratio. Inferred branch-length ratios, t^*/t₀^*, are shown as a function of true divergence time t measured in units of the average number of substitutions per neutral site. Values t^*/t₀^* < 1 indicate sequence conservation, whereas t^*/t₀^* > 1 indicates faster than neutral evolution.

FIG. 3.—

Performance of ML branch-length estimation in the presence of mutational biases, randomly oriented weak constraint, and BGC under the HKY85:ρ and the HKY85:π inference models. The mutational bias is always π_A+T = 0.8 with a transition/transversion ratio of four. (A) and (C) show the results where C is the preferred base, whereas (B) and (D) show the results where A is the preferred base. Different colors indicate different strengths of BGC, which uniformly favors C/G over A/T alleles in both the test and the reference sequence.

FIG. 4.—

Performance of GERP on simulated sequence alignments over a realistic 32 mammalian species tree. The alignment sites were modeled to have evolved under mutational biases, randomly oriented weak constraint, and BGC. (A) Shows the mean RS-scores of all sites where C was the preferred state as a function of the strength of functional selection: γ_func. (B) Shows results for the sites where A was preferred state. The mutational bias was again π_A+T = 0.8 with a transition/transversion ratio of four. Positive RS-scores indicate branch-length reduction as the number of substitutions is lower than expected under neutrality (4.74)—the equivalent of t^*/t₀^*<1. Negative RS-scores indicate more substitutions having occurred than expected—he equivalent of t^*/t₀^*> 1. RS-scores were normalized such that the neutral class (γ_func = 0) has an RS-score of zero. Different colors indicate different strengths of BGC, which uniformly favors C/G over A/T alleles.