A cross-sectional study was conducted in Shihezi City located in the Xinjiang Autonomous Region of China. The study focused on two primary locations: the Shihezi University dairy farm that specializes in breeding dairy cows using advanced technologies and a privately owned dairy farm on the outskirts of the city. Shihezi features a semi-arid climate with extreme seasonal temperatures and fertile plains, and it relies on river irrigation that provides favorable conditions for dairy farming (Yang et al., 2021). The study population consisted exclusively of dairy cows, chosen for the isolation of E. coli from their intestinal tracts to investigate the prevalence and characteristics of the bacteria. Fresh fecal samples were collected immediately after excretion from the middle of the fecal matter to prevent environmental contamination. Samples were placed in 50-ml centrifuge tubes, chilled on ice for preservation, and transported to the laboratory for immediate processing (Lupindu, 2017).

Escherichia coli was isolated from fecal matter using a modified protocol based on the guidelines of Alcock (Ju and Willing, 2018). The isolation was performed via culture techniques on MacConkey, Eosin Methylene Blue, and Luria Bertani agar using the streaking method. Pure isolates were obtained using subculturing that was repeated three times. Colony PCR was then employed for the molecular identification of the 16S rRNA gene using universal primers 27F and 1492R to amplify all nine variable regions (Lane et al., 1985). The PCR was conducted in a 25-μl volume under specific thermal cycling conditions that included an initial denaturation at 95°C followed by 35 cycles of denaturation, annealing, and extension. Gel electrophoresis was performed to confirm the amplification and integrity of the 16S rRNA gene. After validation, 16S rRNA gene sequencing was performed. Raw sequence data were processed using BioEdit (Hall, 1999) and Trimmomatic (Bolger et al., 2014) to remove low-quality reads. This was followed by sequence alignment and analysis in MEGA11 (Tamura et al., 2021) to identify single nucleotide polymorphisms (SNPs), conserved sites, and variable regions, facilitating the identification of the target bacteria.

Based on published data regarding drugs to which E. coli is known to be susceptible (Matuschek et al., 2014; Hombach et al., 2013), we selected 14 different antibiotics for testing. Antimicrobial susceptibility testing was conducted using the Kirby-Bauer disk diffusion method, which involves placing a 6-mm filter paper disk impregnated with a specific amount of antibiotic on a culture plate pre-inoculated with a bacterial suspension (Hudzicki, 2009). Bacterial strains were initially prepared in suspensions, and the optical density was checked at 625 nm to achieve a 0.5 McFarland standard (Benkova et al., 2020). Absorbance values within the range of 0.08–0.13 were considered acceptable according to Hudzicki’s guidelines (Hudzicki, 2009).

After incubation, the diameters of the zones of inhibition around each antimicrobial disk were measured using a ruler and recorded in millimeters. These measurements were used to classify each isolate as susceptible, intermediate, or resistant to the antibiotics tested (Matuschek et al., 2014; Hombach et al., 2013).

Python scripts were used to analyze the data generated from the antimicrobial tests. The Python libraries utilized in this research included Pandas for data manipulation (McKinney, 2010); NumPy for numerical computation and array operations (Harris et al., 2020); Matplotlib for creating static, interactive, and animated visualizations (Hunter, 2007); and Seaborn for its high-level interface for drawing good and informative statistical graphics (Waskom, 2021).

Two multidrug-resistant isolates, designated as E.coli_30 and E.coli_45, were selected for genome sequencing. These two strains were chosen based on their distinctive resistance patterns observed during the susceptibility testing. E.coli_30 demonstrated resistance to all tested antibiotics except aztreonam, while E.coli_45 was resistant to all antibiotics except amikacin. This added layer of resistance compared with the other MDR isolates highlights their unique resistance profiles. This made them ideal candidates for whole-genome sequencing (WGS) to further investigate the genetic mechanisms underlying their extreme drug resistance.

DNA samples were prepared for sequencing by generating indexed libraries for each isolate. The genomic DNA of the two MDR E. coli was fragmented using sonication to an average size of 350 base pairs (bp). The fragmented DNA underwent a series of modifications, including end repair, A-tailing, and ligation with Illumina full-length adapters. This was followed by PCR amplification and purification of the resulting products. The library quality was assessed using qPCR to ensure integrity and concentration. Libraries that met the quality standards were pooled according to their effective concentrations and the desired data output. Sequencing was then performed on an Illumina platform using a paired-end 150 bp (PE150) strategy (Caporaso et al., 2012).

