A New Comparative Framework for Estimating Selection on Synonymous Substitutions

Abstract

Selection on synonymous codon usage is a well-known and widespread phenomenon, yet existing models often do not account for it or its effect on synonymous substitution rates. In this article, we develop and expand the capabilities of multiclass synonymous substitution (MSS) models, which account for such selection by partitioning synonymous substitutions into 2 or more classes and estimating a relative substitution rate for each class, while accounting for important confounders like mutation bias. We identify extensive heterogeneity among relative synonymous substitution rates in an empirical dataset of ∼12,000 gene alignments from 12 Drosophila species. We validate model performance using data simulated under a forward population genetic simulation, demonstrating that MSS models are robust to model misspecification. MSS rates are significantly correlated with other covariates of selection on codon usage (population-level polymorphism data and tRNA abundance data), suggesting that models can detect weak signatures of selection on codon usage. With the MSS model, we can now study selection on synonymous substitutions in diverse taxa, independent of any a priori assumptions about the forces driving that selection.

Keywords: codon substitution models; codon usage bias; molecular evolution; synonymous substitutions.

Fig. 1.

Left: the average (per alignment) number of SNSS for each of the 67 codon pairs in the universal genetic code (gray dots) inferred from the Drosophila gene alignments. Plotted horizontal lines demarcate averages over individual codon pairs for each amino acid. Right (top): the distribution of the mean number of SNSS codon pairs per alignments versus the excess kurtosis shows that all of the MSS rates have a large number of observations in the tail end of the rate distributions. Right (bottom): the distribution of mean SNSS codon pairs per alignment versus the number of alignments where >0 of substitutions for the specific codon-pair are present.

Fig. 2.

Maximum likelihood point estimates for relative rates of each of the 67 pairs of synonymous substitutions in the SynREVCodon model. The rates are sorted by the median of the corresponding distribution (black tick). Rates that are exactly zero or those that are ≥3 are shown as circles, whose size and color reflect the number of alignments contributing to the corresponding estimates. The heatmap for each pair is generated using only the rates that are in (0,3).

Fig. 3.

Maximum likelihood point estimates for relative rates of each of the 18 amino acid level synonymous substitutions in the SynREV model. The rates are sorted by the median of the corresponding distribution (black tick). Rates that are exactly zero or those that are ≥3 are shown as circles, whose size and color reflect the number of alignments contributing to the corresponding estimates. The heatmap for each pair is generated using only the rates that are in (0,3).

Fig. 4.

Median (over genes) estimates of relative substitution rates from SynREV (red ticks) and SynREVCodon (blue dots) models. For amino acids with multiple pairs of synonymous substitutions, individual codon pairs are dispersed around the corresponding amino acid level estimates.

Fig. 5.

There is a positive correlation (Spearman's ρ = 0.48) between the median inferred synonymous substitution rates and the mean SNP density observed in a collection of 197 Drosophila genomes from Zambia. SNP density is an independent measure of selection pressure, as selection pressure to maintain a preferred codon decreases the SNP density for that codon. Correspondingly, the synonymous substitution rate for that codon will also be lower.

Fig. 6.

There is a negative correlation (rank-based Spearman's ρ = −0.41, P = 0.005) between the median inferred synonymous substitution rates (SynREVCodon), and the tRNA abundance ratio for the corresponding codon pair.

Fig. 7.

The relationship between ω (dN/dS) under MG94xREV, SynREV, and SynREVCodon models. Linear regression with the 0 intercept was run using the lm (y∼x + 0) function in R.

Fig. 8.

The relationship between ω under MG94xREV, SynREV, and SynREVCodon models for primate gene alignments. Linear regression with the 0 intercept was run using the lm (y∼x + θ) function in R.