Main

Phylogenetics is central to evolutionary biology1. Phylogenetic trees are graphs representing evolutionary histories and ancestry, and are typically inferred from the genomic data. In genomic epidemiology, for example, DNA sequences of the same pathogen are collected from different hosts (for example, SARS-CoV-2 genomes collected from different patients) and compared with one another. Phylogenetic trees inferred from these genomes can reveal the emergence of drug resistance and new variants of concern, transmission between individuals and countries, and many other details of the evolution and spread of the pathogen4.

Most phylogenetic methods that are scalable to large datasets—such as maximum-likelihood, parsimony-based and heuristic approaches—typically estimate a single phylogenetic tree without intrinsically assessing the reliability or uncertainty of inferences. This issue is typically addressed with additional methods such as Felsenstein’s bootstrap2. For a given dataset, Felsenstein’s bootstrap typically creates 100 or 1,000 replicates by randomly resampling the genetic data with replacement. Phylogenetic inference is performed on each to estimate replicate trees, and the support score of a clade (the set of taxa inferred to be all the descendants of one ancestor in the tree) is defined as the proportion of replicate trees containing that clade. This is also considered to be the support of the phylogenetic branch that separates the clade from all other taxa in the tree.

Felsenstein’s bootstrap has been developed in the context of inter-species evolutionary biology. Consequently, it has a number of drawbacks when applied to genomic epidemiological datasets of closely related genetic sequences. First, repeatedly performing phylogenetic estimation on all replicate datasets can be excessively computationally demanding. Although Felsenstein’s bootstrap approximations5,6,7 are more efficient, these are still not feasibly applicable to pandemic-scale phylogenetic analyses involving millions of genomes8,9.

Second, even a small number of ‘rogue taxa’—that is, sequences whose placement in the inferred phylogenetic tree is highly uncertain (such as incomplete sequences or recombinants)—can substantially lower the Felsenstein’s bootstrap support of many internal branches throughout phylogenetic trees10. Third, Felsenstein’s bootstrap does not measure posterior probability, but rather measures repeatability11. In genomic epidemiology, a single mutation is typically sufficient to define a clade with negligible uncertainty. However, in this scenario, Felsenstein’s bootstrap usually requires three mutations supporting any one clade to be able to assign 95% support to it, making it excessively conservative2,12.

Finally, Felsenstein’s bootstrap and most other phylogenetic branch support measures have a ‘topological’ focus, in that they aim at assessing the reliability of the inferred tree topology via its constituent clades. Although clades are important in taxonomy, they are not as relevant in genomic epidemiology, where the focus is typically on mutation and transmission histories and lineage assignment3,13.

Existing local branch support measures14,15,16,17 are considerably more computationally efficient than Felsenstein’s bootstrap, but also rely on a topological interpretation. They usually compare the likelihood of the inferred phylogenetic tree against the likelihoods of similar alternative trees (Methods, ‘SPRTA and aBayes’). Local branch support measures are particularly appealing because of their computational efficiency and their robustness to rogue taxa.

Here, drawing on concepts from local branch support measures and the approximate Bayes approach (aBayes)15, we present a new measure of branch support that is robust to rogue taxa and can scale to pandemic-scale trees. Subtree pruning and regrafting-based tree assessment (SPRTA) shifts the focus from a topological point of view to one centred around the evolutionary origin of lineages (‘mutational’ or ‘placement’ focus).

Principles of SPRTA

