The isolate L1PEag1 was isolated from fruit peel (exocarp) of cape gooseberry, purchased from the local market of Ibarra city, northern Ecuador, and following a standard procedure (Tenea et al., 2023). The purified isolate was grown overnight on Luria Bertani (LB) agar (Oxoid, UK) at 37°C. Genomic DNA was extracted using a high molecular weight DNA extraction kit (Wizard Genomic DNA purification kit), and following the manufacturer indications (Promega, USA). The DNA was quantified, using a NanoDrop™ 2000 Spectrophotometer (Thermo Fisher, USA). The library was prepared using a Ligation Sequencing kit V14 (SQK-LSK114) (Oxford Nanopore Technologies Ltd., UK). High molecular weight DNA was nick–repaired and end–prepared with the NEBNext® Companion Module for Oxford Nanopore Technologies® Ligation Sequencing (E7180S). The end–preparation protocol phosphorylated 5’ends and dA-tailed 3′ ends of gDNA fragments. A magnetic bead clean-up of enzymes was performed before ligation. Sequencing adapter ligation was performed with the NEBNext® Quick Ligation Module (E6056). After ligation, a second magnetic bead clean–up was performed before sequencing. The library was sequenced on a MinION Mk1C sequencer and a Flow Cell version R10.4.1 (FLO-MIN114, Oxford Nanopore Technologies Ltd., UK). Sequencing run took 48 h and generated 998.45 k reads at an estimated coverage of 400x.

Raw sequence data was high-accuracy basecalled into fastq data (HAC) under the dna_r10.4.1_e8.2_400bps_hac@v4.2.0 model with the Dorado 0.2.1 package (Hariharan et al., 2020). Quality control was performed with NanoPlot version 1.41.3 (Coster et al., 2018) and FastQC version 0.12.1¹. To reduce errors in assembly, we applied filtering and trimming with NanoFilt 2.8.0 (Coster et al., 2018). Reads with less than a quality value of Q14 were filtered out for the assembly. To allow for plasmid assembly, short reads were not filtered out. Trimming was performed on the first 20 bps for each end of the reads to eliminate adapters. A subsequent round of quality control was performed after the filtering and trimming. De novo assembly was made according to the long–read assembly pipeline and the Trycycler tool by Wick et al. (2021), these consist in estimating a final consensus assembly from multiple independent assemblies as input. The reduction of error during the final assembly process, assisted by the Trycycler pipeline and its algorithms, relies on the fact that multiple independent subset assemblies provide alternative sources of evidence that compensate previous assembly biases. In brief: (1) By randomly subsampling the original read set, sequence data was subsampled into 12 independent read sets. To maximize the independence of assembly results, each subset was independently assembled with different assembling tools. The assemblers used for this step were (a) the repeat graph assembler FLYE 2.9.2-b1786 (Kolmogorov et al., 2019), (b) Minipolish (Wick and Holt, 2020), and (c) Raven (Vaser and Šikić, 2021). (2) The different assemblies were clustered based on their k-mer content. For the case of the L1PEag1 isolate, clusters recovered two replicons that consisted of the main circular genome and one accompanying plasmid. (3) The clustered sequences (contigs) were then reconciled by normalizing their directionality or orientation and aligned to each other to repair circularization issues such as missing bases or overlapping bases. (4) Each cluster of sequences was rotated to a common starting sequence, preparing them for multiple sequence alignment. (5) A global multiple sequence alignment was produced for each cluster of sequences. (6) A definitive consensus sequence for each cluster was obtained based on the minimum total Hamming distance to optional or variable regions in the consensus and the best total alignment score. (7) Finally, the consensus sequence was polished according to the Medaka v.1.1.1.3 algorithms².

