In this analysis, we searched all genomic studies investigating multidrug-resistant E. coli isolated from clinical samples in Vietnam. We eliminated studies that did not perform whole-genome sequencing using Illumina platforms or did not publish raw sequencing data, leaving 3 studies for inclusion. Raw sequencing data of 751 multidrug-resistant E. coli strains previously published by Roberts et al. (2022) (ENA Bioproject: PRJEB29424 and PRJEB28400), Hirabayashi et al. (2021) (provided by the authors), and Tuan-Anh et al. (2020) (ENA Bioproject PRJEB39354) were downloaded from the European Nucleotide Archive (ENA) database. Sixteen samples were from environmental surfaces (environmental isolates) and were collected at the National Hospital for Tropical Diseases, Hanoi, Vietnam. Among the other 735 clinical samples, 720 of the 721 isolates from the study by Roberts et al. (2022) were collected in two ICUs in Hanoi: 275 from the Bach Mai Hospital and 445 from the National Hospital for Tropical Diseases. One duplicate isolate (NHP1391) was excluded from the analysis. The two isolates (MH13 and MH17) from the study by Hirabayashi et al. (2021) were from the Military Hospital 103 in Hanoi, and the 27 isolates from the study by Tuan-Anh et al. (2020) were from the Children’s Hospital 1 in Ho Chi Minh city. Two E. coli isolates (TN1393 and XP817) were from our laboratory (NCBI SRA BioProject: PRJNA857185) and were from two patients hospitalized in the Saint Paul Hospital and the 108 Military Central Hospital, respectively. Most of the studied samples were collected between 2017 and 2019, except three samples (MH13, TN1393 and XP817) collected in 2012 and 2013. Each isolate is described in Supplementary Table S1.

The E. coli isolates from the study by Roberts et al. (2022) were classified into two phenotypes, ESBL-producing isolates and CR isolates, in function of their growth on selective media. The phenotype classification of the other isolates was based (i) on the antimicrobial susceptibility test results extracted from the studies by Hirabayashi et al. (2021) and Tuan-Anh et al. (2020), and (ii) for the two isolates from our laboratory (TN1393 and XP817), on antimicrobial susceptibility testing performed using ampicillin, amoxicillin-clavulanic acid, amikacin, piperacillin, ticarcillin, ticarcillin-clavulanic acid, cefalexin, cefoxitin, ceftazidime, colistin, fosfomycin, ofloxacin, imipenem, penicillin, and trimethoprim-sulfamethoxazole and the Kirby Bauer disk diffusion method, as previously described (EUCAST, 2020). In total, 4/751 strains were classified as non-ESBL-producing strains (see also Supplementary Table S1).

Genomic DNA (gDNA) of the two isolates from our laboratory was extracted with the Bacterial Genomic DNA Isolation Kit (Norgen Biotek Corp., Thorold, Ontario, Canada) following the manufacturer’s instructions. The gDNA samples were sheared randomly into small fragments with Covaris S/E210 or Bioruptor to prepare paired-end fragment libraries that were sequenced using an Illumina HiSeq 4,000 system at the Beijing Genomics Institute (Shenzhen, China). According to the information in the articles, gDNA was extracted with the QIACube and QIAamp 96 DNA QIACube HT kit (Qiagen, Hilden, Germany) and sequenced on a Illumina HiSeq X10 machine (Roberts et al., 2022); gDNA was extracted with the Wizard Kit (Promega) and sequenced on a Illumina MiSeq platform (Tuan-Anh et al., 2020); gDNA was extracted with the phenol-chloroform method and sequenced on a MiSeq (MH13 isolate) and on a HiSeq 4000 (MH17 isolate) apparatus (Hirabayashi et al., 2021).

Raw whole-genome sequencing data were analyzed using the bioinformatics workflow described in Figure 1.

