Reversions mask the contribution of adaptive evolution in microbiomes

Abstract
Editor's evaluation
Introduction
Results
Discussion
Methods
Appendix 1
Data availability
References
Article and author information

Abstract

When examining bacterial genomes for evidence of past selection, the results depend heavily on the mutational distance between chosen genomes. Even within a bacterial species, genomes separated by larger mutational distances exhibit stronger evidence of purifying selection as assessed by d_N/d_S, the normalized ratio of nonsynonymous to synonymous mutations. Here, we show that the classical interpretation of this scale dependence, weak purifying selection, leads to problematic mutation accumulation when applied to available gut microbiome data. We propose an alternative, adaptive reversion model with opposite implications for dynamical intuition and applications of d_N/d_S. Reversions that occur and sweep within-host populations are nearly guaranteed in microbiomes due to large population sizes, short generation times, and variable environments. Using analytical and simulation approaches, we show that adaptive reversion can explain the d_N/d_S decay given only dozens of locally fluctuating selective pressures, which is realistic in the context of Bacteroides genomes. The success of the adaptive reversion model argues for interpreting low values of d_N/d_S obtained from long timescales with caution as they may emerge even when adaptive sweeps are frequent. Our work thus inverts the interpretation of an old observation in bacterial evolution, illustrates the potential of mutational reversions to shape genomic landscapes over time, and highlights the importance of studying bacterial genomic evolution on short timescales.

Editor's evaluation

This valuable study addresses the interpretation of patterns of synonymous and nonsynonymous diversity in microbial genomes. The authors present solid theoretical and computational evidence that adaptive mutations that revert the amino acids to an earlier state can significantly impact the observed ratios of synonymous and nonsynonymous mutations in human commensal bacteria. This article will be of interest to microbiologists with a background in evolution and to researchers studying the human microbiome.

https://doi.org/10.7554/eLife.93146.sa0

Introduction

Understanding evolutionary pressures acting upon bacterial populations is crucial for predicting the emergence and future virulence of pathogens (Culyba and Van Tyne, 2021), modeling strategies to combat antimicrobial resistance (Davies and Davies, 2010), and designing genetically modified organisms (Castle et al., 2021). Bacteria can adapt at rapid rates due to their short generation times and large population sizes. Indeed, the rapid evolutionary potential of the microbiome has been proposed to assist in the dietary transitions of mammals (Kolodny and Schulenburg, 2020). However, the vast majority of possible mutations do not increase bacterial fitness and instead result in a neutral or deleterious effect (Kimura, 1977; Davies et al., 1999; Jolley et al., 2000; Dingle et al., 2001). Metrics that estimate the directionality and intensity of past selection at genomic loci of interest have thus become critical tools in modern microbiology and biology more generally.

The normalized ratio of nonsynonymous (N) to synonymous (S) substitutions, known as d_N/d_S or the K_A/K_S ratio, is a widely used indicator of past selection (Jukes and Cantor, 1969; Kryazhimskiy and Plotkin, 2008; Barber and Elde, 2014). Nonsynonymous substitutions change the encoded amino acid and thus are considered likely to impact a protein’s function, while synonymous substitutions do not affect the encoded amino acid and are therefore considered effectively neutral, with limited exceptions (Nowick et al., 2019). To account for the fact that nonsynonymous mutations are more likely than synonymous mutations based on the genomic code (~3× on average; Yang and Nielsen, 2000), the values ‘d_N’ and ‘d_S’ normalize mutation counts to available sites on the genome. The d_N/d_S ratio therefore summarizes past selection on a genetic sequence, which could be a whole genome, pathway, gene, functional domain, or nucleotide; notably values of d_N/d_S averaged genome-wide can obscure signatures of adaptive evolution on other portions of the genome (Loo et al., 2020; Ho et al., 2005; Peterson and Masel, 2009). A d_N/d_S ratio of >1 indicates the dominance of past adaptive evolution (i.e., directional selection) while a ratio of <1 traditionally implies past selection against amino acid change (purifying selection).

Early sequencing work comparing bacterial genomes of the same species reported relatively low d_N/d_S values across the whole genome (<0.15) (Jolley et al., 2000; Dingle et al., 2001). These observations, obtained from comparing distant bacteria within each species, indicated a strong predominance of purifying selection. However, as it became economically feasible to sequence organisms separated by fewer mutations and therefore less evolutionary time, a contrasting pattern emerged in which high d_N/d_S values (~1) were found between closely related strains (Feil et al., 2003; Baker et al., 2004). Recent work in the human microbiome has confirmed such results and furthered the contrast between timescales by finding values of d_N/d_S > 1 (Garud et al., 2019; Lieberman et al., 2011; Shoemaker et al., 2022). The timescale dependence of d_N/d_S has been mainly attributed to the ongoing action of purifying selection (Garud et al., 2019), a model first proposed by Rocha et al., 2006. According to this model, weak purifying selection (or locally inactive purifying selection; Loo et al., 2020) allows for an initially inflated d_N/d_S ratio as deleterious mutations that will eventually be purged remain in the population. As time progresses and purifying selection continuously operates, the d_N/d_S ratio decreases (Loo et al., 2020; Garud et al., 2019; Rocha et al., 2006). However, multiple studies have observed genome-wide values of d_N/d_S > 1 in these same microbial systems, with values substantially >1 in key genes, which are simply unaccounted for in the purifying model (Garud et al., 2019; Lieberman et al., 2011; Marvig et al., 2015; Zhao et al., 2019; Zhao et al., 2020).

Here, we demonstrate fundamental flaws in the purifying selection model in the context of the large within-person population sizes typical to the human microbiome and many bacterial infections (>10¹² bacteria/person). We use analytical, simulation-based, and genomic approaches to support a contrasting model for the timescale dependence of d_N/d_S, in which adaptive evolution predominates but is not apparent on long-timescales due to adaptive reversion. The comparative success of the reversion model suggests that the study of closely related bacteria is needed to fully understand evolutionary dynamics.

Results

A model of purifying selection that fits the data reveals unrealistic parameters

Explaining the timescale dependence of d_N/d_S through an exclusively purifying selection model poses several challenges. Firstly, fitting observed data with purifying selection requires a preponderance of mutations with extraordinarily small effects on fitness (selective coefficients, s), which are challenging to eliminate effectively (Haigh, 1978). Secondly, the occurrence of an adaptive event during the extensive time required to purge weakly deleterious mutations interrupts the purification of such mutations. Lastly, neutral bottlenecking processes, such as those observed during host-to-host transmission, exacerbate the accumulation of deleterious mutations. For most of this section, we will disregard these last two complications and focus on the problem of small s. To provide clarity, we first detail the classic purifying selection model.

Mutations can be divided into three classes, the first two of which accumulate at a constant rate per unit of time: synonymous mutations ( $S$ ), neutral nonsynonymous mutations ( $N_{n e u t}$ ), and non-neutral, transient, nonsynonymous mutations ( $N_{t r a n s i e n t}$ ). We restate the timescale dependence of d_N/d_S as the observation that, in a population starting from a single wild type (WT) cell, the average number of non-neutral nonsynonymous mutations per cell in the population ( ${\bar{N}}_{t r a n s i e n t}$ ) increases and then asymptotes. Assuming an infinitely large population size and an infinite genome size (to circumvent saturation of mutations), the exclusive purifying selection model (Garud et al., 2019; Rocha et al., 2006) can thus be written as

\frac{d_{N}}{d_{S}} = \frac{{\bar{N}}_{n e u t} + {\bar{N}}_{t r a n s i e n t}}{3 \bar{S}}

and

d {\bar{N}}_{t r a n s i e n t} = U_{N} d t - s {\bar{N}}_{t r a n s i e n t} d t .

Here, $U_{N}$ is the non-neutral mutation rate per core genome per generation, $s$ is the selective disadvantage of a non-neutral nonsynonymous mutation (or the harmonic mean of such mutations; see Appendix 1, Section 1.1), and $t$ is the number of generations. The 3 in the denominator of Equation 1 accounts for the discrepancy in the number of potential nonsynonymous and synonymous sites (Yang and Nielsen, 2000). We solve for ${\bar{N}}_{t r a n s i e n t}$ by assuming ${\bar{N}}_{t r a n s i e n t} (t = 0) = 0$ to obtain:

{\bar{N}}_{t r a n s i e n t} (t) = \frac{U_{N} (1 - e^{- s t})}{s} .

We further simplify and combine these equations to create an equation for d_N/d_S with only two parameters as previously done (Garud et al., 2019). First, since d_N/d_S plateaus with time (Figure 1a), we have $\underset{t \to \infty}{l i m} \frac{{\bar{N}}_{n e u t} + {\bar{N}}_{t r a n s i e n t}}{3 \bar{S}} = \frac{{\bar{N}}_{n e u t}}{3 \bar{S}} = α$ . Conveniently, $α$ represents both the asymptote of d_N/d_S and the proportion of nonsynonymous mutations that are neutral. This allows us to leave only $s$ as the other free parameter, obtaining (see Appendix 1, Section 1.1)

\frac{d_{N}}{d_{S}} = α + (1 - α) \frac{1 - e^{- s t}}{s t} .

Figure 1 with 3 supplements see all

Download asset Open asset

The previously proposed explanation for the time dependence of d_N/d_S is weak purifying selection.

(a) Time signature of d_N/d_S as depicted by data points derived from the studies by Garud et al., 2019 and Zhao et al., 2019. Each dot represents a pairwise comparison between the consensus sequence from two gut microbiomes as computed by Garud et al., 2019, using only the top 10 species based on the quality of data points (see ‘Methods’). Where the high initial value of d_N/d_S begins to become the low asymptotic value of d_N/d_S occurs at approximately $d_{S} = \frac{μ_{S}}{s}$ . Fit lines were derived from these points using Equation 4 to depict the trend. The median R² is 0.81 (range 0.54–0.94). Corresponding data from Zhao et al., 2019 confirms these observed trends, demonstrating high levels of d_N/d_S at short timescales and low levels at longer timescales. Adaptive genes are from Zhao et al., 2019 and are defined as those that have high d_N/d_S values in multiple lineages. Insets: d_N vs. d_S on a linear scale. Note that the data was fit to minimize variance in the logarithmic scale, not the linear scale, so the fit is not expected to be as good for the inset. See Figure 1—figure supplement 1 for minimizing variance on a linear scale. See Figure 1—figure supplement 2 for all species on separate panels. (b) Values of s from the output of 999 standard bootstrap iterations of curve fitting, conducted with replacement, demonstrate that only small values of the average selective coefficient can fit the data.

As sequence analysis is not privy to the actual number of generations, we approximate $t$ assuming that synonymous mutations accumulate according to a molecular clock $(t = \frac{d_{S}}{2 (1 / 4) μ})$ , where μ is the mutation rate per generation per base pair, ¼ represents the proportion of random mutations that are synonymous (Yang and Nielsen, 2000), and 2 accounts for the fact that divergence is a measure between a pair of genomes. As selection and mutation are both in units per time, any change in μ results in a corresponding change in s. Both model fits and consequences are largely dependent on the ratio of these two variables (more on this below), and thus are not sensitive to the choice of μ. We use a relatively high mutation rate of 10^-9 per base pair per generation (Drake, 1991; Barrick and Lenski, 2013) as lower rate would imply even weaker purifying selection.

Fitting the data from Garud et al., 2019, we infer median values of α ≈ 0.10 (0.09–0.14) and s ≈ 3.5 × 10^-5 (2.6 × 10^-5-6.5 × 10^–5) across all species (‘Methods’, Figure 1a). Aggregating all of the data at once results in a similar optimal fit of α ≈ 0.11 and s ≈ 2.8 × 10^–5. The similarity across the 10 species is perhaps not surprising, given that all are human gut residents of the order Bacteroidales; these values are also in line with the values obtained previously from aggregating across all species (Garud et al., 2019). These values indicate a model in which only ~10% of nonsynonymous mutations are neutral and the remaining ~90% are so weakly deleterious that they are beyond the limit of detection of any experimental method to date (s ≳ 10^–3) (Gallet et al., 2012). Higher values of s that better reflect experimental observations (Kibota and Lynch, 1996; Trindade et al., 2010; Robert et al., 2018) result in poor fits to the data (Figure 1—figure supplement 3). While the implied proportion of deleterious mutations may seem high, deep mutational scanning experiments have revealed that most amino acid-changing mutations in essential genes are deleterious enough to be measured in the lab (Kelsic et al., 2016; Dewachter et al., 2023); complex real-world environments are expected to constrain an even larger fraction of the genome.

In finite populations, the presence of so many weakly deleterious mutations becomes quickly problematic. When s is smaller than $U_{N}$ , organisms without any deleterious mutations (or with the fewest number of deleterious mutations, the ‘least-loaded class’; Haigh, 1978) can be easily lost from a finite population before they outcompete less fit organisms and fitness decay begins to occur. The likelihood of loss depends on the population size and mutation-selection balance ( $U_{N} / s$ ), a parameter that estimates the average number of deleterious mutations per cell relative to the least-loaded class. Given a core genome of L = 1.5 ×10⁶ bp that can acquire deleterious mutations, we then expect 0.001 new deleterious mutations per genome per generation ( $U_{N} = \frac{3}{4} (1 - α) L μ$ ). Thus, the value of $U_{N} / s$ for the above fits is ~29, indicating that most cells in the population contain dozens of deleterious mutations (see Appendix 1, Section 1.2). With this value of the mutation-selection balance parameter, the frequency of mutation-free organisms in a population is extremely small, even for a population that starts without any deleterious mutations (<10^–12 after 100,000 generations). If the flexible genome also contains deleterious mutations, the least-loaded class is pushed down even further. Simulations substantiate this prediction of mutation accumulation and decrease in frequency of the wild type (Figure 2a, ‘Methods’).

