The ratio of nonsynonymous to synonymous substitutions (dN/dS) encodes important information about the selection pressures acting on protein-coding genes. In bacterial populations, dN/dS often declines with the sequence divergence between strains, but the mechanisms responsible for this broad empirical trend are still debated. Existing models have primarily focused on de novo mutations, overlooking the older genetic variants that are continually introduced through horizontal gene transfer and recombination. Here we introduce a phenomenological model of dN/dS in recombining populations of bacteria, which allows us to disentangle the effects of recombination among pairs of closely related strains. We find that clonally inherited regions of the genome exhibit consistently higher dN/dS ratios, and that the accumulation of recombined segments can quantitatively explain the majority of the decline in dN/dS. We use these observations to re-examine models of purifying selection and adaptive reversion in human gut bacteria, and uncover evidence for widespread weak selection at a large fraction of protein-coding sites. Our findings show that horizontal gene transfer can be an important factor in shaping genome-wide patterns of selective constraint, and raise new questions about the effectiveness of natural selection in complex bacterial populations.

Quantifying the impact of homologous recombination on dN/dS dynamics

To illustrate how horizontal gene transfer and recombination can alter the temporal dynamics of dN/dS, we consider a simple conceptual model based on the genealogies of the underlying strains ([33,34,36,37]; Fig 1). We imagine a genome-wide comparison between two conspecific strains that descended from a single cell division event exactly T generations in the past (Fig 1A). The genetic differences between these bacterial genomes will arise from two primary sources: (i) de novo mutations that are inherited through clonal descent, and (ii) pre-existing variants that are introduced via homologous recombination with other strains (Fig 1B).

Under the assumption that synonymous mutations are effectively neutral, de novo mutations will generate single nucleotide differences (SNVs) at a characteristic rate of , where is the corresponding mutation rate in units of mutations per site per generation. In contrast, the synonymous divergence introduced by recombination will depend on the realized rate of recombination per site, R (which can in principle depend on T), as well as the genetic divergence of the imported segments. Previous work has shown that many of these recombined DNA fragments originate from unrelated donor strains [33,34,36,37], which carry an average synonymous divergence of , where is the average coalescence time within the species (Fig 1A). When , the introduction of SNVs through recombination can substantially inflate relative to the baseline expectation from de novo mutations alone. For closely related strains (), these recombination events will be visible as blocks of elevated sequence divergence against a nearly clonal backbone (Fig 1B). This hybrid divergence process continues until the majority of the clonal backbone has been overwritten by homologous recombination (), after which rapidly converges to the steady-state value (Fig 1C).

This phenomenological model can also be extended to incorporate non-synonymous SNVs. In clonal regions of the genome, non-synonymous divergence accumulates at a time-dependent rate , whose detailed shape will depend on the underlying model of selection (described in more detail in the following section). Recombined segments, by contrast, will carry an average non-synonymous divergence of , which can be estimated from the standing genetic diversity between unrelated pairs of strains (i.e., those with ). When is a non-linear function of time, the value of dN/dS in the clonal regions of the genome will generally differ from the recombined segments. This reflects the fact that mutations on recombined segments have already been circulating in the population for a longer period of time, providing more opportunities for natural selection to act (Figs 1A and 1B). As genomes accumulate more recombination events, the growing contribution from the recombined segments can alter the genome-wide average of dN/dS, similar to how they inflate the overall genetic divergence (Fig 1C). Notably, this recombination process can generate a smoothly declining trend in the genome-wide dN/dS over time, even if the dN/dS ratios are largely constant within both the clonal and recombined regions.

To test whether this hypothesis could explain the dN/dS dynamics in empirical data, we leveraged a previously curated dataset of commensal human gut bacteria [28], containing 5,416 reconstructed strain genotypes from 43 prevalent species (Methods Section 1). The genetic distances between these conspecific strains span nearly three orders of magnitude, from a typical value of to just a handful of detected SNVs (); the corresponding levels vary over an order of magnitude as well (Fig 2B). In a recent study [34], we used a hidden Markov model (HMM) approach to label individual recombination events and clonally inherited regions within the most closely related genomes in this dataset (Fig 2A). Since our HMM method relies on the spatial distribution of synonymous variants, these labels should be independent of the corresponding levels. The conceptual model in Fig 1 then leads to a clear prediction for two different sources of data: if recombination is the primary driver of the dN/dS decay in Fig 2B, then stratifying sites by their recombination status should reveal distinct dN/dS trajectories.

Consistent with this prediction, we found that the dN/dS dynamics in our dataset dramatically differ between clonal and recombined segments (Figs 2C and 2D). Within recombined regions, dN/dS is consistently low (), closely matching the values observed between unrelated pairs of strains. In contrast, dN/dS remains close to one within the clonal regions, even as genome-wide divergence increases by several orders of magnitude.

To demonstrate how these local patterns give rise to the full genome-wide trend in Fig 2B, we formulated a simple mathematical model of dN/dS accumulation that accounts for the growing contribution of recombined segments over time (Methods Section 3). This phenomenological model simply averages the observed dN/dS values in the clonal and recombined regions, and weights them by the empirically determined fraction of the genome covered by recombination at each divergence level dS (Fig A in S1 Text). As a result, once these synonymous divergence dynamics are specified, the model has no additional free parameters that require fitting to the observed dN/dS trend.

Despite its simplicity, we found that this mixture model quantitatively reproduces the observed dN/dS decay across the full range of genome-wide divergence levels in Fig 2B. In particular, it captures the characteristic transition from high to low dN/dS values at intermediate levels of synonymous divergence. In previous dN/dS models, the location of this transition was related to the underlying timescales of the natural selection (e.g., the selection coefficient in the purifying selection model of Ref. [28], or the reversion rate in the adaptive reversion model in Ref. [30]). In our present model, this transition naturally emerges from the increasing genomic footprint of recombination, and occurs around the point when recombined segments begin to dominate the genome-wide values of dN and dS (Fig 1C).

We verified that these trends were not specific to metagenomic reconstructions in Ref. [28] by repeating the same analysis on a collection of isolate genomes from two of the most prevalent species in Fig 2, Phocaeicola vulgatus and Alistipes putredinis (Methods Section 9). We observed a similar stratification of dN/dS values on the clonal and recombined segments of these isolate genomes as well (Figs B and C in S1 Text), confirming that the trends in Fig 2 were not caused by mapping errors or other metagenomic artifacts.

We also explored whether these results extended beyond human gut bacteria by re-analyzing the dN/dS dynamics in a collection of 110 Staphylococcus aureus strains (Methods Section 8). S. aureus is a common bacterial pathogen with a distinct ecology and recombination dynamics compared to many gut commensals [38–42]. Although clinically isolated strains of S. aureus are largely clonal in their recent evolutionary history, three large recombination events have been previously reported in sequence types ST34, ST239, and ST582 [42–44]. We used these previously identified events to stratify the genome into clonal and recombined regions as above, and analyzed the corresponding dN/dS patterns in an analogous manner (Fig D in S1 Text). Once again, we found that dN/dS dynamics separate into distinct trends by genomic region, consistent with the predictions of our basic recombination model. These results suggest that recombination may be a general mechanism underlying the frequently observed temporal trend in dN/dS across diverse bacterial species.

Estimating the strength of purifying selection from recombination-corrected dN/dS curves

While the mixture model in Fig 2B can recapitulate the overall shape of the dN/dS curve, it still requires external estimates of the dN/dS levels in each region (Figs 2C and 2D). These levels are set by the selective forces that act on nonsynonymous variants as they circulate through the broader population, which ultimately produce the low dN/dS values on a typical recombined segment. To probe these dynamics without the confounding influence of recombination, we can use the labels in Fig 2 to construct a “corrected” dN/dS trajectory, using the local values of and in each genomic region.

