Antimicrobial resistance (AMR) is a substantial and growing global health burden. Understanding, and predicting, its evolution in specific pathogens will help responses across scales from individual patient cases to large-scale policy. Here, we use global data on AMR features, predicted from 47k Klebsiella pneumoniae genomes, with hypercubic transition path sampling to infer the evolutionary pathways by which AMR features in K. pneumoniae (KpAMR) are acquired across 102 countries, territories, and areas. We identify “globally consistent” evolutionary behaviors that hold across countries, and “globally divergent” behaviors including carbapenem and fluoroquinolone resistance that vary across countries. We show how these divergent dynamics covary both with public health superregion and drug use policy, and reveal competing evolutionary pathways within and between countries. Using newly sequenced data across several decades from sub-Saharan Africa, we show that this inferred global roadmap of KpAMR evolution successfully predicts prospective evolutionary dynamics. Together, we hope that the ability to characterize and predict evolutionary dynamics of AMR acquisition, connected to socio-economic and drug policy predictors, will help strengthen our understanding of AMR evolution worldwide.

Local and global inferred pathways of AMR evolution in Klebsiella

To learn and compare the likely pathways by which KpAMR is acquired in different countries, we obtained a dataset of 47,721 genomes from 102 different countries, territories and areas, and extracted details of the presence or absence of genes corresponding to each of 22 AMR classes, using PathogenWatch and Kleborate [5,26] (Methods; Fig 1A and 1B). We refer to presence/absence of resistance genes for each of these classes as binary “resistance characters” (“character” referring to a particular property of a species in evolutionary biology). The set of characters we consider, following Kleborate, is given in the Fig 1 caption; in this report, we will use the shorthand defined there, grouping genes by drug resistance class (and Lahey class for β-lactamases) [26–28].

For each country in our dataset, we next used HyperTraPS, a machine learning approach, to infer the “pathways” of KpAMR character acquisition: that is, the ordered sequences with which characters are acquired over time (see Methods). Fig 1C illustrates this pipeline for a single country example (Romania, chosen for its moderate number of samples). HyperTraPS takes presence-absence patterns of characters on a phylogeny as input (Fig 1Ci) and outputs an evolutionary model describing the pathways of character acquisition most supported by the data [9,21–23]. The full output of the inference process is a parameterized transition network describing the probability of different sequences of transitions through the state space of possible resistance characters (Fig 1Cii). This model can be summarized, for example, as the set of posterior probabilities P_ij with which character i is acquired at ordering j in a putative evolutionary process (Fig 1Ciii). We refer to a plot where P_ij is visualized with point size, as in Fig 1Ciii, a “bubble plot”.

In Methods, we describe some nuances involved in using HyperTraPS with these data. KpAMR characters may be acquired reversibly (for example, via plasmid gain and loss), while the model underlying HyperTraPS pictures acquisitions as irreversible. This irreversibility picture reduces a potentially vast (or infinite) intractable set of possible evolutionary pathways to a tractable set (Methods; [29,30]). This increases the uncertainty of inference and challenges the inference of character interactions, but evolutionary orderings—which characters are likely present when another is acquired—can generally be robustly inferred (demonstrated in S1A and S1C Fig and more deeply in [31]). The approach does not assume an “end point” of all acquired KpAMR characters exists, that acquisitions occur one at a time, nor that populations are perfectly sampled, homogeneous, or perfectly representative of the labeled country of origin with no mixing [12,32]. The LINcode taxonomy linking observations (see Methods) is used to guard against potential pseudoreplication due to relatedness of samples [18–20,33], but in practise (given the deep-branching nature of many relationships) has a negligible effect on inference outputs [31,34] (S1B Fig). Finally, HyperTraPS does not assume the presence, absence, or form of interactions between features—evolutionary dynamics can be reported for a “null model” of independent characters (with prevalence reflecting individual rates) and/or interacting characters (where rates depend on acquired characters) from the same model framework. Our focus is on characterizing likely ordering dynamics, and not directly on the interactions between characters.

While Fig 1C demonstrates our approach for a single country, Fig 2A gives the average acquisition probabilities across all countries in our dataset (including grouping by health region from the Global Burden of Disease initiative [2]), and S2 Fig gives the corresponding summary plots for every country in the dataset. Resistance to some antibiotic groups is consistently observed across all countries: following very common presence of Bla-chr and SHV-m, resistance characters against betalactams and extended-spectrum betalactams (ESBLs) (Bla-a, Bla_ESBL-a) are often acquired earlier, while resistance to colistin (drug of last resort for multi-drug resistant Kp; Col-a and Col-m) is acquired later. Country-specific differences exist in the evolutionary acquisition of other characters, including resistance to carbapenems (Bla_Carb-a).