Statistical analyses and data visualization were done with R version 4.0.4 (R Core Team, 2016). The numbers of putative plasmid contigs and of AMR genes were compared in resistance phenotype groups with the Kruskal-Wallis test followed by the Dunn’s post-hoc tests for multiple pair comparisons using the kruskal.test function and the dunn.test function. The Fisher’s exact test was used to compare AMR genes between CR isolates and ESBL-producing isolates through fisher.test function. All p-values from multiple testing were corrected with the Benjamini-Hochberg method, and adjusted p-values <0.05 were considered significant using the p.adjust function. Correlation plots, box plots and bar plots were generated using the ggplot function from the ggplot2 package version 3.4.0. Binary heatmaps were constructed with the Heatmap function from the ComplexHeatmap package version 2.9.4 (Gu et al., 2016).

The number of assembled plasmid contigs, using plasmidSPAdes, hugely varied among the 751 isolates, from 0 to 126 (median = 22, mean ± SD = 26 ± 19). In 8/751 isolates, no plasmid contig was identified. As expected, a moderate positive correlation was found between plasmid contig number and total assembly length (Pearson’s correlation coefficient, r = 0.55) (Figure 2A; Supplementary Table S1). Conversely, no correlation was found between plasmid contig number and sequencing coverage (r = 0.1) (Figure 2B; Supplementary Table S1), suggesting that the number of plasmid contigs in an isolate was linked to the presence of plasmids rather than to the sequencing coverage. Moreover, the number of plasmid contigs was higher in CR isolates (n = 197) than in ESBL-producing isolates (n = 550) (p = 0.008) and non-ESBL-producing isolates (n = 4) (p = 0.041) (Supplementary Figure S1). However, the highest number of plasmid contigs was found in an ESBL-producing strain (n = 126).

Among the 19,894 plasmid contigs identified with plasmidSPAdes in the studied isolates, 18,098 (90.97%) were matched in the plasmid database (putative plasmid contigs) and 1,565 (7.87%) were matched to the E. coli chromosome with at least 95% sequence identity (Figure 3A). The other 231 contigs (1.16%) could not be matched (i.e., unassigned contigs). However, for 223/231 unassigned contigs, at least one homolog was found in the plasmid database or in the E. coli genome database, but with a low BLAST identity (between 70 and 93%). For the other unassigned contigs, a BLASTn homology search against the non-redundant (nr) nucleotide database found that four contigs shared high similarity with the A. baumannii chromosome (sequence identity >98%), one contig shared 96% of identity with an Escherichia phage, and three contigs could not be aligned to any sequence.

Next, the closest match of the putative plasmid contigs in the plasmid database, which was considered as the plasmid homolog, was investigated. Plasmid homologs were mostly from the Escherichia genus (75.79%), followed by the Klebsiella/Raoultella, Salmonella, Vibrio, Enterobacter, Shigella, and Acinetobacter genera. These genera mostly belong to the Enterobacteriaceae family, with the exception of Vibrio that belongs to the Vibrionaceae family and of Acinetobacter that belongs to the Moraxellaceae family. Each genus was dominated by a major species (between 66.37 and 99.64% of plasmid homologs for that genus) (Figure 3B). This indicates that plasmids could be shared between genera, but the transmissions of plasmids are more frequently observed within the same species than those between distinct species. A set of 2,331 plasmid types that represented all plasmid homologs was used for the ANI analysis based on whole plasmid alignments to identify the genetic relatedness among plasmid homologs. This analysis revealed that 1,955/2,331 plasmids from different genera were clustered together with an ANI value >95%, indicating that these plasmid homologs were very closely related. Conversely, 120/2,331 plasmids were in smaller clusters, and 256/2,331 plasmids did not share any similarity with the other plasmids (Figure 3C).