The short-term dynamics of dN/dS can be assessed using the clonally inherited regions that are preserved between closely related genome pairs. Fig 2D shows that the dN/dS levels on the clonal backbone remain close to one even as the nominal dS values span more than three orders of magnitude (; Fig 2D). However, because these genome-wide dS levels are inflated by recombination (Fig A in S1 Text), they overestimate the true divergence time between closely related strains.

To more accurately capture the temporal dynamics within the clonal regions, it is helpful to re-plot the clonal dN/dS values as a function of their clonal synonymous divergence (Fig 3), which provides a better proxy for the elapsed divergence time between closely related pairs of genomes. The observed clonal divergences span a much narrower range (), which can help explain the relatively modest decay of dN/dS within these regions. However, the rate of the decay is still about 10-fold slower than the original data in Fig 2, highlighting the importance of correcting for recombination.

This global trend is also recapitulated at the level of individual species (Fig E in S1 Text), with most species showing only a modest decay from as increases. A notable exception is Alistipes putredinis, which has elevated values () at the shortest genetic distances and declines below one toward the end of the observable range. A closer examination of this species revealed that the elevated dN/dS values are driven by recombination events whose edges are dominated by nonsynonymous variants, which are prematurely truncated by our synonymous HMM detection method (Fig F in S1 Text). Modifying our algorithm to capture these missing edges (Methods Section 10) lowers the clonal dN/dS of A. putredinis back to (Fig G in S1 Text), in line with the other species. Apart from this example in A. putredinis, none of the other species showed consistently elevated dN/dS levels within the subset of sites that we are considering.

For more distantly related strains, recombination events eventually cover the entire genome (Fig 1C), erasing clonal segments and precluding direct measurement of the clonal dN/dS curve over intermediate timescales. However, we can still estimate the long-term dN/dS levels by examining typical pairs of strains, whose genomes consist entirely of recombined segments with an average synonymous divergence . While the spread in this distribution is no longer informative of the differences in the underlying divergence times, the average values still provide a point estimate of the dN/dS levels when . Combining this estimate with the clonal regions above, we arrive at a corrected dN/dS curve with a gap at intermediate times (Fig 3), corresponding to the transition from clonal to quasi-sexual evolution.

Having corrected for the confounding effects of recombination, we can use the residual dN/dS curves in Fig 3 to re-examine the support for different models of selection. We begin by considering the standard purifying selection null model, which assumes that nonsynonymous variants are progressively removed over time [22,28,30]. In this model, the clonal dN/dS curve follows the approximate form,

(1)

where is the underlying distribution of fitness costs of nonsynonymous variants (Methods Section 4). The key features of the data can be captured by a two-class version of this model, where a fraction of the sites are truly neutral (s = 0) and the rest have the same characteristic fitness cost s₁.

Fitting this model to the data in Fig 3 yields a best-fit value of and a corresponding neutral fraction of (Methods Section 4; Figs H–J in S1 Text). Similar values are obtained when A. putredinis is excluded (Fig K in S1 Text), or when the recombined regions are extended to account for nonsynonymous-rich edges (Fig L in S1 Text). These inferred fitness costs are about 10-fold smaller than previous estimates [28–30] obtained from the uncorrected dN/dS curves in Fig 2, highlighting the importance of accounting for recombination. For typical bacterial mutation rates ( per site per generation), the associated selection coefficients are remarkably weak (), and would be difficult to observe on normal experimental timescales. Nevertheless, Fig 3 shows they are still sufficiently costly that they are eliminated by the time that T approaches .

These weak fitness costs are not an artifact of the two-class model in Fig 3. Adding a third class of mutations only marginally improves the fit (Fig M in S1 Text), and the estimated target size of these mutations is comparatively small (<10%). A more general bound on the strength of selection can be obtained by leveraging the functional form of Eq. 1: the fact that dN/dS remains close to one for clonal divergences as large as implies that more than 90% of all nonsynonymous mutations must have fitness costs less than (Methods Section 4; Fig N in S1 Text). At the same time, the fact that dN/dS eventually decays to implies that most of these mutations are eventually purged, which requires that . Taken together, these twin constraints imply that the purifying selection model in Eq. 1 can only explain the data if the fitness costs are weak but pervasive, with only a small minority of sites experiencing very strong () or nearly neutral () levels of selective constraint.

Mutation-type-specific dN/dS dynamics reveal stronger selection against nonsense variants

In the analysis above, all nonsynonymous mutations were aggregated together to compute the temporal dynamics of dN/dS. However, nonsense mutations, which introduce premature stop codons in the gene and often result in a loss of function, are generally expected to be subject to stronger selective pressures than missense mutations. If this is the case, then the model in Eq. 1 predicts that their clonal dN/dS trajectories should diverge from their missense mutation counterparts.

We tested this hypothesis by computing separate dN/dS curves for nonsense and missense mutations, respectively, and examined how each curve varies as a function of the clonal synonymous divergence (Fig 4 and Figs P and Q in S1 Text). The dN/dS decay of missense mutations closely resembled the aggregate results above, as expected given that they constitute the vast majority of nonsynonymous sites. In contrast, nonsense mutations showed a noticeably steeper decline in dN/dS with increasing divergence , reaching a long-term value of . This indicates that nonsense mutations are under stronger purifying selection, with a long-term neutral fraction of only .

To quantify the underlying fitness costs, we fit these data to the same two- and three-class models that we used for the aggregate data above. The missense fits recovered nearly the same estimate as the aggregate data in Fig 3. By contrast, we found that the nonsense mutations were better captured by the three-class model, which can account for at short times (Fig 4; Methods Section 4). The inferred parameters suggest that roughly 30% of nonsense sites fall into a strongly deleterious class with , which are effectively purged before the earliest observable divergence times in Fig 4 (). Given that the total target size for nonsense mutations is about 20-fold lower than missense mutations, these strongly deleterious variants correspond to only a few percent of all nonsynonymous sites—consistent with our earlier estimates from the aggregate dynamics in Fig 3. Interestingly, for the remaining 60% of nonsense sites, the inferred selection coefficients are only about 4-fold stronger than for missense sites, suggesting that many nonsense mutations are only weakly deleterious as well. These results demonstrate that stratifying dN/dS curves by mutation type can help resolve the heterogeneity in selection pressures across the genome: while most nonsynonymous mutations experience weak purifying selection, a subset of nonsense mutations are subject to strong constraint and are rapidly eliminated from the population.

Comparisons to alternative models of adaptive reversion

While the purifying selection model captures much of the observed dN/dS dynamics within the clonal regions, it is not the only possible explanation. One intriguing alternative was recently proposed in Ref. [30], in which the nonsynonymous differences could instead arise from locally adaptive substitutions that revert on longer timescales when the local environment changes. Previous work has shown that this adaptive reversion model can produce clonal dN/dS trajectories with the same functional form as Eq. 1, despite arising from a different set of biological assumptions and parameter meanings ([30]; see Methods Section 6 for a full derivation). To understand the implications of this alternative scenario, we sought to determine how the new recombination-corrected measurements in Figs 3 and 4 further constrain the underlying parameters of the model.

For example, the speed of the decay in the adaptive reversion model is controlled by the reversion timescale , which is dominated by the time that it takes for the environment to shift to a state where the reverse mutation is preferred ([30]; Methods Section 6). Using this mapping, the speed of the decay in Fig 3 yields an associated reversion time of , which translates to years given current estimates of the mutation clock ( per site per year; [8,45]). This revised timescale is significantly longer than a human lifetime, which places important constraints on the environmental processes that could be driving the shifts in the adaptive landscape.