Figure 2 with 4 supplements see all

Download asset Open asset

Models of extremely weak purifying selection that can fit the data suffer from mutation accumulation and fitness decay.

(a) The temporal dynamics of the least-loaded class in a large population under the purifying selection model. The black line represents the predicted frequency of the wild-type (mutation-free) class over time. The simulation curve shows simulation results assuming constant purifying selection in an exceptionally large effective population size (N_e = 10¹⁸; see text for a discussion of population size) under a slightly modified Wright–Fisher model (‘Methods’). (b) As a consequence of the loss of the least-loaded class, fitness declines in finite populations over time. Colored lines indicate simulations from various effective population sizes with mutations of constant selective effect. The deleterious mutation rate in the simulation is 1.01 × 10⁻³ per genome per generation. (c) Using the same simulations as in panel (b), we see that realistic global effective population sizes fail to fit the d_N/d_S curve, with different asymptotes. The black line denotes the infinite population theoretical model, and the colored lines indicate increasing effective population sizes, which change the strength of genetic drift in the simulations. Larger values of s and models in which all mutations are deleterious cannot fit the data (Figure 2—figure supplement 1, Figure 2—figure supplement 2). Generations are assumed to occur once every day.

The time until the least-loaded class is completely lost from the population depends on the strength of genetic drift. The strength of genetic drift is inversely proportional to population size in well-mixed populations (Gillespie, 2004), and in less well-mixed or otherwise nonideal populations, is inversely proportional to a smaller parameter, the effective population size, N_e. N_e is often estimated by assessing polymorphisms in a population (Gillespie, 2004) but is hard to estimate from data following a recent bottleneck. Because each individual’s gut microbiome is thought to be well mixed (census size = 10¹³) (Sender et al., 2016), it has been recently argued that N_e ≈ 10¹¹ reflects drift processes for dominant gut species (Ghosh and Good, 2022; Labavić et al., 2022). On the other hand, lower values of N_e ≈ 10⁹ or less have been estimated for global populations of bacteria (Bobay and Ochman, 2018) because of the slow rates of bacterial transmission across people. While this decrease in N_e when increasing scales may seem paradoxical, we note this use of N_e only reflects the magnitude of the force of drift; for other calculations in nonideal populations, census population size or other parameters should be used.

Without extremely large values of N_e, the least-loaded class will be lost recurrently, rapidly lowering the fitness of the population (i.e., Muller’s ratchet; Haigh, 1978). Assuming $s$ is small and thus approximately additive, this recurrent process of fitness decay occurs roughly when the following inequality is satisfied (see Appendix 1, Section 1.2; Neher and Shraiman, 2012):

2 s N_{e} e^{- \frac{U_{N}}{s}} << 1.

Given $U_{N} / s$ = 29 as derived above, N_e > 10¹⁵ is required to avoid continual deleterious mutation accumulation and fitness decline (Figure 2b). Thus, the purifying model requires levels of drift unrealistic at the within-person or across-globe scales. Simulations confirm that deleterious mutations accumulate and compromise the ability of the purifying model to explain empirical d_N/d_S decay in reasonably finite populations (Figure 2c). Moreover, continuous accumulation of mutations in such populations decreases fitness so much that the average genome contains a sizable fraction (~10%) of deleterious alleles after 1 million years (Figure 2c), assuming N_e = 10⁹ and one generation a day (Korem et al., 2015). Even if this decreased fitness was biologically maintainable, the accumulation of so many deleterious mutations would lead to many potential adaptive back mutations, complicating the efficiency of purifying selection. Consequently, this value of $U_{N} / s$ is simply incompatible with a model where a vast majority of alleles are already optimal.

Lastly, the intolerance of the purifying model to adaptation and transmission is particularly problematic. Within-host adaptive sweeps have been observed in Bacteroides fragilis (Zhao et al., 2019) and other Bacteroides (Garud et al., 2019). Such adaptation interferes with inefficient purifying selection; deleterious mutations are likely to hitchhike (Desai et al., 2013) to fixation on the genomic background of adaptive mutations. Any given weakly deleterious mutation with s = 3.5 × 10^–5 cannot be purged from a within-host population on the timescale of human lifetime (assuming ~1 generation per day), and thus if any adaptive sweep occurred within that host, it would either hitchhike to fixation or be completely removed from the population. Similarly, deleterious mutations can also hitchhike to fixation during neutral transmission bottlenecks, thereby raising the average number of deleterious mutations per cell in the population, furthering mutation accumulation, and hampering the efficiency of purifying selection. Simulations confirm that even infrequent adaptive sweeps and bottlenecks have tangible impacts on d_N/d_S, including raising the asymptote (Figure 2—figure supplement 3).

Neither recombination nor differential selection at transmission can easily rescue a model of weak purifying selection

Homologous recombination, which occurs at detectable rates within human gut microbiomes and within the Bacteroidales order (Liu and Good, 2024), cannot rescue a population from Muller’s ratchet when such weakly deleterious mutations are so frequent. If we assume a generously high rate of recombination, such that a mutated nucleotide is 500 times more likely to be reverted via recombination than mutation (r/m = 500) (Torrance et al., 2024; Liu and Good, 2024) and brings along a single linked synonymous mutation during each recombination event, the decay of d_N/d_S still cannot be recreated in a population of size 10⁹ and fitness will still decay (Figure 2—figure supplement 4). The inability of recombination to suppress mutation accumulation in this regime arises because the selective advantages themselves are still too small to sweep faster than the rate at which mutations accumulate. While recombination does allow d_N/d_S to eventually decay, the rate of decay is much slower than observed, resulting in a poor fit to the data (Figure 2—figure supplement 4). While higher values of N_e or a higher recombination rate could theoretically approximate the absence of linkage and escape of Muller’s ratchet, we note that the maximum r/m across bacteria is estimated to be <50 (Torrance et al., 2024) and our simulations are therefore conservative.

Our presentation so far has implicitly assumed that weak purifying selection has been acting continuously and that values of s are constant for any given allele over time. However, apparently weak purifying selection might theoretically emerge from mutations that spend periods under neutral selection (or even local positive selection) and larger periods under strong negative selection, with the estimated value of s reflecting the harmonic mean (Culyba and Van Tyne, 2021; Loo et al., 2020). However, such models will have a hard time overcoming mutation accumulation. For example, a model in which purifying selection acts only during transmission still cannot prevent mutation accumulation without unrealistic assumptions. In particular, the selection-at-transmission model would still require ~29 non-neutral mutations in the average adult population, which implies a very low frequency of the least-loaded class. Assuming each host’s population gets replaced once every 10,000 bacterial generations (~26 years), such a model would require the least-loaded class to be 6000× more likely to colonize them than the average genotype in the population ( ${(1 + 10,000 s)}^{29}$ ). The presence of rare cells with strong selective advantages would suggest super-spreading across human microbiomes, which has yet to be reported in the human microbiome (Faith et al., 2013). More importantly, Muller’s ratchet would still click because of the low frequency of this least-loaded class.

Adaptive reversions can explain the decay of d_N/d_S

If purifying selection cannot explain the decay in neutral mutations, what can? One particularly attractive process that removes nonsynonymous mutations over time is strong adaptive mutation and subsequent strong adaptive reversion of the same nucleotide when conditions change. Such reversions are likely to sweep in large populations when mutations are adaptive locally but deleterious in other environments (Ascensao et al., 2023). In the gut microbiome, these alternative environments could represent different hosts (Figure 3a) or environmental changes within a single host (e.g., diet, medication, other microbes). As an illustrative example, the presence of a bacteriophage in one gut microbiome might select for a loss-of-function mutation (premature stop codon or otherwise) in a phage receptor, driving this mutation to fixation in its host, but reverting to the wild-type receptor when transmitted to a phage-free host. Reversions are most likely when compensatory mutations that counteract a mutation’s deleterious effects are either scarce or not as beneficial as direct reversion (Levin et al., 2000) (i.e., provided a premature stop codon); we discuss models that include compensatory mutations later in this section.

Figure 3 with 2 supplements see all

Download asset Open asset

Locally adaptive mutations and subsequent reversions can explain the decay of nonsynonymous mutations.

(a) Cartoon schematic depicting a potential reversion event within a single transmitted lineage of bacteria. The color of each individual indicates a different local adaptive pressure. Closed arrows represent mutation while open arrows indicate transmission. (b) Reversions become increasingly likely at larger population sizes and are nearly guaranteed to occur and fix within 1–10 years when strongly beneficial in gut microbiomes. The probability of revertant arising and fixing (top panel) is calculated as $1 - (1 - 2 s_{b e n} U_{r e v})^{N_{e}}$ , and the expected time to fixation of reversion (bottom panel) is calculated as $\frac{1}{2 N_{e} s_{b e n} U_{r e v}} + \frac{l n (N_{e})}{s_{b e n}}$ ( $U_{r e v} = 4.5 \times 10^{- 10}$ per generation). Generation times are assumed to be 1 day. Note that the mutation rate does not affect time to fixation much when *N_e* is large. Here, we assume no clonal interference or bottlenecks, though simulations do take these processes into account. See Appendix 1, Section 2.2 for derivation. Each line type displays a different selective advantage coefficient. (c) The adaptive reversion model can fit the data. Each colored line shows the fit for a different species. The median R² = 0.82 (range 0.54–0.94). Fit minimizes logarithmic variance. See Figure 3—figure supplement 1 for alternative fitting linear variance. See Figure 3—figure supplement 2 for species individually. Insets: fit parameters for $τ_{f l i p}$ , the average number of generations for *a given* environmental pressure to switch directions and $n_{l o c i}$ , the average number of sites under different fluctuating environmental pressures. The scale of the y-axis is linear. Confidence intervals are from 999 bootstrapped resamples.

Adaptive nonsense mutations have been observed to emerge frequently within individual people in both pathogens (Culyba and Van Tyne, 2021; Lieberman et al., 2011; Key et al., 2023; Shopsin et al., 2008) and commensals (Zhao et al., 2019; Barreto et al., 2023). Identifying reversions in vivo requires both high temporal resolution and deep surveillance such that the probability of persistence of ancestral genotype is removed (Snitkin et al., 2013) despite this difficulty, reversions of stop codons have been observed in mouse models (Sousa et al., 2017) and during an outbreak of a pathogen infecting the lungs of people with cystic fibrosis (Poret et al., 2024). While direct reversion has not yet been observed in gut microbiomes, premature stop codons are frequently observed. Among the 325 observed nonsynonymous de novo mutations in a study of within-host B. fragilis adaptation (Zhao et al., 2019), 28 were premature stop codons. This frequency is significantly higher than expected by chance (p=0.015; ‘Methods’). Moreover, 4 of the 44 mutations in 16 genes shown to be under adaptive evolution on this short timescale were stop codons. These same 16 genes show a signature of purifying selection on long timescales (Figure 1a).

Traditionally, mutational reversions of stop codons and other mutations have been considered exceedingly unlikely and have been ignored in population genetics (Tajima, 1996), with a few exceptions (Charlesworth and Eyre-Walker, 2007). However, for a bacterial population within a human gut microbiome, the likelihood of a mutational reversion is quite high. A single species within the gut microbiome can have a census population size of 10¹³, with generation rates ranging from 1 to 10 per day (Sender et al., 2016; Korem et al., 2015). Taking a conservative estimate of one generation per day and a within-person N_e of 10¹⁰ (e.g., bacteria at the end of the colon may not contribute much to the next generation; Labavić et al., 2022), reversions become highly probable (Figure 3b; see Appendix 1, Section 2.2). Given a mutation rate of 10^–9 per site per generation, we anticipate 10 mutants at any given site each generation. In the large population sizes relevant for the gut microbiome, a beneficial mutation will then take substantially longer to sweep the population than occur, with values of s_ben > 1% generally sweeping within 10 years (Figure 3b). Consequently, if selection strongly benefits a reverting mutation, a genotype with a beneficial mutation is essentially guaranteed to emerge within days to weeks and replace its ancestors within the host within months to years.

Given its plausibility, we now consider if the reversion model can explain the observed decay of d_N/d_S. The dynamics of the reversion model can be given by

\frac{d {\bar{N}}_{t r a n s i e n t}}{d t} = \frac{1}{τ_{f l i p}} (n_{l o c i} - {\bar{N}}_{t r a n s i e n t}) - \frac{1}{τ_{f l i p}} {\bar{N}}_{t r a n s i e n t}

With the corresponding solution for ${\bar{N}}_{t r a n s i e n t}$ being (see Appendix 1, Section 3.1)

{\bar{N}}_{t r a n s i e n t} (t) = \frac{n_{l o c i}}{2} (1 - e^{- \frac{2 t}{τ_{f l i p}}}) .

