We analyzed a collection of 53 STEC O139:H1 strains isolated from pigs affected by ED in Italy, maintained in the National Reference Laboratory for E. coli collections at the Istituto Superiore di Sanità. Additionally, we included in this study 83 more STEC O139:H1 genomes retrieved from the GenBank and the European Nucleotide Archive (ENA) databases, from strains isolated from pigs or other sources ( Supplementary Tables 1 , 2 ).

To conduct whole genome sequencing (WGS), total DNA was extracted from a 2 mL overnight TSB culture of each strain grown at 37°C using the GRS Genomic DNA Kit Bacteria (GRISP Research Solutions, Porto, Portugal). The majority of sequences were generated using Ion Torrent sequencing technology (Thermo Fisher Scientific, MA, USA). Sequencing libraries of approximately 400 bp were prepared from 100 ng of total DNA using the NEBNext Fast DNA Fragmentation & Library Prep Set for Ion Torrent (New England BioLabs, MA, USA). These libraries were then processed through emulsion PCR and enrichment on the Ion OneTouch 2 System, followed by sequencing on an Ion Torrent S5 platform (Thermo Fisher Scientific, MA, USA) using the ION 520/530 KIT-OT2 (Thermo Fisher Scientific, MA, USA) according to the manufacturer’s instructions. All genomic sequences are accessible at the GenBank (BioProject: PRJNA1152229).