Quality control of raw sequencing data was conducted using FastQC v0.12.1 (Andrews, 2010) followed by de novo genome assembly using SPAdes genome assembler v3.15.5 (Bankevich et al., 2012). To detect and assemble plasmid DNA, PlasmidSpades pipeline was used, and plasmidFinder-2.1.6 (Carattoli et al., 2014) was employed to confirm the presence of known plasmid replicon sequences in the genomes. The assembled genomes were further evaluated using QUAST v5.2.0 (Gurevich et al., 2013), to ensure the selection of high-quality sequences for downstream analysis. Gene annotations were performed using Prokka v1.14.6 (Seemann, 2014) to identify coding sequences and gene structures. A graphical map of the two genomes and the plasmid was generated using nucleotide sequences of the annotated contigs (fna file) in fasta format. This was achieved using one of the CARD web tools at https://proksee.ca/tools/card (Alcock et al., 2023).

To identify and characterize antimicrobial resistance genes and associated resistance mechanisms, we utilized the Resistance Gene Identifier (RGI) tool v6.0.3 (Alcock et al., 2023; Jia et al., 2016), which analyzes the sequences against the CARD database (Alcock et al., 2020). To effectively visualize and interpret the data, we generated RGI wheels by uploading the RGI results with json extension (JSON file) to the CARD website (Jia et al., 2016).¹ To identify and characterize the virulence genes, we used Abricate 1.0.1 (Seemann, 2024a), which analyzed the sequences against the Virulence Factor Database (VFDB) (Chen et al., 2016).

In addition to the plasmid identified in one of the genomes, other mobile genetic elements, such as IS, prophage regions, and more, were investigated using MGEfinder (Durrant et al., 2020), MobileOG-db (Brown et al., 2022), Virsorter (Roux et al., 2015), Alien Hunter (Vernikos and Parkhill, 2006), PhiGARo (Starikova et al., 2020), and ISfinder (Siguier et al., 2006). This investigation was essential, as many AMR and virulence genes are disseminated via mobile elements. Their presence and roles helped to explain the significant genetic diversity observed between our isolates.

The susceptibility test revealed that 51.1% of the isolates were susceptible to the tested antibiotics, 26% showed an intermediate response, and 22.9% were resistant (Figure 1A). By direct count, 26 out of the 53 isolates tested were multidrug-resistant. The highest level of resistance was observed against imipenem (IMP) (60%). In contrast, amikacin (AMK), ofloxacin (OFX), levofloxacin (LEV), ceftriaxone (CTR), and ciprofloxacin (CIP) demonstrated a more balanced distribution across susceptibility, intermediate, and resistance categories. Trimethoprim/sulfamethoxazole (SXT), ceftazidime (CAZ), chloramphenicol (C), aztreonam (ATM), streptomycin (S), gentamicin (GM), ampicillin (AMP), and tetracycline (TET) were the most effective drugs, with GM and TET showing the highest efficacy against over 65% of the E. coli isolates (Figures 1B,C).

Furthermore, we calculated the correlation coefficients for all the antibiotics based on the isolates’ responses to assess potential relationships between them. Strong positive correlations were observed in the resistance profiles between certain antibiotics, such as STX and AMP, C and AMP, and TET and STX, with correlation coefficients of 0.85, 0.84, and 0.83, respectively (Figure 1D). These strong positive correlations suggest that isolates resistant to one antibiotic are often resistant to the other, possibly due to similar resistance mechanisms, such as those found within the same antibiotic class. Conversely, some antibiotics, such as GM and LEV, ATM and TET, and CAZ and TET, showed little to no correlation, indicating distinct mechanisms of action (MOA) or resistance. In some cases, we also observed negative correlations, suggesting that the isolates exhibited different MOA.

Finally, we clustered the isolates based on their median score/response to antibiotics and identified distinct response patterns (Figure 1E). Nineteen strains were classified as susceptible, 18 as intermediate, and 16 as resistant. A key observation was that isolates in the first cluster (resistant strains) exhibited resistance to multiple antibiotics. This result indicates that these strains may require special attention to manage the spread of resistance effectively. Additionally, some isolates classified as MDR negative or positive appeared to have intermediate responses due to their median scores. Finally, we also observed that some MDR isolates appeared in the susceptible category. This was because, despite being MDR, they were susceptible to most of the antibiotics tested and this resulted in a susceptible median score.