Here, $n_{l o c i}$ denotes the number of loci that experience distinct sources of fluctuating selection. The parameter $τ_{f l i p}$ represents the average number of generations required for the sign of selection at a chosen locus to flip and determines the key point in the d_N/d_S decay curve where ${\bar{N}}_{t r a n s i e n t} (t)$ begins to drop. We note that a locus here could be a nucleotide, gene, or gene set – any contiguous or noncontiguous stretch of DNA in which two knockout mutations would be just as beneficial or harmful as one mutation. We again use α to represent the proportion of nonsynonymous mutations that are neutral. Using Equation 7, we obtain a formula for d_N/d_S that has only three free parameters when a single value for $μ_{S}$ is chosen:

\frac{d_{N}}{d_{S}} = α + \frac{n_{l o c i}}{6 μ_{S} t} (1 - e^{- \frac{2 t}{τ_{f l i p}}}) .

When fitting the d_N/d_S curve, the values obtained are reasonable in the context of bacterial genomics, with median best-fit values across species of $τ_{f l i p} = 46, 000$ bacterial generations (range 25,000–105,000) and $n_{l o c i} = 55$ (range 34–80). Given daily bacterial generations, this value of $τ_{f l i p}$ suggests the sign of selection on a given allele would flip approximately every 110 years. The average time for any pressure to flip would thus be approximately every 2 years, or less frequently if adaptive events occur in bursts (e.g., upon transmission to a new host). While 55 loci under distinct selective pressures may seem high, Bacteroidetes genomes are known to have dozens of invertible promoters (up to 47 in B. fragilis; Jiang et al., 2019). Invertible promoters are restricted out of the genome and re-ligated in the opposite direction to turn gene expression on or off. The number of invertible promoters in a given genome approximates a lower bound on the number of fluctuating selective pressures that these genomes frequently experience. Interestingly, adaptive loss-of-function mutations reported in B. fragilis affect the same genes regulated by invertible promoters (Zhao et al., 2019). The plausibility of these fit parameters lends support to a model in which d_N/d_S decays solely based on strong and recurrent local adaptations.

To ensure a reversion model is robust to finite populations, we performed simulations using fit parameters. These simulations capture the dynamics of a single population evolving as it transmits across a series of hosts through random bottlenecks (Figure 4a; ‘Methods’); these simulations allow for clonal interference between adaptive mutations. We allow new pressures to arise independent of bottlenecks as new selective forces (phage migration [Koskella and Brockhurst, 2014]; immune pressures [Barroso-Batista et al., 2015]; dietary changes [Carmody et al., 2019]) can emerge throughout the lifespan and independent of migration; forcing transmission and bottlenecks to coincide gives similar results (Figure 4—figure supplement 1). As in the purifying selection simulations, the per base pair mutation rate is 10^–9, and 90% of nonsynonymous substitutions are deleterious, but this time they have a larger s of 0.003 (Robert et al., 2018) and are thus purged more quickly from the population. Notably, while some of these deleterious mutations hitchhike to fixation during bottlenecks and adaptive sweeps, fitness does not decay because these mutations are subsequently reverted with adaptive sweeps (Figure 4—figure supplement 2). If deleterious mutations had significantly smaller s, they would be unable to be reverted due to the long time needed to reach fixation, even if bottlenecks and adaptive events are less frequent (Figure 2—figure supplement 3).

Figure 4 with 3 supplements see all

Download asset Open asset

Under a model of reversion, the apparent d_N/d_S on long timescales underestimates the extent of adaptive evolution.

(a) The reversion model successfully fits the data in simulations. We simulate a population of size 10¹⁰ that has a bottleneck to size 10 on average every 10,000 generations (~27 years or a human generation [Wang et al., 2023] given a bacterial generation a day), with one adaptive pressure ( $s_{b e n} = 0.03$ ) occurring on average every 840 generations independently of bottlenecks (see Figure 4—figure supplement 1 for an alternative where bottlenecks and selection are correlated). New pressures either require forward mutations (which can be acquired at a rate of 1.1 × 10^-8 per available locus per generation) or reversions (which can be acquired at a rate of 4.5 × 10^-10 per available locus per generation), the balance of which depends on the history of pressures on the tracked genome (i.e., more past forward pressures implies more potential future reverse pressures). Based on the best fit to the data, we use $n_{l o c i} = 55$ . Deleterious mutations occur at a rate of 1.01 × 10^-3 mutations per genome and have $s = 0.003$ and can themselves be reverted. More details on the simulation can be found in ‘Methods’. Each curve represents the average of 10 runs; the blue line shows the observed pairwise d_N/d_S while the red line includes adaptive mutations and reversions. The theory line is the result of Equation 8. Observable d_N/d_S decays because of reversion, while the actual d_N/d_S of mutations that occurred is >1 when taking into account both forward and reverse mutations. (b) PAML (Yang, 2007) cannot detect true d_N/d_S in a given gene in the presence of adaptive reversions. Both lines are d_N/d_S as calculated by PAML on a simulated gene phylogeny. In the permanent adaptations simulation (blue), adaptive mutations are acquired simply and permanently. In the transient adaptations simulation (red), only more recent mutations will be visible while older mutations are obscured (‘Methods’). Line is the average of 10 simulated phylogenies and shaded regions show the range. (c) Categories of genes in the *Bacteroides fragilis* genome (NCTC_9343) enriched for stop-codon adjacent codons (TTA and TCA) relative to the expectation from the rest of the genome (‘Methods’). The use of these codons suggests these sequences may have recently had premature stop codon mutations. p-Values are displayed above bars and were calculated using a one-proportion Z-test with Bonferroni correction. See ‘Methods’, Supplementary file 1, and Supplementary file 2 for more details.

We note that other complex models that include reversion and other processes are also possible. For example, a model with a very large number of loci with selective tradeoffs and pressures that act only transiently (nonfluctuating) could potentially fit the data. However, the agreement between $n_{l o c i}$ and the number of invertible promoters, and the finding of parallel evolution in vivo, suggests the fluctuating selection model is more realistic than a very many-sites model.

So far, we have assumed that only exact reversions are selected upon when the sign of selection returns to its original state. However, the reversion model can also accommodate compensatory mutations that exclude any selective advantage for reversion; these compensatory mutations can also be subject to reversion themselves. We conceptualize this as a random walk, in which a locus at a nonancestral state acquires a compensatory mutation with probability p or obtains a true reversion with probability 1 – p (see Appendix 1, Section 3.2). As long as $p \leq 0.5$ , d_N/d_S will decay to the same asymptote despite adaptive dynamics occurring. While compensatory mutations shift the timing of d_N/d_S decay to the right, it can be shifted backward by decreasing $n_{l o c i}$ (Figure 4—figure supplement 3a). The condition $p \leq 0.5$ is easily met when $s_{b e n}$ = 0.03, until excluding compensatory mutations are 10 times more likely than true reversion and provide 95% of the selective advantage of the true reversion (Figure 4—figure supplement 3b). If selective pressures are stronger (as they might be in the presence of phage), true reversions will outcompete compensatory mutations even if the supply of compensatory mutations is greater or such mutations provide better relative compensation.

A critical consequence of the reversion model is that apparent and actual d_N/d_S values diverge quickly. Even when the true genome-wide d_N/d_S exceeds 1 – meaning that adaptive sweeps have been a dominant force in shaping genomes – the observed value can be close to 0.1 on long timescales. This disparity complicates the interpretation of d_N/d_S as it becomes challenging to determine whether a genome or gene lacks nonsynonymous mutations due to reversions or negative selection. We confirmed the inability to detect adaptive selection on a gene when reversion is rampant by simulating protein phylogenies; even the advanced software PAML (Phylogenetic Analysis by Maximum Likelihood) (Yang, 2007) significantly underestimates actual d_N/d_S (Figure 4b; ‘Methods’). Without sufficient temporal sampling, no software can realistically estimate these hidden, adaptively driven nonsynonymous mutations.

Lastly, we sought to find evidence of past reversions of stop codons in certain genes by analyzing codon usage. Both leucine and serine have the property that they can be encoded by six codons, only one of which is highly stop codon adjacent (TTA for leucine and TCA for serine). Across the B. fragilis genome, these codons are depleted overall (13.48% usage rather than the neutral expectation of 16.67%). However, specific Clusters of Orthologous Genes (COG) categories are enriched in TTA and TCA codons relative to this baseline, including genes associated with transcription and cell envelope biogenesis (Figure 4c, ‘Methods’, Supplementary file 1). Further, when functionally annotated genes are further categorized by cellular localization, more gene categories exhibit enrichment (Figure 4c, ‘Methods’, Supplementary file 1), most notably genes involved in inorganic ion transport and metabolism that are localized to the outer membrane. Genes implicated in within-host B. fragilis adaptation (Zhao et al., 2019) are also found disproportionately in this category of outer membrane transporters (p=1.22 × 10^–4; ‘Methods’, Supplementary file 1). Both the cell envelope and membrane-bound transporters are known to mediate interactions with the immune system and phage (Sukhithasri et al., 2013; Ongenae et al., 2022), and are therefore expected to experience fluctuating selective pressures. The enrichment of stop-codon adjacent codons in pathways associated with environment-dependent costs further supports a model in which adaptive mutational reversions are frequent.

Discussion

In this study, we present a new interpretation of the time-dependent changes in d_N/d_S for bacterial populations. We show that the traditional weak purifying selection model struggles to replicate theoretical results in realistic population sizes and propose an alternative model with opposite implications that are supported by analytical, simulation, and genomic results. Together, these results challenge the conventional view that high d_N/d_S values on short timescales are an artifact and should not be trusted. Instead, the success of the reversion model suggests that adaptive dynamics are underestimated on long timescales because of the saturation of d_N.

It is perhaps not surprising that reversions have been relatively overlooked in previous literature. First, most population genetics theory focuses on eukaryotic organisms with smaller population sizes and longer generation times, for which reversion is less likely. The low likelihood of reversion in these populations has inspired the use of the convenient infinite-site model (Tajima, 1996), which assumes that reversions never occur and simplifies derivations. While smaller values of N_e can be appropriate for modeling global bacterial dynamics – because bottlenecks and geography limit how many organisms effectively compete – they are inappropriate for within-gut populations, which are less structured. While gut microbiomes do have a spatial structure that reduces competition, theoretical work modeling this biogeography suggests that the census and active population sizes differ only approximately tenfold (Ghosh and Good, 2022; Labavić et al., 2022). This brings the within-gut microbiome N_e to substantially larger than the per-nucleotide mutation rate, invalidating the infinite sites model. Secondly, while bacterial geneticists have long observed adaptive loss-of-function mutations, two common misinterpretations of population genetic parameters can underestimate the probability of reversion: molecular clock rates (μ), which are generally low, can easily be confused with the supply of potential mutations (μN_e) (Lieberman, 2022); and classical approaches that assess N_e from genetic diversity vastly underestimate the currently active population size, particularly if a bottleneck recently occurred (e.g., during transmission). Lastly, simulating large populations, even when appropriate, is computationally difficult. As a consequence, population genetics simulations, including those of bacteria, have used relatively small population sizes (≤10⁶ organisms). We overcome computational limitations by tracking genetic classes rather than individual genotypes (‘Methods’). While our approach does not allow explicit comparison between individuals within a population, we believe this framework represents a powerful method to simulate large population sizes when applicable.

Whether or not a reversion model can be applied beyond host-associated microbial populations remains to be explored. We only analyze microbiome data here, but we anticipate that analyses of highly curated d_N/d_S decay curves from microbial pathogens could yield similarly plausible parameter fits for the reversion model given past observations of d_N/d_S decay (Rocha et al., 2006). When effective population sizes are smaller than 10⁹, reversions are relatively unlikely. For example, while adaptive reversions can sweep individual gut microbiomes, we do not propose that reversions sweep the global bacterial population. Regardless, theoretical work on animal populations has shown that adaptive reversions are possible after local population bottlenecks (Charlesworth and Eyre-Walker, 2007). Similarly, environmental variations that change more rapidly than the timescale required for a local selective sweep (e.g., those imposed by daily dietary changes in the gut; or imposed by light-dark cycles in the environment) would be less likely to drive fixation and subsequent reversion than the less rapid changes considered here (e.g., phage migration) (Cvijović et al., 2015). On the other hand, adaptive reversions may be particularly relevant for viral populations, which are known to undergo within-host adaptation, have very large population sizes, and experience frequent bottlenecks (Feder et al., 2017). Reversions have commonly been observed in certain regions of the HIV genome and have been postulated to diminish measured substitution rates in those regions (Druelle and Neher, 2023).

Despite the success of a model of reversion alone in explaining d_N/d_S decay, it remains possible that other forces could additionally contribute. While we have shown that purifying selection alone, either continuously or during transmission, cannot explain d_N/d_S decay alone, it is possible that some degree of purifying selection could act alongside a reversion model. Similarly, directional selection could be incorporated into the reversion model by adjusting the parameter $α$ . While the true contribution of adaptive evolution to $α$ is likely nonzero, it is difficult to fit with available data and it is therefore left for future work.

