Plasmids shape the diverse accessory resistomes of Escherichia coli ST131.

The human gut microbiome includes beneficial, commensal and pathogenic bacteria that possess antimicrobial resistance (AMR) genes and exchange these predominantly through conjugative plasmids. Escherichia coli is a significant component of the gastrointestinal microbiome and is typically non-pathogenic in this niche. In contrast, extra-intestinal pathogenic E. coli (ExPEC) including ST131 may occupy other environments like the urinary tract or bloodstream where they express genes enabling AMR and host cell adhesion like type 1 fimbriae. The extent to which commensal E. coli and uropathogenic ExPEC ST131 share AMR genes remains understudied at a genomic level, and we examined this here using a preterm infant resistome. We found that individual ST131 had small differences in AMR gene content relative to a larger shared resistome. Comparisons with a range of plasmids common in ST131 showed that AMR gene composition was driven by conjugation, recombination and mobile genetic elements. Plasmid pEK499 had extended regions in most ST131 Clade C isolates, and it had evidence of a co-evolutionary signal based on protein-level interactions with chromosomal gene products, as did pEK204 that had a type IV fimbrial pil operon. ST131 possessed extensive diversity of selective type 1, type IV, P and F17-like fimbriae genes that was highest in subclade C2. The structure and composition of AMR genes, plasmids and fimbriae vary widely in ST131 Clade C and this may mediate pathogenicity and infection outcomes.

Data Summary

The following files are available on the FigShare project ‘Plasmids_ST131_resistome_2020’:

The set of 794 AMR genes derived from [1] are available (with their protein sequence translation) at FigShare at doi: 10.6084/m9.figshare.11961402.

The AMR gene profiles per sample determined by their blast sequence similarity results against CARD are available at FigShare at doi: 10.6084/m9.figshare.11961612. This dataset includes the PlasmidFinder results. It also includes other AMR database comparisons (ARG-ANNOT, ResFinder, MegaRes, VFDB and VirulenceFinder).

The blast sequence similarity results for the fim, pil, pap and ucl operons’ genes versus 4071 ST131 assemblies from Decano & Downing (2019) are available at FigShare at doi: 10.6084/m9.figshare.11961711.

The genome sequences and annotation files for reference genomes NCTC13441, EC958 and SE15, along with the assembled contigs for 83 972 and 3_2_53FAA are available at FigShare at doi: 10.6084/m9.figshare.11961813.

The 4071 ST131 genome assemblies from Decano & Downing (2019) are available at FigShare at doi: 10.6084/m9.figshare.11962278 (the first 1680 assemblies) and at doi: 10.6084/m9.figshare.11962557 (the second 2391 assemblies).

Introduction

Extra-intestinal pathogenic (ExPEC) cause extensive infections outside the gut, from which they can originate. ExPEC have a wide range of virulence factor (VF) [2-5] and antimicrobial resistance (AMR) genes [6], especially sequence type (ST) 131 from serotypes O25b:H4 and O16:H5 in phylogroup B2 [7]. ST131 causes a substantial fraction of ExPEC and extended-spectrum beta-lactamase (ESBL)-producing cases [8]. ST131’s acquisition of key AMR genes encoded on plasmids [5, 9] has coincided with its adaptation to new environments [10, 11].

ExPEC adherence factors allowing colonisation of different host niches are mainly produced by the Chaperone-Usher secretion pathway, including type 1, P, F1C/S and AFA fimbriae [12-15]. Within ExPEC ST131, fluoroquinolone-resistant Clade C is the main cause of human infection globally [16, 17]. Its d-mannose-specific type 1 fimbriae encoded by the fim operon binds the host mucosal epithelium to cause urinary tract and kidney infection [18-21]. The genetic diversity of fimbrial operons beyond fim in Clade C has not yet been examined in large collections of isolates.

ExPEC infection can be associated with a changed microbiome composition [22]. is the most common initial coloniser of infant intestines [23], where most are commensal [24] and some protect against pathogen invasion [25]. Nonetheless, AMR is prevalent in neonatal units [26] and within-host gene exchange between commensal and pathogenic bacteria may occur [27]. ST131 spreads between mother-infant pairs [28] likely via an oro-faecal transmission route [29], and such colonisation of infants can last for long periods [30]. Consequently, AMR gene screening can inform on treatment strategies [31] and the effect of antibiotics on microbiomes [32].

The ST131 resistome (the set of AMR genes) includes ESBL genes [33, 34] allowing third-generation cephalosporin-resistance [35] and are associated with three main cefotaximase (CTX-M) resistance alleles: bla CTX-M-14/15/27 [36]. Like most AMR genes, these are sandwiched by mobile genetic elements (MGEs) on plasmids and thus can be gained by horizontal gene transfer (HGT) [37-39] or lost if not beneficial [40, 41]. Of nine common bacterial pathogens, has the most MGEs, including phage-associated ones and transpose (tnpA) genes [42]. Diverse pathogenic bacteria share MGEs [43-45] and HGT may worsen nosocomial outbreaks [46-48]. Plasmids common to ST131 often have ISEcp1, IS903D and IS26 elements [49] encoding tnpA flanked by short inverted repeats, typically with ESBL genes. Such ESBL genes can be chromosomally transferred and may form part of the core resistome [50].