E.coli_30 and E.coli_45 were selected for genome sequencing after susceptibility testing, as they demonstrated resistance to most of the antibiotics tested, in contrast to the other isolates. After sequencing, de novo assembly of E.coli_30 (Figure 2A) revealed the number of scaffolds (75) was slightly lower than the number of contigs (84). The N50 value was significantly higher for scaffolds (344,767 bp vs. 246,478 bp), indicating better assembly continuity. The total lengths of contigs and scaffolds are similar, at approximately 5.12 million bp, with an identical GC content of 50.92%. However, the scaffolding process introduced gaps, as indicated by an increase in ‘N’s per 100 kbp (17.8 in scaffolds compared to 0.02 in contigs). Annotation of this genome revealed 4,767 genes, including 4,663 coding sequences (CDS), eight rRNA genes, 95 tRNA genes, three repeat regions, and one tmRNA.

For the plasmid identified in E.coli_30 (Plasmid30), the total length of the assembly was approximately 218,143 bp, with the largest contig measuring 58,599 bp. The N50 value was 23,145 bp, while the N90 value was 3,131 bp for contigs and 3,150 bp for scaffolds. The GC content remained consistent at 50.58%, and no mismatches were found in the contigs. In contrast, the scaffolds exhibited 47.05 ‘N’s per 100 kbp and a total of 100 ‘N’s, indicating gaps introduced during the scaffolding process. In general, we can say that the scaffolding process slightly reduced the number of contigs, maintained the assembly’s length and continuity, and introduced a few gaps. Figure 2B shows the map of this plasmid, and it highlights various genes and elements that contribute to its functionality and antibiotic resistance.

The second genome, depicted in Figure 2C, had a total assembly length of approximately 4.99 million bp, with a GC content of 50.57%. The number of scaffolds (176) was slightly lower than the number of contigs (185), and the N50 value was higher for scaffolds (151,645 bp vs. 150,606 bp), indicating better assembly continuity. However, the scaffolding process introduced gaps, as reflected by the significant increase in ‘N’s per 100 kbp (15.83 in scaffolds compared to 0.02 in contigs). Overall, the scaffolding effectively linked contigs into longer sequences and enhanced the assembly’s utility for downstream analyses despite the expected introduction of gaps. Annotation of the scaffolds identified 5,076 genes (several were hypothetical proteins), including 4,953 CDS, eight rRNA genes, 114 tRNA genes, four repeat regions, and one tmRNA.

PlasmidFinder-2.1.6 (Carattoli et al., 2014) was used to confirm the presence of known plasmid replicon sequences within the plasmid genome identified in E.coli_30. Four plasmid replicon sequences were detected: Col440I with 94.74% identity over 114 bp in the 14th locus, IncFIB (AP001918) with 98.39% identity over 682 bp in the 7th contig, IncFII (pCoo) with 96.56% identity over 262 bp in the first contig (contig node 1), and p0111 with 98.64% identity over 885 bp in contig node 5. The presence of these plasmid replicon sequences confirmed that E.coli_30 harbors plasmid within its genome. No hits were detected for plasmid types Inc18, NT_Rep, Rep1, Rep2, Rep3, RepA_N, RepL, and Rep_trans, indicating the absence of these plasmid types in the isolate. A nucleotide BLAST search was conducted to identify the origins of these replicons and the result revealed that all the sequences were from E. coli, except for Col440I, which originated from Klebsiella pneumoniae.

RGI CARD investigations of the two genomes revealed the presence of multiple AMR genes (Figure 3), each contributing to antibiotic resistance through various mechanisms. Both E. coli genomes exhibited MDR, but the plasmid in E.coli_30 (Plasmid30) adds a layer of complexity to its resistance mechanisms. While both genomes rely heavily on efflux pumps as a primary resistance mechanism, E.coli_30 was found to have additional resistance factors from its plasmid.

Both genomes contained resistance-nodulation-cell division (RND) efflux pump systems (acrA, acrB, AcrF, mdtA, mdtB, and mdtC) that confer resistance to multiple classes of antibiotics, including fluoroquinolones, cephalosporins, tetracyclines, and macrolides. The major facilitator superfamily (MFS) pumps (e.g., emrK, emrY, mdfA, and mdtM) in both genomes offer resistance to tetracyclines, phenicols, and other antibiotics. The genomes also include ATP-binding cassette (ABC) transporters like msbA and TolC that contribute to resistance to nitroimidazoles and peptide antibiotics. Plasmid30 significantly contributes to E.coli_30 resistance, particularly to macrolides, phenicols, and fluoroquinolones (Figure 3C). It carries macrolide phosphotransferase resistance gene (mphA) that is not found in E.coli_45 and is involved in antibiotic inactivation.