A second key feature of the adaptive reversion model in Ref. [30] is that it predicts that the vast majority (~90%) of the clonal nonsynonymous SNVs in Fig 3 must be reversible adaptive mutations. This large fraction of adaptive variants introduces an additional constraint linking the total number of adaptive sites and the degree of parallel evolution in other strains (Methods Section 6). For example, if local adaptation was driven by a modest number of reversible adaptive sites, one would expect to see the same mutations independently arise in multiple pairs of strains. We quantified this signature of parallel evolution by identifying all purely clonal strain pairs (i.e., those with no detected recombination events) and asking how often the observed SNV differences between them were also found in other strains within our larger cohort (Methods Section 7). Figure 5A shows an example of this calculation for Phocaeicola vulgatus, one of the most prevalent gut commensals in our dataset. The data show that nearly 90% of all clonal nonsynonymous SNVs were private to a single clonal lineage and were never observed elsewhere among the 120 sampled hosts. Similar patterns are observed across the other species in our dataset (Fig Q in S1 Text), with clonal SNVs rarely recurring in multiple strains (and at comparable rates between synonymous and nonsynonymous mutations). Stratifying SNVs by their recombination status is crucial for this calculation, since much higher rates of “parallelism” are observed when recombined sites are included (Fig 5B and Fig Q in S1 Text).

By extending the calculations in Ref. [30], we show that these low rates of parallel evolution can explain the dN/dS curves in Fig 3 only if more than 1% of all nonsynonymous sites in the genome—corresponding to tens of thousands of sites—are targets of reversible adaptive evolution (Methods Section 6). These numbers could be consistent with a smaller number of effective loci if mutations in the same locus interact with each other epistatically (e.g., loss-of-function mutations) [30]. However, a similar calculation shows that at least ~500 adaptive loci, each containing tens to hundreds of sites, are required to match the observed dN/dS curves in Fig 3 (Methods Section 6). This new estimate is an order of magnitude larger than the ~50 adaptive loci previously estimated in Ref. [30], and reflects the ~10-fold longer reversion timescale implied by the recombination-corrected data.

Finally, while the purifying selection model mechanically predicts that at short times, the adaptive reversion model can only match this feature of the data if the fraction of adaptive sites and their corresponding substitution rates are precisely tuned so that their product is of order (Methods Section 6). Moreover, this condition must separately hold within each species in our dataset, as well as for different classes of sites (e.g., the missense and nonsense mutations in Fig 4), where both the decay rates and fractions of adaptive sites can differ. The mechanisms that could give rise to such fine-tuning currently remain unclear.

Interestingly, the parallelism analysis in Fig 5 also reveals a small subset of clonal mutations that are observed at higher prevalence within the broader population. Some of these mutations likely correspond to shorter recombination events that were not detected by our HMM algorithm. However, we also note that the dN/dS ratios of these recurrent clonal SNVs are somewhat elevated relative to the population average (Fig 5B), which captures the long-term expectation for recombined sites. Although the small number of prevalent clonal SNVs limits statistical power, this relative enrichment of nonsynonymous variants raises the possibility that a small minority of clonal mutations may still correspond to adaptive reversions, consistent with prior direct observations [28,46,47]. Nevertheless, the small overall fraction of these recurrent clonal SNVs (Fig 5 and Fig Q in S1 Text) suggests they cannot be driving the large-scale temporal trends in Fig 3.

1 Metagenomic analysis of human gut commensals