Fluctuating antibiotic exposure, host type, anatomical niche and immunity mean that resistomes vary [51]. Conjugation, recombination and MGEs help shape ST131 resistome dynamics [2, 11, 52–55]. ST131’s plasmids typically are from incompatibility (Inc) groups F, I1, N and A/C [56-58]. Some plasmids in ST131 encode genes for post-segregation killing and stable inheritance to ensure their propagation, but these genes can be lost or may recombine with other plasmids [58-60]. As a result of this mixing and their extensive array of MGEs, plasmids may rearrange extensively even within a clonal radiation [61, 62]. Plasmids may also impair cell reproduction due to the energetic cost of their replication and maintenance, so conjugation and recombination could allow gene dosage optimisation and gene expression coordination.

There is extensive evidence of AMR gene conjugation in the gut between species including [63-66], and also between its STs, such as: transfer of a bla CTX-M-1-positive Incl1 plasmid among ST1640, ST2144 and ST6331 in a single patient’s gut [67]; a bla OXA-48 gene on a IncL/M-type plasmid from ST14 to ST666 [68] and from ST648 to ST866 [69]; a 113.6 kb IncF plasmid with a mercury detoxification operon bla KPC-2, bla OXA-9 and bla TEM-1A genes from ST69 to ST131 [70]; and transfer of a range of sulphonamide- (sul2) and ampicillin-resistance genes (bla TEM-1b) on a pNK29 plasmid within subtypes [71, 72]. Moreover, the frequency of conjugation was ten times higher and more stable in bla CTX-M-15-producing ST131 than in [73].

This study focused on uropathogenic ST131 Clade C (fimH30) because of its high plasmid and AMR gene diversity [50, 61] where the mobile resistome and how plasmids mediate AMR gene composition needs deeper scrutiny. Here, we found a large shared core preterm infant resistome and a smaller accessory one that varied between closely related isolates. This was caused by plasmid turnover and rearrangements, affecting fimbrial genes too. We hypothesise that certain plasmids are more compatible with ST131 genomes, and applied topological data analysis (TDA) to understand plasmid-chromosome co-evolution.

Results

A large core preterm infant resistome in

We examined the preterm infant resistome and plasmid composition of seven ST131 genomes in the context of related reference genomes and non-pathogenic isolates. This resistome of 794 genes was assembled from 2004 contigs originally derived from DNA sequencing of 21 preterm infant faecal samples experimentally tested in vitro for resistance to 16 antibiotics (see Data Access, Table S1, available in the online version of this article) [1]. This resistome provided a culture-unbiased perspective on antibiotic resistance genes [1]. We found that seven ST131 Clade C genome assemblies from adult urinary tract infections (UTIs) together with subclade C2 reference genomes EC958 and NCTC13441 had 262 to 291 unique contigs in this resistome (Table 1). Commensal ST131 Clade A reference SE15 (O150:H5, fimH41) had 251 AMR contigs, and intestinal Human Microbiome Project (HMP) assembly 3_2_53FAA had 269. In contrast, well-studied asymptomatic urinary tract HMP sample 83 972 from ST73 [74] had 213 unique contigs, suggesting a subset of AMR genes may be associated with symptomatic ST131 UTIs.

Table 1.

genome assembly collection and metadata. The table shows the sample ID, SRA accession ID, ST131 Clade, bla CTX-M allele(s), country and year of isolation, and numbers of unique AMR contigs in the preterm infant resistome. All seven ST131 assemblies, EC958 and NCTC13441 came from pathogenic UTIs. SE15 was a faecal commensal isolate. EC958 had 18 AMR contigs on plasmid pEC958A (pEC958B had none). NCTC13441 had 27 AMR contigs on its plasmid. The seven ST131 and two HMP samples were assembled from Illumina reads, and the three references were assembled from PacBio reads

Sample ID	SRA ID	ST131 clade or ST	bla CTX-M allele(s)	Country, year	Unique AMR contigs
8289_1#91	ERR191724	C1	14	Ireland, 2007	268
8289_1#35	ERR191668	C1	14	Ireland, 2008	270
8289_1#3	ERR191636	C2	15	Ireland, 2005	274
8289_1#34	ERR191667	C2	15	Ireland, 2010	274
8289_1#24	ERR191657	C2	14 and 15	Ireland, 2009	283
8289_1#60	ERR191693	C2	15	Ireland, 2010	262
8289_1#27	ERR191660	C2	15	Ireland, 2010	262
SE15	AP009378	A	none	Japan, 2010	251
NCTC13441	ERS530440	C2	15	UK, 2003	280
EC958	HG941718	C2	15	UK, 2005	291
3_2_53FAA	–	ST95	none	Canada, 2007	269
83 972	–	ST73	none	Sweden, 1978	213

–

–

Comparing the AMR contigs’ gene product functions identified a core resistome of 210 AMR contigs shared by all (Fig. 1), highlighting that these may be essential for non-pathogenic isolates. This core resistome differed from the accessory resistome of 109 contigs, of which 50 in ST131 Clade C alone corresponded to 18 unique AMR genes (Fig. 1). These accessory AMR genes were identified with ARIBA v2.11.1 [75] and the Comprehensive Antibiotic Resistance Database (CARD) [76], where matches were verified by read mapping to the reference resistome using GROOT [77].

Fig. 1.

