Bacterial phylogenetic reconstruction from whole genomes is robust to recombination but demographic inference is not.

UNLABELLED: Phylogenetic inference in bacterial genomics is fundamental to understanding problems such as population history, antimicrobial resistance, and transmission dynamics. The field has been plagued by an apparent state of contradiction since the distorting effects of recombination on phylogeny were discovered more than a decade ago. Researchers persist with detailed phylogenetic analyses while simultaneously acknowledging that recombination seriously misleads inference of population dynamics and selection. Here we resolve this paradox by showing that phylogenetic tree topologies based on whole genomes robustly reconstruct the clonal frame topology but that branch lengths are badly skewed. Surprisingly, removing recombining sites can exacerbate branch length distortion caused by recombination. IMPORTANCE: Phylogenetic tree reconstruction is a popular approach for understanding the relatedness of bacteria in a population from differences in their genome sequences. However, bacteria frequently exchange regions of their genomes by a process called homologous recombination, which violates a fundamental assumption of phylogenetic methods. Since many researchers continue to use phylogenetics for recombining bacteria, it is important to understand how recombination affects the conclusions drawn from these analyses. We find that whole-genome sequences afford great accuracy in reconstructing evolutionary relationships despite concerns surrounding the presence of recombination, but the branch lengths of the phylogenetic tree are indeed badly distorted. Surprisingly, methods to reduce the impact of recombination on branch lengths can exacerbate the problem.

Observation

Phylogenetic methods are powerful and widely used tools for reconstructing the ancestral history of pathogen populations. These methods have been used extensively in evolutionary contexts and are increasingly applied to bacterial populations in clinical settings for strain classification and outbreak detection (1). Such applications require accurate estimation of the phylogenetic tree, but this can be problematic for bacteria due to recombination, in which DNA is exchanged via transformation, transduction, or conjugation (2). In the early 2000s, several authors demonstrated that recombination distorts phylogenetic inference, leading to biased estimates of branch lengths, artifactual signals of population expansion (3), false inference of positive selection (4, 5), and unreliable reconstruction of the tree topology (6, 7). Recombination causes tree topology and branch lengths to change along the genome, preventing a single tree from adequately explaining the reticulated ancestry of recombining sequences.

With the advent of accessible whole-genome sequencing, phylogenetic approaches are increasingly being used to reconstruct the evolutionary history of bacterial populations from their genome sequences (1, 8). The prevalence of phylogenetic analyses despite their demonstrable problems raises difficult questions concerning the credibility of conclusions drawn from phylogenetic inference. The esthetic appeal of phylogenetic trees partly explains their continued popularity, but the lack of viable alternatives is also an important factor. Several sophisticated methods attempt to model reticulated ancestries, but their practical application has been limited by computational demands (9–14). However, we contend that phylogenetic approaches have endured because biologists have found they convey meaningful information about the structure and relatedness of bacterial populations that fits with other evidence.

Milkman and Bridges (15) introduced the concept of the clonal frame to describe the phylogeny of sites in the bacterial genome that have not experienced recombination. Since a bacterial recombination event typically affects only a fraction of the genome, continual assault by recombination throughout the genome would be required to obliterate the signal of the clonal frame. Despite the attention given to the effect of recombination on phylogenetic inference, investigation into the accuracy of topological reconstruction has been limited to analyses of single or concatenated gene sequences and small sample sizes (6, 16). Therefore, we reasoned that phylogenetic inference might be reliably recovering the signal of the clonal frame from bacterial genomes, which could explain the continued faith placed in phylogenetic inference despite the problem of recombination.

We set out to test this idea through simulation. We simulated 1,000 populations of 100 bacterial genomes, each 1 Mb long with moderate mutation (substitution rate [θ] = 1%) under three scenarios: high, low, and no recombination (recombination rate [ρ] = 1%, 0.1%, and 0%, respectively). For each simulation, we recorded the clonal frame and estimated the phylogeny using neighbor joining (NJ) (17), unweighted-pair group method with arithmetic means (UPGMA) (18), maximum likelihood (ML) (19), and BEAST (20) (full details in Text S1 in the supplemental material). We quantified accuracy as the percentage of branches in the clonal frame correctly reconstructed. We found that the clonal frame topology was reconstructed remarkably accurately even when recombination was present (>97% [Fig. 1b]). Increasing ρ only modestly reduced accuracy, which appeared to be driven by the shorter branches (see Fig. S1 in the supplemental material). In a model of stable population size, branches nearer the tips tend to be shorter, whereas in an exponentially growing population, the tendency for tips to be shorter than deep branches is reduced, and at high growth rates, it is reversed (21). As such, branches closer to the root are less accurate at high recombination rates for exponentially growing populations (Fig. S2). In contrast, we found that bootstrap values (NJ, UPGMA, and ML) and posterior probabilities (BEAST) were upwardly biased by recombination (Fig. S3). Our results indicate that the accuracy of the tree topology decays progressively with increasing recombination rate. It follows that at very high recombination rates, it would no longer be sensible to pursue tree-based inference, although even at ρ = 8%, we found that topological accuracy remained high (93% based on 100 simulations with constant population size).

FIG 1 