The dN/dS measurements in the main text were obtained from collection of shotgun metagenomic data that have been analyzed in several previous works [28–30,34]. This dataset comprises of 932 fecal metagenomes from 693 individuals from North America, Europe, and China. All samples were processed using the same reference-based pipeline described in Ref. [28]. In brief, we used the MIDAS software package [62] to align the raw sequencing reads to a panel of reference genomes representing common gut bacterial species, and identified single nucleotide variants (SNVs) based on the raw read pileups in each sample. In this work, we restricted our attention to SNVs in core genes (those that were present in more than 90% of samples in which the species was reliably detected). We used the quasi-phasing approach in Ref. [28] to identify intra-host populations in which the genotype of the dominant lineage could be inferred with a high-degree of confidence. Each quasi-phased genotype is therefore a within-host consensus sequence rather than a collection of individual reads: an allele is assigned to a site only when it constitutes a strong majority (≥80%) of the reads covering that site in a given host, while sites at intermediate within-host frequency (or insufficient coverage) are treated as missing data. Applying this method to the entire dataset yielded a collection of 5,416 strain genotypes from 43 species. Filtered genotype and coverage matrices for these quasi-phased genomes are available at Zenodo (https://doi.org/10.5281/zenodo.14853785).

To identify synonymous and nonsynonymous variants, we focused on the subset of SNVs that were biallelic across the entire collection of quasi-phased genomes, ensuring that each site had a single putative ancestral allele and a single derived allele. Under the assumption that the major allele is ancestral, we used the reference genome sequence to determine site degeneracy (e.g., non-degenerate or 4-fold) and to classify each mutation as synonymous, missense, or nonsense. Since every genotype is defined on this common reference coordinate system (rather than through de novo genome alignments) each site is assigned to a fixed reference codon and reading frame, so that frame-shifting indels or alignment ambiguities elsewhere in the gene should not alter the nonsynonymous mutation count. These site-specific annotations were then used to tally the number of observed differences and opportunities in each mutation class.

Note on species names. The analyses in this study were performed using version 1.2 of the MIDAS reference genome database [62]. Subsequent taxonomic updates have divided the Bacteroides genus into Bacteroides and Phocaeicola [63], leading to changes in species names relative to the original MIDAS database. For consistency, we retain the MIDAS species names in supplementary figures, while adopting updated nomenclature in the main text where appropriate. The corresponding name translations are listed below for reference:

2 Estimating dN/dS ratios from pairwise comparisons

To estimate the synonymous and nonsynonymous substitution rates between a given pair of genomes, we used the observed sequence divergence at 4-fold‐degenerate (4D; synonymous) and non-degenerate (1D; nonsynonymous) sites. These divergences were estimated from the ratio of the total number of SNV differences () and the total number of covered sites () for each class of mutations. Since our goal is to examine how varies as a function of , it is important to avoid the mechanical correlation that arises when the sampling noise in is shared between both the input and output variables. Following the Poisson thinning approach in Garud and colleagues [28], we randomly partitioned the synonymous sites into two disjoint sets of sizes and , with corresponding difference counts and . This yields a pair of synonymous divergence estimates for each a pair of strains,

(2)

which should be statistically independent of each other after conditioning on the underlying divergence time. We therefore computed using in the denominator,

(3)

and plotted this ratio as a function of the independent variable .

Recombination-corrected dN/dS values. When comparing closely related genome pairs from the gut commensal dataset, we further stratified sites by recombination status to compare between clonal and recombined regions. We used the recombination tracts that were previously inferred for these same pairs of strains in Ref. [34], using a pairwise hidden Markov model that leverages the spatial profile of synonymous sequence divergence across the genome. Closely related genome pairs were defined empirically as those sharing more than 50% identical genome blocks, where each block corresponds to 1,000 4-fold‐degenerate sites. For each such pair, we used the inferred tract boundaries to assign each site to either the clonal or recombined class, and then computed , , , and separately for each class before applying the Poisson thinning procedure described above.

Computing average dN/dS curves. While the above procedure provides for individual genome pairs, we also sought to compute the average value for all pairs of strains with a similar estimated divergence time. To do so, we first performed our Poisson thinning procedure on the clonal synonymous sites of each genome pair and used to assign pairs into logarithmically spaced divergence bins. Within each bin, we pooled all synonymous and nonsynonymous sites across pairs and computed a “ratio‐of‐means” estimate of the average within the bin:

(4)

Confidence intervals were obtained from 600 bootstrap replicates, resampling both genome pairs and Poisson partitions. The same procedure was used to estimate the long‐term for each species by analyzing a sample of 20 unrelated, fully recombined genome pairs.

Estimating mutation-type specific dN/dS curves. Finally, to compare missense and nonsense substitution rates within clonal regions, we annotated the mutation class of each 1D SNV difference and calculated the total number of missense () and nonsense () opportunities, where . We then repeated the Poisson thinning and bootstrap procedures separately for each mutation class. In all cases, was calculated using the synonymous substitution rate per synonymous opportunity, while the x-axis in the relevant figures (e.g., Fig 4) shows the synonymous substitution rate per synonymous site.

We note that partitioning 1D differences into missense and nonsense classes slightly changes the neutral per–opportunity mutation rate because the underlying transition:transversion (Ti:Tv) mix differs across these opportunity sets. In the standard genetic code, the Ti:Tv composition of missense opportunities is nearly identical to that of 4D sites, and therefore the expected neutral mutation rate per opportunity is effectively the same. Nonsense opportunities, however, are modestly enriched for transversions, with a Ti fraction of roughly 0.27 based on the genetic code. Assuming a typical bacterial rate bias of [64], this composition would yield a neutral mutation rate per nonsense opportunity that is about 10% lower than at 4D sites. Thus, while this mechanism can produce a nonsense‐only modestly less than one even under neutrality, values much lower than this expectation are indicative of the effect of selection.

3 Modeling the dynamics of genome-wide dN/dS with recombination

The phenomenological model in Fig 2 assumes that each pair of bacterial genomes contains a mixture of clonally inherited and recently recombined genomic regions. If we let T denote the time since the two strains last shared a common ancestor, then their synonymous and nonsynonymous divergence can be written as a weighted average,

(5a)(5b)

where denotes the fraction of the genome covered by recombination events, and and are the synonymous and nonsynonymous divergences in the clonal (i = c) and recombined (i = r) regions of the genome, respectively. Rearranging these equations yields an implicit expression for the dN/dS ratio,

(6)

as a function of T and .

We assume that in the clonal regions, sequence divergence arises solely through mutation, so the expected synonymous divergence at time T is given by

(7a)

where is the per-site mutation rate. In the simplest version of the model, we assume that the nonsynonymous clonal divergence is effectively neutral, so that

(7b)

though we also consider more complex generalizations in Sections 4 and 6 below. In contrast, the recombined regions of the pairwise alignment are assumed to originate from unrelated donors in the population. As a result, the sequence divergence within the recombined segments will reflect the long-term evolutionary distance between unrelated strains. In the simplest version of the model, we assume that the recombined divergences and are drawn from a fixed distribution with means and , respectively, which can be independently measured from the larger dataset using the typical divergence between unrelated pairs of genomes.

Finally, the fraction of the genome that is covered by recombination events () will depend on divergence time T between each pair of strains. In our previous work [34], we showed that recombination dynamics in gut bacteria can be both heterogeneous across strains and nonlinear in time. To account for this complexity, we adopted an empirical approach to estimate the functional form of . For example, the expected value can be written as

(8)

where R(T) is the average effective recombination rate since the two strains last shared a common ancestor. In this work, we do not make any assumptions about the functional form of R(T), but note that it can be estimated empirically from the synonymous divergence data by fitting as a function of .

However, the measured dN/dS curve in Fig 2B does not correspond to the conditional distribution of dN/dS at a fixed divergence time T, but rather to the related conditional distribution,

(9)

at a fixed value of dS. These two conditional distributions will generally differ from each other because is a random variable; this implies that the same value of dS can be represented by different combinations of and T.

To calculate the dS-conditioned distribution in Eq. 9, we can leverage the implicit relation in Eq. 6. We will assume that the variation in this expression will be nominated by the noise in T and , which should be a good approximation whenever the genome is sufficiently long. If we also assume that the recombined fraction of the genome is small (), which is the relevant regime for Fig 2B, then the expected value of Eq. 6 reduces to

(10)

where the dependence on is encapsulated in the conditional average

(11)

This expression shows that the functional form of will generally depend on the noise in , in addition to its mean value in Eq. 8. It will also generally depend on the underlying variation in T for the close pairs in the focal species of interest.

Since the full distributions of T and are poorly parameterized by existing data, we used an nonparametric approach to fit empirically from the synonymous divergence dynamics inferred by our HMM. In particular, using the inferred values of and for each pair of closely related strains, we estimated and in sliding windows of dS, and fit the resulting trajectory to a spline curve (Fig 3). We then substituted this function into the analogous version of Eq. 10 to plot the theoretical predictions in Fig 2B:

(12)

We used this median estimator rather than the mean in Eq. 10 to better capture the central tendency of the dN/dS curve on a logarithmic scale.

In contrast to other models [28,30] that have free parameters that must be estimated from the observed curve itself (e.g., mutation rates, selection coefficients, and/or reversion rates), the phenomenological model in Eq. 12 only depends on quantities that can be independently measured from other features of the data: once , , and are specified, there are no additional free parameters required to produce the theoretical curve in Fig 2B.

However, while this phenomenological recombination model provides an explanation for the observed decay of genome-wide , it does not explicitly model the dynamics of within the clonal regions, or the values of and on the recombined segments. In the following two sections, we consider two classes of models that generate these predictions from different underlying mechanisms.

4 Purifying selection model

Multiple formulations of the purifying selection model describe how changes as a function of synonymous divergence () at a locus evolving without recombination [22,28,30]. These models converge on the common functional form in Eq. 1, which captures the impact of selection against nonsynonymous mutations with a distribution of deleterious fitness effects (DFE) . The shape of this curve is governed by the ratio , which sets the timescale over which a mutation of effect s is purged. For an arbitrary distribution , sites with different selection strengths are purged at different ranges of , and the curve reflects their relative contributions. Since the actual DFE is unknown, we used two different strategies to constrain this parameter from the observed data.

Estimating parameters for multi-class DFEs. Following previous work, we first considered discretized versions of the DFE, where nonsynonymous mutations fall into one of K discrete fitness classes with characteristic costs and relative target sizes . In the simplest version of this model with K = 2 classes, a fraction of nonsynonymous sites are truly neutral (s₀ = 0), and the remaining fraction share the same characteristic fitness cost s₁. In this case, the curve reduces to

(13)

which depends on the two compound parameters and . At long times (), asymptotically approaches , which can be interpreted as the equilibrium fraction of neutral sites. At short timescales (), the second term approaches , and thus regardless of s₁, reflecting the fact that deleterious mutations have not yet been purged.

We also considered three-class versions of the model,

(14)

which include an additional fitness class with . This three-class formulation can capture cases like Fig 4B where the short-time is significantly below 1, indicating that a substantial fraction of nonsynonymous mutations are strongly deleterious.

We fit Eqs. 13 and 14 to the binned data using weighted nonlinear least squares. As described in Methods Section 2, we used the Poisson–thinning procedure to partition the synonymous sites into two disjoint sets, “A” and “B”, with used for bin assignment and for the calculation. For each divergence bin i, let and denote the total numbers of nonsynonymous and synonymous differences, respectively, obtained by pooling counts across all genome pairs in the bin for the “B” set of sites, and and the corresponding pooled site counts. The binned estimate is then

(15)

To account for varying statistical uncertainty across bins, we used inverse–variance weights,

(16)

on a log-transformed version of , which down-weights bins with few observed differences. We excluded divergence bins containing fewer than three genome pairs to ensure that there was sufficient data for fitting. The resulting objective function is then defined by

(17)

where is the predicted from the chosen model (two-, or three-class) with parameters . For all model variants, we fixed the neutral-site fraction to the long-term value for each species, estimated from unrelated genome pairs (Methods Section 1).

We quantified the uncertainty in our estimates by bootstrap resampling the genome pairs, repeating the Poisson partitioning of the synonymous sites, and re-fitting the model to each bootstrap replicate. The 2.5%, 50%, and 97.5% quantiles of the bootstrap parameter distribution were used to define the corresponding confidence interval.

General constraints on the DFE from the asymptotic behavior of dN/dS. For an arbitrary DFE, the functional form in Eq. 1 still imposes strict constraints on the fraction of strongly deleterious mutations from the short-time behavior of . A robust empirical observation across most of the species in our dataset is that the clonal curve remains close to one for as large as . Intuitively, this implies that strongly deleterious sites cannot be too abundant: if many mutations were purged almost immediately, would already have dropped well below one at these short times.

We can formalize this argument by partitioning the integral in Eq. 1 at a user-specified threshold , and then bounding the integrand in each term by its value at the lower limit of integration:

(18)

Rearranging this expression yields an upper bound on the fraction of mutations that are stronger than :

(19)

As a concrete example, Fig 3 shows that for . Eq. 19 then implies that more than 90% of nonsynonymous mutations must have and more than half must have , independent of the detailed form of the DFE.

An analogous argument for the large limit of the dN/dS curve yields a lower bound on the fraction of mutations that are more costly than a given value:

(20)

Applying this bound to the right-hand side of the curve in Fig 3 shows that more than 80% of all nonsynonymous mutations must have , and more than half must have . Together, these two conditions provide a general constraint linking the early- and late-time behavior of the curve to the underlying scale of the DFE.

Assessing the spread of clonal dN/dS values. To investigate whether the scatter of the clonal values around the fitted purifying-selection curve exceeds what finite-site sampling alone would produce, we compared the observed spread to a matched null simulation in Figs I and J in S1 Text. For each strain pair, we retained the observed covered 4D and 1D site counts and synonymous differences, and reused the same Poisson-thinning draw that generated the original figure (so the axis, the denominator, and the thinning noise are identical to the data). The only simulated quantity was the nonsynonymous numerator, , which we drew from a Poisson distribution with mean , with taken from the global fitted curve and no additional pair-to-pair biological variation. Pooling across pairs, the observed spread of within bins is only – wider than this sampling null, and the excess within individual species is similarly modest (–) and present across essentially all species. This suggests that most of the apparent scatter can be attributed to finite-site sampling of the small synonymous denominator. The residual overdispersion (largest at higher , where sampling noise is smallest) may reflect a small contribution from undetected recombination events or other sources of heterogeneity across sampled pairs.

5 Forward-time simulations of purifying selection and genetic linkage

The purifying selection model in Eq. 1 is best viewed as a constitutive relation, which connects the expected values of to the observed values of for different pairs of strains. Deriving this relation from first principles—along with the corresponding distributions of and s—requires additional assumptions about the underlying mechanics of intra- and inter-host competition. While a thorough investigation of this connection is beyond the scope of our present study, we investigated the behavior of our fitted models in two commonly considered limits.

Pure intra-host selection. Following Ref [30]., we first considered the extreme limit where there is no inter-host selection, such that each host’s microbiome can be treated as its own independently evolving population. In this case, the values of and in Fig 3 correspond to fixed differences between pairs of independently evolving populations that diverged from each other roughly generations ago. The fact that and correspond to fixed differences (rather than differences between randomly sampled clones in each host) is an important feature of the metagenomic data underlying Fig 3. As we will see below, this empirical constraint will make it difficult for any purely intra-host model of purifying selection to generate a decaying clonal dN/dS curve.

We simulated intra-host evolution using a minimal version of the purifying selection model in Fig 3. Specifically, we considered a Wright-Fisher population of N competing cells with binary genomes of length . We assumed that deleterious mutations could occur at of these sites, with a constant additive fitness cost s drawn from the two-class DFE in Fig 3:

(21)

We also allowed for recombination, where contiguous stretches of sites were transferred into the recipient genome from a randomly selected donor cell. We considered two different sources of recombination:

1. Intra-population transfers, which occur at a total rate per genome per generation, and sample a donor from one of the other cells in the same host.
2. Inter-population transfers, which occur at a total rate , and involve a donor cell from another randomly sampled host.

Based on the empirical distribution of recombination tracts detected in Methods Section 2, we simulated a characteristic transfer length of [34]. The corresponding rates of recombination are more difficult to constrain from the available empirical data. In the absence of selection, the analysis in Ref. [34] suggests a characteristic value of , based on the accumulation of successful transfers between hosts. However, when either the donor or recipient DNA sequences have the potential to carry selected mutations, the link between the observed rate of accumulation and becomes more complex (see the discussion in Ref. [34]). In this case, the more general constraint is that the realized rates of transfer should match the empirically observed rates in Ref. [34], rather than the underlying rate . Moreover, the external transfers in Ref. [34] and Methods Section 2 are not directly informative of , which could occur at a much higher rate than while remaining invisible to our hidden Markov approach. For this reason, we simulated a range of possible values loosely guided by the empirical estimates above (). For these parameters, the relative rate that a given site is involved in a recombination event versus a de novo mutation is given by the ratio , which ranges from 10³ to for the parameter values above; this includes cases where recombination is both much slower and much faster than selection ().

We simulated the resulting model across a range of scaled population sizes () informed by empirical estimates of the rates of genetic drift in the human gut microbiome (; [65–67]). To enhance computational efficiency, we used standard diffusion rescaling [68] to simulate smaller populations with the same compound parameters , N s₁, and . For each set of parameters, we measured the average pairwise intra-host , based on the diversity at a collection of linked neutral sites. We also measured scaled rate of fixation of deleterious mutations over an extended interval of length following an initial burn-in period of length .

Our simulation results partially agree with those previously reported in Ref. [30]. In the absence of homologous recombination (), the large number of weakly selected mutations in Fig 3 generate substantial selective interference [55,56], which reduces the efficiency of natural selection and allows many deleterious mutations to fix [; Fig R in S1 Text]. Since we have assumed that there is no inter-host selection, these deleterious hitchhikers become invisible to natural selection after they fix, and can only be repaired by subsequent reversions. The dN/dS curves in this scenario can therefore be described by the reversion model in Eq. 28 below, which reduces to

(22)

in the special case where , and . This implies that when , dN/dS curve will start at an value, and then decay on a divergence scale that is dictated by the pace of successful reversions []. Previous work [56] has shown that the reversion probability in this class of models is bounded above by , so that . When , this divergence scale is much larger than the values that are observed in the empirical data in Fig 3, which implies that dN/dS will remain close to its initial value and not decay much with time.

Quantitatively similar results arise when one allows for modest rates of internal recombination (). However, for larger internal recombination rates (), the scaling with N qualitatively changes, such that deleterious mutations can no longer appreciably fix at the largest intra-host population sizes []. Assuming that these intra-host transfers are still classified as clonal (since Fig R in S1 Text shows that they are much less diverged than the typical inter-host distances), Eq. 22 implies that will effectively start near its long term value , which is now much less than one.

We note that in both of these cases, the clonal dN/dS curves do not appreciably decay with time. This is ultimately a consequence of the assumption that and are driven by intra-host fixations, which implies that the decay in can only be driven by the fixation of successful reversions. However, in the absence of fluctuating selection pressures (Methods Section 6), the rates of successful reversions are strongly constrained by the corresponding deleterious fixation rate [56]. Nearly neutral hitchhiking [] implies a correspondingly low reversion rate [], while rapid reversions require that the original rate of fixation is low (). These constraints combine to produce a roughly constant dN/dS curve, regardless of the strength of selection (). This suggests that a model of purifying selection in isolated hosts cannot capture the empirical trends in Fig 3.

Allowing for realistic rates of external recombination () can alter this picture slightly, since reversions can now be seeded by recombination events, rather than the slower process of mutation. However, since these external transfers by definition do not contribute to and , they cannot produce a decaying curve. The plausibility of this model therefore strongly depends on how well one can detect the relevant external recombination tracts. If most external transfers involve longer segments of typical divergence (), then they should be efficiently detected by our Hidden Markov method, and the decaying curves in Figs 3 and 4 would be difficult to explain in this model. On the other hand, if we assume that there is also a comparably large influx of shorter or less diverged transfers (), then these would be missed by our current methods, and could contribute to the declines in Figs 3 and 4. Further efforts to investigate this possibility would be an interesting avenue for future work.

Pure inter-host selection. Following Ref. [58], we also considered the opposite limit of the one above, where each host is represented by a single dominant strain and all competition occurs at the between-host scale. In this case, and correspond to polymorphic mutations separating two individuals in the same population. We simulated this limit using an analogous version of the model above, incorporating deleterious mutations, recombination, and genetic drift.

In this limit, the population size N has a different interpretation than in the intra-host model above. It now refers to the global population of competing lineages, which should be proportional to the total the number of effective hosts (). However, we must also account for differences in the effective generation time. If lineages are replaced every divisions, then the rates of mutation, recombination, and selection should be scaled up by a corresponding factor of . Using diffusion rescaling, this is equivalent to assuming a larger effective population size . For the human population, plausible values of might range from , while might range from from 10 to 100 years. Thus, plausible values of the compound parameter might range from 10^–2 to 10⁴.

In this model, the average synonymous diversity among circulating strains () places additional constraints on via the relation . However, this is only an indirect constraint, since may be much smaller than when selective interference is widespread [55]. We therefore simulated a range of plausible values to see which—if any—could produce values that match the empirical data in Fig 3 (). With the appropriate parameter rescalings, these inter-host simulations fall in the same regimes as those performed for the intra-host model above (Fig R in S1 Text).

These results confirm that a pure inter-host competition model also leads to strong interference selection, with most of the simulated parameters producing values far below the empirically observed levels in Fig 3. While it was computationally infeasible to simulate the upper range of plausible values above, extrapolating from the observed trends suggests that enormous values of are likely necessary to approach the empirically observed levels—particularly if the underlying rate of recombination is assumed to be equal to the rate of apparent accumulation (). These results suggest that further work is critically needed to understand how the realized rates of transfer constrain the plausible values of R in the presence of negative selection. If the underlying rates remain close to , then this simple model of inter-host competition will struggle to explain the observed data. However, this simple model also neglects many factors (e.g., spatial structure, epistasis, and spatiotemporally varying selection pressures) that are potentially relevant for the empirical populations in Fig 3. Understanding how these additional complexities might alter the picture above remains an interesting avenue for future work.

6 Adaptive reversion model

We also compared the data in Figs 3 and 4 to a generalized version of the adaptive reversion model in Ref. [30]. This metapopulation model assumes that a subset of nonsynonymous sites undergo recurrent locally adaptive substitutions, which are later reversed when environmental conditions change.

Derivation of model. Following Ref. [30], we assume that these locally adaptive sites exhibit a form of modular epistasis. Each site is part of a larger module (or “locus” in Ref. [30]) with two states that we denote as “mutant” (M) and “wildtype” (WT). We assume that mutations at any of the sites in the module will switch it from the wildtype to the mutant state, but compensatory mutations (MWT) can only occur through a reversion of the original mutation. We also assume that there is strong intra-module epistasis, such that two or more mutations in the same module have the same phenotypic effect as a single mutation.

Substitutions and reversions within a single person’s microbiome are modeled as stochastic transitions with total rates and , respectively. These two effective parameters encapsulate both the rate of environmental switching and the time that it takes for the population to reach the new locally adapted state. In many cases of interest, this latter timescale may be negligible, and and can be interpreted as switching rate between environmental states where either the mutant or wildtype version of the module is preferred. Ref. [30] focused on a symmetric version of the model where ; for completeness, we consider general values of and to allow for cases where the mutant state is preferred in only a small fraction of microbiomes ().

Under these assumptions, the probability that a given module is in its mutant state at time t follows a simple Markov process,

(23)

where is the steady-state probability of finding the module in the mutant state. However, since dN/dS is a site-based metric, we must also understand the probability that a given site within the module will be found in its mutant or wildtype state. Under similar assumptions as above, we can model the dynamics at a given focal site as a three-state Markov process, corresponding to cases where (i) neither the module or the focal site are mutated (WT,WT); (ii) both the focal site and module are mutated (M,M); and (iii) cases where the module is mutated at another site within the module, while the focal site retains its wildtype allele (WT,M). The probabilities of finding the site in each of the three states evolve as

(24)

where is an additional parameter that describes the relative probability that a locally adapted module is driven by a mutation at the focal site. The solution to Eq. 24 is given by

(25)

where

(26)

are the steady-state probabilities of finding the site in a given state. This expression reduces to Eq. 23 when each site comprises its own module (), while allowing for more complex dynamics when .

To calculate the expected genetic divergence between two independently evolving microbiomes, we assume that fixation within each microbiome is fast relative to the divergence time and that the shared ancestral allele is drawn from the steady-state multinomial distribution in Eq. 26. The expected probability of observing a difference at a given site is given by

(27)(28)

In our parameter regime of interest where each site makes a small contribution to its module (), this further reduces to

(29)

Summing over all nonsynonymous sites within a given focal set yields a corresponding prediction for the ratio,

(30)

where is the aggregate distribution of , , and across all of the sites within the focal set. This generalizes the derivation in Ref. [30] to allow for variable on and off rates, and applies to arbitrary subsets of nonsynonymous sites, in addition to the total genome-wide value. Notably, the functional form of Eq. 30 is identical to the purifying selection model in Eq. 1, even though the two expressions emerge from distinct evolutionary mechanisms.

Three-class model. To understand the constraints that this equivalence places on the underlying parameters of Eq. 30, we first considered a simplified version of the reversion model, similar to the one proposed in Ref. [30], where the nonsynonymous sites are divided into three discrete categories: (i) a set of effectively neutral sites with , (ii) a set of reversable adaptive loci with characteristic rates and , and (iii) a set of strong, permanently deleterious sites with and . Denoting the fraction of sites in each class by , , and , respectively, the curve in Eq. 30 reduces to

(31)

where

(32)

is the average value of across all of the sites within class (ii).

To more directly compare this expression with the one provided in Ref. [30], it is necessary to invoke some additional assumptions about the functional form of . We will assume that the probability that a given site will be the driver mutation of its module is proportional to its relative target size—a term we use in a generalized sense to refer to both its underlying mutation rate and any differences in its within-host fixation probability. Using this notation, we can express the value of for a given site i as

(33)

where denotes the relative target size for mutations at site i, m(i) denotes the module that site i belongs to, and denotes the total target size for adaptive mutations in module m. This last term includes the contributions from (i) other nonsynonymous sites in the focal set, (ii) any nonsynonymous sites outside the focal set (if is a subset of the genome), and (iii) any non-coding or non-point mutations (e.g., invertible promoters; [69]) that can contribute to the overall adaptation of the module but are not included in dN/dS.

Summing over all of the nonsynonymous sites in class (ii), the average value of can then be expressed as

(34)

where is the total number of modules, denotes the total target size of module m within the focal set, and the average in the angle brackets is now taken over all modules m that have at least one site in the focal set. Substituting this result into Eq. 31 yields a related expression for dN/dS:

(35)

which reduces to the expression in Ref. [30] in the special case where and .

In contrast to the purifying selection model in Eq. 1, all of the time-dependence in Eqs. 31 and 35 is generated by reversions of adaptive sites. This places strong constraints on the underlying parameters, which we can start to understand by considering the large and small limits of this expression.

Comparisons to data. Like the purifying selection model in Eq. 1, the long-time behavior of Eqs. 31 and 35 is governed by the fraction of effectively neutral nonsynonymous sites . In order to match the empirical observations of between unrelated strains, we must continue to have . Similarly, the timescale of the decay of suggests that the typical reversion rate must be on the order of , analogous to the effective selection coefficient inferred under the purifying selection model in Fig 3.

One important difference between the two pictures is that while the purifying selection model in Eq. 1 is mechanistically constrained to approach one at short times, the adaptive reversion model possesses a different short-time limit,

(36)

which functions as an additional fitting parameter. This extra degree of freedom allows the adaptive reversion model to produce dN/dS values that are substantially greater than one, which can be useful in cases where adaptive evolution is observed [8]. However, for empirical curves like Figs 3 and 4, where dN/dS remains comparatively close to one at short times, this additional degree of freedom must be precisely tuned in order to match this feature of the data. In particular, given the constraint that , matching the empirical values of at short times requires that the reversible component of Eq. 36 must also be order one, yielding a corresponding condition

(37)

The left-hand side of Eq. 37 is a product of several independent parameters (, , , and ) which can differ from each other by several orders of magnitude in real biological systems. Achieving a dN/dS value of order one requires that these different terms must precisely balance one another. Since there is not an a priori reason why these different evolutionary parameters should be related to each other in this way, we argue that the adaptive reversion model requires an additional degree of fine-tuning to match this feature of the data.

Notably, this fine-tuning constraint must hold for any collection of sites within the genome that exhibit at short times. For example, Fig 4 shows that Eq. 37 must separately hold within the subsets of missense and nonsense mutations, even though the underlying parameters could differ between them. Such agreement is not automatically guaranteed by the model and would require an additional degree of fine-tuning.

To illustrate this idea, we recall that the speed of the dN/dS decline in Fig 4 was 4x faster for nonsense versus missense mutations, which we interpreted as a 4-fold higher fitness cost of nonsense mutations under the purifying selection model. In the adaptive reversion picture, this same difference in the dN/dS curves implies that the environmental switching rates must be 4-fold faster for nonsense versus missense mutations, which would require in turn that the nonsense mutations must arise from a systematically different set of modules than their missense counterparts. However, in order for the condition in Eq. 37 to remain satisfied in this new setting, the larger value of must be balanced by a 4-fold reduction in the remaining terms (e.g., the number of modules per focal site or the relative on-vs-off rates) to ensure that at short times. At present, the mechanisms that would allow for such fine-tuning remain unclear.

Even stronger tests of the adaptive reversion model are possible if one can identify the subset of reversible adaptive loci. In this case, the remaining, non-adaptive portion of the genome should provide a direct readout of the effects of purifying selection, which is assumed to be very strong in the version of the model in Eqs. 31 and 35. While it is difficult to fully enumerate the set of adaptive loci, several previous studies have performed conceptually related analyses by leveraging recurrently mutated genes in collections of isolates sampled from individual hosts [8,70].

In one study in Bacteroides fragilis [8] (one of the species included in Fig 2), this subset of recurrently mutated genes exhibited a dN/dS ratio of 6.0, indicating strong positive selection [8]. However, the remaining portion of the genome retained a dN/dS indistinguishable from one, similar to the genome-wide patterns in Fig 3. A similar analysis in the pathogenic bacterium Burkholderia dolosa [70] observed a dN/dS ratio of 18 in the most recurrently mutated genes, and a residual dN/dS that was measurably less than one (dN/dS0.63), but still larger than the long-term value of estimated from Ref. [71]. In both species, these dN/dS values imply that either the vast majority of mutations in the non-parallel genes () must also be locally adaptive—and subject to the same fine-tuning constraints above—or that they must be eliminated on longer timescales through a purifying-selection-like process similar to Eq. 1. By contrast, no studies (to our knowledge) have identified a suitably large collection of sites with a short-term dN/dS 0.1, as predicted by the strong form of the adaptive reversion model in Eq. 35.

Together, these analysis show how comparing dN/dS curves for different classes of sites can provide additional information about the underlying selective forces that operate within a population. Extending this approach to other subsets of the genome would be an interesting avenue for future work.

Constraints from parallel evolution. An additional set of parameter constraints can be obtained by considering signatures of parallel evolution—i.e., the probability that the same adaptive mutations are independently sampled in unrelated lineages. In the model defined above, the degree of parallel evolution depends on the equilibrium frequencies in Eq. 26, and how they vary across the different sites in our focal set.

To make this intuition more concrete, we recall that the short-time limit of Eq. 36 implies that of all nonsynonymous substitutions at short times must be locally adaptive mutations. This suggests that any signatures of parallel evolution observed among these recent clonal mutations will provide additional information about the realized values of . We implemented this idea using clonal strain pairs in P. vulgatus in Fig 5 and found that the vast majority (>90%) of nonsynonymous differences are private to a single strain pair (Methods Section 7). For the adaptive reversion model to be consistent with this observation, the average steady-state probability must be less than , where n is the total size of the cohort. For our cohort size of , this implies that

(38)

This constraint introduces a tension: if adaptive substitutions are rare because is small, then the fraction of adaptive reversible sites () must be large to maintain the necessary nonsynonymous substitution rate. Combining this result with Eq. 36 and the estimated value of from Fig 3, we find that must satisfy the condition

(39)

In other words, at least one percent of all nonsynonymous sites in the genome—corresponding to tens of thousands of sites—would need to participate in reversible adaptive dynamics to generate the observed pattern in Fig 3. If recurrent clonal SNVs continue to remain as rare in cohorts of increasing sizes, the upper bound on would become more stringent, requiring an even larger fraction of nonsynonymous sites to be reversible in order to explain the observed trend.

Similar arguments hold for more general distributions of . Echoing the purifying selection model in Eq. 1, the timescale of the curve requires that the vast majority of adaptive reversions have and . These constraints suggest that the adaptive reversion model can account for the empirical patterns only if a substantial fraction of nonsynonymous sites experience reversible selection on environmental timescales much longer than the human generation time.

Limitations and future directions. Our analysis above has focused on a particularly simple form of the adaptive reversion model (Eq. 35), similar to the one employed in Ref. [30]. Extensions of this model that incorporate more elaborate forms of epistasis, compensatory mutations [30], or mixtures of adaptation and purifying selection [61] could potentially some of the constraints identified above. Efforts to explore this space of models—similar to Ref. [30]—remain an important avenue for future work.

We have also neglected the effects of recombination. While this was motivated by our ability to focus on the clonal regions of the genome where the recombined segments have been removed, recombination can still affect this clonal diversity in an indirect manner by removing portions of the genome (and any variants they might contain) from both the numerator and denominator of dN/dS. If recombination events affect synonymous and nonsynonymous variants in an unbiased manner, their effects will cancel out when computing the clonal dN/dS ratio in Figs 3 and 4. However, recombination could also be an important mechanism for generating the adaptive reversion events considered above [30], which could in principle lead to systematic biases in the realized recombination rates at synonymous versus nonsynonymous sites. In the simple models considered above, these biases would effectively contribute to the overall reversion rate . However, recombination can also be an important source of forward mutations, which would reduce the effective value of in the presence of modular epistasis, and further limit the contribution of reversible adaptive mutations to the clonal dN/dS curves in Figs 3 and 4. Accounting for this behavior with more explicit recombination-aware models is another important avenue for future work.

7 Analysis of clonal versus population-wide SNVs

To examine the distribution of clonal SNVs across the broader host population, we focused on the subset of strain pairs that were predicted to be fully clonal (i.e., those with no inferred recombination events). Because multiple pairs of strains can all be closely related to each other, clonal SNVs can occur in more than one pair by common descent rather than by homoplasy. To avoid such redundancy, we filtered the initial strain set so that the remaining genomes were either clonal relatives or nearly fully recombined with respect to one another, defined as sharing less than 50% of their genome in identical blocks (Methods Section 2). Within each group of strains connected by clonal relationships, we retained a single representative clonal pair. This procedure ensured that clonal variants were assessed only against otherwise unrelated genomic backgrounds.

In the case of P. vulgatus, our metagenomic dataset contains two major clades that are sometimes classified as separate species (P. vulgatus and P. dorei; [72]). We restricted our analysis to the major P. vulgatus clade. After filtering, this procedure yielded a total of 120 strains, including 17 independent pairs of clonal strains. We catalogued all SNV differences at non-degenerate (1D) and 4-fold degenerate (4D) coding sites between clonal pairs and computed the frequency with which each such variant appeared across the full set of 120 strains. For comparison, we also tabulated the prevalence of all SNVs present in at least one of the 120 strains, irrespective of whether they were originally identified in clonal pairs.

We repeated this procedure for each species with at least five clonal pairs, applying the same strain filtering and allele prevalence analysis. The resulting distributions of clonal versus population-wide SNVs, stratified by degeneracy class and allele prevalence, are shown in Fig Q in S1 Text.

8 Analysis of Staphylococcus aureus

To assess whether our results extend beyond human gut commensals, we also analyzed a collection of 110 publicly available Staphylococcus aureus genomes obtained from Ref. [42] and originally sequenced by Ref. [73]. These assemblies had already been aligned to the reference genome of strain MRSA252 in Ref. [42], and we annotated site degeneracy using the same reference–based in Methods Section 1. When estimating between genome pairs, we utilized all sites jointly covered between the two genomes, rather than restricting to a predefined core genome set.

Unlike most gut bacteria, the recent evolutionary history of S. aureus is dominated by clonal expansion, with relatively few recombination events [42]. However, three large hybridization events have been documented previously, corresponding to the formation of sequence types (STs) 34, 239, and 582. These strains arose through the replacement of 200–600 kb regions in ST30, ST8, and ST25 genetic backgrounds, respectively [42–44]. We used these previously identified tracts to stratify sites into “clonal” versus “recombined” categories and then estimated separately within each class using the same Poisson–thinning procedures described above (Methods Section 2).

Pairs of genomes that did not belong to these hybrid sequence types were included in the genome-wide analysis, but were omitted from the stratified clonal versus recombined comparison. Despite the limited number of recombination events, the resulting dynamics again revealed distinct trajectories in between clonal and recombined regions, consistent with the predictions of our recombination model.

9 Analysis of isolate genomes from P. vulgatus and A. putredinis

Isolate genome selection. Cultured-isolate genome assemblies were retrieved from NCBI (RefSeq/GenBank) via the NCBI Datasets v2 API (accessed June 1, 2026) for two gut species with sufficient publicly available isolates: Phocaeicola vulgatus (MIDAS reference Bacteroides_vulgatus_57955) and Alistipes putredinis (Alistipes_putredinis_61533). For each species, paired GenBank/RefSeq records were collapsed to unique physical assemblies (retaining RefSeq), and metagenome-assembled genomes (MAGs) were excluded both by the NCBI assembly_type field (“derived from metagenome”) and by a keyword screen of the BioSample, BioProject, and organism metadata (e.g., “metagenome”, “metaSPAdes”, “MEGAHIT”, “bin”). The remaining cultured isolates were then restricted to high-quality assemblies (CheckM completeness 95%, contamination 5%, or a closed genome). For P. vulgatus this funnel proceeded from 3,041 assemblies to 2,265 unique assemblies, 1,139 cultured isolates (1,126 MAGs removed), and 369 high-quality isolates; for A. putredinis, from 1,816 assemblies to 1,048 unique assemblies, 69 cultured isolates (979 MAGs removed), and 41 high-quality isolates. Species identity was confirmed for every retained assembly by NCBI average-nucleotide-identity assignment to the corresponding type strain (taxonomy check “OK”, species match). For P. vulgatus, where the assembler was reported it was in every case a single-isolate assembler (e.g. SPAdes, Unicycler, Flye), with no metagenomic assembler detected. The resulting 369 P. vulgatus and 41 A. putredinis assemblies were used for downstream SNV catalog construction and dN/dS analysis.

SNV catalog construction and downstream analysis. For each isolate assembly, whole-genome alignments to the corresponding MIDAS reference genome [62] were computed using NUCmer (MUMmer 4) [74], and per-site base calls and coverage masks were derived from the resulting delta files after filtering with delta-filter and parsing with show-snps. Sites were annotated as synonymous (4D) or nonsynonymous (1D) using the reference gene features, and SNV tables were stored in the same sample site format used for the metagenomic analysis. The resulting SNV tables for both isolate and metagenomic datasets therefore share the same coordinate system and site-type annotations, and all downstream analyses—recombination detection with the CP-HMM, stratification of mutation counts into clonal and recombined regions, and dN/dS estimation—were applied identically to both.

10 Correcting for missing recombination calls at nonsynonymous-rich tract boundaries

Our HMM procedure [34] detects recombination events from the local density of synonymous (4D) differences. In some species (A. putredinis), we found that the detected transfers often carried flanking regions with anomalously large numbers of 1D variants (Fig F in S1 Text); these flanking regions would be misclassified as clonal sequences by our synonymous-only HMM approach, since they lack a corresponding enrichment of 4D variants. To correct for this issue, we implemented a heuristic post-processing step in our HMM package that extends the boundaries of each detected transfer to cover adjacent nonsynonymous-rich sequence, without altering the underlying 4D HMM or its priors. Starting from each originally detected boundary, we walked outward over successive sites carrying a 1D difference between the pair, accepting a step to the next such site if (i) its distance from the current frontier did not exceed a gap threshold (1,000 bp; a longer 1D-free run was taken to mark background clonal sequence and terminated the extension), and (ii) the running density of 1D differences in the extended flank exceeded max(3perkb,5λ0), where is the pair’s genome-wide clonal 1D-difference density. Extensions were capped at 5 kb per boundary, and tracts that came to abut after extension were merged. To remain conservative, we retained only extensions that absorbed at least two 1D differences, as a diagnostic analysis showed that single-difference extensions, though individually tight, produced a large aggregate change driven by many marginal boundary adjustments. Because this procedure extends existing tracts, it cannot recover transfers that were entirely undetected; pairs with no detected recombination are unchanged.

Masked clonal dN/dS and refitting. For each strain pair, the nonsynonymous and synonymous differences absorbed into the extended tracts were reassigned from the clonal to the recombined region of the genome, while the clonal opportunity counts (denominators) were held at their original values, since the corresponding change in flank length is negligible. We then recomputed clonal per species and re-pooled the counts across the cohort (Figs G and L in S1 Text). The aggregate purifying-selection fit was repeated exactly as in the main analysis above.

This correction reduced the elevated clonal dN/dS of A. putredinis from 2 to a nearly-neutral value of 1, while leaving the pooled selection-strength estimate essentially unchanged ( varied by only about 2-fold across masking scenarios, preserving its order of magnitude; Figs G and L in S1 Text). The magnitude of this correction was markedly larger in Alistipes than in other genera like Bacteroides, though appreciable adjustments also occurred in a few other species, such as Akkermansia muciniphila and Dialister invisus). Why this signal is so pronounced in Alistipes is currently unclear. Further investigations of these nonsynonymous-rich transfers would be an interesting direction for future work, which could potentially shed light on the functional consequences of the pervasive recombination observed in these bacterial populations.