As typical in molecular phylogenetics, we assume that genetic data are represented as a multiple sequence alignment D: a matrix in which each row corresponds to the genetic sequence of a considered taxon (in genomic epidemiology this is often the genome sequence of one of the considered samples, aligned to the reference genome) and each column contains all homologous nucleotides (those descending from the same ancestral nucleotide). We also assume we have a rooted phylogenetic tree T (Methods, ‘SPRTA and aBayes’) inferred from D. Our aim is to assign confidence scores to branches b of T (see Fig. 1). Most branch support methods assign scores that represent the confidence that the sequences contained in subtree S**b (black triangle in Fig. 1), containing all descendants of b, indeed form a clade within T—we refer to this as the topological focus. Instead, we want to assess the probability that b correctly represents the evolutionary origin of S**b—we call this a mutational focus. To clarify this, we describe branch b as having immediate ancestor A and descendant B (the root of subtree S**b), and dividing T into S**b and its complement T*Sb* (Fig. 1). A and B might correspond to individual genomes if, as typical in genomic epidemiology, uncertainty regarding ancestral genomes is low; otherwise, they might represent sets of probable genomes. Our aim is to efficiently approximate the probability that B evolved from A through mutations along branch b (that is, B and *Sb* are ‘placed’ at A), as opposed to the alternative hypothesis that B descended from some other node of T*Sb* (‘alternative placements’ of *Sb* (blue arrows in Fig. 1)). In other words, SPRTA(b) represents the confidence that branch b is the evolutionary origin of B and S**b.

Fig. 1: Graphical representation of SPRTA branch support measurement.

Here we represent SPRTA branch support calculation for an example tree branch b (highlighted as a solid red line). The subtree S**b below b is represented with a black triangle and is not affected by any of the SPR moves considered to evaluate the support of b. A and B represent the ancestors of b and S**b, respectively. Solid black lines represent the branches of T*Sb*, which are assessed as alternative placements for *Sb*, and the dashed black line represents the rest of T (which is an arbitrarily large tree), not shown here. Some possible SPR moves are highlighted with dotted blue arrows, and cause hypothetical new branches (shown as shaded red dashed lines), leading to copies of S**b (grey triangles) descending from nodes A2–A4. The relative likelihood of the original tree is (\Pr (D| {T}_{1}^{{b})) while the likelihoods of alternative topologies are represented by (\Pr (D| {T}_{2}}{b})), (\Pr (D| {T}_{3}^{{b})) and (\Pr (D| {T}_{4}}{b})). In this case, the support for branch b will be (\Pr (D| {T}_{1}^{{b})/{\sum }_{i=1}}{4}\Pr (D| {T}_{i}^{b})). In practice, a large number I**b of alternative topologies might be considered for each branch b (Methods, ‘SPRTA and aBayes’).

To achieve this, for each branch b, SPRTA considers a number I**b of alternative topologies ({T}_{i}^{{b}) (1 ⩽ i ⩽ I**b) of T obtained by performing single subtree pruning and regrafting (SPR) moves—changes to T that relocate S**b as a descendant of other parts of T*Sb* (Fig. 1; *Ib* is defined in Methods, ‘SPRTA and aBayes’). Each such SPR move represents a possible different origin of B, as a direct descendant of a node A**i other than A. For notational convenience, we assume that ({T}_{1}}{b}=T) is the original tree topology. The likelihood (\Pr (D| {T}_{i}^{b})) of such topologies is routinely and efficiently calculated by MAPLE9 (up to a constant common multiplicative factor) and we use it to define the SPRTA support score:

$${\rm{SPRTA}}(b)=\mbox{Pr}(b| D,T\backslash b)=\frac{\mbox{Pr}(D| T)}{{\sum }_{1\leqslant i\leqslant {I}_{b}}\mbox{Pr}(D| {T}_{i}^{b})}$$

(1)

where T*b* represents the union of subtree S**b and its complement subtree T*S**b*. SPRTA is described in more detail in Methods ‘SPRTA and aBayes’.

SPRTA mutational branch support scores need to be interpreted differently from those of existing topological branch support methods. SPR moves involving branch b preserve the clade defined by b (Fig. 1), and so SPRTA scores do not represent an assessment of this clade. Rather, SPRTA scores are to be interpreted as an approximate probability that B evolved directly from A. These scores are a particularly useful assessment of tree reliability in the context of genomic epidemiology, where the placement of individual incomplete sequences (rogue taxa) can be highly uncertain. SPRTA scores are expected to be robust to rogue taxa, since their placement is expected to have negligible effect on relative likelihood scores at internal tree nodes.