Most of the bioinformatic analyses to characterize the genomes were carried out using the tools available on the Galaxy public server ARIES (Istituto Superiore di Sanità, https://www.iss.it/site/aries) (Knijn et al., 2020).

Single-end reads from the Ion Torrent S5 platform were assembled using SPADES version 3.12.0 with default parameters (Bankevich et al., 2012), and filtered using the Filter SPAdes repeats tool (https://github.com/phac-nml/galaxy_tools) with default settings to eliminate repeated contigs or those <1,000 bases in length. Paired-end reads were trimmed, filtered using the Extended Randomized Numerical alignEr–filter (Del Fabbro et al., 2013), and de novo assembled using SPAdes version 3.10.0 (Bankevich et al., 2012).

Multilocus sequence typing was performed using the MentaLiST tool version 0.2.3 (Feijao et al., 2018), following the scheme developed by Wirth et al (Wirth et al., 2006). The assembled contigs were analyzed using BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) and the blastn algorithm version 2.7.1 using different databases. Serotyping was determined by aligning the contigs with reference sequences for the O and H antigen genes (Joensen et al., 2015). Subtyping of stx was performed using the Shiga toxin typer tool v2.0 (https://github.com/aknijn/shigatoxin-galaxy), which conducts an optimized blastn search against the sequence database of stx subtypes developed by the Statens Serum Institut (https://bitbucket.org/genomicepidemiology/virulencefinder_db/src/master/stx.fsa).

BLAST + blastn was also used to detect the presence of plasmid-related sequences and their replicon types within the whole genome of the O139:H1 strains using the PlasmidFinder and PLSDB tools, with a minimum identity of 95% and minimum coverage of 60%. Additionally, we retrieved the most closely related plasmid sequences available online from the National Center for Biotechnology Information (NCBI) nucleotide collection database (accessed June 1, 2024). The focus of this analysis was on determining the similarity of these identified sequences to known reference plasmid sequences.

The identification of resistance genes was conducted by submitting the complete reference plasmid nucleotide sequence to the ResFinder web server with default parameters, which required a minimum identity of 90% and minimum coverage of 60% (http://genepi.food.dtu.dk/resfinder) (Zankari et al., 2012). Similarly, the reference plasmid was submitted to the VFDB web server via the BLAST sequence-similarity (setB), also with default parameters, to identify virulence factors (http://www.mgc.ac.cn/VFs/search_VFs.htm) (Liu et al., 2022).

The Prokka tool (Galaxy Version 1.14.5) (Seemann, 2014) was used to perform functional annotation on the assembled sequences of the O139:H1 strains and reference plasmids, utilizing the E. coli-specific gene database and default parameters. Additionally, the Blast Ring Image Generator (BRIG) software v0.95 (Alikhan et al., 2011) was used with default parameters to compare and visually represent these plasmid-like sequences in relation to reference plasmids, highlighting regions of similarity and potential functional significance.

The associations and clustering between the virulence genes identified on the reference plasmids were assessed using the STRING database (https://string-db.org/). The analyses were carried out following the updated instructions (The STRING database in 2023) (Szklarczyk et al., 2023).

We performed the analysis using core genome multilocus sequence typing (cgMLST) with the chewBBACA tool and the INNUENDO project’s scheme, available on the Galaxy public server ARIES, which includes 2,360 loci (Llarena et al., 2018; Silva et al., 2018). Pairwise comparisons were deemed reliable when more than 80% of loci for each sample were assigned to an allele. The distances between strains were calculated by comparing allelic profiles pairwise, using the chewTree tool on the ARIES webserver. For each sample pair, alleles that were missing, partially identified, or incorrectly assigned to any locus were excluded. The resulting dendrogram was visualized using Newick Display on Galaxy Version 1.6 (Junier and Zdobnov, 2010).

Among the genomes of the 53 Italian strains examined, 50 were classified as ST1, and 3 as ST955. All strains were identified as O139:H1 serotype, with all the 53 strains carrying the stx2 gene only, specifically subtype stx2e.

Among the additional 83 genomes retrieved from the public domain, 79 were classified as ST1, two as ST955, one as ST10859, and one as ST114. Of these, 79 strains carried the stx2 gene, subtype stx2e, while the stx2 subtypes in four isolates were not identified ( Supplementary Tables 1 , 2 ).

Three different reference plasmids were identified in the genome of studied isolates based on the similarity between plasmid-like sequences in the O139:H1 strains and the identified reference plasmid sequences according to the PLSDB ( Figure 1 ): Escherichia coli O139:H1 strain W13-16 plasmid pW1316-2 (Accession number: NZ_CP080237.1), Shigella flexneri 2a strain ATCC 29903 (Accession number: CP026790.1), and Salmonella enterica subsp. enterica serovar Typhimurium strain 21G7 isolate B71 plasmid pB71 (Accession number: NZ_KP899806.1). Moreover, some genome returned hits against Escherichia coli O111:H- str. 11128 plasmid pO111_1 (Accession number: NC_013365.1), Escherichia coli strain 15OD0495 plasmid p15ODAR (Accession number: NZ_MG904995.1), Escherichia coli strain ESBL3153 plasmid pESBL3153-IncX4 (Accession number: NZ_MW390521.1), and Escherichia coli strain 20Ec-P-124 plasmid pMRY16-002_3 (Accession number: NZ_AP017613.1) plasmids ( Figure 1 ).

We identified seven different replicon types among the plasmid-related sequences in the O139:H1 genomes based on the PlasmidFinder which included: IncI1-I(Alpha), IncI2, IncFIA(HI1), IncHI1B(R27), IncFII, IncX1, IncX4 ( Figures 2 , 3 ).

As far as the additional genomes included in the study are concerned, the plasmids Escherichia coli O139:H1 strain W13-16 plasmid pW1316-2 (Accession number: NZ_CP080237.1), Shigella flexneri 2a strain ATCC 29903 (Accession number: CP026790.1), and Salmonella enterica subsp. enterica serovar Kentucky plasmid pCS0010A_9 (Accession number: NC_019104.1) were the most frequently identified reference plasmids to which the genomes demonstrated similarity. Moreover, the replicon types IncFII, IncX1, IncI1-I(Alpha), and IncHI1B(R27) were the most frequently detected replicon sequences ( Supplementary Tables 1 , 2 ).

Overall, the most frequent plasmid signatures identified in all strains (O139:H1 STEC isolates of our and other studies), were to pW1316-2, Shigella flexneri 2a, and pB71 plasmids.

Our analysis using BRIG software revealed the presence of hha, yhcR, finO, aidA-I, tibC, rhaR_1, tpx, rhaR_2, faeE, elfC, bin3, pir, dnaT, hlyD, hlyB, hlyA, hlyC, and topB genes on the pW1316-2 plasmid ( Figure 4 ).

Shigella flexneri 2a plasmid, harbored speE, higB_1, RepB_1, RepB_2, dam, noc, traC, virB, parM, uvrD, tus, hha, repE, hns, smc, dcm, umuD, umuC, dsbC, and resA genes ( Figure 5 ).

Ultimately, in pB71 plasmid, the genes uvrD, tus, hha, repE, hns, smc, dcm, yhcR, cynR, gltS, yqjZ, tetR, tetA, tetC, folP, emrE, ant1, xerC, hin, cat, catM, ifcA, proP, umuD, umuC, dsbC, resA, corA, higB-1, repB_2, dam, repB_1, traC, virB, and parM were present ( Figure 6 ). Most of the genes detected were present in all the reference plasmid sequences identified in the studied genomes ( Figures 4 - 6 ).

The virulence factors hlyA, hlyB, aidA-I, and faeE were identified on pW1316-2 plasmid and hns, traC and hha were detected both on the Shigella flexneri 2a and pB71 plasmids according to the VFDB web server ( Table 1 ).

Based on the ResFinder analysis on the reference plasmids, the pW1316-2 and Shigella flexneri 2a plasmids had no known antibiotic resistance genes; and in the pB71 plasmid, we detected tet genes (tetA, tetC, and tetR) responsible for Tetracycline antibiotic resistance ( Table 1 ).

Amid the identified plasmids, we observed functional and regulatory associations among the genes encoded on the reference plasmid derived from the Shigella flexneri 2a ( Figure 7 ). The core of these clustering is centered around the Hha-H-NS complex, the functions of which ultimately lead to the production of α-hemolysin, an important virulence factor in ED (Details of associations are provided in tabular form in Supplementary Table 3 ).

To explore the phylogenetic relationships and their association with plasmid characteristics, we conducted a cluster analysis using cgMLST for comparative purposes. In addition, we included 28 STEC O139:H1 genomes from various global sources, retrieved from GenBank and ENA databases ( Supplementary Table 2 ). We then calculated the number of allelic differences between strains ( Supplementary Table 4 ). The analysis grouped the strains into nine main clades ( Figure 8 ). Clades 1-2 and 6-9 exhibited significantly lower allelic distances (AD) compared to the other clades, all showing fewer than 90 AD. The majority of the 28 additional strains from external studies were incorporated into clade 5, which was phylogenetically close to clade 3, with most strains showing fewer than 30 AD relative to clade 3 ( Figure 8 ). The most genetically distant strains were found in clade 4 ( Figure 8 ). Notably, a remarkable correlation was observed between allelic differences in the cgMLST analysis and plasmid characteristics, with strains carrying specific plasmids and replicon types clustering into particular phylogenetic clades ( Figure 8 ).