While we have presented evidence that recombination alone is unlikely to rescue a model of weak purifying selection, it remains possible that recombination could be included in a model that includes adaptation and, notably, could drive adaptive reversions. Microbial geneticists have frequently observed that recombined regions exhibit lower d_N/d_S values compared to non-recombined regions (Castillo-Ramírez et al., 2011), a signature consistent with having already experienced reversion or purifying selection. Recombination could potentially revert multiple mutations at specific loci simultaneously, which might be particularly beneficial in the presence of genomic epistasis. Thus, despite the success of the mutation-driven model, it is likely that recombination plays some role in the decay of d_N/d_S.

While more direct observation of adaptive reversions is currently lacking, we propose that this paucity is simply an artifact of lacking samples along a line of descent with sufficient genomic resolution. Despite this challenge in observation, a recent study tracking de novo mutations between mothers and infants revealed several cases of apparent reversion, with elevated values of d_N/d_S > 1, though not significantly so (Chen and Garud, 2022). Moreover, many short-term studies in the gut microbiome and beyond have revealed strong evidence of within-person adaptation, including parallel evolution (Lieberman et al., 2011; Marvig et al., 2015; Zhao et al., 2019; Cooper and Lenski, 2000) and loss-of-function changes like premature stop codons (Key et al., 2023) – with low long-term d_N/d_S values in these same short-term genes (Vigué and Tenaillon, 2023). We note that adaptation and reversion do not result in parallel evolution in the genomic record if various initial mutations result in the same phenotype (i.e., loss-of-function mutations); however, it would result in changes in codon usage bias we have shown (Figure 4c).

The shortcomings of the purifying model and the success of the reversion model under realistic assumptions highlight the importance of studying evolution in real time for understanding evolutionary dynamics. In addition, our results emphasize the importance of simulating large population sizes for explaining observations in bacterial population genomics, spotlight the potential for strong adaptation in bacterial populations, and underscore the need for continued development of population genetics theory for microbial populations.

Methods

Data and parameter estimation

Data was obtained from Shoemaker et al., 2022 and was initially generated by Garud et al., 2019. Pairwise d_N/d_S values can be found in the GitHub repository. The parameters are estimated using scipy.optimize.curve_fit. The fit minimizes the RMSD of the logarithmic d_N/d_S. If we fit the data by minimizing just d_N/d_S on a linear scale, we get s ≈ 2.0 × 10^–5, which suggests an even weaker purifying selection (Figure 1—figure supplement 1). We analyzed the 10 species with the most data points reflecting short divergence times (d_S < 0.0005), which is critical for data fit.

Population simulations overview

The majority of the computational simulations performed are built upon the idea of the Wright–Fisher model with selection (Tataru et al., 2017) that population generations can be determined from a multinomial distribution. However, we have made some changes to generalize this model for our purposes.

First, the simulations do not necessarily assume a constant population size but rather assume the population grows via a logistic growth model with a capacity to allow for the implementation of bottlenecks. Specifically, if $P [t]$ is the population on generation $t$ and $K$ is the population capacity, then

\bar{P [t + 1]} = P [t] + P [t] (1 - \frac{P [t]}{K})

And

P [t + 1] = P o i s s (\bar{P [t + 1]}) .

The population size is a Poisson random variable as we choose to determine the offspring of individual genetic classes as a Poisson random variable. We note that except for the very first few generations and after bottlenecks, the population size only has small fluctuations around a fixed capacity.

To speed the simulation up and enable the simulation of very large population sizes, we implemented a variety of genotype classes, rather than tracking each genotype individually. Genotype classes are similar to the practice of simulating fitness classes (Desai and Fisher, 2007), though we manage the number of unique classes via Poisson merging and splitting.

For all simulations, we start from a single organism that begins with 500,000 neutral alleles, representing a core genome size of this many codons that has yet to receive any mutations. When a mutation occurs, one allele may change types or stay the same, depending on the mutation received and the state of the randomly chosen codon. For example, a deleterious mutation occurring at a codon already in a deleterious state does not change the genotype class of the organism.

The specific implementations and additional parameters used for this model are provided in the following sections. Here, we outline the theory that ensures that genotype classes accurately represent such a population and enable the calculation of fitness. Consider the total population of size, $P [t]$ , at generation $t$ , as a composite of multiple different classes. The number of individuals in class $j$ on generation $t$ will be $A_{j} [t]$ . We have

P [t] = \sum_{j \geq 0} A_{j} [t] .

Within class $j$ , we store several variables that provide information about the genotype of members of $A_{j} [t]$ . Specifically, we store a number $j_{k}$ that specifies the number of alleles of type $k$ in the class $j$ . Examples of potential types that are used in our work include deleterious alleles, adaptive alleles, and alleles that result from reversion. Each type of allele is associated with a specific selective advantage $s_{k}$ . We can now write a formula to calculate the absolute fitness $F_{j}$ of class $j$ :

F_{j} = \prod_{k \geq 0} {(1 + s_{k})}^{j_{k}} .

From the absolute fitness $F_{j}$ , we calculate the average absolute fitness of the population on generation $t$ via

E [F] = \frac{\sum_{j \geq 0} F_{j} A_{j} [t]}{P [t]} .

We now calculate the relative fitness of class $j$ on generation $t$ as

f_{j} = \frac{F_{j}}{E [F]} .

Next, we calculate the expected size of class $j$ in the next generation with

A_{j} [t + 1] = f_{j} (A_{j} [t] + A_{j} [t] (1 - \frac{P [t]}{K})) .

Note that

P [t + 1] = \sum_{j \geq 0} A_{j} [t + 1] = \sum_{j \geq 0} f_{j} (A_{j} [t] + A_{j} [t] (1 - \frac{P [t]}{K})) =

\sum_{j \geq 0} \frac{F_{j}}{E [F]} (A_{j} [t] + A_{j} [t] (1 - \frac{P [t]}{K})) = P [t] + P [t] (1 - \frac{P [t]}{K}) .

This allows us to use a logistic model of growth to represent population size rather than being constrained to fixing it, which is useful for simulating bottlenecks.

To account for genetic drift through random fluctuations, we rewrite the above equations to be

A_{j} [t + 1] = P o i s s (f_{j} (A_{j} [t] + A_{j} [t] (1 - \frac{P [t]}{K}))) .

which also implies

P [t + 1] = P o i s s (P [t] + P [t] (1 - \frac{P [t]}{K})) .

Note that this simulation still has equivalent dynamics of the frequencies of classes as a Wright–Fisher model with selection in the case of a fixed population size due to the ability to split Poisson processes, that is,

P r (A_{j} [t + 1] = X | P [t + 1]) = P r (B i n (P [t + 1], \frac{f_{j} A_{j} [t + 1]}{P [t]}) = X) .

Mutations are added in every generation depending on the mutation rate. Only single mutants are generated per generation, and an organism cannot get more than one mutation per generation. The number of new mutants is determined by the binomial distribution. New mutants are added then to their appropriate class $j$ . For example, if a deleterious mutation is gained in a class with 10 deleterious alleles (and nothing else), this new mutant will increase the population size of the class with 11 deleterious alleles (and nothing else) while decreasing the population size of the class with 10 deleterious alleles (and nothing else).

By grouping individuals in classes rather than by genotype, computational costs can be greatly cut down. Grouping individuals does not affect the dynamics of the simulation because the merging of Poisson processes is still Poisson. The downside to this approach is information loss, though by designing custom alleles, we can track specific mutational histories like reversions.

Purifying selection simulations

The purifying selection simulations (Figure 2) utilized the base framework as mentioned above. Effective population sizes (N_e) varied from 10⁶ to 10¹⁸ depending on the simulation. The simulation begins with an initial organism with 500,000 neutral alleles (representing a WT core genome). The population quickly grows to the carrying capacity (N_e) and follows logistic growth (see ‘Population simulations overview’). The deleterious nonsynonymous mutation rate per genome per generation is 1.01 × 10^–3. These deleterious mutations have a selective disadvantage of s ≈ 3.5 × 10^–5. Synonymous and neutral nonsynonymous mutations are not simulated directly as they are neutral and are instead assumed to accumulate in the population with an average rate of 3.75 × 10^–4 and 1.12 × 10^–4 per genome per generation, respectively.

We estimate the average d_N/d_S of the population by taking the average number of codon differences between two individuals in the population to be twice the average number of mutations in the population. This approximation is valid due to the lack of selective sweeps, bottlenecks, and large effective population size, which results in expected coalescent time between random individuals being 10⁶–10¹⁸ generations (far longer than our simulations).

Variations of this basic purifying selection model are performed as described in the article, including increasing the mutation rate to 1.13 × 10^–3 mutations per genome per generation (Figure 2—figure supplement 2) and the simulation of recombination (Figure 2—figure supplement 3). For simulations of recombination, we assume that transitions to the ancestral state (purging the deleterious allele) occur at a rate of 2.5 × 10^–7 per codon per genome per generation and bring along a synonymous mutation (tracked via the number of recombinations to the ancestral state). This procedure does not allow for recombination to purge multiple deleterious alleles at a time; such events are unlikely given that deleterious alleles are rare and randomly distributed. We also include a model in which synonymous mutations are not included in this reversion event.

The simulation of purifying selection through bottlenecks and infrequent adaptive sweeps was performed as in the modified version of the reversion model (see ‘Reversion simulations’), though with less frequent adaptation and larger and less frequent bottlenecks.

Stop codon enrichment in the Zhao and Lieberman et al. dataset

Table S7 in Zhao et al., 2019 provides an Excel sheet detailing all observed mutations. There were 325 observed nonsynonymous mutations of which 28 were stop codons. Under a null model, there are 415 possible permutations of initial codon and codon one mutation (see ‘Code availability’) away that result in a nonsynonymous substitution of which 23 lead to a stop codon. Assuming no preference for specific mutation or initial codon, we would expect roughly 18 stop codons in this data. Under a null binomial distribution, the p-value for obtaining 28 or more is 0.015.

Reversion simulations

We simulated gut bacterial populations using a modified Wright–Fisher model (see ‘Population simulations overview’) to monitor mutation acquisition over time compared to an ancestor. Like all simulations, we begin with a single organism with 500,000 neutral alleles to represent the WT core genome. The population can grow to a capacity of 10¹⁰ via a logistic growth model. Environmental changes occur with a probability of $\frac{n_{l o c i}}{τ_{f l i p}}$ per generation, triggering an average of one selective pressure per environmental change, modeled by a Poisson distribution. Population bottlenecks to 10 individuals occur independently of environmental changes with a probability of 10^–4 per generation (see Figure 4—figure supplement 1 for an alternative in which bottlenecks and environmental change are correlated).

Both adaptive selective pressures and adaptive mutations are categorized into two allele types: forward and reverse. These classes are designed to enable tracking of complete mutational history and therefore recorded relative to the ancestral state rather than the current state. Thus, actual d_N/d_S is calculated as the sum of these mutations and observed d_N/d_S using their difference (plus asymptomatic d_N/d_S). When releasing beneficial selective pressures, their classification as forward or reverse is based on the balance of previously released selective pressures: the probability of an adaptive pressure being classified as a reverse adaptation increases as the number of forward pressures increases and is equal to the difference between the forward pressures previously released ( $q_{F}$ ) and reverse pressures previously released ( $q_{R}$ ) divided by the number of loci (i.e., $\frac{q_{F} - q_{R}}{n_{l o c i}}$ ). All beneficial mutations have a selective advantage of s_ben = 0.03 (for forward or reverse). The rate at which mutations occur given an available pressure depends on whether the mutation is adapting to a forward or reverse pressure: the reversion rate is set at one-fifth the rate of the nonsynonymous per codon mutation rate (4.5 × 10^–10 per generation per cell), while forward mutations are set at a rate five times higher than the nonsynonymous mutation rate because they can happen at multiple sites (25× the reversion rate; 1.1 × 10^–8 per generation per cell).

This simulation treats each adaptive mutation as occupying a unique codon in the core genome for simplicity. This assumes that the ancestral allele at a given locus has been purged before the next environmental change affecting that locus (or selective pressure); as theory suggests that a beneficial mutation takes 768 generations to fix (see Appendix 1, Section 2.2), compared to 46,200 generations for pressure shifts at any locus we believe this assumption is reasonable. To confirm this theory still holds in the presence of bottlenecks and clonal interference, we tracked the average number of beneficial mutations in the population relative to the number of selective pressures released at any generation in simulations; we found only a 2% deviation between the average and expected total beneficial mutations over 2 × 10⁶ generations.

Throughout the simulation, deleterious mutations occur at a rate of 1.01 × 10^–3 mutations per genome per generation, with a selective disadvantage of s = 0.003, and can be reverted to a deleterious reversion allele class (separate from the adaptive reversion allele class).

To calculate d_N/d_S, we assume the simulated population could be compared to an equivalent population but with distinct mutations, allowing us to calculate the d_N/d_S as using double the current observed substitutions.

We assume during the reversion simulations that the ancestor has no initial transient mutations. We make this assumption for computational simplicity but the theoretical curve is equivalent whether starting from no revertible mutations or the equilibrium where half of the loci currently have forward mutations (assuming forward and reverse mutations occur at equal rates; see Appendix 1, Section 3.3).

Testing standard d_N/d_S software