In total, 5,554 plasmid-carried AMR genes were found in 83.98% (624/743) of isolates with plasmids. Among these genes, 99.06% (5,502/5,554) were detected in contigs that were matched in the plasmid database (Supplementary Table S2) and the other 52 AMR genes were detected in contigs that matched to chromosomes and in unassigned contigs. More than 60% of AMR genes identified in putative plasmid contigs belonged to the major facilitator superfamily antibiotic efflux pump, sulfonamide resistance, trimethoprim resistant dihydrofolate reductase, aminoglycoside acetyltransferase ANT(3″), resistance-nodulation-cell division antibiotic efflux pump, macrolide phosphotransferase, CTX-M β-lactamase, and TEM β-lactamase families (Supplementary Table S2 and Supplementary Table S2). The mean number of AMR genes detected in plasmid contigs for each isolate was nine. Overall, the number of AMR genes in plasmid contigs was higher in CR isolates than ESBL-producing isolates (p < 0.001). As expected, non-ESBL-producing isolates, which had a lower resistance phenotype, had fewer AMR genes than CR and ESBL-producing isolates (Figure 4A). Both CR and ESBL-producing isolates harbored many different plasmid-carried AMR genes belonging to >40 different AMR gene families. Conversely, only five AMR gene families [i.e., major facilitator superfamily antibiotic efflux pump, sulfonamide resistance, TEM β-lactamase, aminoglycoside 6-phosphotransferase APH (6), aminoglycoside 3″-phosphotransferase were found in non-ESBL-producing isolates]. Although CR and ESBL-producing isolates showed similar frequencies of AMR gene families across the putative plasmid contigs, the frequency of genes encoding OXA, NDM and KPC β-lactamases, which are implicated in carbapenem resistance, varied between CR and ESBL-producing isolates (Figure 4B and Table 1). The presence/absence of AMR genes classed into AMR gene families was further investigated in each isolate. AMR genes belonging to the major facilitator superfamily, sulfonamide resistance, trimethoprim resistant dihydrofolate reductase, and aminoglycoside acetyltransferase ANT(3″) families were highly represented in the studied isolates. While there was no clear separation between groups on the binary heatmap of AMR gene family across the studied strains, the presence of AMR genes in isolates was not likely associated with their resistance phenotype, sequence type (MLST), location, and study. However, isolates belonging to the sequence type 617 (ST617) could be grouped in clusters on the basis of their AMR gene profile, and these clusters included mainly CR isolates (Figure 4C).

Then, the analysis focused on 22 AMR genes that are associated with resistance to carbapenems. AMR genes belonging to the efflux pump and β-lactamase families were found in both CR and ESBL-producing isolates. Efflux pump genes (TolC, adeN, adeI, adeJ and adeK) were found in <1% of the studied strains. MarA and its homolog ramA, which encode bacterial transcriptional activators involved in regulating the AcrAB/TolC multidrug efflux pump, and porin genes were detected in ~11% of isolates. Carbapenemase-encoding genes (bla_CTX-M, bla_KPC, bla_NDM and bla_OXA) were detected in CR isolates and also in ESBL-producing isolates. Specifically, bla_CTX-M-27, a β-lactamase-encoding gene that leads to resistance to carbapenems, cephalosporins, and penams, was the most frequently detected CR gene in putative plasmids (9.64% of CR and 15.82% of ESBL-producing isolates). Its presence in plasmids could be responsible for high resistance to antibiotics in CR and ESBL-producing isolates. Moreover, the frequency of bla_KPC-2, bla_NDM-5, bla_OXA-1, bla_OXA-48, and bla_OXA-181 (β-lactamase-encoding genes) on putative plasmid contigs was higher in CR than ESBL-producing isolates (Fisher’s exact test: adjusted p-values <0.05). These genes are probably strongly associated with carbapenem resistance (Figure 5; Table 1). Additionally, >90% (n = 178/197) of the studied CR isolates carried at least one copy of the bla_KPC, bla_NDM, or bla_OXA gene in plasmid or chromosome contigs. Specifically, 67 of these CR isolates carried these genes only on plasmid contigs, 60 only on chromosome contigs, and 51 on both plasmid and chromosome contigs. Together, these findings indicated that carbapenem hydrolysis by KPC, NDM, and OXA carbapenemases is the major mechanism underlying carbapenem resistance in E. coli isolates with a significant contribution from plasmids.

