Longitudinal genomic surveillance of multidrug-resistant Escherichia coli carriage in critical care patients

A prospective observational cohort study was conducted over a 6-month period from June 2021 to December 2021 at Makassed General Hospital in Lebanon. This is a 200-bed hospital located in a heavily populated area of Beirut with a medically underserved population. It serves nearly 15,000 in-patients per year that are mostly from middle- to low socioeconomic status groups.

A standard protocol for infection and prevention control to screen for any patient admitted to the intensive care unit and to the hematology oncology unit has previously been established at Makassed General Hospital. As described below, a perianal swab was collected from each patient and screened for extended-spectrum beta-lactamases and carbapenem-resistant Enterobacterales (ESBL and CRE, respectively). In total, 144 patients were admitted to the intensive care and oncology units; they were recruited for this study and followed during their hospitalization. Only patients that were hospitalized for more than 3 days were included. In total, 97 patients were dropped because of their early discharge, death, or transfer to another hospital ward, leaving 47 remaining patients. Perianal swabs were collected by healthcare workers at admission and during hospitalization. A minimum of 5 days was required between the first, second, or third collection time. A standardized questionnaire for each patient was filled by the resident consultee. The questionnaire included information about age, gender, cause of hospitalization, duration of hospitalization, antibiotics exposure (past and during hospitalization), types of antibiotics used, and bacterial infection during hospitalization.

Anal swabs were cultured on three different culture media plates: MacConkey agar without antibiotics to confirm the absence of susceptible Gram-negative strains, MacConkey agar with ertapenem (0.5  mg/L) to isolate carbapenamase-producing bacteria, and MacConkey agar with ceftriaxone (4  µg/mL) for isolation of ESBL-producing bacteria. Resistant strains were collected from the selective plates that includes antibiotics for further analysis. All plates were incubated at 37°C for 24 h. Following bacterial growth, two colonies of the same color, morphology, and shape were picked and pooled and then taxonomically classified using matrix-assisted laser desorption/ionization time-of-flight mass spectrometry (MALDI-TOF). Antimicrobial susceptibility using the disk diffusion method was determined for all E. coli isolates for the following antibiotics: amoxicillin (20 µg), amoxicillin/clavulanic acid (20/10 µg), cefepime (30 µg), ceftriaxone (30 µg), piperacillin-tazobactam (100/10 µg), ceftazidime-avibactam (20/10 µg), cephalothin (30 µg), ertapenem (10 µg), fosfomycin (200 µg), trimethoprim/sulfamethoxazole (1.25/23.75 µg), ciprofloxacin (5 µg), colistin (10 µg), imipenem (10 µg), nitrofurantoin (300 µg), amikacin (30 µg), tetracyclin (30 µg), and gentamicin (10 µg).

Genomic DNA was extracted from the E. coli isolates using the QIAamp DNA Mini Kit (Qiagen) following the manufacturer’s guidelines. DNA concentrations were measured using Qubit (Thermo Fisher Scientific) with the Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). Sequencing libraries were made with Nextera XT DNA Library Preparation Kit (Illumina) according to the manufacturer’s instructions. Illumina sequencing was performed on a Novaseq-6000 producing paired-end 2 × 150-bp reads, with an average base quality (Phred score) of 36 and G + C content of 54.1%.