We simulated gene sequences with selective pressures acting at specific sites for Figure 4b. For each genomic distance investigated (every 500,000 bacterial generations), we ran 10 simulations as described below, with each simulation resulting in 10 sequences derived from a branching process. In the permanent adaptations simulation (blue), adaptive mutations in the phylogeny are acquired simply and permanently. In the transient adaptations simulation (red), only more recent mutations in the same phylogeny will be visible while older mutations are obscured by reversion. Both sets of sequences were then fed to PAML v4.8 (Yang, 2007) for the estimation of d_N/d_S values. PAML uses maximum likelihood analysis to estimate the rate of substitution that best explains a given phylogenetic tree.

For each simulation, we generated a random 1500 bp open-reading frame and designated 10% of codons as neutral, 10% under positive selection, and 80% under purifying selection. We introduced mutations and branches across several cycles, with each cycle representing 100,000 generations. For each cycle, we assigned mutations at random according to the following probabilities: 67.5% that a nonsynonymous mutation occurred at a codon under positive selection, 12.5% for synonymous mutation at any codon, 3.75% for nonsynonymous mutation at a nonselected site (neutral), and 16.25% for no mutation. These values were selected to give a d_N/d_S of around 2 and to match the general ratios in the reversion model.

Two phylogenies were constructed from each simulation: both received identical mutations, but they differed in how nonsynonymous mutations at selective codon sites were visible at the end of the simulation. In the transient adaptations version of the phylogeny, nonsynonymous mutations at selective codon sites were reverted at the end of the simulation, except those that occurred within the last 500,000 generations. Reverted sites were converted to a synonymous substitution at a frequency based on the codon table (assuming an equal probability of all nucleotide mutations). Both versions of the sequences underwent multiple sequence alignment and neighbor-joining tree construction (Biopython; Cock et al., 2009). We calculated treewide d_N/d_S ratios using PAML v4.8’s codeML feature (Álvarez-Carretero et al., 2023), employing the M2a model to analyze site-specific selection.

Closeness to stop codons

To evaluate possible enrichment for stop codon adjacency, we focused on TTA and TCA codons. TTA and TCA are ideal for measuring the likelihood of nonsense mutations because each has two point mutations that yield a stop codon, unlike the five other redundant codons encoding for the same amino acids (for both leucine and serine, one codon is singly stop codon adjacent and the last four are not stop codon adjacent). In B. fragilis, these codons have a codon usage rate of 13% for leucine and 14% for serine.

We annotated the reference genome NCTC_9343 with Bakta v1.9 (Schwengers et al., 2021) and obtained COG categories for each gene using eggNOG v5.0 (Cantalapiedra et al., 2021). Genes that did not have a functional COG category (35%) were removed. To control for unusual outlier genes skewing results, only the 15 COG groups that had at least 50 genes were considered for enrichment analyses. For each COG category, we calculated a null codon usage proportion based on the proportion of leucine and serine codons and compared this to the actual proportion using a one-proportion Z test. To address the fact that genes in the same functional category but localized to different parts of the cell may be under different selective pressures, we analyzed cellular location classifications from PSORTb v3.02 (Yu et al., 2010) and categorized genes by the combination of function and localization. We analyzed the 15 function-location combinations with more than 50 genes. After identifying those categories that were significantly enriched, we cross-referenced which categories of genes shown to be under adaptive within-person evolution in a previous study of B. fragilis within-person evolution were in Zhao et al., 2019. Of the 16 genes reported in that paper, 8 were assigned a functional COG category/cellular location and in the reference genome (NCTC_9343). Four of these were in outer membrane inorganic ion transport and metabolism, a significant enrichment ( $p = 1.22 \times 10^{- 4}$ ; binomial test) (Supplementary file 1, Supplementary file 2).

Code availability

Code and simulation results are available at https://github.com/PaulTorrillo/Microbiome_Reversions (copy archived at Torrillo, 2023).

Appendix 1

Supporting information

1.1 d_N/d_S theory

The following is meant to provide a more in-depth walkthrough of how one can build up and interpret the purifying selection model and its effect of time dependence on d_N/d_S. To be accessible to a wide audience and self-contained, we have included enough detail that most sections should be followable with basic knowledge of calculus.

Assume an infinite population of organisms. Consider the existence of $m$ classes of nonsynonymous mutations. The number of mutations of the ith class in the population is represented by the variable ${\bar{N}}_{i}$ . Each class ${\bar{N}}_{i}$ has an associated mutation rate $U_{i}$ (mutations per genome per generation per unit time) and an associated selective disadvantage $s_{i}$ (mutation purification per unit time). In the purifying selection model, we assume that $s_{0} = 0$ and $s_{i > 0} > 0$ . We assume that both the mutation rate and the selective disadvantage of each class remain constant throughout time. We assume this is the global population and hence no migration. We then have

d {\bar{N}}_{i} = U_{i} d t - s_{i} {\bar{N}}_{i} d t

Assuming $i > 0$ , we can integrate

\frac{d N_{i}}{U_{i} - s_{i} {\bar{N}}_{i}} = d t .

Using u-substitution, we let $u = μ_{i} - s_{i} {\bar{N}}_{i_{i}}$ so that $\frac{d u}{d {\bar{N}}_{i}} = - s_{i}$ which implies $d {\bar{N}}_{i} = \frac{d u}{- s_{i}}$ so that

\frac{d u}{- s_{i} u} = d t .

Integrating both sides, we get

\frac{- \ln (u)}{s_{i}} = t + C .

where $C$ is the constant of integration. Continuing we have

u = e^{- s_{i} t - s_{i} C}

μ_{i} - s_{i} {\bar{N}}_{i} = e^{- s_{i} t - s_{i} C}

{\bar{N}}_{i} = \frac{μ_{i}}{s_{i}} - \frac{e^{- s_{i} t - s_{i} C}}{s_{i}} .

We can remove the constant of integration and instead replace it with the initial condition ${\bar{N}}_{i} (0)$

{\bar{N}}_{i} (0) = \frac{U_{i}}{s_{i}} - \frac{e^{- s_{i} C}}{s_{i}} .

So, then we have that

C = \frac{l n (U_{i} - s_{i} {\bar{N}}_{i} (0))}{- s_{i}} .

So that we get

{\bar{N}}_{i} = \frac{U_{i}}{s_{i}} (1 - e^{- s_{i} t}) + {\bar{N}}_{i} (0) e^{- s_{i} t} .

Since ${\bar{N}}_{i} (0)$ should occur in a homogeneous population that is the most recent common ancestor of the individuals in the population we are observing, we assume ${\bar{N}}_{i} (0)$ = 0, so we have

{\bar{N}}_{i} = \frac{U_{i}}{s_{i}} (1 - e^{- s_{i} t}) .

If $i = 0$ , we have

{\bar{N}}_{0} = U_{0} t .

We also assume that there are synonymous mutations $\bar{S}$ that are neutral and occur with new mutations per unit time $U_{S}$ so that

\bar{S} = U_{S} t .

So, if we want to find $\bar{N}$ (the total number of nonsynonymous mutations), that will be given by

\bar{N} = U_{0} t + \sum_{i \geq 1} \frac{U_{i}}{s_{i}} (1 - e^{- s_{i} t}) .

Now observe the following with G as the number of base pairs in the core genome, then $d_{N} = 2 \bar{N} / (G \times 3 / 4)$ and $d_{S} = 2 \bar{S} / (G \times 1 / 4)$ . The 2 comes from the fact that there are two diverged lineages when calculating d_N/d_S.

\frac{d_{N}}{d_{S}} = \frac{1}{3} \frac{\bar{N} / G}{\bar{S} / G} = \frac{U_{0} t + \sum_{i \geq 1} \frac{U_{i}}{s_{i}} (1 - e^{- s_{i} t})}{U_{S} t} = \frac{1}{3} (\frac{U_{0}}{U_{S}} + \sum_{i \geq 1} \frac{U_{i} / U_{S}}{s_{i} t} (1 - e^{- s_{i} t})) .

The ⅓ is for normalization. While the above form is easier to analyze in terms of actual time, it should also be noted that data is given in terms of d_S rather than $t$ so the following equivalent form can also be helpful when discussing fitting of the timescale dependence of d_N/d_S:

\frac{d_{N}}{d_{S}} = \frac{1}{3} (\frac{U_{0}}{U_{S}}) + \sum_{i \geq 1} \frac{2 U_{i}}{G s_{i} d_{S}} (1 - e^{- \frac{G s_{i}}{2 U_{S}} d_{S}}) .

Furthermore, we can rewrite $U_{i} = {3 α}_{i} U_{S}$ so that we have

\frac{d_{N}}{d_{S}} = α_{0} + \sum_{i \geq 1} \frac{2 α_{i} U_{S}}{G s_{i} d_{S}} (1 - e^{- \frac{G s_{i}}{2 U_{S}} d_{S}}) .

Note $α_{0}$ is equivalent to $α$ in the main text. Next, we can simply rewrite $β_{i} = \frac{G s_{i}}{{2 U}_{S}}$ so that we have

\frac{d_{N}}{d_{S}} = α_{0} + \sum_{i \geq 1} \frac{α_{i}}{β_{i} d_{S}} (1 - e^{- β_{i} d_{S}}) .

To begin analyzing (S5), we consider the asymptotic behavior. First, we note that

\underset{d S \to + \infty}{l i m} α_{0} + \sum_{i \geq 1} \frac{α_{i}}{β_{i} d_{S}} (1 - e^{- β_{i} d_{S}}) = α_{0} = \frac{U_{0}}{U_{S}} .

For initial behavior, we can use L’Hopital’s rule

\underset{d S \to + 0}{l i m} α_{0} + \sum_{i \geq 1} \frac{α_{i}}{β_{i} d S} (1 - e^{- β_{i} d S}) = α_{0} + \sum_{i \geq 1} α_{i} = \frac{U_{0} + \sum_{i \geq 1} U_{i}}{U_{S}} .

Finally, we are interested in when d_N/d_S ends up being in between these initial and final values. Each class $i$ has a different midpoint at which its contribution to d_N /d_S is a half. In mathematical terms, this can be summarized as

\frac{α_{i}}{2} = \frac{α_{i}}{β_{i} d_{S}} (1 - e^{- β_{i} d_{S}}) .

which implies

\frac{β_{i} d_{S}}{2} = (1 - e^{- β_{i} d_{S}}) .

Thus, we see that for any class of mutation, $\frac{1}{s_{i}}$ will determine the midpoint (assuming some constant $U_{S}$ and $G$ ). We can then imagine the curve as similar to a step function where the location of each step is determined by the corresponding $\frac{1}{s_{i}}$ and how far the function steps down will be determined by the size of the mutation rate $U_{i}$ . Finally, we can surmise that the average step-down occurs at approximately the harmonic average

s^{- 1} = \sum_{i \geq 1} \frac{U_{i}}{U_{N} s_{i}} .

where $U_{N} = \sum_{i \geq 1} U_{i}$ , which represents the total non-neutral nonsynonymous mutation rate. Taking this all into account, here is exactly what is being fit in the d_N/d_S curve. We are fitting the equation

\frac{d_{N}}{d_{S}} = α + \frac{1 - α}{β d_{S}} (1 - e^{- β d_{S}}) .

$α$ is the fraction of nonsynonymous mutations that are neutral. Now $β$ is a compound parameter and we can fit $β = \frac{G s}{{2 U}_{S}}$ which is essentially half of the harmonic average of selective disadvantages in ratio with the synonymous mutation rate per site per generation. Now we can estimate the synonymous mutation per site per generation with the following:

10^{- 9} \frac{m u t a t i o n s}{b a s e p a i r \times g e n e r a t i o n} \times 3 \frac{b a s e p a i r s}{s i t e} \times \frac{1}{4} \frac{s y n o n y m o u s m u t a t i o n}{m u t a t i o n s} = \frac{3}{4} \times 10^{- 9} \times \frac{s y n o n y m o u s m u t a t i o n s}{s i t e \times g e n e r a t i o n} .

Thus, we are truly fitting

\frac{d_{N}}{d_{S}} = α + \frac{\frac{3}{4} \times (1 - α) \times 10^{- 9}}{2 s d_{S}} (1 - e^{- \frac{s}{\frac{3}{2} \times 10^{- 9}} d_{S}}) .

1.2 Mutation accumulation

To analyze genetic drift and Muller’s ratchet (Haigh, 1978; Neher and Shraiman, 2012), we will provide a brief overview of the approach well suited for our work. For any nonsynonymous mutation class $N_{i}$ with $i > 0$ , we can track the size of a population with 0 mutations of class $N_{i}$ , denoted by variable $W_{i}$ , via

d W_{i} = s_{i} {\bar{N}}_{i} (t) W_{i} d t - U_{i} W_{i} d t .

Here, ${\bar{N}}_{i} (t)$ is the average number of mutations of class ${\bar{N}}_{i}$ in the population and $s_{i} {\bar{N}}_{i} (t) W_{i}$ is equivalent to mean fitness. This equation reflects how every unit of time, the mutation-free class should increase by its selective advantage relative to the population though also loses members of the population to the mutation rate.

If we substitute in ${\bar{N}}_{i} (t)$ , we get

d W_{i} = s_{i} \frac{U_{i}}{s_{i}} (1 - e^{- s_{i} t}) d t - U_{i} W_{i} d t .