The genome FASTA sequence was uploaded to the Type (Strain) Genome Server (TYGS) to conduct a complementary genome-based taxonomic study and a phylogenetic relationship inference (Meier-Kolthoff and Göker, 2019). The closest strain genome type was determined by comparing the L1PEag1 genome to all accessible strain genomes in the TYGS database using the MASH algorithm (Meier-Kolthoff and Göker, 2019). Strains with the shortest MASH distance were selected automatically. Thus, the precise distance was calculated using the Genome BLAST distance phylogeny (GBDP) approach under the “coverage” algorithm and distance formula (Meier-Kolthoff et al., 2013). Phylogroups were determined using the in silico Clermont Phylotyper EzClermont³ (Waters et al., 2020). In addition, multi-locus sequence typing (MLST)⁴ (Larsen et al., 2012) was used for accurate subtyping, with the E. coli # 1 (Wirth et al., 2006), and E. coli # 2 (Jaureguy et al., 2008) schemes. Further, SerotypeFinder 2.0⁵ was used to identify the serotype (Joensen et al., 2015), using 95% sequence identity and 60% sequence coverage setup. To enrich the phylogenetic contrast, an additional genome search for similar E. coli sequences from Ecuador was made at NCBI’s databases and guided by previous published work (Montero et al., 2021; Rothstein et al., 2023). While genome sequences from earlier research have been made available at NCBI, these are either assembly-level contigs or raw sequences (SRA experiments in NCBI terminology), which do not meet the requirements for phylogenetic comparisons of the type carried out in this work on the whole assembled genome.

The L1PEag1 genome was annotated using the PROKKA tool (Seemann, 2014). Plasmid annotation was performed with pLannotate (McGuffie and Barrick, 2021). The predictions of CDS, rRNA, tRNA/tmRNA, signal leader peptide, and noncoding RNA were performed using Prodigal (Hyatt et al., 2010), RNAmmer⁶, Aragorn⁷, Signal IP⁸ and Infernal (Nawrocki and Eddy, 2013). Additionally, the predicted genes were annotated through the Global Catalog of Microorganisms genome annotation project pipeline⁹ and based ondatabases SwissProt¹⁰, MetaCyc¹¹ (this is a database that contains pathways responsible for both primary and secondary metabolism) (Caspi et al., 2018), VFDB (a virulence factor database¹²) (Liu B. et al., 2019), PHI: The Pathogen-Host Interaction database¹³, KEGG (Kyoto Encyclopedia of Genes and Genomes), Orthology Database of Prokaryotes, COG (Clusters of Orthologous Groups of proteins), and pFam¹⁴ (Wu and Ma, 2019). A circular map was generated by importing the FASTA sequence into the PROKSEE server¹⁵ (Grant et al., 2023).

CRISPRFinder¹⁶ and PHAge Search Tool Enhanced Release (PHASTER)¹⁷ (Arndt et al., 2016) were used to detect CRISPR, Cas sequences, truncated Cas sequences, and prophage sequences. Only predicted prophages with intact completeness in PHASTER, as defined by Arndt et al. (2017), were considered for report. In addition, mobileOG-DB was used to annotate the mobile elements (Brown et al., 2022).

The Comprehensive Antibiotic Resistance Database (CARD)¹⁸ (Jia et al., 2017) and the Resistance Gene Identifier tool (RGI) were used to detect the antibiotic resistance genes by importing the contig files in FASTA format to the database (Zankari et al., 2012). The ResFinder 4.3.3 server¹⁹ was used to identify acquired antimicrobial resistance genes (Bortolaia et al., 2020). The putative virulence factors were predicted using the VirulenceFinder 2.0 web server²⁰ (Cosentino et al., 2013). The bacterial pathogenicity was predicted using the PathogenFinder web server²¹ (Cosentino et al., 2013). Whereas FimTyper 1.0²² was used to type de fimH alleles (Roer et al., 2017). The detection standard parameters were set at 90% sequence identity and 60% sequence coverage for VirulenceFinder and ResFinder.

GIs were predicted with the webserver IslandViewer 4, using E. coli SMS-3-5 as a reference strain (Bertelli et al., 2017). To search for ISs, the ISfinder tool (ISsaga V.2.0) was used (Siguier et al., 2006).