The overlap of preterm infant AMR contigs across four ST131 subclade C2 assemblies (8289_1#3, 8289_1#24, 8289_1#27, 8289_1#34), two ST131 reference genomes (EC958 and SE15) and two HMP assemblies (83 972 and 3_2_53FAA). Top: The intersection sizes (y-axis) and the numbers of AMR contigs per set showed that most (210) AMR contigs were shared across all isolates. The non-unique AMR contigs indicated smaller numbers of unique AMR genes, so EC958’s 17 AMR contigs (red) corresponded to nine unique bla CMY genes (Table S2); and 8289_1#24’s nine contigs (green) represented four aminoglycoside and tetracycline resistance genes (Table S3). All Clade C bar 8289_1#27 had 12 contigs (blue) denoting to five bla TEM genes (Table S4). Bottom: The numbers of AMR contigs per sample with the corresponding sets (coloured circles).

A diverse accessory preterm infant resistome among closely related ST131 isolates

Bla CMY encoding ESBL AmpC is typically associated with cephalosporin-resistance [78] and here nine variants (alleles CMY-37, 51, 66, 67, 85, 98, 101, 105) resolved from 17 AMR contig matches were unique to EC958 (Fig. 1). All Subclade C2 isolates bar 8289_1#27 had five bla TEM alleles (TEM-57, 104, 116, 215, 220) encoding an Ambler class A β-lactamase. The sole bla CTX-M-14- and bla CTX-M-15-positive isolate (8289_1#24) had unique four genes linked to aminoglycoside and tetracycline resistance (Table S2), a combination that may neutralise the effect of aminoglycoside susceptibility associated with tetracycline efflux pump TetA expression [79, 80].

Additional comparisons using the NCTC13441 genome showed that it had five bla OXY alleles (OXY-1–1, 1–2, 2–3, 6–2, 6–4) encoding an Ambler class A ESBL (Text S1). Ten AMR contigs in all C2 corresponded to a single helix-turn-helix-like transcriptional regulator gene (Fig. 1). Two contigs present in all ST131 encoded araC (a transcriptional regulator) and sugE (associated with resistance to quaternary ammonium salts). The 39 AMR contigs in all bar 83 972 included 28 copies of mdfA, encoding a multi-drug efflux pump linked to broad-spectrum AMR [81]. The 3_2_53FAA and 83 972 shared ftsI encoding a transpeptidase catalysing peptidoglycan crosslinking [82] with the synonym penicillin-binding protein 3 (PBP3) [83, 84], and is sensitive to ESBLs [85].

For a wider context on ST131 dynamics, comparing this preterm infant resistome with showed that an interspecies shared resistome had 58 genes and that the larger core resistome of 308 genes included bla SHV. Nine of ten genes in pathogenic isolate PMK1 that were absent from a urinary tract microbiome assembly (WGLW1) were bla OXY, NCTC13441. Repeating this with the more divergent and lugdunensis found 14 shared AMR genes on plasmids only.

Extensive plasmid rearrangements in closely related ST131 Clade C genomes

To evaluate the origin and relationships of the AMR genes in these seven Clade C (Table 1), we initially found with PlasmidFinder [86] that all had IncF1A and some C2 had lost IncF1B, suggesting potential plasmid changes (Table S3). We aligned these seven assemblies to plasmids common in ST131 (pEK499, pEK516, pEK204, pCA14, pV130) and repeated this via read mapping. Plasmids pEK499 (from England), pEK516 (England) and pEK204 (Belfast) were geographically near these ST131 (all from Ireland). SE15’s pECSF1 122.3 Kb IncF2A/F1B plasmid was used as a control because it had no known AMR genes [87]. The 8289_1#35 and 8289_1#91 (both from C1) had identical results, and vice versa for 8289_1#60 and 8289_1#27 (both from C2) and so are not discussed below.

Clade C ST131 had more similarity to IncF2/F1A bla CTX-M-15-positive plasmid pEK499 than SE15 and the two HMP samples, but to varying degrees (Fig. 2). This plasmid has 185 genes, lacks traX for conjugation, is stably inherited (Table S4), and may have been gained early in the origins of Clade C as a fusion of type F2 and F1A replicons [88]. Nine of pEK499’s ten key AMR genes are in a 25 Kb segment spanning bla TEM, bla OXA-1 and bla CTX-M-15 at 40, 58 and 63 Kb (respectively). The 8289_1#3 and 8289_1#24 (both C2) had all three, 8289_1#27, 8289_1#60 and 8289_1#34 (all C2) had bla OXA-1 and bla CTX-M-15 but not bla TEM, whereas C1 (8289_1#35 and 8289_1#91) were bla TEM- and bla CTX-M-14-positive. Conjugation (tra) genes were in all ST131 except 8289_1#34. IncF2A plasmid pEK516 has ~75 % similarity to pEK499 but is shorter (64.6 Kb) [59] and is usually non-conjugative [89]. This plasmid has 103 genes including bla TEM, bla OXA-1 and bla CTX-M-15, and a type I partitioning locus for stable inheritance absent in pEK499 (Table S4). Read mapping to pEK516 showed that Clade C had high similarity to it, unlike SE15 and the HMP samples (Fig. S1).

Fig. 2.

Comparison of pEK499 (top, 117 536 bp) and pEK204 (bottom, 93 732 bp) with five ST131 C1 and C2 isolates, ST131 Clade A reference SE15 (navy) and two HMP assemblies (3_2_53FAA in dark green and 83 972 in cyan). Top: Normalised read coverage showed high copy numbers of IS1 (at 41 and 103 Kb of pEK499) and IS66 (at 75–77 Kb of pEK204, and at 113 Kb of pEK499). Bottom: blast alignments showed limited similarity for the HMP assemblies and SE15 relative to higher levels for Clade C for pEK499 : 8289_1#27 (C2, mauve), 8289_1#24 (C2, pink), 8289_1#35 (C1, grey), 8289_1#3 (C2, light green), 8289_1#34 (C2, beige). For pEK204, only 8289_1#27 had many regions of similarity. Genes encoding bla TEM, bla OXA-1 and bla CTX-M-15 are at 40, 58 and 63 Kb on pEK499 (respectively). The bla CTX-M gene was at 8 Kb, bla TEM was at 13 Kb, followed by mixed conjugation and segregation genes at 36–70 Kb on pEK204. For pEK204, the bla CTX-M-3 gene differs from bla CTX-M-15 by a single R240G substitution, and so bla CTX-M genes detected here were bla CTX-M-15. The annotation was modified from [59]. Matches spanning >300 bp are shown.

Contrasting sharply with pEK499 and pEK516, 8289_1#27 (C2) reads mapped to conjugative IncI1 plasmid pEK204 found 72.8 Kb regions of similarity, unlike the others that had <10 Kb (Fig. 2). PEK204 is similar to IncI1 plasmid R64, has 112 genes, a complete tra region (Table S4), and a type I partitioning locus [59]. 8289_1#27 had a complete I1 replicon, oriT, tra, shufflons subject to site-specific recombinase activity, and a pil cluster encoding a type IV fimbriae associated with enhanced cell adherence and biofilm formation in entero-aggregative and Shiga toxin-producing [90]. A 9.3 Kb region on pEK204 contained bla TEM-1b, inactive transposase Tn3-tnpA, ISEcp1, bla CTX-M and a 5′ orf477-tnpA-tnpR region. The 14 bp IRL 5′ of ISEcp1 and IRR at the distal end of the inverted orf477 element can mobilise bla CTX-M, and an additional IRR at impB 3′ of the bla TEM-1b gene (7.4 Kb further away) can mobilise the 9.3 Kb unit [91], which has been found on IncFIA, IncFIA-FIB, IncN and IncY plasmids after originating on a pCOL1b-P9-like plasmid [59]. ImpB and impA encode an error-prone DNA repair protein (like UmuDC) [92].

Further investigation of the ST131 plasmid-matching regions showed variable plasmid similarity within closely related isolates. Read mapping to conjugative plasmid pCA14 showed that all Clade C bar 8289_1#3 had Mrx and mph(A) genes associated with macrolide resistance (Fig. S3). This was also found by comparing with contigs from non-conjugative plasmids pV130a and pV130b from sewage treatment plant water in India [93], where 8289_1#27 and 8289_1#60 (both C2) had similarity spanning all the pV130a contigs (Fig. S2).

High rates of pEK499 and pEK204 protein-protein interactions with chromosomal Proteins

The extensive diversity of Clade C AMR genes and plasmids raised the question of how the plasmids’ gene products interact with chromosomal ones. Protein-protein interaction (PPI) networks can examine plasmid-chromosome coevolution based on gene products’ topological proximity [94, 95], which could be higher for plasmid and chromosomal proteins that have co-existed and so may interact more [96]. We used topological data analysis (TDA) to measure the number of non-trivial (indirect) loops where a loop is a chain of at least four PPIs ending at the same protein where it started. We focus on a type of loops called ‘non-trivial’ (see Methods). The number of non-trivial loops was scaled by the number of PPIs per dataset (Table 2).

Table 2.

The numbers of protein-coding genes per plasmid, unique genes per plasmid, length, and numbers of non-redundant PPIs within the plasmid and with the chromosome, and the rate of indirect interactions (non-trivial loops) per PPI. Plasmid pEK516 had the same results as pEK499

Plasmid	No. of genes	No. of unique genes	Length (Kb)	No. of unique PPIs		Non-trivial loops per PPI
				Within plasmid	With chromosome
pEK204	112	87	93.7	21	548	0.096
pEK499	185	87	117.5	8	758	0.117
pEK516	103	55	64.6	8	758	0.117
pJIE186-2	138	78	137.7	2	32	0.015
pCA14	181	93	155.4	5	0	0
pEC958A	142	75	135.6	5	0	0
pECSF1	150	70	122.3	0	0	0

PPI data for 4146 protein-coding genes with 105 457 PPIs (all combined scores >400) from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database v10 [97] was used to get the numbers of non-redundant PPIs and loops per PPI for a plasmid among its own proteins, and then between the plasmid’s proteins and the 4146 chromosomal proteins. We assessed the numbers of PPIs from a combined score threshold of 400 to 900 with a step size of 25. We tested our approach using a published clustering of 60 genes [98] that had a rate of 0.244 non-trivial loops per PPI, which was lower than the rate of 0.518 obtained for the 4146 chromosomal proteins (Table S5, Fig. S4). There were lower rates of non-trivial loops per PPI for pEK499 (0.117) and pEK204 (0.096), but the plasmids’ rates were almost constant across the combined score thresholds (Fig. S5). This held even though the number of PPIs per protein was negatively correlated with the combined score (Fig. S6), suggesting that the rate of non-trivial loops per PPI was robust to changes in the combined score threshold. We found that pEK499 and pEK204 interacted with chromosomally encoded proteins, but pCA14, pEC958A and pECSF1 did not (Table S6), and pJIE186-2 had a small number of PPIs (Table 2) that may be VFs [99].

ST131 genomes had an ancestral pEK499-like plasmid but some gained a pEK204-like one

The seven ST131 assemblies above were from a set of 4071 ST131 from Clades A, B and C [61] that were aligned here to pEK499 and pEK204 with blast. Of these 837 (20.6 %) had >10 Kb of regions similar to pEK204 (Fig. S7), whereas 3108 (76.3 %) had >10 Kb like pEK499 (Fig. S8). All 193 assemblies with >40 Kb of pEK204-like segments were from Clade C: 17 from C0 (out of 52, 33 %), 82 from C1 (out of 1119, 7 %) and 94 from C2 (out of 2051, 5 %) (Fig. 3). All 17 from C0 had an I1 replicon, bla CTX-M gene, bla TEM gene (bar one isolate) and at least a partial tra region, and 13 had the pil operon (76 %). The 82 from C1 and 94 from C2 had lower rates of I1 replicon presence (C1 77, C2 81 %), partial tra regions (C1 94, C2 86 %) and pil clusters (C1 n=60 or 73 %, C2 n=71 or 76 %) but differed in the rates of gene presence for bla CTX-M (C1 28 % vs C2 83 %, odds ratio=12.5, 95 % CI 6.1–25.8, P=e-11) and bla TEM (C1 59 % vs C2 30 %, odds ratio=3.3, 95 % CI 1.8–6.2, P=2.2e-5). This implied bla CTX-M (via orf477) was common in C2 and bla TEM (via impB) in C1 due to different ancestral transposition of the 9.3 Kb region.

Fig. 3.

The distributions of pEK204-like regions in ST131 Clade C genome assemblies with >1 Kb (left, n=3234) and >40 Kb (right, n=193) of matching segments. The regions of similarity were examined in (left) Subclades B0 (n=14 in black), C0 (n=51 in red), C1 (n=1119 in green) and C2 (n=2051 in blue), and (right) C0 (n=17 in red), C1 (n=82 in green) and C2 (n=94 in blue). The 8289_1#27 (c2) is highlighted in yellow. Similarity across pEK204 (93 732 bp) was based on blast alignment. The bla CTX-M gene was at 8 Kb, bla TEM was at 13 Kb, followed by a mix of conjugation and segregation genes at 36–70 Kb. This suggested independent integrations and rearrangements of pEK204-like plasmids across the subclades. The annotation was modified from [59].

To resolve the origin of pEK204-related pil HGT, we searched for the 14 pil operon genes (pilIJKMNOPQRSTUVI) individually in the 4071 assemblies [61]. The entire pil operon was conserved in 376 (9 %), including 61 from A (out of 414, 15 %), 97 from B (out of 433, 22 %), 17 from C0 (out of 52, 33 %), 95 from C1 (out of 1,121, 8.5 %) and 106 from C2 (out of 2051, 5.2 %) (Table 3) (see Data Access). Linked with the pEK204 matches, this indicated potential non-pEK204 pil ancestral acquisition in Clades A and B, and the putative Clade C ancestor may have had a pEK204-like plasmid with no pil that subsequently gained by recombination in some isolates.

Table 3.

Fim, pil, pap and ucl operon genes detected in 4071 ST131 genome assembles. The numbers and percentages of isolates per subclade with complete or absent (-) pil (pilK-pilV) or pap or ucl (uclABCD) operon genes where a subset also had three pap genes (papIBA). This focused on isolates allocated to the eleven main allelic options, isolates with some level of partial operon presence or absence were too sparse to be informative. UclBCD but not uclA was in an additional 2 % (n=77), including seven from A (2 %), 28 from B (6 %), none from C0, 24 from C1 (2 %), and 18 from C2 (1 %)

Operon			Subclade
pil	pap	ucl	A	%	B	%	C0	%	C1	%	C2	%	Totals
–	papIBA	–	279	67 %	107	25 %	3	6 %	772	69 %	1192	58 %	2353
–	all	all	21	5 %	13	3 %	23	44 %	82	7 %	481	24 %	620
all	papIBA	–	57	14 %	2	<1 %			75	7 %	66	3 %	200
–	all	–	1	<1 %	94	22 %			8	<1 %	14	<1 %	117
all	–	–	3	<1 %	71	16 %			5	<1 %	3	<1 %	82
all	all	all			16	4 %	1	2 %	10	<1 %	33	2 %	60
all	all	–	1	<1 %	6	1 %	16	31 %			4	<1 %	27
all	–	all							1	<1 %			1
–	papIBA	all									1	<1 %	1
all	papIBA	all											0
–	–	all											0
Partial matches			52	13 %	124	29 %	9	17 %	168	15 %	257	13 %	610
	Totals		414		433		52		1121		2051		4071

%

%

%

%

%

–

–

–

–

–

–

–

–

–

–

–

–

–

We extended this to two other important fimbrial operons: pap (pyelonephritis-associated pili) and ucl (uroepithelial cell adhesin-like). All four ucl genes (uclABCD) were in 17.5 % (712) of isolates, including 21 from A (5 %), 30 from B (7 %), 25 from C0 (48 %), 94 from C1 (8 %) and 542 from C2 (26 %) (Table 3). All ucl-positive C2 genomes (bar one sample) had a complete pap operon. The pap operon was in a minority (26 %, 1049) of isolates, including 28 from A (7 %), 163 from B (38 %), 49 from C0 (94 %), 127 from C1 (11 %) and 682 from C2 (33 %) (Table 3). Most isolates had papIBA but not the genes 3′ of them (papHCDJKEF) (65 %, n=2661): 359 from A (87 %), 118 from B (27 %), 3 from C0 (6 %), 888 from C1 (79 %) and 1293 from C2 (63 %). PapI (234 bp) is associated with higher P fimbriae expression in CFT073 [100] and here was amplified to two, three or four copies in 22 % ST131 : 26 in Clade A (6 %), 51 in B (12 %), 45 in C0 (87 %), 119 in C1 (11 %) and 658 in C2 (32 %). Only 45 assemblies had an amplified papB, and 265 had an amplified papA. These two genes (papBA) share a promoter pBA, and papI has a separate one (pI).

Most isolates from A, C1 and C2 and a minority of B had no pil nor ucl operons and lost the papH-F segment (Table 3). Clade B differed from A because some from this diverse clade have either pil or pap, but not both. The 83 ST131 with no detected pap genes all had a pil operon, and most (71 of 83) were from Clade B. Given the fim operon’s key role in ST131 evolution, we verified that fim was intact in >99 % of the 4071. As expected, the ISEc55 insertion at fimB [101] that delays fim expression in Clade C (like EC958 [20]) was in >99 % of C1 (1117, 99.6 %) and C2 (2041, 99.5%), but rarer in A (236, 57 %), B (93, 21 %) and C0 (18, 35 %) (32 samples had additional fimB rearrangements). During th e period prior to the divergence of C0 from C1 and C2 [50], the Clade C ancestor likely had a complete pap operon (but not pil nor ucl). The C1/C2 ancestor gained the fimB ISEc55 insertion and a minority of C1 gained either pil (80, 7 %) or ucl (82, 7 %) or both (11, 1 %), whereas 24 % of C2 gained ucl (482), 4 % pil alone (73) and 2 % both ucl and pil (33).

Unlike ucl gene products that are functionally independent [102], regulatory protein PapB reduces fim operon expression by inhibiting FimB and activating FimE, both tyrosine site-specific recombinases that invert fimS, including in CFT073 and 536 [103-105]. PapI and PapB regulate pap expression [104-106] depending on their concentrations and protein binding at the 416 bp regulatory region between papI and papB. Thus, retention of papIBA in >65 % of ST131 could be linked to fim transcriptional regulation by PapB. Using the STRING PPI network data above for CFT073 and 536 that have fim, pap and ucl [107, 108] (no pil gene products were found), there were ten fim-pap inter-operon PPIs in both 536 and CFT073, including pilus rod subunit PapA with FimD and FimF (as well as FimC), potentially matching the PPIs of FimA [109]. If only papIBA was present, the functional effects (if any) of papI and papA remain uncertain because it is unclear if FimC (replacing PapD at the inner cell membrane) could mediate periplasmic transport of PapA (like FimA) subunits via the chaperone FimI (instead of PapH) to usher FimD (rather than PapC) for pilus rod assembly at the outer cell membrane, extended later by FimF (in lieu of PapK) for the base of the tip. Nonetheless, presence of these genes can inform pilicide and coilicide design, such as antibodies targeting VF PapA [109].

Discussion

Commensal and environmental bacteria are major reservoirs of AMR genes [110], which are driven by HGT and recombination events [111], including in [112-115] and microbiome species [116, 117]. Previous work screened contigs from preterm infants for resistance to 16 antibiotics [1]. Here, we showed that this resistome was shared extensively between ExPEC ST131, commensal ST131 and microbiome , indicating likely transfer of these genes across commensal and pathogenic bacteria inhabiting the human gut and urinary tract, as expected [118].

83 972 does not express functional fimbrial adhesins and is used for therapeutic urinary bladder colonisation in patients in which it protects against super-infection [119-121]. Although 83 972 lost virulence during its adaptation to commensalism [122], here it had AMR genes like PBP3 that it shared with an gut microbiome (3_2_53FAA). This retention of certain AMR genes in asymptomatic specimens is important when assessing the AMR gene evolution in ST131 [123, 124].

Within ST131, Clade C core genomes are highly conserved but accessory genomes have extensive differences in AMR gene content [50, 52, 61]. Here, this was supported by NCTC13441’s repertoire of bla TEM and bla OXY genes, contrasting with EC958’s variety of bla CMY ones. This reinforced the view that ST131’s accessory resistome is shaped by the environment mediated by plasmids, rather than population structure or geographic factors [125], so perhaps tracking plasmids and MGEs in addition to AMR genes [42, 126] could assist treatment diagnostics [127].

By combining seven Clade C genome assemblies, then these with 4064 Clade C assemblies, we found that most (78 %) ST131 had >10 Kb of regions similar to plasmid pEK499 and some (9 %) had regions with high similarity to IncI1 plasmid pEK204, suggesting either discrete gains of pEK204-like plasmids in each C subclade or its presence in the Clade C ancestral lineage. IncI1 plasmids are common in ExPEC and are associated with different ESBL genes [67, 128–130]. Given that backbone plasmid genes may determine fitness effects more than ESBL genes [131] and long-term IncF plasmid persistence in ST131 [62], plasmids pEK204, (pEK516,) and pEK499 proteins’ higher interaction rates with chromosomally-encoded proteins relative to other relevant plasmids could be due to co-evolution. Measuring the indirect connectivity as PPI network loops per interaction showed evidence of this long-term retention, and may help identify plasmids compatible with chromosomes.

ST131’s fitness advantage is tightly correlated with type 1 fimbriae variant fimH30 and delayed fim operon expression due to an insertion at fimB [17, 132, 133]. Here, a minority of ST131 had type IV pilus biosynthesis (pil) genes and most Clade C with pil had >40 Kb of regions similar to pEK204, whereas Clade A and B isolates with pil did not have pEK204-like regions. Pil may allow different epithelial cell adhesion from fim via a thinner pilus and biofilm formation [90]. Additionally, F17-like ucl operon was more common in C2 (24 %): UclD (like FimH) is a two-domain tip adhesin that binds intestinal epithelial cells via O-glycans (FimH uses N-glycans) [102]. The ucl operon is associated with uroepithelial cell adhesion, biofilm formation and could be mobilised by its flanking MGEs [134].

The pap operon encodes a P fimbriae with high specificity for kidney epithelial cell and erythrocyte receptor glycolipids (α-d-galactopyranosol(1-4)-α-d-galactopyranoside) [135]. Delayed fim expression in Clade C [20, 101] may be reduced further by PapB: papIBA alone was retained in >65 % of Clade C here. Delayed fim expression [136] and having multiple P fimbrial gene clusters [106] associates more with pyleonephritis than cystitis. Given that express single fimbriae at the cell surface [103, 137] and isogenic cell populations can express distinct fimbrial types [138], type 1 fimbriae may allow bladder colonisation followed by reduced expression in favour of P fimbriae when at the kidney [137]. Our results on plasmid-linked changes across the ST131 resistome and diversity of fimbrial gene composition could inform on potential infection mechanisms [139] and biofilm formation [20].

Methods

E. coli genome isolate collection

Of the 12 assemblies examined (Table 1), three were ST131 references: SE15, NCTC13441 and EC958. Seven were ST131 from Ireland in 2005–2010, of which two from C1 were bla CTX-M-14-positive, as was 8289_1#24 from C2, and all five from C2 were bla CTX-M-15-positive. All seven were from urine except 8289_1#34, which was a rectal swab. The two HMP samples were: 3_2_53FAA (aka EC3_2_53FAA) from a 52-year-old male Canadian with Crohn’s disease’s colon biopsy [140]. The 83 972 (aka EC83972) sample was from the urine of a Swedish girl with a 3 year history of asymptomatic bacteriuria and stable renal function [141]. The 83972’s have common ancestry with virulent pyelonephritis-causing ExPEC CFT073 [142]. ST131 C2 reference genomes NCTC13441 and EC958 were from UTIs in the UK, and had bla CTX-M-15-positive plasmids [10, 59, 132]. NCTC13441 has 4983 predicted protein-coding genes and EC958 has 4982. EC958 has numerous virulence-associated genes that encode adhesins, autotransporter proteins and siderophore receptors, and can cause impairment of uterine contractility in mice [143]. Plasmid pEC958A has 85 % similarity with pEK499 but lacks the latter’s second tra region due to an IS26-mediated bla TEM-1 insertion [10, 60, 132]. ST131 Clade A reference SE15 was examined as a genetic outgroup and a commensal control because it lacks many virulence-associated genes. It has a 4 717 338 bp chromosome with 4338 protein-coding genes and a 122 Kb plasmid pSE15 with 150 protein-coding genes [87].

Illumina library quality control and read mapping

The paired-end Illumina HiSeq libraries for each sample were screened for low quality (Phred score <30) and short (<50 bp) reads using Trimmomatic v0.36 [144] and corrected using BayesHammer from SPAdes v3.9 [144] (Table S7). These corrected read libraries were mapped to references with SMALT v7.6 (www.sanger.ac.uk/resources/software/smalt/), and the SAM files were converted to BAM format, sorted and PCR duplicates removed using SAMtools v1.19 [145].

Homology-based resistome screening and comparison

The read mapping to the reference resistome (Table S1, see Methods of [1]) used GROOT [77] where the reads were indexed using the median read length (Table S7). Contig and protein domain annotation was derived from the Pfam v27.0 and ProSite databases using InterProScan v5.22–61 [146]. The protein homolog dataset in the CARD (2239 genes) [76] was aligned with the genomes to annotate the resistomes with blast v2.2.31, where matches with a bit score >500 and >99 % homology were considered valid. Alignments were visualised using Artemis and the Artemis Comparison Tool (ACT) [147], and also R packages VennDiagram v1.6.1, Seqinr v3.4–5, UpSetR v1.4.0 and WriteXLS v5.0.0.

Assessment of plasmids prevalent in ST131

Sequence and annotation files for pEK499 (NC_013122.1, EU935739), pEK516 (NC_013121.1, EU935738), pEK204 (NC_013120.1, EU935740) and all pV130 contigs (LC056314.1 to LC056328.1) [93] were aligned with Clustal Omega v1.2.1 [148]. Replicon typing used PlasmidFinder [56] and each plasmid was compared to the CARD. Each read library was mapped to each plasmid to verify local genetic features and quantify copy number levels, visualised with Artemis and R v3.4.2’s Reshape2 v1.4.3, Ggridges v0.5.1, Ggplot2 v2.3.2.1, Readr v1.3.1, Dplyr v0.8.3 and Ape v5.3 packages. Sequence similarity was calculated using the Sequence Identity and Similarity (SIAS) tool (http://imed.med.ucm.es/Tools/sias.html). Homology searches against pEK204 and pEK499 in the 4071 genomes [60] examined matches >300 bp length with >98 % sequence identity.

Protein-protein interaction network construction and topological data analysis

TDA has been used to investigate complex and high-dimensional datasets, like breast cancer genomes [149]. We applied it to PPI networks with a focus on direct and indirect connectivity across different network topologies by quantifying rates of absent PPIs (non-trivial loops). The numbers of non-trivial loops per PPI was used as a metric for indirect connectivity because it was consistent across parameters and was independent of network size. We assessed pCA14, pEK204, pEK499, pEK516 and pEC958A, along with two plasmids with no known AMR genes as negative controls: pECSF1 and pJIE186-2 [99]. We examined unique genes with PPI network information only. Results for pEK516 were identical to those for pEK499.

We extracted K12 MG1655 PPI data from the STRING database v10 [97] with R v3.5.2 packages BiocManager v1.30.4, Dplyr v0.8.0.1, Genbankr v1.10.0, Rentrez v1.2.1, STRINGdb v1.22.0, tidyverse v1.2.1 and VennDiagram v1.6.20. STRING’s combined scores were used because they integrate multiple types of evidence while controlling for random PPIs [97]. Using a score threshold of 400, K12 had 4146 protein-coding genes, of which 4121 had interactions, resulting in 105 457 PPIs - 14 555 PPIs were present for a score threshold of 900. For each plasmid or set of genes, we obtained the unique genes and numbers of pairwise PPIs within that set and with the chromosome. We tested our TDA-based approach using previous work [98]. For 536, FooB was used as an equivalent for PapB, and likewise for F7-2 that matched VF PapA in CFT073, and FooG for PapG.

We constructed the Vietoris-Rips complex [150] where the proteins were the vertices with PPIs as edges, so that proteins A and B with a PPI would be joined by a single edge (two proteins in one dimension, 1-D). Proteins A, B and C joined by three PPIs would have a filled 2-D triangle. For four proteins A, B, C and D joined by all six possible pairwise PPIs, their four filled 2-D triangles A/B/C, B/C/D, A/C/D and A/B/D constitute a tetrahedral surface filled with a 3-D tetrahedron. This can be extended to m+1 pairwise connected proteins, which get a m-D polytope filled into their skeleton of edges, triangles, tetrahedra, and so on. For m ≥4, a loop connects m proteins by m PPIs, such that each protein is involved in precisely two of these PPIs. For instance, four proteins A, B, C, D can be joined to a loop by four PPIs A-B, B-C, C-D, D-A, where A-B is a PPI between the proteins A and B (etc). Any loop that can be filled with triangles, along existing PPIs, is called trivial. So if there is an interaction between A and C, the loop is filled by the two triangles AB-BC-CA and CD-DA-CA and hence trivial. If there are no interactions between A and C, and none between B and D, then the loop is non-trivial.

Specialised software [151] was used to optimise the Betti number computations instead of using general purpose implementations of the Vietoris-Rips complex (such as SAGE) due to the large K12 dataset size (4146 proteins with 105 457 interactions). This software used the sparsity of the boundary matrices to process the rank computations efficiently in LinBox [152], with a hard-coded dimensional truncation of the Vietoris-Rips complex to avoid the large number of high-dimension simplices that would be obtained for the full dataset. For each analysis across STRING combined scores 400 to 900 with a step of 25, we computed the first Betti number of the Vietoris-Rips complex that counted the numbers of non-trivial loops (with missing PPIs inside), which was adjusted for the numbers of PPIs above the score threshold (i.e. loops per PPI). The number of PPIs for the complete K12 chromosomal dataset was negatively correlated with the combined score threshold (r2=0.964), as was the loops per PPI (r2=0.905), so the chromosomal loops per PPI was used as a baseline for different score thresholds.

Operon gene homology search approach

Homology searches for each fim operon gene were implemented using the NCTC13441 annotation coordinates to extract the corresponding sequence with SAMtools v1.9 and align each gene to the 4071 ST131 assemblies [60] with blast v2.2.31, processed with R packages Tidyverse v1.2.1 and Dplyr v0.8.3. For all operons, minor individual gene partial matches or losses were not examined due to the large number of samples and lack of consistent patterns for rare combinations. Results for amplified genes were restricted to those with prevalence >1 %. The CDS of fimB is 600 bp, and with the ISEc55 insertion the region typically spanned 2493 bp. NCTC13441 has only three pap genes (papIBA) and lacks any ucl or pil genes, so the pil genes were from pEK204. PilI and pilJ encoding IncI1 conjugal transfer proteins were mainly absent with no clear association with the other pil genes, and so were not examined here. PilV encoding a pilus tip adhesin had >1 copy in 292 assemblies, as expected for a locus rearranged to change pilus binding specificity [153]. The sequences for the four genes in the 5 Kb ucl operon were determined from 83 972 (CP001671): major subunit (uclA), chaperone (uclB), usher (uclC) and adhesin (uclD). UTI89’s chromosome (CP000243) and plasmid (CP000244) [111] were used to get the pap genes.

Click here for additional data file.