Although both genomes had MDR profiles, E.coli_45 displayed a broader range of efflux pump systems and additional resistance mechanisms (Table 1), such as protein overexpression and regulatory mutations (e.g., soxR, soxS, and marA). These mutations enhance antibiotic efflux and reduce β-lactam permeability. Moreover, E.coli_45 includes more entries related to aminoglycoside resistance through inactivation mechanisms (e.g., APH(3′)-IIa, and APH(3″)-Ib). Supplementary Tables S1, S2 summarize the detailed RGI findings for both genomes.

Abricate analysis of the two genomes revealed not only similarities in their pathogenic mechanisms; but also some differences in specific gene clusters (Supplementary Tables S3, S4). Both E.coli_30 and E.coli_45 rely on several core virulence factors, particularly in adhesion, nutrient acquisition, and secretion systems. For instance, both strains possess genes involved in type 1 fimbriae (fimA, fimB, and fimC), E. coli common pilus (ykgK/ecpR and yagZ/ecpA), and colonization factor antigens. Both strains also possess the curli fiber system (csgA, csgB, and csgC) that is important for biofilm formation. They share similar enterobactin siderophore systems for iron acquisition, with genes such as entA, entB, entC, fepB, fepC, fepD, fepE, and fepG. However, the salmochelin siderophore system was found only in E.coli_30 and its plasmid, giving it an advantage in nutrient acquisition over E.coli_45.

Type III secretion systems (TTSS) were detected in both genomes, with genes such as espX1, espX4, and espX5 encoding effector proteins that manipulate host cell signaling and immune defenses. However, E.coli_45 showed a slightly broader range of TTSS-related genes, such as espL1 and espR1, indicating a wider repertoire of effector proteins. Furthermore, we detected a type VI secretion system (T6SS) in E.coli_45, with genes such as tssA and hcp1/tssD1 that were absent in E.coli_30. T6SS gives E.coli_45 an advantage in bacterial competition and host interaction by enabling it to deliver toxic effectors to other bacteria or host cells.

Both strains relied on the enterobactin siderophore system for iron/nutrient acquisition. However, E.coli_45 was found to have additional enterobactin genes, such as entS, fepD, and fepE, making its siderophore system potentially more efficient than that of E.coli_30. We also identified an unknown protein related to the TraJ family that is encoded by AAA92657 in E.coli_45 and is involved in the regulation of invasin.

Several MGEs were identified in the two MDR E. coli genomes. MGEfinder was used to detect several insertion sequences and their insertion sites in the two genomes. Some insertions were very short sequences, approximately 20 bp, while others were larger, ranging from 1,000 to 3,000 bp. Annotation of the larger regions revealed the presence of IS110 family transposases, IS3, and several passenger genes, including invasin, electron transfer flavoprotein-ubiquinone oxidoreductase, ygcP, and fixB, among others, in E.coli_45. We also found IS110 and IS3 transposases in E.coli_30, with passenger genes including rhm1, algC, and btuD. We also found 20 inserted genes in one larger insertion site in E.coli_30.

Furthermore, the genomes were then annotated using MobileOG-db, which tracks mobile elements by organizing them into orthologous groups. We tracked a wide range of MGEs, including transposons, integrons, IS elements, phage related genes, and integrative conjugative elements, which explains the large number of MGEs found within the two genomes (Figure 8A). Among the identified MGEs, we found genes involved in critical processes such as replication, recombination, repair, integration/excision, and phage-related functions (Figures 8C,D). Notable genes include sdiA, associated with conjugation and relaxosome activity, and ruvA, ruvC, and holE_2, which play key roles in DNA repair, recombination, genome stability, and horizontal gene transfer (HGT). Additionally, phage-related genes such as dicB, essD, and quuQ_2 that contribute to lysis, lysogeny, and infection regulation, were detected. Transposon-associated genes such as tnpA and tnpB facilitate integration and excision, while stability/transfer/defense genes like hipA and cnu help maintain genetic integrity. We observed that some entries for both genomes were manually curated, while others were annotated through automated homology-based methods using several databases.