An estimated genome size of 4825.72 Kbp (400X coverage) and a plasmid of 3.561 Kbp (400X coverage) were generated by the sequence analysis. Species matching resulted in a 100% hit on the entire contig for Escherichia coli. The results of the whole genome analysis placed the L1PEag1 strain on the same lineage as E. coli DSM30083 (Figure 1). Based on the BLASTN analysis, the plasmid pL1PEag1 showed 99.96% sequence identity with plasmid pECQ4552_IHU08 from Klebsiella quasipneumoniae subsp. similipneumoniae, strain IF3SW-P1 (CP092122.1); 99.89% identity with an unnamed plasmid from Staphylococcus aureus strain, Alexandria 2020–19 (CP113245.1); and 99.81% sequence identity with plasmid pECQ4552_IHU08 from E. coli strain Q4552 (CP077071.1).

Based on the Clermont Escherichia-typing method (Waters et al., 2020), the L1PEag1 isolate belongs to phylogroup B2 (TspE4 +, chu +, yjaA +, and arpA –). Generally, B2 phylogroup strains are associated with Extraintestinal pathogenic E. coli (ExPEC) (Xia et al., 2022). A comparison of several E. coli strains prevalent in inflammatory bowel disease indicated that the B2 phylogroup strains share a specific metabolism and carry virulence genes enabling them to colonize intestinal mucosa (Fang et al., 2018). In E. coli, MLST requires assessment of 7 housekeeping genes (adk, fumC, gyrB, icd, mdh, purA and, recA) to assign sequence types (ST) genes (Wirth et al., 2006); thus, by using the E. coli #1 scheme, L1PEag1 was identified as the ST1170 sequence type. By using the E. coli #2 scheme (by the assessment of genes dinB, icdA, pabB, polB, putP, trpA, trpB, and uidA), no sequence type was discovered; nonetheless, the closest ST518 and ST966 were recommended (Supplementary Table S2). Several studies indicated the presence of specific STs in foods, whether of plant (ST2847, ST973) (Liu and Song, 2019; Song et al., 2020; Montero et al., 2021) or animal origin (ST1170) (Castellanos et al., 2017). E. coli isolates with phylogroup B2 and sequence type ST1170 were detected in wastewater from pig transport trucks (Savin et al., 2021) and from nosocomial and ambulatory infections in Germany (Pietsch et al., 2017). The ST966 type was identified in Klebsiella pneumoniae isolated from animal feed (Wu et al., 2022), suggesting environmental contamination. Cilantro-isolated E. coli samples, from a market in Germany, were classified as phylogroup B1 with multi-locus sequence types ST6186, ST165, ST58, and ST641 (Reid et al., 2020). Several E. coli strains, isolated from strawberries and banana obtained from agricultural farms in Ecuador, were classified into phylogroup B1, with the most frequent multi-locus sequence types being ST453, ST847, and ST6598 (Montero et al., 2021). Observed patterns of sequence types and phylogroups may be related with country- or producer-specific production methods, such as those pertaining to antimicrobial use, hygiene, and irrigation water systems (Molina et al., 2024). Based on the fimH typing, the L1PEag1 serotype ST1170 was classified as fimH20, with 100% sequence identity. Moreover, the subtyping results indicated the presence of the fliC gene, which matches the H-type H4 (100% identity with E. coli U1-41, AJ605764.1) (Beutin et al., 2005) and wzy genes (98.92% identity with the E. coli O1,H12 variant and 99.14% identity with the O1,H20 variant, KY115223.1). The presence of these variant alleles provided evidence for classifying the L1PEag1 isolate into the O1:H4 serotype. A single structural subunit (flagellin) encoded by the fliC gene defines the H antigen of E. coli (Reid et al., 1998), whereas the wzy gene encodes the O-antigen polymerase, which plays an important role in the synthesis of the lipopolysaccharide of bacteria (Zuo et al., 2019). The E. coli O1:H12 variant was detected in pig feces (Delannoy et al., 2017).