d W_{i} = U_{i} (1 - e^{- s_{i} t}) W_{i} d t - U_{i} W_{i} d t .

d W_{i} = - U_{i} e^{- s_{i} t} W_{i} d t .

So then

\ln (W_{i}) = \frac{U_{i}}{s_{i}} e^{- s_{i} t} + C .

Assuming that $W_{i}$ at time 0 is given by $W_{0}$ (representing the initial population size which under our assumptions is always free of mutations and hence wild type), then

\ln (W_{0}) = \frac{U_{i}}{s_{i}} + C .

C = \ln (W_{0}) - \frac{U_{i}}{s_{i}} .

W_{i} (t) = W_{0} e^{- \frac{U_{i} (1 - e^{- s_{i} t})}{s_{i}}} .

We see that if we set $W_{0} = 1$ , then everything can be given in terms of frequencies within the population. Asymptotically, we have that

\underset{t \to \infty}{l i m} W_{0} e^{- \frac{U_{i} (1 - e^{- s_{i} t})}{s_{i}}} = W_{0} e^{- \frac{U_{i}}{s_{i}}} .

Importantly, we also note the following: let $W$ be the frequency of the wild type (mutation-free class) and $W_{0} = 1$ . Then,

W (t) = \prod_{i > 0} W_{i} (t) = e^{- \sum_{i > 0} \frac{U_{i} (1 - e^{- s_{i} t})}{s_{i}}} .

which implies

\underset{t \to \infty}{l i m} e^{- \sum_{i > 0} \frac{U_{i} (1 - e^{- s_{i} t})}{s_{i}}} = e^{- \sum_{i > 0} \frac{U_{i}}{s_{i}}} = e^{- \frac{U_{N}}{s}} .

Once again using

U_{N} = \sum_{i \geq 1} U_{i} .

And

s^{- 1} = \sum_{i \geq 1} \frac{U_{i}}{U_{N} s_{i}} .

Furthermore, the average time for the frequency in the population without mutations of class i to reach one-half of the logarithm of the asymptotic frequency is the same when using the above simplifications. In other words, the difference in how much and how fast the wild type is lost should not depend too heavily on the distribution of selective coefficients.

Finally, we can predict if the least-loaded class will be lost to drift. We can form this prediction via the following. We assume that if the least-loaded class $W$ drops below its steady-state frequency, $e^{- \frac{U_{N}}{s}}$ , it will have advantage $s$ every generation. If the least-loaded class has advantage $s$ , then it has a 2s probability of extinction (see Appendix 1, Section 2.1). Hence, if we expect there to be $N_{e} e^{- \frac{U_{N}}{s}}$ individuals, we can estimate that mutation accumulation will occur when

2 s N_{e} e^{- \frac{U_{N}}{s}} << 1.

2.1. Extinction and fixation probability

Here, we will derive how the fixation probability of a mutation with selective advantage $s_{b e n}$ is approximately $2 s_{b e n}$ . This is a standard result that can be found in classic population genetics textbooks. It is usually derived via differential equations but can also be obtained more classically from the study of branching processes (Haldane, 1927) ,which we will use here. First, we assume that individuals reproduce via a Galton–Watson branching process with mean of 1 + $s_{b e n}$ . We also assume there are no other mutations in the population that can interfere with fixation or extinction. A fundamental result from the study of branching processes is that the extinction probability is given by the smallest non-negative root of the branching processes corresponding probability generating function, f(z). Being a Galton–Watson branching process, the probability generating function is the probability generating function of a Poisson process so

f (z) = e^{(1 + s_{b e n}) (1 - z)} .

Thus, we need to find the smallest non-negative solution, z, to

e^{(1 + s_{b e n}) (1 - z)} - z = 0.

We can use a second-order Taylor approximation to approximate the exponential so we have

1 + (1 + s_{b e n}) (z - 1) + \frac{{((1 + s_{b e n}) (z - 1))}^{2}}{2} - z = 0.

which has the smallest non-negative solution of

z = \frac{s_{b e n}^{2} + 1}{{(s_{b e n} + 1)}^{2}} .

If $s_{b e n}$ is small, then we have

z \approx \frac{1}{1 + 2 s_{b e n}} .

which we can further approximate by doing

z \approx \frac{1}{1 + 2 s_{b e n}} \frac{1 - 2 s_{b e n}}{1 - 2 s_{b e n}} \approx 1 - 2 s_{b e n} .

If $1 - 2 s_{b e n}$ is the extinction probability, then the fixation probability will be $2 s_{b e n}$ .

2.2. Likelihood of reversion

We calculate the time for mutations to occur and fix in the population for demonstration in Figure 3b. For the sake of simplicity, we consider the weak-mutation strong selection regime (no clonal interference) in our theory but do include such dynamics in our simulations. Regardless, clonal interference will only minorly change the frequency of a revertant with high $s_{b e n}$ in the population as multiple backgrounds will find the same reversion when $N_{e}$ is sufficiently high. The expected time to fixation given a mutation with selective coefficient $s_{b e n}$ is estimated using

{(1 + s_{b e n})}^{t} = N_{e} .

where $N_{e}$ is the census population size. The above implies

t l n (1 + s_{b e n}) = l n (N_{e}) .

And if $s_{b e n}$ is small, then we have

t = \frac{l n (N_{e})}{s_{b e n}} .

Now we need to know the expected time for a fixing mutant to arise. First, we want the probability the mutation arises and fixes, which will be given by

1 - {(1 - 2 s_{b e n} U_{b e n})}^{N_{e}} \approx 2 N_{e} s_{b e n} U_{b e n} .

This implies the expected time to arrive is approximately

\frac{1}{2 N_{e} s_{b e n} U_{b e n}} .

So, therefore, in the absence of clonal interference, the time for reversion to fix in the population depends more on the time for it to fix over time than the time for the reversion to arise (and is, therefore, less dependent on the mutation rate) if

\frac{1}{2 N_{e} s_{b e n} U_{b e n}} < \frac{l n (N_{e})}{s_{b e n}} .

Alternatively written as

1 < 2 N_{e} U_{b e n} l n (N_{e}) .

Finally, note that the expected time to fixation is given by

\frac{1}{2 N_{e} s_{b e n} U_{b e n}} + \frac{l n (N_{e})}{s_{b e n}} .

The second term can be larger because of bottlenecks and clonal interference.

3.1 Reversion model with fluctuating loci

The reversion model can be derived in the following way. First, one can rewrite Equation 2 as a generic negative feedback model in which mutations emerge and are purged from the population at a rate proportional to how many mutations have accumulated:

d {\bar{N}}_{t r a n s i e n t} = R_{i n} d t - P_{o u t} {\bar{N}}_{t r a n s i e n t} d t .

In this formulation, R_in is simply the rate of accumulation of transient non-neutral nonsynonymous mutations per genome per unit time and P_out is the loss rate of these mutations per unit time.

From this more general form, we can develop the reversion model. Similar to the original purifying selection model, it is possible to assume a variety of classes of mutations, but for simplicity, we only assume 1. If we make the following definitions for

R_{i n} = \frac{n_{l o c i}}{τ_{f l i p}} .

P_{o u t} = \frac{2}{τ_{f l i p}} .

we can link Equation S21 to a fluctuating loci model via

\begin{array}{cc} \frac{d {\bar{N}}_{t r a n s i e n t}}{d t} & = R_{i n} - P_{o u t} {\bar{N}}_{t r a n s i e n t} = \frac{n_{l o c i}}{τ_{f l i p}} - \frac{2}{τ_{f l i p}} {\bar{N}}_{t r a n s i e n t} \\ = \frac{1}{τ_{f l i p}} (n_{l o c i} - {\bar{N}}_{t r a n s i e n t}) - \frac{1}{τ_{f l i p}} {\bar{N}}_{t r a n s i e n t} \end{array} .

Here, $τ_{f l i p}$ represents the average number of generations it takes for a given selective pressure on a locus to switch direction, and $n_{l o c i}$ , the number of loci under fluctuating selective pressures. Equation S22 can be more directly linked to a fluctuating loci model as we see the rate of nonsynonymous mutations is proportional to $n_{l o c i} - {\bar{N}}_{t r a n s i e n t}$ (the number of loci unmutated) and the rate out is proportional to ${\bar{N}}_{t r a n s i e n t}$ the number of loci mutated.

Solving similar to the purifying selection model, we have

{\bar{N}}_{t r a n s i e n t} (t) = \frac{n_{l o c i}}{2} (1 - e^{- \frac{2 t}{τ_{f l i p}}}) .

This can then be used to obtain d_N/d_S

\frac{d_{N}}{d_{S}} = \frac{1}{3} (\frac{U_{0}}{U_{S}} + \frac{n_{l o c i}}{2 U_{S} t} (1 - e^{- \frac{2 t}{τ_{f l i p}}})) .

Or equivalently

\frac{d_{N}}{d_{S}} = \frac{1}{3} (\frac{U_{0}}{U_{S}} + \frac{n_{l o c i}}{G d S} (1 - e^{- \frac{G d_{S}}{U_{S} τ_{f l i p}}})) .

If we set $α = \frac{U_{0}}{3 U_{S}}, β = \frac{G}{U_{S} τ_{f l i p}}$ , and $γ = \frac{n_{l o c i}}{3 U_{S} τ_{f l i p}}$ so that

\frac{d_{N}}{d_{S}} = α + \frac{γ}{β d S} (1 - e^{- β d_{S}}) .

which is equivalent to Equation S9 with one more free parameter.

3.2 Effect of compensatory mutations

We can further extend the theory to include compensatory mutations by calculating the expected value of ${\bar{N}}_{t r a n s i e n t}$ as a Markov chain. First, let $v$ be the state vector for a locus under selection where the index of each row corresponds to the number of observed mutations currently at that locus and $A$ be the corresponding stochastic matrix with rows i and columns $j$ . First, we need to consider the probability the state does not change on a given generation (i.e., the diagonal of $A$ ). This will be $\frac{1}{τ_{f l i p}}$ . For every element on the diagonal of $A$ , we thus have $A_{i = j} = 1 - \frac{1}{τ_{f l i p}}$ .

Supposing there is a state change, let us define $p$ as the probability there is an increase in observed mutations (a forward or compensatory mutation) and $1 - p$ , the probability there is a decrease in observed mutations (a reversion). Then, $A_{j = i - 1} = \frac{1}{τ_{f l i p}} p$ and $A_{j = i + 1} = \frac{1}{τ_{f l i p}} (1 - p)$ (with the exception of $A_{i = 1, j = 0} = \frac{1}{τ_{f l i p}}$ ). Finally, with $Δ t$ being the number of time steps and assuming $v_{0} = 1$ and $v_{i > 0} = 0$ , then

{\bar{N}}_{t r a n s i e n t} (Δ t) = n_{l o c i} \sum_{i \geq 0} {(A^{Δ t} v)}_{i} \cdot i .

If $p > 1 - p$ , compensatory mutations are gained at a linear rate and the compensatory mutation rate will factor into the asymptotic d_N/d_S value. Conversely, if $p < 1 - p$ , the number of compensatory mutations is almost surely finite by the central limit theorem and hence will not factor into the asymptotic d_N/d_S. Finally, if $p = 1 - p$ , this is the classic elementary random walk well known to deviate from the origin with $O (\sqrt{n})$ . Interestingly, this is also sublinear and will not factor into asymptotic d_N/d_S.

3.3 Reversion model starting from equilibrium conditions

Our simulations and theory assume that the initial population starts with no forward mutations (i.e., WT) for simplicity. However, starting at equilibrium conditions does not impact the shape of the curve. The intuition here is that while starting from equilibrium enables the identification of reversions of initially transient mutations, these will be subsequently hidden by parallel evolution. Noting Equation S23, we have that for our simulations and base theory:

d_{N} = n_{l o c i} (1 - e^{- \frac{2 t}{τ_{f l i p}}}) .

We can show the same result occurs starting at equilibrium conditions. A Python script confirming the algebra is available on the GitHub repository. First, let there be three states a given locus can be in. The first state (i = 1) will be initial transient mutations, the second state (i = 2) will be ancestral alleles, and the third state (i = 3) will be subsequent transient mutations. We can then build a transition matrix:

A = (\begin{array}{ccc} - \frac{1}{τ_{f l i p}} & 0 & 0 \\ \frac{1}{τ_{f l i p}} & - \frac{1}{τ_{f l i p}} & \frac{1}{τ_{f l i p}} \\ 0 & \frac{1}{τ_{f l i p}} & - \frac{1}{τ_{f l i p}} \end{array}) .

We can then use the transition matrix to find the probability of being in a state at any given time via the matrix exponential, which is

e^{A t} = (\begin{array}{ccc} e^{- \frac{t}{τ_{f l i p}}} & 0 & 0 \\ \frac{1 - e^{- \frac{2 t}{τ_{f l i p}}}}{2} & \frac{1 + e^{- \frac{2 t}{τ_{f l i p}}}}{2} & \frac{1 - e^{- \frac{2 t}{τ_{f l i p}}}}{2} \\ \frac{1 + e^{- \frac{2 t}{τ_{f l i p}}}}{2} - e^{- \frac{t}{τ_{f l i p}}} & \frac{1 - e^{- \frac{2 t}{τ_{f l i p}}}}{2} & \frac{1 + e^{- \frac{2 t}{τ_{f l i p}}}}{2} \end{array}) .