SPRTA scores for the terminal branches of a tree evaluate the placement probability of individual observed sequences, and in fact correspond closely to the probabilistic support measure used by some tools that map query sequences onto a pre-estimated phylogeny18,19. By contrast, topological support methods cannot assess terminal branches and sequence placements.

Computational demand

The SPR search required by SPRTA is typically performed as part of the phylogenetic tree search in many maximum-likelihood phylogenetic methods such as RaxML20 and MAPLE9, and so it is not expected to lead to significant additional runtime when executed in conjunction with phylogenetic inference. For comparison, we assessed the computational demand of SPRTA against the main existing measures of branch support (Felsenstein’s bootstrap2, local bootstrap probability (LBP)14, approximate likelihood ratio test (aLRT)15, aLRT with Shimodaira–Hasegawa test (aLRT-SH)16, Bayesian-like transformation of aLRT (aBayes)17, transfer bootstrap expectation (TBE)10 and ultrafast bootstrap approximation (UFBoot)7). SPRTA reduces runtime and memory demands by at least two orders of magnitude compared with all these methods, with the difference growing as dataset size increases (Fig. 2 and Extended Data Fig. 1a,b). Note also the premature termination of the lines in Fig. 2 and Extended Data Fig. 1 for the other methods, indicating cases where they could not be run successfully. Among these other methods, those based on Felsenstein’s bootstrap (Felsenstein’s bootstrap, UFBoot and TBE) have substantially higher computational demand than local branch support measures. This comparison shows that SPRTA can assess much larger phylogenetic trees than existing approaches.

Fig. 2: Computational demand.

a,b, Time (a) and maximum memory usage (b) of different branch support methods. The x axis shows the number of simulated SARS-CoV-2 genomes (Methods, ‘Simulated genomes’) included in each replicate. For each dataset size considered, we ran 20 replicates. Dots show the mean across replicates and violin plots show variation between replicates. All analyses were run on one core of an Intel Xeon Gold 6252 processor at 2.10 GHz. Each method was run on the same replicates for the results in a,b. In b, the lines for all methods other than UFBoot and SPRTA overlap. FB, Felsenstein’s bootstrap.

Accuracy

To benchmark different branch support methods, we simulated SARS-CoV-2-like genome data (Methods, ‘Simulated genomes’) for which we know the true tree and the true mutational history along the tree. Because genomes in genomic epidemiological datasets are typically very similar to each other, we usually infer ancestral genomes and mutation events conditional on an inferred phylogenetic tree T with negligible uncertainty. In this case the task of assessing whether B evolved from A along branch b is equivalent to assessing the correctness of the mutation events implied by T on b, since alternative placements of S**b often result in alternative mutational histories leading to B. Therefore, in this benchmark we consider a mutational interpretation of branch support, and interpret branch support scores as an estimate of the posterior probability of the mutation events implied by T on the considered branch b. Each individual branch of a tree inferred from simulated data is considered correct if the mutations inferred on that branch actually happened in our simulations, and so if B at the lower end of b actually evolved from A. An advantage of this approach is that, whereas S**b or the clade defined by it are often uncertain and might not exist in the real tree, B will typically not be affected by rogue taxa or other topological uncertainty within S**b. More details on how we define accuracy of branch support are provided in Methods, ‘Benchmarking of SPRTA support’. This is not the typical interpretation of branch support measures other than SPRTA, and so results here are not meant as an evaluation of other branch support measures in their classical topological interpretation. Below, all branch support measures are evaluated based on their ability to assign higher support to correctly inferred mutation events, and lower support to erroneously inferred mutation events.