The genome contains 4,473 genes, 88 tRNAs, and 22 rRNA (8 copies of 5S rRNA, 7 copies of 16S rRNA, and 7 copies of 23S rRNA). A physical genomic map is shown in Figure 2A. Prokka was used to predict the location, while BLAST was used to infer the function and identification of assembled sequences against nucleotide and protein sequence databases. Predicted genes in the previous step were aligned with several databases to obtain their corresponding annotations with the aligners (Table 1). The number of genes associated with COG (3,723 genes) and KEGG (4,439 genes) functional annotation categories are shown in Supplementary Figures S1A,B. Plasmid annotation (Figure 2B) consists of the following features: (1) R–region: an endolysin with trans glycosylase activity that degrades host peptidoglycans and participates with the holin and spanin proteins in the sequential events which lead to the programmed host cell lysis releasing the mature viral particles; (2) Rz: Component of the spanin complex that disrupts the host outer membrane and participates in cell lysis during virus exit (Berry et al., 2013; Rajaure et al., 2015); (3) bor-not known function, is expressed during lysogeny in E. coli (bacteriophage); (4) ydfO-putative protein (COG5562), 74.9% similarity with E. coli K-12 strain; (5) S- [Isoform Antiholin]: Counteracts the aggregation of the holin molecules and therefore pore formation; (6) cos: lambda cos site; allows packaging into phage lambda particles; (7) ybcW: Protein inferred from homology (E. coli O157:H7); (8) ylcI: unknown function, protein predicted in the reference genome of E. coli K12.

Two sequences with CRISPR were found within the L1PEag1 genome (beginning at position165680 and ending at 1656883 or at position 4,818,020 and ending at 4818143, respectively). Both encompass a short spacer sequence, spanning between degenerate repeats (consensus DRs) with 95.65 and 96.87% conserved repeats respectively, and 100% spacer conservation. No Cas elements were detected. There are seven identified prophage regions; of which three are complete, three are incomplete, and one is questionable. The Enterobacteria phages mEp460, SfV, and cdtI were the most frequently occurring intact phages (Table 2). Previous studies indicated the presence of Enterobacteria phage mEp460 in three E. coli strains isolated from seawater and marine sediment, indicating its high adaptability to different environments (Lai et al., 2017). Phage SfV is a temperate serotype-converting phage of Shigella flexneri and encodes the factors involved in type V O-antigen modification (Allison et al., 2002). Shigella phage SfIV and Enterobacteria phage cdtI (cytolethal distending toxins) were detected in the human gut microbiome and were associated with diarrhea (Federici et al., 2023). Additionally, a total of 23 GIs of 238,133 bp were predicted with IslandViewer using as a reference the multi-drug resistant E. coli SMS-3-5 genome (Fricke et al., 2008). A total of 229 hypothetical proteins, virulence genes, and antibiotic resistance genes were annotated in the GIs. Due to their frequent rearrangements, excisions, transfers, and further acquisition of additional DNA, these genomic regions play a significant role in the rapid evolution, diversification, and adaptation of E. coli. Most of them encode several proteins involved in the defense mechanism, virulence factors, iron-acquisition systems, and toxins (Supplementary Table S3). Our study included both in vitro and in silico approaches to analyze the pathogenicity of L1PEag1. The virulence factors detected in silico in the L1PEag1 genome are shown in Table 3. Among 25 genes from the virulence database, 8 genes showed perfect identity match (ID % = 100 match for the given gene, covering the entire length of the virulence genes in the database, 100% alignment), 12 were non-perfect matches (ID% < 100%, covering the entire length of the virulence genes in the database, 100 alignment), and 5 were non-perfect matches (ID < 85%, input sequence length is shorter than the virulence gene length, 60% alignment). It has been shown that virulence-associated regions are specific to carriage patterns, generally corresponding with specific E. coli STs. Out of 120 selected E. coli strains, 96% showed the presence of fimbrial adhesin gene fimH, whereas 63% contained glutamate decarboxylase gene gad (Reid et al., 2020). In our study, both virulence factors, fimH and gad, were detected with 100% identity. fimH is a type 1 fimbriae detected in more than 80% of E. coli strains and plays a significant role in the virulence of E. coli extraintestinal enteropathogenic strains (Poole et al., 2017). The enzyme GAD (glutamate decarboxylase) has been reported to be limited to pathogenic E. coli, being a virulence marker for contaminated food and water (Grant et al., 2001). In addition, the yfcV gene encodes a fimbrinal protein, the OmpT gene encodes an outer membrane protease (protein protease 7), the fyuA encodes a siderophore receptor, and sitA encodes an iron transport protein. In the L1PEag1 genome, all three genes, OmpT, fyuA, and fimH, were found with 100% match identity, within the island picks and were previously detected in urinary infections (Spurbeck et al., 2012; He et al., 2015). Although not found with a 100% match, the capsular polysaccharide gene kpsM, serine protease autotransporter vat, and increased serum survival iss, which has been previously associated with ExPEC infection (Ravan and Amandadi, 2015), were detected in the target L1PEag1 genome. In addition, the virulence fimA, fimC, fimD, fimH, fimG, focC, kpsM, and papC genes were annotated with Prokka. Complementary PCR analysis confirmed that L1PEag1 harbored a gene that encodes for type 1 fimbriae (fimH) and pap (pili associated with pyelonephritis) but not for S fimbriae (sfa) (Supplementary Table S1). This virulence gene profile was previously detected in E. coli isolated from urine culture of a patient with urine infection, and they belonged to phylogroup B2 (Rahdar et al., 2015).