Reads were processed with Trimmomatic v0.39 (43) to remove adaptor sequences and quality filter the reads; we specified the TruSeq3 adapters and the quality filtering parameters of LEADING:3, TRAILING:3, SLIDINGWINDOW:4:15, and MINLEN:30. The assembly was performed using SPAdes Genome Assembler v3.15.4 with default parameters (44). Assemblies were annotated with Prokka v1.14.5 (45) with the kingdom specified as bacteria. The resulting annotated genes were processed by Panaroo v1.3.0 to create the core genome alignment and pangenome (46). This was run with the “strict” clean mode parameter and with the “clustal” option specified for the core genome alignment. Plots were generated using Roary v3.12.0 including the pangenome frequency plot, a presence and absence matrix against a tree and a pie chart of the pangenome, breaking down the core, soft core, shell, and cloud (47). The phylogenetic tree along with the pangenome of the isolates was visualized with Phandango v1.1.0 (48). The genomes of this study were aligned with 32 previously published E. coli genomes identified from clinical settings to create first an extended core genome alignment from Panaroo and second to generate a maximum-likelihood phylogeny from the combined alignment using RAxML (46, 49). The resulting tree was visualized with the Interactive Tree of Life (iTOL v6) (50). SRST2 (v0.2.0) was used to identify multi-locus sequence types (MLST) and serotypes (51). Phylogroups were classified based on the ClermonTyping method (52). FastANI v1.33 was used to calculate the average nucleotide identity (ANI) of orthologous genomic regions between pairs of genome assemblies (53). The presence of antimicrobial resistance genes, putative virulence factors, and plasmid replicons were studied using ABRicate with the ResFinder database, Virulence Factor database, and PlasmidFinder database (54). In silico plasmid replicon typing was performed using PlasmidFinder v2.1 (https://cge.food.dtu.dk/services/PlasmidFinder/) and IncF replicon sequence types were characterized through the pMLST v2.0 tool (https://cge.food.dtu.dk/services/pMLST/) (55). We used Snippy v4.6.0 to compare E. coli genomes from the same patient and identify single nucleotide variants (SNVs). The reference genome used was the E. coli strain identified at T0 (at hospital admission) or at T1 (the first collection time during hospitalization) (56). The major allele frequency distribution was also based on the identified variants, including SNVs and insertions/deletions.

In total, 20 (43%) patients were eligible for this study, all of whom were colonized by antibiotic-resistant isolates. Throughout their hospitalization, 14 (70%) patients were colonized with E. coli, one with Raoultella ornithinolytica and one with Klebsiella pneumoniae. Four patients were colonized by multiple bacterial species during their hospitalization. Here, we focus on the E. coli-colonized patients. The median age of patients infected with E. coli was 67 years (range = 17–88 years), eight were female (57%) and six (43%) were male. These 14 patients were admitted to different hospital wards: the intensive care unit (n = 7), the hematology and oncology ward (n = 3), and the isolation ward (n = 4). The median duration of hospitalization for these patients was 29 days (range = 9–73 days). The cause of admission was either pneumonia, chemotherapy, bleeding, or sepsis (Table S1). In total, 11 out of 14 patients had an indwelling catheter during hospitalization and 13 patients were under corticosteroid therapy. The patients were exposed mostly to broad-spectrum antibiotics during their hospitalization, including the following: meropenem which was used in 50% of patients; colistin and piperacillin/tazobactam in 50%; amikacin, ceftriaxone, and ceftazidime/avibactam in 50%; and levofloxacin in 29% of patients. Vancomycin was used in 43% of patients; trimethoprim-sulfamethoxazole, teicoplanin, gentamycin, and clarithromycin were used in less than 20% of patients. The median time between first sample collection (T0) and second collection (T1) was 9 days (range = 5–28 days). During hospitalization, 7/14 (50%) of patients were documented to have infectious episodes (Table S1), including respiratory tract infections (n = 5), followed by bacteremia (n = 5), urinary tract (n = 3), and skin and soft tissues (n = 3).

In total, 12 out of 14 patients infected with E. coli were selected for further analysis, of which 8 patients had isolates from multiple time points, yielding 22 whole-genome sequences in total. Two patients (MGH 3 and MGH 12) and two isolates (MGH 2A and MGH 6B) were dropped since their assembly sizes exceeded 8 Mbp, much larger than expected for an E. coli genome. Two of these assemblies also had median SNV major allele frequencies in the 0.93–0.97 range, whereas all others were very close to 1 (Fig. S1). This suggests that some of the unexpectedly large assemblies could have been due to contamination or due to pooling colonies of two distinct E. coli strains. The remaining 22 genomes considered for further analysis can be considered single strains.

These 22 E. coli genomes were genetically diverse, with 11 STs (ST10, ST44, ST69, ST131, ST224, ST227, ST540, ST648, ST1286, ST1431, and ST1491) identified (Fig. 1). The most common was ST131 which was carried by 5/12 patients (41%).

We inferred the maximum-likelihood phylogeny of these isolates using a core-gene alignment. We included 32 E. coli genomes previously sequenced from clinical samples collected in Lebanon (31, 32). Based on this phylogeny, we identified five phylogroups. Five isolates were phylogroup C, eight were phylogroup A, and phylogroups F and D were each represented by one isolate. The ST131 strains (within phylogroup B2) were the most frequent (n = 7). In total, 11 serotypes were identified with the most common being O25:H4 and O9:H30 (Fig. 2).

Of the six patients sampled at two time points 5–11 days apart, five patients (MGH 1, MGH 8, MGH 9, MGH 10, and MGH 11) had pairs of isolates that clustered very closely on the phylogeny (Fig. 2) and had identical profiles of AMR genes (Fig. 1). This is consistent with E. coli lineages persisting in patients between sampling points, although we cannot exclude transmission events of near-identical genomes among patients, for example, the two isolates from patient MGH 8 were very closely related to isolates from patients MGH 4, 6, and 14, suggesting possibly recent transmission of ST131 strains. In contrast, the two isolates from patient MGH 13 sampled 7 days apart came from entirely different phylogroups (C and D; Fig. 2) and had distinct AMR gene and plasmid profiles (Fig. 1). This suggests either a mixed infection in patient MGH 13 (with phylogroup C sampled first and D sampled second) or a strain replacement that occurred between time points.

Among the 12 patients with sequenced isolates, five harbored the ST131 type. Of these, four strains belonged to the O25:H4 serotype and fimH30 variant and one strain to the O16:H5 serotype and fimH41 variant. ST131 formed a distinct clade on the phylogenetic tree, classified as phylogroup B2 (Fig. 2). All ST131 genomes encode bla_CTX-M-27, bla_EC-5, and the IncF plasmid group F1:A2:B20 plasmid; 80% of them also encode bla_OXA-181. They included type 1 fimbriae, P fimbriae, and enterotoxin TieB. All ST131 genomes were highly genetically similar, with average nucleotide identity values ranging from 99.5% to 99.9% (Fig. S2). A pangenome analysis revealed 3,308 core genes present in all isolates, 85 “soft core” genes, 2,520 “shell genes,” and 2,653 “cloud” genes. The gene presence/absence patterns are shown alongside the phylogeny and strain typing information (Fig. 3). Notably, ST131 strains with the O25:H4 serotype have a unique pattern of gene presence/absence, which is distinct from other ST131 strains and from other phylogroups.

Among the 14 patients carrying E. coli, 13 carried MDR strains. Only one patient (MGH 11) carried a strain that was resistant to cephalosporins, yet it was susceptible to all other antibiotics (Table S2).

Susceptibility testing was performed on 29 E. coli isolates collected from patients at different time points. Identical strains (defined as ANI of 99.9% or more) isolated from the same patient at different time points were omitted from the analysis (n = 7). Among the remaining 22 E. coli strains, 55% (n = 12/22) were ESBL producing, 5% (n = 1/22) CRE, and 41% (n = 9/22) ESBL/CRE. All patients enrolled in this study carried resistant strains upon admission. A total of 95% (n = 21/22) were resistant or intermediate to cefotaxime, 86% to ceftriaxone (n = 19/22), 54% (n = 12/22) to ceftazidime, 63% (n = 14/22) to ciprofloxacin, 50% (n = 11/22) to ertapenem, 32% (n = 7/22) to imipenem and meropenem, and 10% (n = 2/22) to ceftazidime/avibactam (Table S2).

The genomes of these isolates were screened for known genetic determinants of AMR using the ABRicate pipeline. More than 46 different acquired AMR genes were detected; β-lactamases (n = 19 genes) and aminoglycosidases (n = 9 genes) were among the most common, and patients contained 2–16 resistance genes, with a median of 11 per patient (Fig. 1). Among the 12 patients with non-identical E. coli strains, Class A, C, and D beta-lactamases were mostly identified, including bla_CTX-M-15 (n = 5/12, 41%), bla_CTX-M-27 (n = 5/12, 41%), bla_EC-5 (n = 5/12, 41%), bla_EC-15 (n = 4/12, 33%), and bla_OXA-181(n = 4/12, 33%). The beta-lactamases bla_NDM-5 and bla_OXA-244 were seen in one and two patients, respectively. Several other AMR genes were also detected (Fig. 2), with qnrS1 (encoding quinolone resistance) being quite prevalent (n = 9/12, 75%). There was a correlation between the observed antibiotic resistance phenotypic profile and the specific genes identified through sequencing.

During hospitalization, beta-lactam antibiotics were used to treat critically ill patients, all of whom were colonized with cephalosporinase-producing E. coli. Notably, in addition to this treatment regimen, six patients were infected with E. coli encoding carbapenemase genes and were administered meropenem. This approach raises critical questions about the selection for resistance genes and the effectiveness of treatment in individuals colonized with these resistant isolates. Evaluating whether patients with such colonization could still respond positively to antibiotic therapy is essential in our ongoing efforts to address antibiotic resistance and improve patient outcomes.

The ABRicate pipeline was also used to identify the incompatibility groups for the recovered plasmid replicons (Fig. S3). IncFI was the most common plasmid (73%; n = 16/22), followed by IncFII (63%; n = 14/22), IncX (41%, n = 9/22), IncI (23%; n = 5/22), IncB (27%; n = 6/22), and IncH (14%; n = 3/22).

We screened the isolates for virulence genes which are often encoded by E. coli and detected a wide variety (Fig. S4). Some virulence genes were commonly detected, including ompA (100%) which is an outer membrane protein required for conjugation. Iron regulatory proteins (IRP1 and IRP2) were detected in 86% (n = 19/22) of isolates.

We recovered some of the known markers of intestinal and extraintestinal virulence in some isolates, such as the invasion and evasion factors (kpsM and kpsD [1%, n = 9/22]) and adherence factors (fim [86%, n = 19/22]) and pap genes (41%, n = 9/22) that are associated with intestinal and extraintestinal infections in humans. One patient (MGH 9) harbored an enteroaggregative E. coli strain with an anti-aggregation gene and dispersin-encoding gene (aap) and another patient (MGH 5) an enterotoxigenic E. coli with heat-labile (eltA and eltB) enterotoxin genes. Other virulence genes detected in the study included chu, which codes for an outer membrane hemin receptor (41%, n = 9/22), and the enterotoxin gene, senB (32%, n = 7/22).

A comparative genomic analysis was performed to compare the strains isolated from patients at different time points. Among the eight patients who were followed over time, five (MGH 1, MGH 8, MGH 9, MGH 10, and MGH 11) had very closely related E. coli isolates at T0 and T1 (Fig. 1 and 2). The duration of hospitalization for these patients ranged between 5 and 11 days, and genomes differed by an average of 16 SNVs (range 10–22 SNVs). In contrast, three patients (MGH 4, MGH 7, and MGH 13) were colonized by a highly divergent strain (different by >1000 SNVs or <97% ANI; Fig. S2) between time points, leading to marked changes in AMR gene profiles (Fig. 1). In two patients, MGH 4 and MGH 7, there was a putative strain replacement between T0 and T1, with the same strain persisting between T1 and T2, differing by a few SNVs (Fig. 4). Comparing only closely related pairs of genomes from the same patient yielded no evidence for a molecular clock (Spearman’s correlation between number of SNVs and number of days separating isolates, rho = 0.24, P = 0.60). The lack of clock signal could be due to a small sample size, short duration of sampling, or sampling only one isolate per patient per time point.

Of the 120 genetic variants (SNVs or complex variants of multiple mismatch nucleotides or insertions/deletions) observed between closely related isolates from the same patient, half (n = 60/120) were located in non-coding regions and 29% (n = 35/120) were synonymous changes. Most occurred in hypothetical proteins (n = 18/35, 51%) (Table S3). Missense variants or insertions were observed in 21% of the total variants detected (n = 25/120). These were detected in hypothetical proteins (n = 12/25, 48%), transposases (n = 5/25, 20%), a putative protein YjdJ (n = 1/25, 4%), antigen 43 (n = 2/25, 8%), cytoskeleton bundling-enhancing antitoxin (n = 1/25, 4%), ompF outer membrane porin F (n = 1/25, 4%), and D-alanine-D-alanine ligase A (n = 1/25, 4%) (Table S3). We next tested the hypothesis that the frequency of known AMR genes or virulence factors (VFs) might increase within patients over time, due to exposure to antibiotics and other selective pressures in the critical care hospital environment. Although certain patients (e.g., MGH 13) acquired several AMR genes and VFs over time, there was no consistent trend across patients (Table S3; Fig. 5). Changes were most pronounced in patients MGH 4, MGH 7, and MGH 13, who experienced a putative strain replacement between time points (Fig. 4 and 5).