Generally, all methods give high support scores (>80%) to all mutations, both correctly and wrongly inferred ones (Fig. 3 and Extended Data Fig. 1c,d). This is to be expected given the low level of divergence in the considered dataset, and that these mutations have been inferred by maximum likelihood and so by definition will typically have higher likelihood than alternative mutation histories. However, SPRTA is the only method that reliably assigns higher support to correctly inferred mutation events (typically a mean of 98–99%) and lower support scores to wrongly inferred ones (typically a mean of 85–90%) (Fig. 3 and Extended Data Fig. 1e,f). Although some methods like UFBoot, aLRT, and aBayes give higher support to correctly inferred mutations than does SPRTA, they also similarly support wrongly inferred mutations, and so do not discriminate between the two. This is also reflected in higher AUROC and AUPRC values for SPRTA (Extended Data Figs. 2–5). It should be remembered, however, that topological approaches have been developed for a different interpretation of branch support scores, and so these results are not indicative of their performance for that task. What our analysis does show is that SPRTA can usefully and effectively assess the confidence in mutational histories implied by inferred phylogenetic trees.

Fig. 3: Benchmark of branch support methods.

a, Mean support of correctly inferred mutations (supp. correct mut.) within each replicate. b, Mean support of wrongly inferred mutations (supp. wrong mut.) within each replicate. Note that values shown are distributions of mean support values across all mutations in a replicate simulation, and not distributions of support scores of individual branches or mutations. Other details are as in Fig. 2. Branch support scores for methods other than SPRTA were again only calculated when computationally feasible.

Uncertainty of SARS-CoV-2 evolution

We applied SPRTA to a global public dataset of 2,072,111 SARS-CoV-2 genomes (Methods, ‘Viridian genome dataset’) that is too large for the feasible application of existing topological branch support measures. Phylogenetic tree estimation with MAPLE, parallelized over 14 cores of an Intel Xeon Gold 6252 processor at 2.10 GHz, took around 10 days (ref. [21](https://www.nature.com/articles/s41586-025-09567-x#ref-CR21 “De Maio, N. et al. Rate variation and recurrent sequence errors in pandemic-scale phylogenetics. Preprint at bioRxiv https://doi.org/10.1101/2024.07.12.603240

(2024).“)). Post hoc SPRTA mutational support assessment, re-performing SPR move evaluations that had already been performed during tree inference, required 7 h 27 m on a single core and maximum memory 26.93 Gb. This branch support assessment can be also performed simultaneously with tree inference, at negligible additional computational cost compared to the tree inference itself. During tree inference, however, we do not evaluate possible alternative placements of less informative genomes (a sequence is not informative at a position if it contains the ‘N’ character at that position, and a sequence s1 is less informative than s2 if it coincides with s2 except at positions where s1 is not informative), since alternative placements of these genomes do not increase the likelihood of the considered tree[21](https://www.nature.com/articles/s41586-025-09567-x#ref-CR21 “De Maio, N. et al. Rate variation and recurrent sequence errors in pandemic-scale phylogenetics. Preprint at bioRxiv https://doi.org/10.1101/2024.07.12.603240

(2024).“). To evaluate alternative placements of all genomes—for example, to identify a larger number of rogue taxa in our dataset—we also performed a more in-depth and computationally demanding SPRTA assessment that did not collapse genomes less informative than other genomes in the dataset; this required 22 h 42 m on a single core and maximum memory 27.75 Gb.

Phylogenetic uncertainty

Among the 2,072,111 genomes considered here, 636,022 are mutationally informative (the branch separating them from the tree has inferred mutations on it), but for many of these the mutational history is uncertain: 87,406 have SPRTA placement below 90% and 53,365 have SPRTA placement below 50%. From the remaining 1,436,089 non-mutationally informative genomes (those that are separated from the tree by a terminal branch of length 0 and so, conditional on the rest of the tree, are not informative of mutation events), 162,100 have SPRTA placement support below 90% and 115,358 have SPRTA placement below 50%.