The phage genome annotator PhiGARo predicted prophages and other virus-like sequences in conjunction with virus orthologous groups (VOGs). We identified five prophage regions and 138 prophage genes in E.coli_30. In the second genome, we predicted 11 prophage regions and 259 genes, significantly more than those in E.coli_30. The taxonomy of these prophage genes revealed that most belonged to the Myoviridae family, a few were of unknown origin, and one was classified as Podoviridae. The prophage genes identified were associated with roles such as lysis, integration, termination, tail assembly, replication, coat formation, and portal functions. HGT regions and their insertion sites within the two genomes were predicted using Alien Hunter, while viral signals were detected using VirSorter within the two genomic sequences (Figure 8A).

Finally, we further investigated the locus harboring the ESBL genes in E.coli_45 and found that the bla_TEM-1 gene was likely inserted into the genome via the Tn3 insertion sequence. Annotation of this locus revealed a Tn3 family transposase (Tn2) and a site-specific DNA recombinase (tnpR) located just a few bases upstream of the ESBL gene (Figure 8B), which likely facilitated its integration. The closest match we found using ISfinder BLAST for this locus showed the presence of a DDE recombinase (IS110 family transposase), Tn3 resolvase, and a passenger gene (bla_TEM-1), and this result aligns closely with our findings.

The findings from the antibiotic resistance profiles among the E. coli isolates provided vital insights into the response of each antibiotic against the tested bacterial isolates. In this study, GEN exhibited the most efficacy against the E. coli isolates tested, while IMP and CIP, on the other hand, exhibited high levels of resistance. Shams et al. (2018) documented high levels of carbapenem-resistant E. coli, while Shariati et al. (2022) reported CIP resistance in the context of urinary tract infections and other E. coli-associated diseases.

Interestingly, the susceptibility test showed that our isolates displayed much greater levels of resistance to IMP compared to AMP, a first-generation penicillin. This situation could have occurred due to several reasons. One possibility is that when certain antibiotics begin to show reduced efficacy against bacteria, their increased consumption may result in heightened resistance. This explains why IMP-resistant bacteria are often associated with carbapenem intake, possibly corroborating our findings (Yang et al., 2018).

Another plausible reason is that the IMP-resistance isolates in our study may habor carbapenem-specific β-lactamases such as New Delhi Metallo-β-lactamase (NDM) or Klebsiella pneumoniae Carbapenemase (KPC) without ESBLs or AmpC enzymes, thus remaining susceptibility to AMP. For instance, Mojica et al. (2022) reported that some E. coli strains harboring NDM remain susceptible to some penicillins due to the specific substrate profile of NDM enzymes.

Additionally, IMP resistance could be due to alterations (mutations) in porin channels that limit the entry of IMP without affecting AMP. This phenomenon has been reported in clinical isolates where specific mutations in porin proteins confer carbapenem resistance (Poirel et al., 2020). Furthermore, efflux pumps, such as those of the RND family, which confer MDR, may have substrate specificity that includes IMP but not AMP, potentially contributing to the resistance patterns observed (Tambat et al., 2022).

The identification and characterization of AMR genes in our E. coli isolates revealed several key resistance mechanisms. The presence of pmr phosphoethanolamine transferases (e.g., marA and acrD), efflux pumps (e.g., acrAB-tolC), and β-lactamase indicate that these isolates possess vast mechanisms to evade antibiotic treatment. These findings agree with other studies that have reported the widespread presence of AMR genes in some E. coli strains, posing major problems to public health (Yang et al., 2018; Majumder et al., 2021). Efflux pumps play a crucial role in mediating bacterial resistance by expelling antibiotics from the cells. To combat this, Lu et al. (2024) developed a potential efflux pump inhibitor (EPI) derived from a seaweed compound, which displayed efficacy in reducing AMR in drug-resistant E. coli strains. Another study explored the use of a microbe-derived EPI, ethyl 4-bromopyrrole-2-carboxylate (RP1), which shows promise in combination therapies (Tambat et al., 2022). The findings from that study indicate significant reductions in minimum inhibitory concentrations and extended post-antibiotic effect, making RPI a promising candidate for enhancing existing antibiotic treatments.