In several instances, bimodality in the inferred evolutionary dynamics is observed. That is, a given character may be acquired at an earlier stage or a later stage in the evolutionary process, but not at intermediate stages. Such bimodality is characteristic of mutually repressing interactions between characters, and equivalently of several different evolutionary pathways existing in a system. Examples can readily be seen in Fig 1C for Romania and Fig 2A for the global mean behavior, where SHV-m is likely acquired either very early or at more intermediate stages, with little probability in between. These competing dynamics are directly visible in Fig 1Cii, where distinct, canalized pathways exist: one involving SHV-m as the first acquisition, one involving Bla_chr as the first acquisition and SHV-m only following rather later. Across countries (S2 Fig), bimodality is observed in several characters: some examples are in Bla_chr (Tanzania, Zambia); Bla-a (Ireland, S Korea; Bla_ESBL-a (Laos); and SHV-m (Romania and at least 15 other countries across GDB regions), reflecting multiple evolution pathways for both HGT-mediated and other KpAMR character acquisitions. Correspondingly, the country-specific inferred orderings of KpAMR characters generally correlate, but rather weakly (R² = 0.41), with the simple prevalence of characters in that country’s source data (S3 Fig).

Across subsets of countries, interactions between our KpAMR resistance characters are observed (S4 Fig). For example, in around a quarter of cases, the acquisition of tigecycline resistance (Tgc_a) is inferred to promote the acquisition of glycopeptide resistance (Gly_a), and the acquisition of sulfonamide resistance (Sul_a) is inferred to promote the acquisition of trimethoprim resistance (Tmt_a). This latter observation reflects an independently known connection: sul and dfr genes providing Sul_a and Tmt_a resistance, respectively, are often genetically linked, and mechanistically act on connected steps in the folate synthesis pathway [35,36]. In some cases, the acquisition of MLS resistance (MLS_a) is inferred to repress the acquisition of carbapenem resistance (Bla_carb_a), suggesting possible interactions that could be clinically exploited in combination therapies (see Discussion). However, the inference of interactions between reversibly acquired features is challenging using HyperTraPS [31], so these observations must be regarded with some caution.

While S2 Fig gives a detailed “roadmap” of KpAMR evolution across countries, it is hard to immediately extract global-scale insights from this representation. We therefore next considered how to compare these results in a reduced-dimensionality picture.

Global consistency in evolutionary acquisition of resistance to a subset of drug families

To proceed, we used principal components analysis (PCA) to embed the high-dimensional set of posterior probabilities in a 3D space (Fig 2B). Here, the axes of highest variance in the inferred evolutionary pathways are identified and used to naturally lay out the individual country outputs. Comparison of country-by-country behavior along these axes can then reveal major sources of variability (and similarity) in inferred KpAMR evolution across countries. We underline that migration means that samples from one country likely represent a mixture of bacteria that have evolved in, and outside, that country, and that perfect separation between intrinsic dynamics of countries is therefore not possible with these data (see Methods).

Inferred evolutionary pathways for Kp across countries consistently involve early acquisition of characters conferring resistance to aminoglycosides, sulfonamides, trimethoprim, and betalactams (Agly, Sul, Tmt, and Bla characters). Aminoglycosides and sulfonamides were both discovered and entered clinical use early in the mid 20th century. Resistance properties inferred to be acquired at later orderings include extended-spectrum beta-lactam and phenicol resistance, acquired via mobile genetic elements (Bla_ESBL-a, Phe-a). Acquisitions inferred to be acquired late across countries include inhR-driven resistance (Bla_inhR-a, Bla_ESBL_inhR-a), and resistance to later- and final-line drugs including colistin, fosfomycin, tigecyclin, and glycopeptides (Col-a, Fcyn-a, Tgc-a, Gly-a). These form a set of “globally consistent” evolutionary observations.

Country-to-country differences in evolutionary acquisition of resistance to carbapenems, fluoroquinolones, tetracycline, and others

The general patterns above appear across evolutionary dynamics across countries and regions, with the first principal component of inferred Kp variability (accounting for 23% of the overall variance) simply reporting the precision with which these patterns can be characterized (Figs 2Ci and 3). In countries supporting less precise inference (for example, with limited datasets), the expected ordering of these characters is more uniform and closer to the average null hypothesis of all characters behaving equally; in countries with more precise output the expected orderings diverge to early and late values (Figs 2Ci and 3).

Once variability in this precision is accounted for (by considering the first principal component), variability corresponding to other differences in the inferred dynamics can be extracted. Several characters covary comparatively little with the (precision-aligned) PCA1 and instead show stronger covariance with PCA2 or PCA3, suggesting that there is genuine country-to-country variability in these characters after controlling for differences in observational noise (Figs 2Cii, 3, and 5). PCA2 (accounting for 15% of overall variance) displays strong covariance with carbapenem, fluoroquinolone, MLS, and tetracycline resistance and OMP mutation (Bla_Carb-a, Flq-m, MLS-a, Tet-a, Omp-m,). PCA3 (accounting for 8% of overall variance) displays covariance with Flq-a, Rif-a, SHV-m. None of these “globally divergent” characters display the “null-to-extreme” divergence seen in those characters most closely related to PCA1, suggesting genuine country-to-country differences in evolutionary dynamics, rather than observational differences, are responsible for this variability.

Geographical correlates with Kp evolutionary dynamics

To explore the sources of this evolutionary variability in more detail, we asked whether geographical region was linked to different variants of the inferred evolutionary dynamics (after accounting for observational differences). We first asked whether countries from different regions showed systematic differences in their PCA2 and PCA3 projections—corresponding to different orderings of the “globally divergent” country-specific characters above (Fig 4). We found that sub-Saharan Africa was a dramatic outlier in PCA2, involving notably later acquisition of carbapenem and fluoroquinolone resistance and OMP mutation than the other regions (Fig 4A and 4C). This ordering result agrees with the continuous-time history of treatment in the region, where carbapenem and fluoroquinolone use was established later (especially for children) than in other regions. Other regions were broadly consistent in PCA2 but showed substantial differences in PCA3. These differences correspond mainly to acquired resistance to fluoroquinolone and rifampicin, with substantially earlier inferred acquisition in South-Southeast-East Asia and Oceania, and substantially later in Central and Eastern Europe, Central Asia, and sub-Saharan Africa (Fig 4B and 4C).

Policy predictors of variance in AMR evolutionary dynamics

We next asked whether known policy of antibiotic use influence the inferred evolutionary dynamics from our analysis. We gathered statistics on drug use from 2016 to 2021 in countries available from WHO GLASS surveillance programme [37] and looked at correlations between per-capita drug use statistics and the inferred evolutionary dynamics of resistance to those drugs (Fig 4D). We found that for those drug types associated with PCA1 (where variability in inferred dynamics is largely determined by observation uncertainty) there was little correlation between drug use and the ordering of resistance acquisition. However, for carbapenems, associated with PCA2 (more non-observational country-to-country variability), there was a link between drug usage statistics and inferred resistance evolution, with carbepenam resistance (Bla-carb-a) more generally present before other acquisitions in countries with high drug use rates. This observation could arise through increased Bla-carb-a acquisition rate, or (perhaps more likely) a reduced loss rate of the Bla-carb-a feature once acquired (see Methods). The direction of this correlation was also observed in MLS, also associated with PCA2, though this relationship did not show statistical support at the p < 0.05 level. This pattern of correlations is consistent with the picture above: acquisition of some (PCA1) KpAMR characters is relatively constant across countries and circumstances, while other (PCA2) characters are more plastic and depend more tightly on country- and region-specific drivers. We note, however, that the GLASS data only cover a subset of years for which bacterial samples exist in our data (Fig 1A) and that KpAMR evolution is not constrained only within these dates (see Methods), so our drug use measure is only an imperfect estimate of the strength of any selective pressure shaping the dynamics.

Testing predictions of AMR evolutionary dynamics with newly sequenced, cross-decadal data

One motivating application of evolutionary inference to AMR problems is the ability to predict the next steps in evolutionary pathways [9]. A question of potentially direct clinical relevance is: given a clinically observed isolate with a given set of AMR characters, which character(s) will likely evolve next—and can treatment guidelines be adapted accordingly? Previous studies of AMR with HyperTraPS have demonstrated its ability to predict properties of evolutionary transitions in a withheld subset of the overall dataset, following the common training-test machine learning paradigm [21]. Here, we are in the unique position of using independently obtained Kp genome data to test the predictions of the evolutionary models trained above.

We sequenced Kp isolates from patients with bloodstream infections in 2001−2002 and 2017−2018 in Tanzania and 2015−2016 in Zanzibar, collected in previous studies [38–41] (Fig 5; see Methods). Bioinformatic analysis placed these new genomes in a phylogenetic tree with a subset of existing genomes from the training data, including their AMR profiles (Figs 5A and 6). Using only transitions independent of those in the training data (see Methods), we tested the predictions of the HyperTraPS model trained on data from Tanzania. Specifically, for each new transition, we queried the trained model about likely next steps from the ancestral state, and recorded the predicted ranks of the next steps actually present in the transition (Fig 5B). By far the most common outcome was that the true next step was ranked first (predicted to be most likely). For these data, we considered the number of instances where one model outperformed the other, in the sense of assigning a lower predicted rank (closer to the effective rank-1 “truth”) to a given transition. The trained model outperformed a naïve model (based only on the relative prevalence of KpAMR traits) in a limited subset (6/34 = 18%) of cases, and underperformed in 1/34 = 3% of cases. This is because the newly-observed data showed little evidence of higher-order effects (including bimodal ordering distributions and interactions between characters), which HyperTraPS can readily capture (demonstrated in previous applications to AMR data [21,20], but a prevalence-only model cannot (hence the variability in S3 Fig).

This analysis demonstrates that within-country predictions of future evolutionary dynamics are possible given a trained HyperTraPS model. We also asked whether the model could predict evolutionary dynamics at the regional scale. We use Zanzibar as a test case here—a region of Tanzania that was not explicitly represented in the original Kleborate training dataset. Our findings would predict KpAMR in Zanzibar to fall into the sub-Saharan Africa set of behaviors, with high values on PCA2 and associated signatures of late carbapenem resistance and others (Fig 4A). To this end, we applied HyperTraPS inference to newly sequenced isolates from Zanzibar. The resulting inferred dynamics clearly place Zanzibar within the set of sub-Saharan Africa behaviors (Fig 2B), with high values on PCA2 arising from later acquisition of carbapenem resistance, and high values on PCA3 arising from later acquisition of rifampicin resistance (Fig 4A and 4B).

We also explored the representation of different KpAMR characters in the snapshot genome sequences alone, without connecting with evolutionary dynamics. Fig 5C shows the relative prevalence of acquired KpAMR characters in existing and newly-sequenced data, not accounting for phylogenetic relatedness (and without representation of an explicit non-AMR control group). With this sampling, most characters showed an increase from 2001−2002to 2015−2016, particularly in mobile genetic acquisition of fluoroquinolone resistance, while the relative prevalence of acquired phenicol resistance and fluoroquinolone resistance mutations decreased. Proportions of these characters were in broad agreement with known coarser-grained observations from Tanzania [42,43] and existing Tanzanian genomes from Pathogenwatch, but we observed more acquired rifampicin resistance and lower acquired phenicol resistance than in existing samples. Our samples from Zanzibar showed lower relative proportions of aminoglycoside and β-lactam resistance, and higher prevalence of outer membrane protein mutations.

We have embedded the trained HyperTraPS KpAMR models in an online Shiny app https://stochasticbiology.shinyapps.io/amr-predict/. Here, a country of interest can be specified, and the KpAMR profile of a sample given by specifying presence or absence of each of our characters. The app will then give the inferred probabilities of each next possible character acquisition, as well as a transition graph describing the likely subsequent pathways for up to five acquisitions in the future (S7 Fig).

Existing Klebsiella data acquisition

We obtained all 47,721 Klebsiella pneumoniae sensu stricto records from pathogenwatch as of March 2024 [5]. The most common isolation years were 2015–2020, with a peak at 2018 (Fig 1A). We used AMR features (genes and mutations) reported by drug class as reported by Kleborate version 2.3 [26]. Kleborate groups genes or mutations known to confer resistance to clinically relevant antibiotic groups and related phenotypes. Betalactamases are further grouped by enzyme activity [26]. We consider a set of L = 22 drug classes for each genome, dichotomizing Kleborate output by the presence of any gene or mutation related to each group. That is, if there is any gene or mutation present for a resistance phenotype, the isolate is considered resistant, otherwise it is considered susceptible. This assumption holds for most resistance phenotypes considered in this study. However, fluoroquinolone resistance (Flq-m) often requires multiple mutations to cause resistance. Thus, the presence of these characters cannot be interpreted as a resistance phenotype, but rather an increased MIC (minimum inhibitory concentration). The isolates were divided by country based on metadata, and countries were grouped into superregions as defined by the Global Burden of Disease regional classification system [49]. A coarse-grained phylogeny was generated via LIN-codes (Life Identification Numbers: a systematic, stable quantitative system for genomic classification based on nucleotide identities [50], pre-computed for our dataset [51]) using lincoding.py supplied by [51]. 2,102 genomes were excluded due to missing metadata. Any country with less than two genomes available were excluded due to the inability to construct a tree. The resistance profiles and phylogenetic tree were combined into a phylogenetic tree annotated with resistance to drug classes.

Evolutionary pathway inference with hypercubic transition path sampling (HyperTraPS)

HyperTraPS [21] models an “evolutionary space” containing every possible state of a system involving L binary characters, then infers the probabilities of (Markovian) transitions between states in this space (here, involving the ordered accumulation of different AMR characters) that is most compatible with observed data [22,23]. A recent overview of the method in an AMR context is given in [9]. The source data is a collection of length-L “barcodes” labeling presence or absence of our characters—resistance to drug classes as defined in Kleborate—on the tips of a phylogeny. Ancestral state reconstruction (here, assuming that AMR resistance character accumulation is rare and irreversible, so that an ancestor possesses a character if and only if all descendants possess it) is by default used to infer ancestral states, and thereby construct a set of transitions from the data. This process is used solely to guard against pseudoreplication, and for these (largely deeply-branching) datasets, the inclusion of phylogenetic detail (and hence the details of ancestral state reconstruction) has a negligible influence on the inference outputs (S1B Fig). The resulting transitions are the input data for the HyperTraPS algorithm. The target of inference is an L × L matrix , where is the base rate for acquiring character i, and is the influence that acquisition of character j has on the base rate of acquisition of character i (as in mutual hazard networks [52]). Specifically, the probability of acquiring character i from the state corresponding to binary vector s is , where is over all characters absent in s [21]. In this way, positive and negative interactions between the acquisition of different AMR characters is supported, so that (for example) acquisition of resistance to drug X can make acquisition of resistance to drug Y more (or less) likely. The output of the inference process is, more broadly, a transition matrix , with elements describing the probability with which an isolate in state a will next transition to state b. For each of the 102 countries, territories and areas we ran HyperTraPS on three different random seeds to ensure that the results were consistent across runs (S2 Fig). Simulations were run for 10⁶ steps with 10³ random walkers for sampling, employing a penalized likelihood involving the HyperTraPS-estimated likelihood for a given set of transitions given a parameterization , computed with the probability term above calculated over sampled paths, and a unit penalty for each non-zero parameter (see [21]). The US genomes were sub-sampled to a subset of 1,000 genomes; three different subsamples clustered closer together than with any other countries in the PCA plot.

Synthetic control study for inference of reversible characters

To determine the effect of our assumption of irreversible accumulation of AMR characters (given that the loss of plasmids is not uncommon in Kp evolution [7]), we constructed a synthetic control study. A random phylogeny over 128 tips was constructed using a birth-death process with birth rate 1 and death rate 0.5 (chosen empirically to roughly match the phylogenies from our case studies) using phangorn [53]. A synthetic accumulation process reflecting a single evolutionary pathway was simulated on this phylogeny. Different loss processes—including uniform rate over a subset of features, heterogeneous, and random rates—were simulated in different instances (details in S1A Fig caption). The same HyperTraPS pipeline (assuming irreversibility), including ancestral state reconstruction and pathway inference, was run on both the pre-loss and post-loss datasets, and the outputs compared (S1A Fig). We also applied HyperMk [15], an EvAM approach supporting reversible dynamics, to a reduced subset of the KpAMR data (HyperMk is limited in the scale of data it can analyze), showing both a tight agreement with the HyperTraPS output and that model selection criteria favor the irreversible model (S1C Fig). These are a small-scale illustrative examples; the inference of reversible dynamics under an irreversible model is the focus of study in [31].

Global structure in inferred pathways

The outputs of HyperTraPS can be summarized as a matrix P, where element P_ij gives the posterior probability that character i is acquired at ordering j in a putative evolutionary process from an ancestor with no AMR characters towards a final state with all AMR characters. For example, P₃₂ gives the probability that the third character is the second to be acquired in such an evolutionary process. We obtained this matrix for each country in our dataset, then performed PCA on the set of matrices [54]. We observed the structure of matrices with “bubble plots”, where an array of points is plotted with areas proportional to P_ij.

Drug usage data

We collected antibiotic consumption data from 57 countries, territories and areas enrolled in WHO Glass surveillance program in the period 2016–2021 [37]. Overlaying the consumption data with the countries that had 3 or more available genomes reduced the number to 53. We use the median consumption statistics over the available time points for each country.

New Klebsiella data acquisition

The 79 newly sequenced Kp isolates in this article are derived from a collection of studies performed in Tanzania and Zanzibar over the last three decades. These comprise bacterial isolates from (a) blood samples from pediatric patients with bloodstream infections in Dar-es-Salaam, Tanzania in 2001–2002; (b) blood samples from adult patients with bloodstream infections in Dar-es-Salaam, Tanzania in 2017–2018; (c) fecal samples from the patients in (b); and (d) blood samples from adult and pediatric patients with bloodstream infections in Mnazi Mmoja hospital, Zanzibar in 2015–2016. These cultures were sequenced by MicrobesNG (MicrobesNG, Birmingham, UK) using Illumina HiSeq technology. Assembled contigs were provided by the sequencing service, which formed the initial step for the bioinformatics pipeline below.

New Klebsiella data analysis

The analysis pipeline begins with the contigs from the initial sequencing and assembly process. We first used dnadiff, which wraps nucmer from the MUMmer 3.0 package [55], to scaffold the contigs from a given isolate on the Klebsiella reference genome GCA_000240185.2_ASM24018v2 [56] and report statistics of the alignment. We retained only isolates with <20% unaligned bases as instances of Klebsiella. Combining the newly sequenced isolates with the existing Tanzanian genomes from Pathogenwatch, we then ran pairwise dnadiff across the combined set. We extracted the average identity across the aligned regions for each pair and recorded this as the ANI (average nucleotide identity). We then used d = 1 − ANI as a distance measure for phylogenetic tree estimation using the unweighted pair group method with arithmetic mean (UPGMA) [57] implemented in the phangorn [53] package in R (R Core 58]. We confirmed that an alternative method, neighborhood joining [59], did not give qualitatively different results for the predictions of AMR evolution. We used AMR profiles of existing and new isolates to reconstruct ancestral states across this complete tree, then identified the set of subtrees that had only newly-sequenced isolates at their tips. As transitions within such subtrees are independent of transitions in the (Pathogenwatch) training data set, we used this set of transitions as our independent test set.

Testing predicted dynamics from the trained model on new data

We queried the model fitted from the training data to predict the likely next states from each precursor state in this independent set of transitions [9,21]. We ranked the possible next steps from each state from most likely to least likely under the fitted model. We then recorded the ranks corresponding to the actual step(s) corresponding to the observed transition, and compared these to a null model where each possible transition could occur with probability given by its relative prevalence in the training data. For the Shiny app in S7 Fig, we report the probabilities of each possible next acquisition step given a provided current state, and construct the transition network starting from that current state and involving up to five further acquisitions.

Notes on nuances and resolutions in applying HyperTraPS to KpAMR data

The model underlying HyperTraPS pictures acquisitions as occurring irreversibly, while KpAMR characters—particularly those acquired through plasmid exchange—may be gained and lost many times by a lineage [20]. The fundamental reason we use an irreversible model is computational tractability. If reversible transitions are allowed, a potentially infinite set of evolutionary pathways (involving repeated acquisition and loss of characters) can connect any two states, complicating the computation of these pathways and their interpretation. In practise—including in AMR systems—only a limited number of the vast space of possible pathways play important roles in evolutionary dynamics [29,30], and powerful results have been demonstrated for accessibility of pathways sets on hypercubic evolution landscapes [60,61].

While applying irreversible models to systems involving reversibility increases the uncertainty of the inference process and challenges the inference of the detailed interactions between features, inferred evolutionary dynamics and particularly acquisition orderings is robust to this challenge—illustrated with some simple examples in S1A Fig, comparison with an EvAM approach allowing for reversibility [15] in S1C Fig, and the focus of study in [31]. This holds both for comparable loss rates across features and for distinct loss propensities, in which case the approach functions to infer the gain/loss rate ratio—which itself determines evolutionary orderings and prevalence patterns [31]. Generally, the outputs of HyperTraPS infer and report which features are likely present when a given feature is gained—a quantity that remains well-defined across irreversible and reversible cases. The congruence of this inference across reversible and irreversible approaches can be seen, for example, in S1C Fig.

We also note that the model does not assume that a final “end point” with all characters present is ever reached (or indeed possible), instead leaving uncertain all acquisitions that are not observed in real data. HyperTraPS invokes a “emission” picture, comparable to a hidden Markov model, to connect with observations [23], and therefore does not require complete information or homogeneous populations—it can naturally account for limited sampling of a broader, potentially heterogeneous population [32].

HyperTraPS does not a priori assume that the inferred evolutionary dynamics are a consequence of independent character accumulation (where each character’s observed prevalence is determined by a single acquisition rate) or interactions between characters (where the acquisition of a character is influenced by the presence of others) [31]. The “null model” of independence is challenged, for example, by instances of bimodality in inferred orderings – where a character may be acquired before many others or after many others, but not at intermediate stages. But the inferred evolutionary trajectories may be reported and analyzed regardless of which generative model is most supported. (We note, in parallel to this point, that simultaneous acquisition of multiple features is very compatible with our approach; the resulting output simply places equal probabilities over all compatible ordering pathways.)

Our samples are connected by a phylogeny, which we estimate with a taxonomy constructed from LINcodes (see Methods). Accounting for this phylogeny allows us to control against pseudoreplication, where individuals inheriting characters from a common ancestor are interpreted as independent acquisition instances of those characters [18–20,33]. Our approach uses the relative evolutionary timings from the source data phylogeny (for example, Fig 1Ci) to report the orderings of events: which acquisitions precede or follow which others, given that characters may influence each other, with the timescale of these events (see Discussion). In our case, much of the taxonomy corresponds to LINcode elements reflecting a deep timescale, with only separations within the final clusters corresponding to separation timescales under 1,000 years (following a back-of-the-envelope calculation in footnote below). For many KpAMR features, a picture of independent acquisitions is then not unreasonable, as many KpAMR acquisitions may have occurred on a shorter timescale (that of modern human medicine). However, our phylogenetically-embedded approach only biases the inference in the case of strong connection between acquisition patterns and sampling [31], otherwise simply increasing the uncertainty of the inferences [31,34]. In many real-world EvAM cases, whether or not phylogeny is included has an effect on uncertainty but no detectable impact on estimated orderings—demonstrated for two countries in our dataset in S1B Fig. As we control for the different uncertainties across countries in the followup PCA analysis, we use this approach to correct for potential instances of pseudoreplication.

The large majority of genomes with timing information were isolated in the period 2015–2020 (Fig 1A). Our approach here does not directly consider the specific, real-world timings of the evolutionary transitions involved. The precise timescales of the evolutionary process are in principle addressable using a version of our approach called HyperTraPS-CT (continuous time) [21], but this requires a consistent estimation of divergence times between all isolates in the dataset. Importantly, our approach here does not completely neglect a timescale—it is implicit in the phylogenetic relationships included in the source data, and this inclusion both accounts for the temporal ordering of events and “pseudoreplication” as above. However, this approach based on relative evolutionary timings means that we cannot directly tie differences in inferred evolutionary behavior to known real-world historical events shaping AMR profiles: shifts in drug policy, for example, or outbreaks or other events changing the effective size of the evolving population. While our results, focusing on the ordering of acquisition events, already demonstrate descriptive and predictive power in KpAMR evolution across countries, a continuous time picture on a smaller scale in the future (perhaps a single country or small set) would enable finer-grained linking between model prediction and the real-world details of policy and disease history.

Our samples are labeled by country of collection, which we treat as a covariate in our analysis. However, of course, bacteria and their hosts are not physically limited to remain within a given country. The phylogeography of Kp has been studied in depth, for example, in [6,62,63] and certainly involves mixing of strains across borders. Our analysis and interpretation do not assume that observations from a country all and exclusively reflect evolution in that country alone; our regional and correlative downstream analysis accepts that mixing can readily take place and attempts to identify variability that can be observed despite this mixing (which, in the limiting case, would lead to homogeneous inferred dynamics worldwide).

The terminal clusters in our LIN-derived taxonomy reflect ANI > 99.84%, giving a per-site substitution proportion of at most 0.0016. Estimates for Kp substitution rates are around 1.5 × 10⁻⁶ substitutions per site per year [64]. A substitution proportion of 0.0016 might then be expected over a timescale of 0.0016/1.5 × 10⁻⁶ ≈ 10³ years, reflecting an estimate for the longest timescale of separation within the terminal clusters.

S1 Fig. Limited issues with modeling reversible dynamics.

Plots involve a dataset (left) and “bubble plot” summary of inference as in the text (right). (A) HyperTraPS recovers accumulation orderings from reversible processes. Synthetic datasets involving the progressive accumulation, and loss, of characters X1 to X10, with different accumulation and loss rates (i)–(vi). For (i) accumulation rate = 1 and no loss take place. For (ii)–(vi) accumulation rate = 2 and (ii) loss rate is 0.1 for features 1–5 and 0 otherwise; (iii) loss rate increases linearly from 0.05 for feature 1 to 0.5 for feature 10; (iv) loss rate increases linearly from 0.2 for feature 1–2 for feature 10; (v) loss rate for each feature is drawn from U(0, 0.5); (vi) loss rate for each feature is drawn from U(0, 1). The output inferred from HyperTraPS accurately reports the accumulation ordering of characters, with differences in uncertainty, and ambiguity about characters which are not observed across all samples. This small-scale case study reflects a specific instance of the picture explored more generally in [31]. In all cases, inferred ordering can be interpreted as reporting which features are present when another is gained (see Methods). (B) Limited impact of neglecting phylogenies. Datasets from (i) Romania; (ii) Senegal. Inferences shown ignoring phylogeny (cross-sectional picture) vs. including phylogenetic information for ancestral state reconstruction and considering transitions down lineages. Here (and generally), many samples are related only via deep branches and are effectively independent; more recent relatedness within clades has a marginal effect on the uncertainty of posteriors but no noticeable effect on inferred orderings, in agreement with [34,31]. (C) Agreement of irreversible and reversible inference. (i) A reduced subset of the full KpAMR dataset, involving 5 features chosen to reflect a spread of prevalences: Bla_ESBL_a (1), Sul_a (2), Bla_chr (3), Tet_a (4), Rif_a (5). With this reduced set, the HyperMk approach [15] can be used to perform inference allowing for reversibility. Both the bubble plot summary (ii) and the inferred transition network (iii) for HyperMk (“Rev”) are very similar to the HyperTraPS outputs assuming an irreversible picture (“Irrev”), with HyperMk supporting reversible versions of all the same transitions identified (also observed in [31]); the Akaike Information Criterion favors the irreversible picture (393 vs. 553). The data and code underlying this figure can be found at https://doi.org/10.5281/zenodo.20408311.

https://doi.org/10.1371/journal.pbio.3003848.s001

(TIF)

S2 Fig. All posterior orderings of AMR character acquisitions.

This “global roadmap” summarizes “bubble plots” describing inferred evolutionary dynamics of KpAMR characters in different countries. Each row is a country, ordered vertically by descending number of Kp samples. Each column is a KpAMR character; the horizontal axis within columns gives evolutionary orderings (from early to late). The size of a point gives the probability that that character is acquired at that ordering in that country. For example, across the vast majority of datasets (including the top row, USA, for example), Bla_chr has a high probability of early acquisition and Gly_acquired has a high probability of late acquisition. In Cameroon, all characters have similar inferred evolutionary orderings, reflecting the limited data available to make more precise statements. Character names given in Fig 1 caption. The data and code underlying this figure can be found at https://doi.org/10.5281/zenodo.20408311.

https://doi.org/10.1371/journal.pbio.3003848.s002

(TIF)

S4 Fig. Inferred interactions between KpAMR characters.

An arrow from character X to character Y reports an inferred interaction between the evolutionary acquisition of those characters. Blue arrow: acquiring X makes acquiring Y more likely (“promoting”). Red arrow: acquiring X makes acquiring Y less likely (“repressing”). The number on each arrow gives the number of countries in which that interaction was inferred; only interactions appearing in >15% of countries with a posterior coefficient of variation under 1 are shown. Character names given in Fig 1 caption. However, as the inference of character interactions with HyperTraPS is challenging for reversible systems [31], these results should be treated with caution. The data and code underlying this figure can be found at https://doi.org/10.5281/zenodo.20408311.

https://doi.org/10.1371/journal.pbio.3003848.s004

(TIF)

S5 Fig. KpAMR characters determining aspects of global variability in evolutionary dynamics (alternative visualization).

Connected to Fig 3. KpAMR characters (horizontal axis) are plotted by their expected acquisition ordering (vertical axis) for each country (points). The country’s projection on PCA1, PCA2, and PCA3, for those characters that most strongly covary with each PCA axis, are given by the color of a point. As in Fig 3, those characters that covary with PCA1 form a consistent spectrum linked to precision: low values of PCA1 correspond to less precise, more uniform timing, while higher values lead to more precise earlier or later timings. Characters covarying with the other PCAs display wider ranges connected to PCA1-independent variability. Character names given in Fig 1 caption. The data and code underlying this figure can be found at https://doi.org/10.5281/zenodo.20408311.

https://doi.org/10.1371/journal.pbio.3003848.s005

(TIF)

S7 Fig. Shiny app for predictions of future KpAMR character acquisitions.

A country of interest is specified, and the KpAMR profile of a sample given by specifying presence or absence of each of our characters. The graphical outputs then give the inferred probabilities of each next possible character acquisition (which can be downloaded in CSV format), as well as a transition graph describing the likely subsequent pathways for up to five acquisitions in the future. Available at https://stochasticbiology.shinyapps.io/amr-predict/.

https://doi.org/10.1371/journal.pbio.3003848.s007

(TIF)