Phylogenetic uncertainty also affects many internal branches of the phylogenetic tree, highlighting substantial uncertainty in the ancestral mutational history, and not only in terminal branches: from 453,976 internal branches with inferred mutations, 59,523 have SPRTA support below 90% and 29,641 have SPRTA support below 50%. In terms of inferred mutations (which are assigned the same SPRTA score as the branch they are inferred to have occurred on), out of a total of 1,827,786 there are 249,092 with SPRTA support below 90% and 126,308 with SPRTA support less than 50%.

Although in our simulations we model site-specific genome evolution and incomplete sequence data to closely match patterns observed in real SARS-CoV-2 data (Methods, ’Simulated genomes’), we still find lower SPRTA scores in real data than in simulations (Extended Data Fig. 1g). This suggests further complexities in real data such as site- and nucleotide-specific mutational patterns22 that are not accounted for in our models of genome evolution.

A prominent examples of uncertainty in the SARS-CoV-2 mutational history is the evolution of the AY.4 Delta sub-lineage, one of the most abundant SARS-CoV-2 lineages, represented here by more than 480,000 genomes. Two mutations appear to be ancestral to most AY.4 genomes: T17040C and C21846T (Fig. 4, middle). After the appearance of mutation T17040C, the reversion C17040T appears to have occurred many times: we infer around 650 reversions. Position 17040 is inferred by MAPLE to have a substitution rate 31.9 times higher than average[21](https://www.nature.com/articles/s41586-025-09567-x#ref-CR21 “De Maio, N. et al. Rate variation and recurrent sequence errors in pandemic-scale phylogenetics. Preprint at bioRxiv https://doi.org/10.1101/2024.07.12.603240

(2024).“), mostly as a result of these reversions within AY.4. Although inferred reversions can be due to reference biases in consensus genome calling methods, our consensus genomes were called with Viridian23, which is not affected by this issue[21](https://www.nature.com/articles/s41586-025-09567-x#ref-CR21 “De Maio, N. et al. Rate variation and recurrent sequence errors in pandemic-scale phylogenetics. Preprint at bioRxiv https://doi.org/10.1101/2024.07.12.603240

(2024).“),23. Furthermore, we did not observe any issues with read data or substitution distribution along the tree that would suggest the presence of recurrent sequence errors at position 17040[21](https://www.nature.com/articles/s41586-025-09567-x#ref-CR21 “De Maio, N. et al. Rate variation and recurrent sequence errors in pandemic-scale phylogenetics. Preprint at bioRxiv https://doi.org/10.1101/2024.07.12.603240

(2024).“). This suggests that 17040 might be a genuinely hyper-mutable genome position in SARS-CoV-2, but only when in the mutated C nucleotide state, and not in the ancestral T. This is in line with the observation that mutation patterns in SARS-CoV-2 are highly position- and nucleotide-specific22.

Fig. 4: Uncertain evolutionary history of SARS-CoV-2 lineage AY.4.

Left, the global SARS-CoV-2 phylogenetic tree inferred by MAPLE and visualized in Taxonium25. Tips are coloured according to the Pango3 lineage assigned by Pangolin35 v.4.3 (with Pangolin-data v.1.21) to the corresponding genomes. We also show the names of some of these lineages. Middle, magnified view of lineage AY.4, showing the locations of inferred T17040C (green arrow) and C21846T (blue arrow) mutations that define the inferred early evolution of the lineage. We highlight with red circles the approximately 650 C17040T reversions within AY.4. One of these, highlighted with a red arrow, is ancestral to more than 163,000 genomes and has only 7.6% SPRTA support. Right, further magnification of the location of an alternative origin of this last sub-lineage (with corresponding SPR move represented by the dark blue arrow) that has 5.0% SPRTA support and entails re-placing the aforementioned sub-lineage as a direct descendant of the phylogenetic node with genome containing the T27484C and C17040T mutations.