A total of two insertion sequences (ISs) and two miniature inverted-repeat transposable elements (MITEs) were grouped into 2 families (ISAs1 and IS630) and were retrieved with the ISfinder web tool. These genomic datasets can provide information on the ecological strength and adaptability of strains as well as their function in various environments. Besides, a total of 399 mobile elements were predicted with mobileOG-db (Brown et al., 2022), from which 50 were insertion/excision, 117 were replication/recombination/repair, 141 belonged to phage, 46 were stability/transfer/defense and 45 were transfer. Based on the 464 matched pathogenic families (10.44% of proteome), the results showed that the L1PEag1 isolate is a human pathogen (0.934 likelihood). Prokka annotation indicated the presence of the hlyE hemolysis gene in the L1PEag1 genome. The BLASTP analysis indicated a 100% identity to the hemolysin E protein family of the E. coli strain AR_451 (a clinical isolate, AWZ82741.1). Hemolysin hlyE is a novel pore-forming toxin and was first discovered in E. coli K-12 (Wyborn et al., 2004). This result was also complemented by the inhibition zone on blood agar media with the effect of hemolytic activity (data not shown). Although gelatinase is a less significant virulence factor in E. coli (Shruthi et al., 2012); in the current investigation, the gelatinase enzyme was produced by both L1PEag1 and UTNEc1. These results concur with earlier research, which showed that many E. coli isolates from urinary infections were gelatinase positive (Sharma et al., 2007).

