A megaplasmid family driving dissemination of multidrug resistance in Pseudomonas

Abstract

Multidrug resistance (MDR) represents a global threat to health. Here, we used whole genome sequencing to characterise Pseudomonas aeruginosa MDR clinical isolates from a hospital in Thailand. Using long-read sequence data we obtained complete sequences of two closely related megaplasmids (>420 kb) carrying large arrays of antibiotic resistance genes located in discrete, complex and dynamic resistance regions, and revealing evidence of extensive duplication and recombination events. A comprehensive pangenomic and phylogenomic analysis indicates that: 1) these large plasmids comprise an emerging family present in different members of the Pseudomonas genus, and associated with multiple sources (geographical, clinical or environmental); 2) the megaplasmids encode diverse niche-adaptive accessory traits, including multidrug resistance; 3) the accessory genome of the megaplasmid family is highly flexible and diverse. The history of the megaplasmid family, inferred from our analysis of the available database, suggests that members carrying multiple resistance genes date back to at least the 1970s.

Fig. 1. Maps of the pBT2436 and pBT2101 AMR regions.

The coordinates of the AMR regions within the megaplasmid sequences are shown flanking the maps. ORFs and repeats, represented as coloured arrows and blocks, respectively, are indicated in the maps along with the location of the detected mercury resistance operons and complete integrons and clusters of attC sites lacking integron-integrases (CALIN) as identified by the Integron Finder tool. Insertion sequences (ISs) recognised with ISfinder are represented as orange blocks in the maps along with their corresponding names. Different contiguous IS identified within a plasmid region are separated by dotted lines. Maps are drawn to scale. Gene names or products of selected ORFs of interest are indicated above their corresponding arrows. Gene duplications occurring within the same AMR region are numbered and distinguished with an “_N” suffix in their names. Darker shade-coloured arrows denote ORFs shared by the two megaplasmids’ RR1 regions. The different identified repeats are numbered according to their first occurrence in the AMR region and a suffix letter (a–e) is included to distinguish between them. Repeat size in base pairs is displayed next to their corresponding blocks. Repeat blocks sharing matching colours between the two megaplasmids’ RR1 regions represent shared regions.

Fig. 2. Genome comparison of megaplasmids of the pBT2436-like family.

Fourteen complete Pseudomonas megaplasmid sequences, indicated at the bottom of the figure and represented as colour rings, were aligned to the pBT2436 (left) and pBT2101 (right) genomes at the nucleotide level. Solid circles denote sequence homology to pBT2436/pBT2101, whereas gaps within the rings correspond to regions lacking sequence similarity. Innermost rings (black) represent the GC content deviation from the average in the reference genomes. The three outermost rings (from innermost to outermost) indicate the location in pBT2436/pBT2101 of the core genes identified from the pangenome analysis of the family (grey), AMR genes (red), and genes encoding integrases or transposases (blue). Names of genes of interest and the position of the pBT2436/pBT2101 AMR (RR) and variable (VR) regions are indicated in the figure.

Fig. 3. Rarefaction curves for the core and pangenome of the pBT2436-like megaplasmid family.

Graphs show the estimated size of the core (left) and pangenome (right) of the pBT2436-like group. Plotted data come from sampling experiments with ten random-seeded replicas from the BDBH- (core) or OMLC-based (pan) protein clustering. Fitted curves follow the Tettelin function.

Fig. 4. Phylogenomic analysis of complete megaplasmids of the pBT2436-like family.

Maximum-likelihood unrooted phylogenetic trees displaying the relationships among the megaplasmids reported in this work (blue branches) and 13 Pseudomonas homologous plasmids detected in GenBank. Branches corresponding to plasmids from non-aeruginosa species are coloured in orange. Pangenome phylogeny (left) was estimated from the pangenome matrix and reflects the evolutionary relationship of the megaplasmids in terms of their gene content, i.e. presence/absence patterns of the 1164 clusters contained in the matrix. The core genome phylogeny (right) was estimated from the concatenated set of 105 top-ranking alignments from core genes selected as phylogenetic markers and depicts the patterns of divergence of the megaplasmid genomic backbone. Taxa clustered together in both trees are highlighted with dotted-line ovals of matching colours. Purple dots represent approximate Bayesian posterior probability values >0.99 or >0.81 for the corresponding nodes in the core genome and pangenome trees, respectively. Scale bars correspond to 0.1 expected substitutions per site under the binary GTR+FO (pangenome tree) or best-fitting GTR+F+ASC+R2 (core genome tree) models. Where known, sources are indicated by nation flags. Asterisks in plasmid names denote cases where the geographical source is putative, based on the information of the submitter’s country in GenBank.

Fig. 5. AMR gene content of the pBT2436-like megaplasmid group.

The heatmap shows the collection of AMR genes identified in 15 members of the pBT2436-like megaplasmid family through BLASTn searches against the Comprehensive Antibiotic Resistance Database (CARD). Percentage of coverage of the AMR genes identified from the searches is indicated. The megaplasmids (Y axis) are hierarchically clustered based on their content of AMR genes (X axis) using the “complete” method with Euclidean distance. The “_N” suffix in gene names denotes cases where the gene occurs more than once in the same megaplasmid. AMR genes are additionally classified based on the drug class they confer resistance to according to CARD.

Fig. 6. pBT2436-like megaplasmids recovered from Pseudomonas genome assemblies deposited in GenBank.

Nucleotide sequences of Pseudomonas genomes from the four GenBank assembly categories were aligned to pBT2436 to identify unreported or overlooked related megaplasmids. Search was performed in November 2018. a The chart shows the pBT2436 estimated coverage and collection date reported in the metadata (when available) of matches covering >60% of the pBT2436 sequence (Supplementary Data 4). Bacterial species of the matches are colour coded and indicated in the figure. b Phylogeny of pBT2436-like megaplasmids recovered from Pseudomonas assemblies. The maximum-likelihood unrooted tree (GTR+F+ASC+R2 model) was inferred from the alignment of 4 core genes selected as phylogenetic markers (see “Methods”) and represents the relationships among 53 sequences from GenBank displaying >80% pBT2436 coverage and 15 complete megaplasmids (bold letters). Dotted lines connect branch tips with the corresponding taxa names and sources (three-letter codes in parentheses, where known). Lines corresponding to non-aeruginosa species’ sequences are colour coded as in a. Purple dots indicate posterior probability values >0.7. The scale bar corresponds to 0.1 expected substitutions per site. No clear correlations between cluster composition and geographical origin or bacterial species of its taxa are observed.

Fig. 7. pBT2436-like megaplasmid conjugation and fitness cost.

a The relative fitness (growth in competition with plasmid-free) of P. fluorescens SBW25 with one of two different pBT2436-like megaplasmids was measured from quadruplicate experiments for three independent transconjugants. Filled points show separate measurements; the open point and bars indicate means calculated for each transconjugant and standard error. pQBR103 and pQBR57 (one transconjugant each), known to be costly, are included as a positive control for plasmid cost. We detected a significant, positive effect of pOZ176 acquisition on fitness (Tukey’s HSD pOZ176 vs no-plasmid comparison padj = 0.009) and no significant effect of p1 acquisition (padj = 0.44). b Plotted are the ratios of transconjugants (T) to recipients (R) from the endpoints of the competition assays (see “Methods”). Note that we retrieved transconjugant data from only 2/4 experiments with the mercury-resistant plasmids (p1, pQBR103, pQBR57). Source data is provided as a Source Data file.