The largest inferred sub-lineage within AY.4 descending from a C17040T reversion contains more than 163,000 genomes (Fig. 4, middle). However, the branch defining this lineage (Fig. 4, middle, red arrow) has only 7.6% SPRTA support. The reason is that, owing to the hyper-mutability of C17040, there are many alternative plausible mutational histories within AY.4. The most likely alternative origin of this sub-lineage has SPRTA support of 5.0% and involves replacing the C17040T reversion defining this lineage with a C27484T reversion within a genomic background that already contains the C17040T reversion (Fig. 4, right, dark blue arrow). Although many equally parsimonious alternative mutation histories exist, this one is inferred by SPRTA to be the most likely alternative origin of this sub-lineage because the background mutation rate from C to T is very high in SARS-CoV-2, and because position 27484 also has an inferred substitution rate 20.8 times above average[21](https://www.nature.com/articles/s41586-025-09567-x#ref-CR21 “De Maio, N. et al. Rate variation and recurrent sequence errors in pandemic-scale phylogenetics. Preprint at bioRxiv https://doi.org/10.1101/2024.07.12.603240

(2024).“). These results show that the evolutionary history of lineage AY.4 is highly uncertain, and illustrate how SPRTA scores can not only highlight uncertain parts of an inferred phylogenetic tree and mutational history, but also identify and probabilistically assess alternative evolutionary origins of considered pathogen variants.

Regarding individual samples with uncertain placement, we show two examples (of many) in Fig. 5. In the top example, uncertainty is caused by the sample having an incomplete genome sequence; in the bottom example, uncertainty is caused by the existence of two mutually plausible mutational histories. This shows how SPRTA can effectively identify the uncertainty related to the placement of rogue taxa, which by definition will have many possible placements but all with similarly low SPRTA support. Other branch support methods, by contrast, cannot provide an evaluation of sample placement, as they only measure support for internal branches of the phylogenetic tree. Also, unlike for Felsenstein’s bootstrap and most other topological branch support measures, rogue taxa mostly do not affect the SPRTA support scores of ancestral nodes, as their placement typically only affects the inference of mutation events on phylogenetic branches near the placement itself.

Fig. 5: Sample placement uncertainty.

Left, the global SARS-CoV-2 tree as in Fig. 4. Right, magnified view of two uncertain sample placements. Top right, within lineage Alpha, the considered sample (marked by the black arrow) has no sequence information at genome position 27281 (ambiguity character N), the same position where a G-to-T mutation occurs nearby in the phylogenetic tree (marked by the blue arrow). This sample has 50% SPRTA support both at the current placement and at the other side of G27281T (red arrow). Bottom right, within the AY.44 lineage, the placement of a sample has 80.5% SPRTA support, since the mutation C10748T implied by the current placement of the sample (top green arrow) also occurs on a nearby branch (bottom green arrow), permitting an alternative placement with 19.5% SPRTA support (red arrow). This alternative placement would require one fewer C10748T mutation on the tree but one additional G5743T mutation.

These are only some examples of all mutations and placements in the SARS-CoV-2 tree that are substantially uncertain (126,308 out of all 1,827,786 inferred mutations have SPRTA support below 50%; see above). Our full maximum-likelihood tree annotated throughout with SPRTA support scores concisely represents the evolution of SARS-CoV-2 during the COVID-19 pandemic and presents the SPRTA probabilistic assessment of each sampled genome placement and each inferred mutation. It also offers a concise summary of plausible alternative placements and mutational histories. Such pandemic-scale probabilistic phylogenetic assessment—and in general the assessment of mutation events, placements and alternative mutational histories—is not possible with existing popular branch support methods. This annotated tree is available on Zenodo[24](https://www.nature.com/articles/s41586-025-09567-x#ref-CR24 “De Maio, N. Global (2M) SARS-CoV-2 genomes dataset, from Viridian, processed with MAPLE0.6.11. Zenodo https://doi.org/10.5281/zenodo.14974813

(2025).“) together with the considered SARS-CoV-2 genome alignment, and can be visualized easily within Taxonium25.

Effect on inferred mutation rates

