Role of mobile genetic elements in the global dissemination of the carbapenem resistance gene blaNDM

Abstract

The mobile resistance gene bla_NDM encodes the NDM enzyme which hydrolyses carbapenems, a class of antibiotics used to treat some of the most severe bacterial infections. The bla_NDM gene is globally distributed across a variety of Gram-negative bacteria on multiple plasmids, typically located within highly recombining and transposon-rich genomic regions, which leads to the dynamics underlying the global dissemination of bla_NDM to remain poorly resolved. Here, we compile a dataset of over 6000 bacterial genomes harbouring the bla_NDM gene, including 104 newly generated PacBio hybrid assemblies from clinical and livestock-associated isolates across China. We develop a computational approach to track structural variants surrounding bla_NDM, which allows us to identify prevalent genomic contexts, mobile genetic elements, and likely events in the gene's global spread. We estimate that bla_NDM emerged on a Tn125 transposon before 1985, but only reached global prevalence around a decade after its first recorded observation in 2005. The Tn125 transposon seems to have played an important role in early plasmid-mediated jumps of bla_NDM, but was overtaken in recent years by other elements including IS26-flanked pseudo-composite transposons and Tn3000. We found a strong association between bla_NDM-carrying plasmid backbones and the sampling location of isolates. This observation suggests that the global dissemination of the bla_NDM gene was primarily driven by successive between-plasmid transposon jumps, with far more restricted subsequent plasmid exchange, possibly due to adaptation of plasmids to their specific bacterial hosts.

Fig. 1. Composition of the global dataset of 6155 NDM-positive samples.

a Geographic distribution of bla_NDM-positive assemblies. Points are coloured by geographic region and their size reflects the number of samples they encompass. The world map was rendered from coordinates provided in rworldmap package in R. b Distribution of host bacterial genera of NDM-positive samples. c Distribution of sample collection years. d Distribution of contigs according to the plasmid backbone.

Fig. 2. A network of bla_NDM-carrying contigs (circles) mapping to the bacterial plasmid reference sequences (diamonds).

The network is visualized using Cytoscape and coloured according to communities identified by the Infomap algorithm. The largest communities are annotated according to the predominant plasmid type or the reference plasmid. Plasmids of <10 Kb in length include China pKP-YQ12450 (n = 610 contigs), Thai pKP-M11 (n = 399), Indian pSE-GMI17-002_2 (n = 354), Swiss pEC-C-F-163_C (n = 324), Peruvian pKP-NDM-1_isoform1-4 (n = 318), Peruvian pKP-NDM-1_isoform5 (n = 226), Thai pKP-M33 (n = 91), Indian pKP-MGR-K8 (n = 66), and 39 other <10 Kb putative plasmids.

Fig. 3. Schematic representation of the tracking algorithm splitting structural variants upstream or downstream of bla_NDM gene.

a A pairwise BLAST search is performed on all NDM-positive contigs. Starting from bla_NDM and continuing downstream or upstream, the inspected region is gradually increased using the ‘splitting threshold’. b At each step, a graph is constructed connecting contigs (nodes) that share a BLAST hit with a minimum length as given by the ‘splitting threshold’. Contigs that have the same structural variant at a certain position of the threshold belong to the same graph component, while the short contigs are singled out. c The splitting is visualized as a tree where branch lengths are scaled to match the position within the sequence, and the thickness and the colour intensity of the branches correspond to the number of sequences carrying the homology. For a more detailed explanation of the algorithm please refer to the Methods section.

Fig. 4. Splitting of structural variants downstream of bla_NDM.

The ‘splitting’ tree for the most common (i.e., ≥10 contigs) structural variants is shown on the left-hand side. The labels on the nodes indicate the number of contigs remaining on each branch. Labels of (yellow) branches with <20 contigs are not shown. The other contigs either belong to other structural variants or were removed due to being too short in length. The number of contigs cutting short is indicated by the area chart at the bottom. Similarly, the number of less common structural variants is indicated by the upper area chart. Genome annotations provided by the Prokka and Roary pipelines of the most common structural variants are shown in the middle of the figure. The homologous regions among structural variants are indicated by the grey shading. Some of the structural variants and branches were intentionally cut short even though their contigs were of sufficient size or longer. This was done to prevent excessive bifurcation and to make the tree easier to interpret. In particular, branches with more than 75% of contigs lost due to variation and short length were truncated. The distribution of genera, plasmid backbones and geographical regions of samples that belong to each of the common structural variants is shown on the right-hand side.

Fig. 5. Global prevalence and genetic context of most frequent putative (pseudo-) composite transposons and insertion sequences containing common regions (ISCRs) capable of mobilizing bla_NDM gene.

Transposons were defined as stretches of DNA flanked by two matching complete or partial ISs <30 Kb apart and enclosing bla_NDM. Putative pseudo-composite transposons were labelled according to their constituent ISs (IS15 and IS26). a Marginal distributions of transposon lengths or distances of ISCRs from bla_NDM start codon. b Overall counts of the frequent transposable elements (i.e., >25 representatives). c–e Bar plots, respectively, indicating proportions of plasmid backbones, bacterial genera and sampling location associated with most frequent transposable elements.

Fig. 6. Molecular dating of bla_NDM mobilization by Tn125 and Tn3000.

Posterior density distributions of ancestral sequence age (i.e., root height) for the Tn125 (a) and Tn3000 (b) transposons. The ancestral sequence emergence was estimated using two Bayesian tip-dating approaches implemented in BactDating and BEAST2. Three different population growth priors were used in case of BEAST2: Coalescent Constant Population, Coalescent Exponential Population, and Coalescent Bayesian Skyline as given by the colour scheme and legend at right. Median estimates with 95% highest density interval (HDI) are provided in the panel legends.