Now note that we can use $e^{A t} (j, i)$ to find the probability of being in a given state j after starting at initial state i. First, let us calculate expected d_N at a locus that started from the ancestral allele and has subsequently diverged for time t. Take note that subsequent transient mutations are assumed to be distinct (i.e., they will always lead to at least one difference when compared to a different lineage).

d_{N} (i n i t a l a n c e s t r a l) = e^{A t} (2, 2) e^{A t} (3, 2) + e^{A t} (3, 2) e^{A t} (2, 2) + 2 e^{A t} (3, 2) e^{A t} (3, 2) .

Now let us calculate $d_{N}$ assuming the allele was initially a transient mutation.

+ e^{A t} (1, 1) e^{A t} (2, 1) + e^{A t} (2, 1) e^{A t} (1, 1) + 2 e^{A t} (1, 1) e^{A t} (3, 1) + 2 e^{A t} (3, 1) e^{A t} (1, 1) .

Noting that we have assumed equal rates of forward and reverse adaptations throughout, the equilibrium would be composed of $n_{l o c i}$ /2 initially ancestral loci and of $n_{l o c i}$ /2 initially transient loci. Working through all of the algebra, we find

d_{N} = \frac{n_{l o c i} (d_{N} (i n i t a l a n c e s t r a l) + d_{N} (i n i t a l t r a n s i e n t))}{2} .

d_{N} = \frac{n_{l o c i} (1 - e^{- \frac{2 t}{τ_{f l i p}}}) + n_{l o c i} (1 - e^{- \frac{2 t}{τ_{f l i p}}})}{2} .

d_{N} = n_{l o c i} (1 - e^{- \frac{2 t}{τ_{f l i p}}}) .

Data availability

Code and results of simulations are available at Github repository https://github.com/PaulTorrillo/Microbiome_Reversions (copy archived at Torrillo, 2023).

References

(2023) Beginner’s guide on the use of paml to detect positive selection
Molecular Biology and Evolution 40:msad041.

https://doi.org/10.1093/molbev/msad041
- PubMed
- Google Scholar
(2023) Quantifying the local adaptive landscape of a nascent bacterial community
Nature Communications 14:248.

https://doi.org/10.1038/s41467-022-35677-5
- PubMed
- Google Scholar
(2004) Silent nucleotide polymorphisms and a phylogeny for Mycobacterium tuberculosis
Emerging Infectious Diseases 10:1568–1577.

https://doi.org/10.3201/eid1009.040046
- PubMed
- Google Scholar
1. Barber MF
2. Elde NC
(2014) Escape from bacterial iron piracy through rapid evolution of transferrin
Science 346:1362–1366.

https://doi.org/10.1126/science.1259329
- PubMed
- Google Scholar
(2023) The landscape of adaptive evolution of a gut commensal bacteria in aging mice
Current Biology 33:5272.

https://doi.org/10.1016/j.cub.2023.11.021
- PubMed
- Google Scholar
1. Barrick JE
2. Lenski RE
(2013) Genome dynamics during experimental evolution
Nature Reviews. Genetics 14:827–839.

https://doi.org/10.1038/nrg3564
- PubMed
- Google Scholar
(2015) Adaptive immunity increases the pace and predictability of evolutionary change in commensal gut bacteria
Nature Communications 6:8945.

https://doi.org/10.1038/ncomms9945
- PubMed
- Google Scholar
1. Bobay LM
2. Ochman H
(2018) Factors driving effective population size and pan-genome evolution in bacteria
BMC Evolutionary Biology 18:153.

https://doi.org/10.1186/s12862-018-1272-4
- PubMed
- Google Scholar
(2021) eggNOG-mapper v2: functional annotation, orthology assignments, and domain prediction at the metagenomic scale
Molecular Biology and Evolution 38:5825–5829.

https://doi.org/10.1093/molbev/msab293
- PubMed
- Google Scholar
1. Carmody RN
2. Bisanz JE
3. Bowen BP
4. Maurice CF
5. Lyalina S
6. Louie KB
7. Treen D
8. Chadaideh KS
9. Maini Rekdal V
10. Bess EN
11. Spanogiannopoulos P
12. Ang QY
13. Bauer KC
14. Balon TW
15. Pollard KS
16. Northen TR
17. Turnbaugh PJ
(2019) Cooking shapes the structure and function of the gut microbiome
Nature Microbiology 4:2052–2063.

https://doi.org/10.1038/s41564-019-0569-4
- PubMed
- Google Scholar
(2011) The impact of recombination on dN/dS within recently emerged bacterial clones
PLOS Pathogens 7:e1002129.

https://doi.org/10.1371/journal.ppat.1002129
- PubMed
- Google Scholar
(2021) Towards an engineering theory of evolution
Nature Communications 12:3326.

https://doi.org/10.1038/s41467-021-23573-3
- PubMed
- Google Scholar
1. Charlesworth J
2. Eyre-Walker A
(2007) The other side of the nearly neutral theory, evidence of slightly advantageous back-mutations
PNAS 104:16992–16997.

https://doi.org/10.1073/pnas.0705456104
- PubMed
- Google Scholar
1. Chen DW
2. Garud NR
(2022) Rapid evolution and strain turnover in the infant gut microbiome
Genome Research 32:1124–1136.

https://doi.org/10.1101/gr.276306.121
- PubMed
- Google Scholar
1. Cock PJA
2. Antao T
3. Chang JT
4. Chapman BA
5. Cox CJ
6. Dalke A
7. Friedberg I
8. Hamelryck T
9. Kauff F
10. Wilczynski B
11. de Hoon MJL
(2009) Biopython: freely available Python tools for computational molecular biology and bioinformatics
Bioinformatics 25:1422–1423.

https://doi.org/10.1093/bioinformatics/btp163
- PubMed
- Google Scholar
1. Cooper VS
2. Lenski RE
(2000) The population genetics of ecological specialization in evolving Escherichia coli populations
Nature 407:736–739.

https://doi.org/10.1038/35037572
- PubMed
- Google Scholar
1. Culyba MJ
2. Van Tyne D
(2021) Bacterial evolution during human infection: Adapt and live or adapt and die
PLOS Pathogens 17:e1009872.

https://doi.org/10.1371/journal.ppat.1009872
- PubMed
- Google Scholar
(2015) Fate of a mutation in a fluctuating environment
PNAS 112:E5021–E5028.

https://doi.org/10.1073/pnas.1505406112
- PubMed
- Google Scholar
(1999) High frequency of cryptic deleterious mutations in Caenorhabditis elegans
Science 285:1748–1751.

https://doi.org/10.1126/science.285.5434.1748
- PubMed
- Google Scholar
1. Davies J
2. Davies D
(2010) Origins and evolution of antibiotic resistance
Microbiology and Molecular Biology Reviews 74:417–433.

https://doi.org/10.1128/MMBR.00016-10
- PubMed
- Google Scholar
1. Desai MM
2. Fisher DS
(2007) Beneficial mutation selection balance and the effect of linkage on positive selection
Genetics 176:1759–1798.

https://doi.org/10.1534/genetics.106.067678
- PubMed
- Google Scholar
(2013) Genetic diversity and the structure of genealogies in rapidly adapting populations
Genetics 193:565–585.

https://doi.org/10.1534/genetics.112.147157
- PubMed
- Google Scholar
1. Dewachter L
2. Brooks AN
3. Noon K
4. Cialek C
5. Clark-ElSayed A
6. Schalck T
7. Krishnamurthy N
8. Versées W
9. Vranken W
10. Michiels J
(2023) Deep mutational scanning of essential bacterial proteins can guide antibiotic development
Nature Communications 14:241.

https://doi.org/10.1038/s41467-023-35940-3
- PubMed
- Google Scholar
1. Dingle KE
2. Colles FM
3. Wareing DR
4. Ure R
5. Fox AJ
6. Bolton FE
7. Bootsma HJ
8. Willems RJ
9. Urwin R
10. Maiden MC
(2001) Multilocus sequence typing system for Campylobacter jejuni
Journal of Clinical Microbiology 39:14–23.

https://doi.org/10.1128/JCM.39.1.14-23.2001
- PubMed
- Google Scholar
1. Drake JW
(1991) A constant rate of spontaneous mutation in DNA-based microbes
PNAS 88:7160–7164.

https://doi.org/10.1073/pnas.88.16.7160
- PubMed
- Google Scholar
1. Druelle V
2. Neher RA
(2023) Reversions to consensus are positively selected in HIV-1 and bias substitution rate estimates
Virus Evolution 9:veac118.

https://doi.org/10.1093/ve/veac118
- PubMed
- Google Scholar
1. Faith JJ
2. Guruge JL
3. Charbonneau M
4. Subramanian S
5. Seedorf H
6. Goodman AL
7. Clemente JC
8. Knight R
9. Heath AC
10. Leibel RL
11. Rosenbaum M
12. Gordon JI
(2013) The long-term stability of the human gut microbiota
Science 341:1237439.

https://doi.org/10.1126/science.1237439
- PubMed
- Google Scholar
1. Feder AF
2. Kline C
3. Polacino P
4. Cottrell M
5. Kashuba ADM
6. Keele BF
7. Hu SL
8. Petrov DA
9. Pennings PS
10. Ambrose Z
(2017) A spatio-temporal assessment of simian/human immunodeficiency virus (SHIV) evolution reveals A highly dynamic process within the host
PLOS Pathogens 13:e1006358.

https://doi.org/10.1371/journal.ppat.1006358
- PubMed
- Google Scholar
1. Feil EJ
2. Cooper JE
3. Grundmann H
4. Robinson DA
5. Enright MC
6. Berendt T
7. Peacock SJ
8. Smith JM
9. Murphy M
10. Spratt BG
11. Moore CE
12. Day NPJ
(2003) How clonal is Staphylococcus aureus?
Journal of Bacteriology 185:3307–3316.

https://doi.org/10.1128/JB.185.11.3307-3316.2003
- PubMed
- Google Scholar
(2012) Measuring selection coefficients below 10(-3): method, questions, and prospects
Genetics 190:175–186.

https://doi.org/10.1534/genetics.111.133454
- PubMed
- Google Scholar
(2019) Evolutionary dynamics of bacteria in the gut microbiome within and across hosts
PLOS Biology 17:e3000102.

https://doi.org/10.1371/journal.pbio.3000102
- PubMed
- Google Scholar
1. Ghosh OM
2. Good BH
(2022) Emergent evolutionary forces in spatial models of luminal growth and their application to the human gut microbiota
PNAS 119:e2114931119.

https://doi.org/10.1073/pnas.2114931119
- PubMed
- Google Scholar
Book
1. Gillespie JH
(2004)
Population Genetics: A Concise Guide

JHU Press.
- Google Scholar
1. Haigh J
(1978) The accumulation of deleterious genes in a population--Muller’s Ratchet
Theoretical Population Biology 14:251–267.

https://doi.org/10.1016/0040-5809(78)90027-8
- PubMed
- Google Scholar
1. Haldane JBS
(1927) A mathematical theory of natural and artificial selection, part v: Selection and mutation
Mathematical Proceedings of the Cambridge Philosophical Society 23:838–844.

https://doi.org/10.1017/S0305004100015644
- Google Scholar
(2005) Time dependency of molecular rate estimates and systematic overestimation of recent divergence times
Molecular Biology and Evolution 22:1561–1568.

https://doi.org/10.1093/molbev/msi145
- PubMed
- Google Scholar
1. Jiang X
2. Hall AB
3. Arthur TD
4. Plichta DR
5. Covington CT
6. Poyet M
7. Crothers J
8. Moses PL
9. Tolonen AC
10. Vlamakis H
11. Alm EJ
12. Xavier RJ
(2019) Invertible promoters mediate bacterial phase variation, antibiotic resistance, and host adaptation in the gut
Science 363:181–187.

https://doi.org/10.1126/science.aau5238
- PubMed
- Google Scholar
1. Jolley KA
2. Kalmusova J
3. Feil EJ
4. Gupta S
5. Musilek M
6. Kriz P
7. Maiden MC
(2000) Carried meningococci in the Czech Republic: a diverse recombining population
Journal of Clinical Microbiology 38:4492–4498.

https://doi.org/10.1128/JCM.38.12.4492-4498.2000
- PubMed
- Google Scholar
1. Jukes TH
2. Cantor CR
(1969) Evolution of protein molecules
Mammalian Protein Metabolism 3:21–132.

https://doi.org/10.1016/B978-1-4832-3211-9.50009-7
- Google Scholar
1. Kelsic ED
2. Chung H
3. Cohen N
4. Park J
5. Wang HH
6. Kishony R
(2016) RNA structural determinants of optimal codons revealed by MAGE-Seq
Cell Systems 3:563–571.

https://doi.org/10.1016/j.cels.2016.11.004
- PubMed
- Google Scholar
1. Key FM
2. Khadka VD
3. Romo-González C
4. Blake KJ
5. Deng L
6. Lynn TC
7. Lee JC
8. Chiu IM
9. García-Romero MT
10. Lieberman TD
(2023) On-person adaptive evolution of Staphylococcus aureus during treatment for atopic dermatitis
Cell Host & Microbe 31:593–603.