Effects of recombination in bacteria on phylogenetic tree topology and growth rate estimates. (a) The true clonal frame (left) and ML phylogenies constructed from all sites (center) and only nonhomoplastic sites (right) representing the evolutionary history of a population of 100 bacterial genomes of 1 million base pairs. The recombination rate (ρ) and substitution rate (θ) were fixed at 1%. The number of homoplasies per branch is shown for the center tree. (b) Estimates of branch accuracy for trees reconstructed using ML, BEAST, NJ, and UPGMA at three different values of ρ. The means and standard errors are based on 1,000 simulations of a demographic model of constant population size. (c) Mean posterior estimates of the exponential growth rate parameter (g) from BEAST, averaged over analyses of 1,000 simulated data sets. Data were simulated under a demographic model of constant population size (gray), low exponential growth (blue), and high exponential growth (red) and at three different values of ρ. Error bars represent the mean 95% confidence intervals. Estimates from analyses using either all sites in the sequence alignment (filled triangles) or only those sites without homoplasies (open circles) are plotted. Black dashed horizontal lines represent the true value of the exponential growth rate parameter used in the simulations.

In contrast to the robustness of the phylogenetic topology, recombination gave rise to a spurious or inflated signal of demographic growth when we fitted a model of exponential growth using BEAST (Fig. 1c). In simulations under high, low, and no growth (exponential growth rate parameter [g] = 10, 1, and 0, respectively), growth rates were systematically overestimated, even though tree topology remained accurate (>98% for ρ = 0.1% and 1%; see Fig. S4 in the supplemental material).

Some authors have recommended the removal of recombining sites to ameliorate their detrimental effect on phylogenetic analysis, in particular the tendency for recombination to produce a spurious signal of exponential growth (22–24). Recombination generates various signatures including homoplasy, in which the same substitution is observed in different parts of the tree. Homoplasy can be generated by repeat and back mutation, but it also results from reshuffling diversity among ancestral lineages by recombination, so that excess homoplasy is indicative of levels of recombination sufficient to cause problems for phylogenetic inference (25). We investigated whether removing homoplastic sites improved the estimation of exponential growth rates by BEAST. We found that removing homoplasies actually exacerbated the spurious signal of demographic growth generated by recombination (Fig. 1c), because older recombination events were more likely to be detected as homoplasies. This led to preferential removal of substitutions from the deep branches of the tree, producing trees that appeared even more star-like (Fig. 1a). The magnitude of the effect increased with higher recombination rates, producing 95% confidence intervals that excluded the true growth rate. The number of homoplastic sites removed due to repeat and back mutation amounted to 0.2% of the genome and had a negligible effect on the estimation of growth rates (observed in the absence of recombination in Fig. 1c). We found that removal of homoplasies followed by reestimation of the phylogeny had limited effect on the accuracy of the topology itself (see Fig. S5 in the supplemental material).

In summary, our results show that the clonal frame topology is robustly reconstructed from bacterial whole genomes by phylogenetic methods even in the presence of recombination, but the branch lengths of the clonal frame are not. Removal of recombining sites exacerbates branch length distortion, because older events are easier to detect than young ones, meaning that phylogenetic-based demographic inference should still be viewed with caution in recombining species.

Supplemental methods. Additional details of the simulations, phylogenetic tree construction, calculation of tree accuracy, and identification of recombining sites are given. Download

Text S1, PDF file, 0.1 MB

Branch accuracy for trees reconstructed using ML, BEAST, NJ, and UPGMA at three different values of the recombination rate (ρ) and growth rate (g). Branches are partitioned into three intervals according to their length, selected in an attempt to keep the number of branches within intervals the same (mean of 65.3 branches). The mean number of branches per interval for each method is displayed above each bar. Means and standard errors are based on analyses of 1,000 simulations under  a demographic model of constant population size (g = 0) (gray), low exponential growth (g = 1) (blue), and high exponential growth (g = 10) (red). Download

Figure S1, PDF file, 0.05 MB

Branch accuracy for trees reconstructed using ML, BEAST, NJ, and UPGMA at three different values of the recombination rate (ρ) and growth rate (g). Branches are partitioned into three intervals according to the distance between the end of the branch and the root node. These intervals differed between growth rates in an attempt to keep the number of branches within intervals the same (mean of 65.3 branches). The mean number of branches per interval for each method is displayed above each bar. Means and standard errors are based on analyses of 1,000 simulations under a demographic model of constant population size (g = 0) (gray), low exponential growth (g = 1) (blue), and high exponential growth (g = 10) (red). Download

Figure S2, PDF file, 0.04 MB

Accuracy of branches in estimated trees partitioned by the accuracy of either the bootstrap value or posterior probability support for each branch. Trees were reconstructed using ML, BEAST, NJ, and UPGMA at three different values of the recombination rate (ρ). The expected linear relationship between support and accuracy is plotted in blue. Means and standard errors are based on analyses of 1,000 simulations under a demographic model of constant population size. Download

Figure S3, PDF file, 0.1 MB

Branch accuracy for trees reconstructed using ML, BEAST, NJ, and UPGMA at three different values of the recombination rate (ρ) and growth rate. Means and standard errors are based on analyses of 1,000 simulations under a demographic model of constant population size (g = 0) (gray), low exponential growth (g = 1) (blue), and high exponential growth (g = 10) (red). Download

Figure S4, PDF file, 0.04 MB

Branch accuracy for trees reconstructed using ML, BEAST, NJ, and UPGMA from genome sequence alignments after the removal of homoplasies. Data were simulated under three different values of the recombination rate (ρ) and growth rate (g). Means and standard errors are based on analyses of 1,000 simulations under a demographic model of constant population size (g = 0) (gray), low exponential growth (g = 1) (blue), and high exponential growth (g = 10) (red). The accuracy of NJ trees is marginally greater than for those constructed from all sites (see Fig. S3 in the supplemental material), due to the alignment being enriched for sites supporting the clonal frame. However, the accuracy of UPGMA trees is lower after the removal of homoplasies when ρ = 1%. Download

Figure S5, PDF file, 0.03 MB