Next, a sequence similarity network was built using 304 putative plasmid contigs that harbored bla_KPC, bla_NDM and bla_OXA genes to obtain a global view of the genetic organization of β-lactamase-encoding genes on plasmids. Plasmid contigs carrying the same β-lactamase-encoding gene shared a high degree of genetic similarity (sequence identity >99% and minimum query coverage of 90%) in pairwise alignments. All contigs carrying different bla_NDM genes were grouped into one cluster. Contigs carrying different bla_OXA genes showed a great diversity of genetic structures and were organized in four clusters. Moreover, plasmid contigs harboring bla_OXA-48 were classified in two sub-clusters. Although the clusters of the plasmid contigs carrying bla_OXA-1 and bla_KPC-2 did not share any sequence similarity, they were connected through a contig that contained both the bla_OXA-1 and bla_KPC-2 genes (Figure 6A). The genetic environment surrounding β-lactamase-encoding genes was similar in CR and in ESBL-producing isolates. In agreement with the result of the sequence similarity network, the structural organization of β-lactamase-encoding genes within each cluster was highly conserved across plasmid contigs (Figure 6B). Transposable elements (TEs; including transposons and insertion sequences), which play important roles in the spread of AMR genes, particularly bla_KPC, bla_NDM and bla_OXA, in healthcare settings (Partridge et al., 2018), were often located close to the β-lactamase-encoding genes.

The BLASTn homology search of contigs (assembled with plasmidSPAdes) against the NCBI RefSeq plasmid database and E. coli complete genomes (to determine the closest relative putative plasmid contigs) showed that apart from E. coli origin, putative plasmid contigs were found to be originated from other Enterobacteriaceae species or from the Vibrionaceae and Moraxellaceae families. More importantly, despite their different species/genus origin, theplasmid homologs shared high genomic similarity, which is broadly in agreement with a previous study on the prokaryotic plasmidome (Redondo-Salvo et al., 2020). These findings suggest that genetic exchange events by plasmid conjugation occur frequently within the Enterobacteriaceae family, especially in E. coli species but also between species. The potential plasmidic horizontal transfers between bacterial species, e.g., A. baumanni and E. coli have been reported (Chatterjee et al., 2017; Cooper et al., 2017). This aspect is fundamental to explain the ecology of AMR and to understand the emergence of drug resistance in the different bacterial species of public health interest.

ESBL-producing and CR isolates harbored many different AMR genes on their plasmids. The most frequently detected genes encoded antibiotic efflux pumps that eliminate antimicrobial agents from host cells (Kumar et al., 2020). Genes encoding CTX-M β-lactamases also were frequently identified on plasmids of ESBL-producing and CR isolates, as already previously observed in clinical ESBL-producing E. coli isolates (Ma et al., 2005; Ghenea et al., 2022). Strikingly, each CTX-M variant is widespread over a particular geographic area, thus providing epidemiological evidence for the transmission of β-lactamase-encoding genes within and between communities (Canton and Coque, 2006; Rossolini et al., 2008). In the present study, the bla_CTX-M-27, bla_CTX-M-15 and bla_CTX-M-55 genes were widely represented on plasmids of ESBL-producing and CR strains. Previous studies showed that the genes encoding these CTX-M variants are predominant in ESBL-producing E. coli isolates from humans and farm animals in Vietnam (Biedenbach et al., 2014; Nguyen et al., 2019). CTX-M-15 is the dominant CTX-M variant worldwide. Other CTX-M variants (e.g., CTX-M-55, CTX-M-27) that were initially specific to some world regions have now spread globally due to international travel (Bevan et al., 2017; Nakayama et al., 2020). Thus, monitoring the emergence of plasmid-carried bla_CTX-M genes, particularly genes encoding new variants, could provide important information for developing strategies to reduce the transmission of AMR genes via plasmids in clinical settings.