In phylogenetic investigations of mutation rates and selective pressures, mutations are usually estimated on a fixed phylogenetic tree, without accounting for phylogenetic uncertainty22. We find that ignoring phylogenetic uncertainty does not substantially affect the inference of genome-wide mean nucleotide substitution rates for SARS-CoV-2 (Extended Data Fig. 6a). However, it can substantially affect the numbers of inferred mutations at individual sites (Extended Data Fig. 6b). One of the sites strongly affected is 17040, owing to previously described uncertain reversions (Fig. 4). Other significantly affected sites also show similar patterns, with many reversions following a substitution (for example, see Extended Data Fig. 7 for position 7926 and Extended Data Fig. 8 for position 21595). The presence of several genome positions with high reversion rate and therefore highly uncertain mutational history might explain why we find lower SPRTA scores in real data than in simulations (Extended Data Fig. 1g). Details of these methods are provided in Methods, ‘Assessing the impact on mutation rates’.

Effect on Pango lineages

The definition of pathogen outbreak lineages (for example, with the Pango3 or Autolin26 systems) and the assignment of samples to these lineages (for example, with the Pangolin tool27) often rely on a fixed phylogenetic tree. Uncertainty in this tree is typically ignored for these tasks. We investigated how much phylogenetic uncertainty might affect the definition of and assignment to Pango lineages, a globally important system for classifying and naming SARS-CoV-2 lineages including variants of concern (detailed methods in ‘Assessing the impact on Pango lineages’).

We found that out of 1,127 Pango lineages in our tree, 26 had more than 5% probability (and 39 had more than 1% probability) of origin from a different Pango lineage than the one implied by the maximum-likelihood tree. Among these, the most noticeable (due to its size and uncertainty) was Pango lineage BA.2.13. This lineage contains more than 500 genomes and is inferred to have descended from Pango lineage BA.2.56 in our maximum-likelihood tree (Fig. 6). However, a likely alternative placement of this lineage suggests that it might have evolved directly from BA.2. The uncertainty appears to be caused by the joint occurrence of the pair of substitutions A23767G and C22792T, both within lineages BA.2 and BA.2.56, giving rise to two almost equally probable evolutionary origins (with approximately 54% and 46% total SPRTA support, respectively) of BA.2.13 (Fig. 6).

Fig. 6: Uncertain origin of Pango lineage BA.2.13.

Left, the global SARS-CoV-2 tree as in Fig. 4. Middle, magnified view of lineage BA.2. Right, magnified view of the two plausible origins of Pango lineage BA.2.13 within BA.2 and its descendants. Top right, the maximum-likelihood placement of BA.2.13 as a descendant of sub-lineage BA.2.56. BA.2.13 has three substitutions that are not part of BA.2: C22916A (red, shared with BA.2.56 and BA.2.9), C22792T (purple, shared with BA.2.9) and A23767G (green). The inferred maximum-likelihood tree suggests that BA.2.13 evolved through a A23767G mutation within a BA.2.56 background already containing substitutions C22916A and C22792T. Bottom right, the distinct occurrence of A23767G and C22792T within the BA.2 lineage. SPRTA shows that lineage BA.2.13 might plausibly have evolved from an additional C22916A substitution within this BA.2 background already containing substitutions A23767G and C22792T, corresponding to a reversed ordering of the three substitutions that define BA.2.13.

We also found that, out of the 2,072,111 genomes in our tree, 3,036 had less than 95% total SPRTA support for being assigned to the Pango lineage onto which they are placed in the maximum-likelihood tree; 754 of these had less than 70% support. We illustrate these findings using sample ERR10476226. This sample has 6 possible placements in the tree, all with the same likelihood and SPRTA support (approximately 16.7%). This is owing to the incomplete consensus sequence of the sample, causing it to be potentially identical to six different genomes in the tree (Extended Data Fig. 9), and so possibly being considered a ‘rogue taxon’. Although the maximum-likelihood placement in the tree of this sample is in Pango lineage BQ.1, all five other placements are on lineage BQ.1.1, which can therefore be interpreted as the most probable lineage assignment for this sample.

Discussion