Previous studies indicated that E. coli is a significant source of resistance genes that could be the cause of treatment failures in veterinary and human medicine. An increasing number of resistance genes have been found in E. coli isolates, and many of these resistance genes were acquired through horizontal gene transfer (Poirel et al., 2018). The putative antibiotic resistance genes identified in the genome of L1PEag1 is shown in Figure 3. A total of 4 perfect and 51 strict hits were detected, and the CARD resistant gene identifier report is shown in Supplementary Table S4. Within the former set of hits, the following antibiotic resistance genes were detected cpxA, mdtH, H-NS, and evgA, these confer resistance to aminoglycoside, aminocoumarin, fluoroquinolone and macrolide, cephalosporin, cephamycin, penam, and tetracycline. Previous research in E. coli indicated that Cpx-regulated genes are centrally involved in cell energetics, mediated envelope stress adaptation, and antibiotic resistance (Raivio et al., 2013). The mdtH (multidrug transport) gene is one of the 35 efflux pumping encoding genes detected in E. coli. Overexpression of cloned mdtD in the E. coli K-12 strain, which lacks the major efflux pump permease AcrB (E. coli KAM3) (Hirakawa et al., 2003a), results in a two-fold increase in resistance to norfloxacin and enoxacin (compared to the mutant parent) but does not impact the resistance to a range of other toxic compounds (dyes, detergents, antibiotics, and others) (Nishino and Yamaguchi, 2001). H-NS is a histone-like protein involved in global gene regulation in Gram-negative bacteria (Nishino and Yamaguchi, 2004). It is a repressor of the membrane fusion protein genes acrE, mdtE, and emrK as well as nearby genes of many intrinsic multidrug exporters (Nishino and Yamaguchi, 2004). EvgA, when phosphorylated, is a positive regulator for efflux protein complexes emrKY and mdtEF (Hirakawa et al., 2003b) which confer resistance to erythromycin, cloxacillin, tetracycline, oxacillin, and norfloxacin. Notably, an SNP mutation (R234F) was found in the E. coli elongation factor EF-Tu, which confers resistance to pulvomycin. The E. coli glpT gene (E448K) confers resistance to fosfomycin. The UhpT (E350Q) carries a mutation that confers resistance to fosfomycin; The PBP3 (D350N, S357N) gene, similar to Haemophilus influenzae, carries a mutation that confers resistance to beta-lactam antibiotics; and E. coli AcrAB-TolC bears MarR mutations (Y137H, G103S), which confer resistance to ciprofloxacin and tetracycline (Supplementary Table S4). The resistance to azithromycin (AZM15), clindamycin (DA2), novobiocin (NV30), amikacin (AMK30), oxacillin (OX1), erythromycin (E15), bacitracin, and other antibiotics was confirmed by the disk assay method; out of the 21 antibiotics tested, the target strain showed resistance to 14 antibiotics belonging to several classes as shown in Figure 4. Antibiotic resistance was previously evaluated in several E. coli strains isolated from lettuce, cabbage, cucumbers, and tomatoes (Rahman et al., 2021). Besides, 100% of the E. coli isolates obtained from fruits and vegetables collected from various regions in Ecuador, exhibited resistance to ampicillin, cefazolin, cefotaxime, and tetracycline (Montero et al., 2021). The strains isolated from fresh produce have demonstrated variable resistance to colistin (CST), AMK, E15, cefotaxime (CTX), ceftazidime (CAZ), AMP, gentamycin (GEN), and amoxicillin (AMX) (Shah et al., 2015; Gómez-Aldapa et al., 2016). In addition, E. coli isolated from fresh vegetables (i.e., pak choi and lettuce) were resistant to kanamycin, levofloxacin, doxycycline, fosfomycin (FOS), CTX, AMP, GEN, NA, TE, and colistin (CST) (Liu B-T et al., 2019; Manageiro et al., 2020). PCR analysis results confirmed the resistance to tetracycline (bla_TEM-1) and carbapenemase (bla_KPC) (Supplementary Table S1). The gene bla_TEM-1 is widely found in animal-derived E. coli and encodes for β-lactamases with a narrow spectrum that are capable of inactivating aminopenicillins and penicillins (Poirel et al., 2018). In the study of Montero et al. (2021) there were differences in a subset of ESBL-E. coli isolates that tested positive for certain genes by PCR but negative by whole genome sequencing (WGS); these were 12 isolates with blaTEM, nine isolates with blaSHV, and one isolate with blaCTX-M. Furthermore, WGS revealed blaSHV and blaTEM in two isolates, but PCR results were negative. A recent study revealed a relatively high prevalence of bla_KPC–carrying, carbapenemase-resistant Enterobacteriaceae in fecal isolates, among food handlers in Kuwait (Moghnia et al., 2021). Nonetheless, non-specific amplification products were detected for bla_VIM (Verona integron-encoded metallo-β-lactamase), oxacillinases (OXA48/181), and cefotaxime (bla_CTXM-9). This might be produced by primers binding to seemingly random locations in the sample DNA. No amplification was detected for the other tested β-lactam genes (Supplementary Table S1).

Based on RGI analysis, the sitABCD genes were detected in the genome of the L1PEag1 isolate (i.e., sitA: a periplasmic binding protein, sitB: an ATP-binding component, sitC: an inner membrane component, and sitD: an inner membrane component, which conferred resistance to hydrogen peroxide). Along with other ion transport systems, SitABCD is an iron and manganese transporter that may confer oxidative stress resistance and bactericidal effects. These genes showed 97.31% identity with a sitABCD gene located on the virulence plasmid pAPEC-1 of an avian pathogenic E. coli strain (AY598030.2) (Sabri et al., 2006).