Plasmid-carried carbapenemase-encoding genes are considered the main cause of the resistance phenotype in CR E. coli. In our analysis, genes encoding the KPC-2, NDM-5, OXA-1, OXA-48, and OXA-181 β-lactamases were more frequently detected on plasmids in CR than ESBL-producing isolates. The bla_KPC gene was first detected on a plasmid of a carbapenem-resistant K. pneumoniae strain in the southern United States in 2001 (Yigit et al., 2001). It has then rapidly spread to other bacterial pathogens, including E. coli, worldwide, including in Vietnam (Petrella et al., 2008; Woodford et al., 2008; Richter et al., 2011; Linh et al., 2021). Similarly, plasmid-carried bla_NDM-5 was identified in patients in the United Kingdom and the United States (patient hospitalized in India) in 2016–2017 (Hornsey et al., 2011; Rojas et al., 2017), and in multidrug-resistant clinical E. coli strains in many countries (Pitart et al., 2015; Soliman et al., 2016; Zou et al., 2020). In Vietnam, plasmid-encoded bla_NDM-1 was detected in clinical Enterobacteriaceae isolates from a surgical hospital in 2010 (Tran et al., 2015). To our knowledge, Roberts et al. (2022) were the first to describe the presence of bla_NDM-5 in clinical isolates in Vietnam. In our study, plasmid-carried bla_NDM-1, but not bla_NDM-5, was detected in the XP817 and MH-13 isolates collected in 2012–2013, whereas plasmid-carried bla_NDM-5 was detected in clinical samples isolated in 2017–2019. This suggests the emergence of plasmid-carried bla_NDM-5 in Vietnam approximately in 2017, close to the year of its first report in other countries.

Plasmid-encoded OXA β-lactamases can hydrolyze β-lactam substrates, including carbapenems (Evans and Amyes, 2014). The KPC, NDM and OXA enzymes can mediate carbapenem resistance. Their hydrolytic activity depends on the expression level of the genes encoding these enzymes that is mostly influenced by their genetic environment and copy number (Roth et al., 2011; Evans and Amyes, 2014; Paul et al., 2017). This could explain why some ESBL-producing isolates that harbor these β-lactamase-encoding genes do not show any phenotypic resistance to carbapenems.

The sequence similarity network based on pairwise alignments of plasmid contigs harboring the bla_KPC, bla_NDM and bla_OXA genes captured a high degree of similarity between contigs carrying the same AMR gene. This confirmed the results by Roberts et al. (2022) showing that AMR genes are transmitted within and between hospital settings in Vietnam. The genetic environment around the bla_KPC-2 and bla_NDM genes is largely consistent with what reported by previous studies (Dortet et al., 2012; Cuzon et al., 2013; Martins et al., 2020). Specifically, our analysis showed that bla_KPC-2 is located between two TEs, whereas bla_NDM genes are close to the bleomycin resistance gene located downstream of a TE and share the same promoter. Conversely, bla_OXA genes were found in several clusters with different environmental structures, as previously reported by Evan & Amyes (Evans and Amyes, 2014). This could be explained by spontaneous DNA rearrangements following acquisition of the bla_OXA by horizontal transfer.

Our study, based on short-read whole-genome sequencing data, provided general and predictive information on plasmid profiles in 751 ESBL-producing and CR E. coli isolates. This work has certain limitations. We could only assemble partially plasmids from short-read data using the existing assembly tool and therefore part of plasmid information remains unknown. Another limitation of our study was sample isolation bias since most of the studied strains were isolated from patients of two hospitals in northern Vietnam. Nevertheless, the current study serves as a good foundation for further analysis on a larger scale to have a more detailed look into resistance plasmid epidemiology at the national and global level.