With the increasing use of genomic epidemiology, pandemic-scale phylogenetics is set to become an essential tool for pandemic preparedness and epidemiology in general. New approaches such as MAPLE9 and UShER8 can be used to infer the huge phylogenetic trees necessary at this scale, but there is considerable uncertainty inherent in these estimates due to complications including low divergence, recurrent mutations and errors, and incomplete genome sequences[21](https://www.nature.com/articles/s41586-025-09567-x#ref-CR21 “De Maio, N. et al. Rate variation and recurrent sequence errors in pandemic-scale phylogenetics. Preprint at bioRxiv https://doi.org/10.1101/2024.07.12.603240

(2024).“),28,29.

Traditional methods to quantify and represent this uncertainty cannot be used with pandemic-scale datasets. To overcome this problem, we present a new approach, SPRTA. In addition to addressing the limitations of existing methods in terms of computational demand, SPRTA also offers a new mutational (or ‘lineage evolution’) interpretation of branch support scores that is particularly useful in genomic epidemiology, replacing the topological focus of previous approaches.

Owing to the excessive computational demand of applying Bayesian phylogenetic methods or Felsenstein’s bootstrap to large genome datasets, pandemic-scale phylogenies are usually taken at face value in downstream analyses such as inference of viral geographic spread30, mutational patterns22, lineage assignment26, recombination31 and variant fitness advantage32. Unchecked errors and uncertainty in phylogenetic trees can therefore propagate in these analyses and affect their accuracy and measures of uncertainty. Our approach allows us to efficiently distinguish between reliable and unreliable parts of an inferred phylogenetic tree, so that downstream analyses can focus on reliable phylogenetic signals or integrate over phylogenetic uncertainty.

Although SPRTA is currently implemented in MAPLE and optimized for applications to datasets with a low level of divergence, it would be possible to implement the same approach within any likelihood-based phylogenetic inference tool. Its mutationally focused scores could therefore be applied to any phylogenetic analysis, without constraints on the number of sequences considered, their divergence or their length. The only limitation to be considered is that at higher divergence, mutational histories conditional on a tree T become more uncertain, and so the mutation-based simulation benchmark that we used here might not be readily interpretable.

SPRTA could also be used to create a phylogenetic network based on a backbone maximum-likelihood phylogenetic tree that is extended to include additional branches with lower but substantial support. In this way, SPRTA could be used to efficiently summarize vast numbers of possible phylogenetic trees (see for example, ref. 33). In future, this might serve as a foundation for developing an efficient and complementary approach to Bayesian phylogenetics, helping to account for tree uncertainty in applications such as phylodynamics34.

In conclusion, SPRTA not only addresses a fundamental outstanding problem in genomic epidemiology, but also offers a new paradigm in evolutionary biology for the interpretation and representation of phylogenetic information and uncertainty.

Methods

SPRTA and aBayes

Given an estimated rooted phylogenetic tree T and data D in the form of a multiple sequence alignment, our aim is to assign confidence scores to branches b of T.

We take inspiration from the approximate Bayes (aBayes) approach17. aBayes assigns to b a probability score Pr(b∣D, T*b*) based on the ratio of the likelihood Pr(D∣T) of the estimated binary tree T versus the likelihoods of the trees ({T}_{i}^{{b}) obtained using nearest neighbour interchange36 (NNI) moves centred around b (which comprise (T={T}_{1}}{b}) itself in addition to two tree topologies not containing the clade in T descending from b):

$${\rm{aBayes}}(b)=\Pr (b| D,T\backslash b)=\frac{\Pr (D| T)}{{\sum }_{i=1}^{{3}\Pr (D| {T}_{i}}{b})}.$$

(2)

These NNI moves perform small changes to T adjacent to branch b. One of the appeals of aBayes is that it can score not only T, but also the two alternative topologies considered, ({T}_{2}^{{b}) and ({T}_{3}}{b}), not containing the clade defined by b:

$$\Pr ({T}_{j}^{b}| D,T\backslash b)=