https://doi.org/10.1016/j.chom.2023.03.009
- PubMed
- Google Scholar
1. Kibota TT
2. Lynch M
(1996) Estimate of the genomic mutation rate deleterious to overall fitness in E. coli
Nature 381:694–696.

https://doi.org/10.1038/381694a0
- PubMed
- Google Scholar
1. Kimura M
(1977) Preponderance of synonymous changes as evidence for the neutral theory of molecular evolution
Nature 267:275–276.

https://doi.org/10.1038/267275a0
- PubMed
- Google Scholar
1. Kolodny O
2. Schulenburg H
(2020) Microbiome-mediated plasticity directs host evolution along several distinct time scales
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 375:20190589.

https://doi.org/10.1098/rstb.2019.0589
- PubMed
- Google Scholar
1. Korem T
2. Zeevi D
3. Suez J
4. Weinberger A
5. Avnit-Sagi T
6. Pompan-Lotan M
7. Matot E
8. Jona G
9. Harmelin A
10. Cohen N
11. Sirota-Madi A
12. Thaiss CA
13. Pevsner-Fischer M
14. Sorek R
15. Xavier R
16. Elinav E
17. Segal E
(2015) Growth dynamics of gut microbiota in health and disease inferred from single metagenomic samples
Science 349:1101–1106.

https://doi.org/10.1126/science.aac4812
- PubMed
- Google Scholar
1. Koskella B
2. Brockhurst MA
(2014) Bacteria-phage coevolution as a driver of ecological and evolutionary processes in microbial communities
FEMS Microbiology Reviews 38:916–931.

https://doi.org/10.1111/1574-6976.12072
- PubMed
- Google Scholar
1. Kryazhimskiy S
2. Plotkin JB
(2008) The population genetics of dN/dS
PLOS Genetics 4:e1000304.

https://doi.org/10.1371/journal.pgen.1000304
- PubMed
- Google Scholar
(2022) Hydrodynamic flow and concentration gradients in the gut enhance neutral bacterial diversity
PNAS 119:e2108671119.

https://doi.org/10.1073/pnas.2108671119
- PubMed
- Google Scholar
(2000) Compensatory mutations, antibiotic resistance and the population genetics of adaptive evolution in bacteria
Genetics 154:985–997.

https://doi.org/10.1093/genetics/154.3.985
- PubMed
- Google Scholar
1. Lieberman TD
2. Michel J-B
3. Aingaran M
4. Potter-Bynoe G
5. Roux D
6. Davis MR Jr
7. Skurnik D
8. Leiby N
9. LiPuma JJ
10. Goldberg JB
11. McAdam AJ
12. Priebe GP
13. Kishony R
(2011) Parallel bacterial evolution within multiple patients identifies candidate pathogenicity genes
Nature Genetics 43:1275–1280.

https://doi.org/10.1038/ng.997
- PubMed
- Google Scholar
1. Lieberman TD
(2022) Detecting bacterial adaptation within individual microbiomes
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 377:20210243.

https://doi.org/10.1098/rstb.2021.0243
- PubMed
- Google Scholar
1. Liu Z
2. Good BH
(2024) Dynamics of bacterial recombination in the human gut microbiome
PLOS Biology 22:e3002472.

https://doi.org/10.1371/journal.pbio.3002472
- PubMed
- Google Scholar
1. Loo SL
2. Ong A
3. Kyaw W
4. Thibaut LM
5. Lan R
6. Tanaka MM
(2020) Nonsynonymous polymorphism counts in bacterial genomes: A comparative examination
Applied and Environmental Microbiology 87:e02002-20.

https://doi.org/10.1128/AEM.02002-20
- PubMed
- Google Scholar
(2015) Convergent evolution and adaptation of Pseudomonas aeruginosa within patients with cystic fibrosis
Nature Genetics 47:57–64.

https://doi.org/10.1038/ng.3148
- PubMed
- Google Scholar
1. Neher RA
2. Shraiman BI
(2012) Fluctuations of fitness distributions and the rate of Muller’s ratchet
Genetics 191:1283–1293.

https://doi.org/10.1534/genetics.112.141325
- PubMed
- Google Scholar
(2019) Selection pressures on RNA sequences and structures
Evolutionary Bioinformatics Online 15:1176934319871919.

https://doi.org/10.1177/1176934319871919
- PubMed
- Google Scholar
(2022) Reversible bacteriophage resistance by shedding the bacterial cell wall
Open Biology 12:210379.

https://doi.org/10.1098/rsob.210379
- PubMed
- Google Scholar
1. Peterson GI
2. Masel J
(2009) Quantitative prediction of molecular clock and ka/ks at short timescales
Molecular Biology and Evolution 26:2595–2603.

https://doi.org/10.1093/molbev/msp175
- PubMed
- Google Scholar
Preprint
1. Poret AJ
2. Schaefers M
3. Merakou C
4. Mansour KE
5. Lagoudas GK
6. Cross AR
7. Goldberg JB
8. Kishony R
9. Uluer AZ
10. McAdam AJ
11. Blainey PC
12. Vargas SO
13. Lieberman TD
14. Priebe GP
(2024) De Novo Mutations Mediate Phenotypic Switching in an Opportunistic Human Lung Pathogen
bioRxiv.

https://doi.org/10.1101/2024.02.06.579193
- Google Scholar
1. Robert L
2. Ollion J
3. Robert J
4. Song X
5. Matic I
6. Elez M
(2018) Mutation dynamics and fitness effects followed in single cells
Science 359:1283–1286.

https://doi.org/10.1126/science.aan0797
- PubMed
- Google Scholar
1. Rocha EPC
2. Smith JM
3. Hurst LD
4. Holden MTG
5. Cooper JE
6. Smith NH
7. Feil EJ
(2006) Comparisons of dN/dS are time dependent for closely related bacterial genomes
Journal of Theoretical Biology 239:226–235.

https://doi.org/10.1016/j.jtbi.2005.08.037
- PubMed
- Google Scholar
(2021) Bakta: rapid and standardized annotation of bacterial genomes via alignment-free sequence identification
Microbial Genomics 7:000685.

https://doi.org/10.1099/mgen.0.000685
- PubMed
- Google Scholar
1. Sender R
2. Fuchs S
3. Milo R
(2016) Revised estimates for the number of human and bacteria cells in the body
PLOS Biology 14:e1002533.

https://doi.org/10.1371/journal.pbio.1002533
- PubMed
- Google Scholar
(2022) Comparative population genetics in the human gut microbiome
Genome Biology and Evolution 14:evab116.

https://doi.org/10.1093/gbe/evab116
- PubMed
- Google Scholar
(2008) Prevalence of agr dysfunction among colonizing Staphylococcus aureus strains
The Journal of Infectious Diseases 198:1171–1174.

https://doi.org/10.1086/592051
- PubMed
- Google Scholar
(2013) Genomic insights into the fate of colistin resistance and Acinetobacter baumannii during patient treatment
Genome Research 23:1155–1162.

https://doi.org/10.1101/gr.154328.112
- PubMed
- Google Scholar
(2017) Recurrent reverse evolution maintains polymorphism after strong bottlenecks in commensal gut bacteria
Molecular Biology and Evolution 34:2879–2892.

https://doi.org/10.1093/molbev/msx221
- PubMed
- Google Scholar
(2013) Innate immune recognition of microbial cell wall components and microbial strategies to evade such recognitions
Microbiological Research 168:396–406.

https://doi.org/10.1016/j.micres.2013.02.005
- PubMed
- Google Scholar
1. Tajima F
(1996) Infinite-allele model and infinite-site model in population genetics
Journal of Genetics 75:27–31.

https://doi.org/10.1007/BF02931749
- Google Scholar
(2017) Statistical inference in the wright-fisher model using allele frequency data
Systematic Biology 66:e30–e46.

https://doi.org/10.1093/sysbio/syw056
- PubMed
- Google Scholar
1. Torrance EL
2. Burton C
3. Diop A
4. Bobay LM
(2024) Evolution of homologous recombination rates across bacteria
PNAS 121:e2316302121.

https://doi.org/10.1073/pnas.2316302121
- PubMed
- Google Scholar
Software
1. Torrillo P
(2023) Microbiome reversions, version swh:1:rev:70cb1d1484e57a0aa5cfdd719c705515d21e818e
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:cc5db861b46ab1ead5c8a80bf8a29598a614021a;origin=https://github.com/PaulTorrillo/Microbiome_Reversions;visit=swh:1:snp:c07f8832389fa1e5285bde66070bc62e19c69633;anchor=swh:1:rev:70cb1d1484e57a0aa5cfdd719c705515d21e818e
(2010) Rate and effects of spontaneous mutations that affect fitness in mutator Escherichia coli
Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 365:1177–1186.

https://doi.org/10.1098/rstb.2009.0287
- PubMed
- Google Scholar
1. Vigué L
2. Tenaillon O
(2023) Predicting the effect of mutations to investigate recent events of selection across 60,472 Escherichia coli strains
PNAS 120:e2304177120.

https://doi.org/10.1073/pnas.2304177120
- PubMed
- Google Scholar
1. Wang RJ
2. Al-Saffar SI
3. Rogers J
4. Hahn MW
(2023) Human generation times across the past 250,000 years
Science Advances 9:eabm7047.

https://doi.org/10.1126/sciadv.abm7047
- PubMed
- Google Scholar
1. Yang Z
2. Nielsen R
(2000) Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models
Molecular Biology and Evolution 17:32–43.

https://doi.org/10.1093/oxfordjournals.molbev.a026236
- PubMed
- Google Scholar
1. Yang Z
(2007) PAML 4: phylogenetic analysis by maximum likelihood
Molecular Biology and Evolution 24:1586–1591.

https://doi.org/10.1093/molbev/msm088
- PubMed
- Google Scholar
1. Yu NY
2. Wagner JR
3. Laird MR
4. Melli G
5. Rey S
6. Lo R
7. Dao P
8. Sahinalp SC
9. Ester M
10. Foster LJ
11. Brinkman FSL
(2010) PSORTb 3.0: improved protein subcellular localization prediction with refined localization subcategories and predictive capabilities for all prokaryotes
Bioinformatics 26:1608–1615.

https://doi.org/10.1093/bioinformatics/btq249
- PubMed
- Google Scholar
1. Zhao S
2. Lieberman TD
3. Poyet M
4. Kauffman KM
5. Gibbons SM
6. Groussin M
7. Xavier RJ
8. Alm EJ
(2019) Adaptive evolution within gut microbiomes of healthy people
Cell Host & Microbe 25:656–667.

https://doi.org/10.1016/j.chom.2019.03.007
- PubMed
- Google Scholar
Preprint
1. Zhao S
2. Dai CL
3. Evans ED
4. Lu Z
5. Alm EJ
(2020) Tracking Strains Predicts Personal Microbiomes and Reveals Recent Adaptive Evolution
bioRxiv.

https://doi.org/10.1101/2020.09.14.296970
- Google Scholar

Article and author information

Author details

Paul A Torrillo
1. Institute for Medical Engineering and Sciences, Massachusetts Institute of Technology, Cambridge, United States
2. Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, United States
Contribution
Conceptualization, Software, Formal analysis, Investigation, Visualization, Methodology, Writing – original draft

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-4618-6061
Tami D Lieberman
1. Institute for Medical Engineering and Sciences, Massachusetts Institute of Technology, Cambridge, United States
2. Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, Cambridge, United States
3. Broad Institute of MIT and Harvard, Cambridge, United States
4. Ragon Institute of MGH, MIT and Harvard, Cambridge, United States
Contribution
Conceptualization, Supervision, Funding acquisition, Investigation, Writing – original draft, Project administration

For correspondence
tami@mit.edu

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-5430-3937

Funding

National Institutes of Health (1DP2GM140922-01)

Tami D Lieberman

National Science Foundation (Graduate Research Fellowship Program)

Paul A Torrillo

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Daniel Fisher, two anonymous reviewers, Benjamin Good, Erik van Nimwegen, and all members of the Lieberman Lab for their thoughtful feedback on this manuscript. We also thank William Shoemaker for making the data used in this work easily accessible and for his feedback on the manuscript. This work was funded by a grant from the National Institutes of Health (1DP2GM140922-01 to TDL) and a fellowship for the National Sciences Foundation (to PAT).

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Paul A Torrillo
Tami D Lieberman

(2024)

Reversions mask the contribution of adaptive evolution in microbiomes

eLife 13:e93146.

https://doi.org/10.7554/eLife.93146

Share this article

Cite this article

The previously proposed explanation for the time dependence of dN/dS is weak purifying selection.

Models of extremely weak purifying selection that can fit the data suffer from mutation accumulation and fitness decay.

Locally adaptive mutations and subsequent reversions can explain the decay of nonsynonymous mutations.

Under a model of reversion, the apparent dN/dS on long timescales underestimates the extent of adaptive evolution.

Author details

Paul A Torrillo

Contribution

Competing interests

Tami D Lieberman

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading

The previously proposed explanation for the time dependence of d_N/d_S is weak purifying selection.

Under a model of reversion, the apparent d_N/d_S on long timescales underestimates the extent of adaptive evolution.