We noted the presence of specific β-lactamase genes (bla_TEM-1, bla_CTX-M-55) in E.coli_45 (Figure 8B), making this isolate an ESBL-producing E. coli (Chen et al., 2019; Papouskova et al., 2020). The presence of ESBL enzyme in E.coli_45 indicates that this strain is resistant to most β-lactam antibiotics, posing a significant public health risk due to limited treatment options (Chen et al., 2019). A recent pharmaceutical strategy to address this issue involves the use of β-lactamase inhibitors. These inhibitors are crucial in restoring the activity of β-lactam antibiotics against bacteria that produce β-lactamase enzymes. Le Terrier et al. (2024) studied taniborbactam, a novel broad-spectrum β-lactamase inhibitor effective against metallo-β-lactamases (MBLs). Their study showed that taniborbactam exhibited significant activity against B1 subclass MBLs, including NDM and VIM (Verona integrin-encoded metallo-β-lactamase) enzymes (Le Terrier et al., 2024). However, certain variants such as NDM-9 and VIM-83 demonstrated resistance, raising concerns regarding the potential emergence of resistant strains (Le Terrier et al., 2024). Taniborbactam’s ability to inhibit a wide range of MBLs, including those resistant to traditional β-lactamase inhibitors, positions it as a critical tool in the fight against MDR pathogens.

In addition to β-lactam resistance, several other AMR genes were detected in both isolates. This underscores the critical need for monitoring and developing novel therapeutic strategies, such as bacteriophage therapy, to combat MDR (Montso et al., 2019). We also identified multiple AMR genes associated with resistance to various classes of antibiotics, including quinolones (qnr genes), sulfonamides (sul genes), and tetracyclines (tet genes). The presence of numerous variants of these genes within the MDR isolates suggests HGT as a major facilitator of AMR (Partridge et al., 2018). This finding aligns with research that highlights the role of plasmids, integrons, and transposable elements, such as IS, in spreading AMR across E. coli and other bacterial species (Zingali et al., 2020; Vinayamohan et al., 2022).

The analysis of virulence-associated genes in the test isolates revealed the presence of genes associated with adhesion (e.g., cfaB and fimH), invasion (e.g., ibeC and ompA), effector delivery systems (TTSS effectors), and metabolic factors or iron acquisition (e.g., siderophores like enterobactin and salmochelin) (Algammal et al., 2020). Notably, salmochelin siderophore virulence genes, such as iroB and iroC, were exclusively found in E.coli_30 and its plasmid. These virulence identifiers (adhesion, invasion, and nutrition) are vital for the pathogenic potential of E. coli, enabling the bacterium to colonize hosts, evade immune responses, and cause disease (Croxen et al., 2013).

Our findings are supported by previous research highlighting the wide range of virulence-associated genes found in pathogenic E. coli isolates, which significantly contribute to the bacteria’s ability to cause different types of infections (Sarowska et al., 2019; Sora et al., 2021). The identification of virulence genes, such as the TTSS effectors (e.g., espX1, espX4, and espL1) and iron uptake systems like salmochelin, in our isolates suggests their potential to cause severe infections (Radics et al., 2014; Magistro et al., 2015). These identified factors, including invasin, adhesin, and metabolic factors, facilitate bacterial invasion, intracellular survival, and nutrient acquisition, respectively, which are essential for causing infections in a host organism (Liu et al., 2022).

The presence of virulence genes in both commensal and pathogenic bacterial strains has been reported and these genes underscore the complex relationship between E. coli and its host, determining the outcome of infections (Dionisio et al., 2023).

The phylogenetic tree generated in this study revealed distinct clades, each representing different evolutionary lineages. The isolates we sequenced, E.coli_30 (ST3579) and E.coli_45 (ST1121), were placed in different clades, highlighting a major genetic divergence and suggesting that these two isolates have been through different evolutionary pressures or they may have come from different ancestral populations. What we observed here agreed with previous studies that reported on the genetic diversity within E. coli strains driven by various evolutionary factors such as mutations, HGT via MGEs, and genetic recombination (Vinayamohan et al., 2022; Touchon et al., 2009). The diverse placement of isolates with different STs in the phylogenetic tree indicate a wide genetic diversity even within smaller groups of E. coli strains. This aligned with the result of a study by Didelot et al. (2012) in which they reported on the diversity of E. coli with dynamic adaptation to many hosts and environmental conditions, leading to major genome plasticity. The 45 isolates from the NCBI database were distributed across various clades in the phylogeny, except for isolates with the same STs (mostly tend to clade together). This highlights the extensive genetic diversity present in the NCBI database for E. coli strains. The varied clustering patterns we observed in the phylogenetic tree suggest that these E. coli isolates were arranged into multiple sub-types with each strain following its unique evolutionary paths (Touchon et al., 2020).

The phylogenetic analyses of bla_TEM-1 and bla_CTX-M genes in E. coli, Salmonella enterica, and Klebsiella pneumoniae from China showed how widely these resistance genes have spread, thanks to MGEs. The discovery of multiple copies of these genes within individual genomes and across species highlights the important role of MGEs, such as IS and transposons, in moving these genes between bacteria. bla_TEM-1 appeared to be highly conserved, meaning it has changed little across species, while bla_CTX-M variants, such as CTX-M-55 and CTX-M-15, showed more diversity, forming distinct clusters that reflect different evolutionary paths. This mirrors the findings from other studies that show plasmids helping to spread resistance genes across bacterial species in China (Wu et al., 2015; Kuang et al., 2018). The presence of small mutations, or SNPs, in some of the genes suggests that they are continuing to evolve. Overall, these results emphasize how MGEs drive the spread of antibiotic resistance across bacteria, agreeing with global trends observed in similar studies (Peng et al., 2024; Jin et al., 2015; Awosile and Agbaje, 2021).

The pangenome analysis conducted in this study included both core and accessory genomes to elaborate extensively on the genetic diversity within 50 E. coli strains. The core genomes has conserved genes that are present in all the isolates, shedding knowledge on essential functions and evolutionary conserved mechanisms. The accessory genome on the other hand, consisting genes present in some but not all isolatesprovided insights into the adaptive potential of E. coli to many different environments and antibiotic treatments (Vernikos et al., 2015). The results of the pangenome analysis revealed a high proportion of unique genes in individual isolates, indicating strain-specific E. coli adaptations. Recent studies on E. coli pangenome (gene-specific) support what we reported that new genes continue to be obtained through HGT, helping in the bacteria’s adaptation and vast ecology (Hall et al., 2021). The identification of genes associated with metabolism (iron uptake), stress response, and adhesion, highlights the adaptive nature of the E. coli genome, facilitating its survival in diverse environments (Liu et al., 2022).

Our investigation into MGEs across the two genomes revealed that the results varied significantly depending on the analytical tools used. Each tool offered distinct methodologies for MGE identification and these helped illuminate different aspects of MGE diversity and behavior. MGEfinder emerged as the most reliable tool in this study, providing detailed mapping of insertion sites and accurately identifying associated insertion sequences. This tool successfully recognized IS elements, whether or not they contained inverted terminal repeats, which are essential features for many MGEs. In particular, sequences with terminal repeats, such as those from the IS3 family transposase (IS3) and ISEc17 were detected in the two genomes and across several other genomes when studying the spread of ESBL genes. These elements are widespread in E. coli genomes and often carry critical passenger genes such as bla_TEM-1 and bla_CTX-M that contribute to antibiotic resistance. This finding aligned with prior studies that emphasized the role of IS elements in disseminating resistance genes (Partridge et al., 2018). Additionally, MGEfinder was able to identify elements that lacked inverted terminal repeats, which are predominantly part of the IS110 family (Durrant et al., 2024). In both the two genomes, we noted the presence of IS621, an essential member of the IS110 recombinase. The integration mechanism for this family involves a recently discovered bridge RNA recombination process (Durrant et al., 2024; Hiraizumi et al., 2024), and this underscores a novel pathway for MGE mobility.

Complementary to MGEfinder, other tools such as MobileOG-db provided MGE characterization by annotating the genomes against a variety of databases. MobileOG-db matched sequences with several key databases, including pVOG, Plasmid RefSeq, GPD, and ICE. This broadened the scope of our analysis and yielded a large number of hits for MGEs in our isolate. In addition, we employed VirSorter, PhiGARo, and Alien Hunter, which further contributed to MGE detection across the genomes by targeting specific MGE types. VirSorter has been widely used in recent virome studies (Roux et al., 2015) to identify viral sequences embedded in bacterial genomes, a feature relevant to our study given the potential role of viruses in HGT. Similarly, PhiGARo and Alien Hunter have previously demonstrated effectiveness in detecting prophages and alien DNA, respectively, complementing each other by expanding the types of MGEs detected (Akhter et al., 2012).