Evolutionary responses of Escherichia coli to phage pressure: insights into mucoidy and colanic acid overexpression

Abstract

Background

Antibiotic resistance is a major issue affecting all spheres of human activity, including agriculture. One significant example is the Avian Pathogenic Escherichia coli (APEC), a bacterium that infects poultry and leads to substantial economic losses in the farming industry. As antibiotics lose efficacity, bacteriophages (phages) —viruses that specifically target bacteria—are emerging as a promising alternative to antibiotics for treating and preventing bacterial infections. However, bacteria can develop resistance to phages through various mechanisms. Studying the coevolution between a phage and its host bacterium is important to gain insight into the phage's potential as a therapeutic agent. This study investigates the evolutionary responses of an APEC strain and a laboratory E. coli strain to a commercial phage originally isolated from APEC.

Results

In most cases, phage resistance resulted in a significant increase in mucoidy. Genomic analysis revealed that this resistance consistently correlated with amino acid changes, particularly in proteins involved in colanic acid production, such as YrfF. Further investigation of a mutation found in the YrfF protein demonstrated that this mutation altered the protein’s structure and its interaction with the membrane. Transcriptomic analysis confirmed that the genes involved in colanic acid production were significantly overexpressed. Although the strains possessed a CRISPR-Cas system, it did not contribute to phage resistance.

Conclusions

This study suggests that specific amino acid changes in key proteins may be a mechanism employed by E. coli, including APEC, to defend against phage infections.

Introduction

Bacteria play an important role in maintaining the health of the digestive system, including preventing colonization by pathogens in the gastrointestinal tract by competitive exclusion [1]. In chickens, the gut microbiota is particularly abundant and complex in the ceca, with approximately 10¹⁰ colony-forming units (CFU) per gram of digesta and approximately 1000 different species [2]. Escherichia coli is a commensal bacterium of the gut in warm-blooded animals, including chickens [1, 3,4,5]. In newly hatched chicks, E. coli becomes the first colonizer and rapidly dominates the ceca during the first week of life [2, 3].

While most E. coli strains are harmless, a specific subset of extraintestinal pathogenic E. coli, known as Avian Pathogenic E. coli (APEC), colonize the gut of chickens asymptomatically [6,7,8] and can trigger invasive infections outside the gastrointestinal tract of avian species. This pathogen is prevalent across all stages of the chicken life cycle [9] and is the etiologic agent of avian colibacillosis, a disease that can increase mortality by up to 53.5% in young chickens [10]. APEC can also lead to increased morbidity and carcass condemnation rates at slaughter [3, 11]. The disease typically begins with the infection of the respiratory tract, likely following the inhalation of contaminated dust particles [3, 7], and rapidly progresses to septicemia [3, 7, 12, 13].

Previous extensive use of antibiotics in livestock farming as a means of disease prevention and as a growth factor has caused the emergence of bacteria that are multi-resistant to antibiotics [14, 15]. Poultry producers now require other means of controlling pathogenic bacteria, such as APEC, to ensure the profitability of their production and animal welfare. The current threat of antibiotic-resistant bacteria has spurred further research into antibiotic alternatives including bacteriophages (phages)—viruses infecting only bacteria—particularly in Western countries [16, 17]. Various studies have shown that phages could be an interesting alternative to antibiotics in poultry farming [5, 18,19,20,21,22,23].

However, the use of phages presents several challenges that can undermine the effectiveness of the treatment, including the development of bacterial resistance to phages [24,25,26,27]. Bacteria can become resistant through numerous mechanisms, as reviewed by several authors [25, 27,28,29], and are referred to as “bacteriophage insensitive mutants” (BIMs). Phage resistance mechanisms include prevention of phage adsorption (production of extracellular matrix, production of competitive inhibitors, and blocking phage receptors), prevention of phage DNA entry (superinfection exclusion systems), degradation of phage nucleic acids that enter the bacterial cell (restriction-modification systems, prokaryotic argonautes, CRISPR-Cas systems), spontaneous mutation of the receptor recognized by the phage, and abortive infection systems [16, 25, 28,29,30,31]. Also, phages tend to have a narrow lytic spectrum, often infecting only a few strains of a bacterial species [32,33,34,35]. Studying the coevolution between a phage and its host bacterium is important to provide information for the eventual therapeutic use of phages in animal husbandry.

The present study examined the interaction between a commercial phage, phage 66, and two bacterial strains, APEC and the laboratory strain E. coli K-12, both susceptible to infection by this phage. BIMs were isolated and characterized, and amino acid changes in proteins were correlated with phage resistance. Further investigation of the consequences of a specific amino acid change in YrfF through molecular dynamics (MD) and transcriptomic analyses demonstrated that genes involved in colanic acid (CA) production play a decisive role in phage protection. This study provides a better understanding of the evolutionary dynamics between E. coli and its phages.

Material & methods

Wild-type bacterial strains, bacteriophage, and growth medium

SyntBioLab Inc. (Lévis, Canada) supplied phage 66 (cat # VR- 21.0056) in lysogeny broth (LB) medium and its host bacterium, E. coli APEC17, as a frozen glycerol stock solution (LB medium, 20% v/v). The APEC strain was originally provided by the Quebec Department of Agriculture, Fisheries, and Food (MAPAQ), while phage 66 was isolated by SyntBioLab Inc. from a poultry farm, both located in the province of Quebec, Canada. The bacterium E. coli K-12 MG1655, also frozen in a cryotube, was available at the Département des sciences animales at Université Laval (Quebec, Canada).

Throughout the study, bacteria were grown at 37 °C in liquid heart infusion broth (HIB) medium (Oxoid Inc, Ottawa, ON, Canada, cat # CM1032B) at 200 rotations per minute (RPM) or on HIB agar plates (1.5% agar). New glycerol stocks (15%–20% glycerol) of these strains were prepared in fresh HIB medium and used for the rest of the study.

The phage 66 was amplified using the host strain APEC17 in fresh HIB medium by adding 100 µL phage to 100 µL exponential phase bacteria in 4 mL fresh media, followed by incubation at 37 °C for 3.5 h. After incubation, the culture was centrifuged at 3,220 g and the lysate was filtered through a 0.22 µm filter. The amplification cycle was repeated three times before it was titered to obtain the plaque-forming units (PFU) per mL.

The final lysate was prepared for phage titration following the method described by Leduc et al. (2021) [36], with minor modifications. Using the double layer agar assay, 100 µL serial dilutions of phage 66 (10⁻⁵ to 10⁻⁸) in phage buffer 1X (50 mM Tris–HCl pH 7.5, 100 mM NaCl, 8 mM MgSO₄) were mixed with 100 µL exponential phase subculture of strain APEC17, plated on HIB agar, and incubated O/N at 37 °C.

One-step growth curve

A one-step growth curve experiment was performed to determine the burst size and latent period of phage 66, with modifications to a previously described method [37]. Briefly, 2 mL of an exponential phase APEC17 subculture (optical density [OD]_{600 nm} 0.6–0.8) was centrifuged at 17,000 g for 1 min at room temperature (RT). The pellet was resuspended in 900 µL fresh HIB medium, and 100 µL phage 66 (10⁸ PFU/mL) was added to achieve a multiplicity of infection (MOI) of 0.01, assuming that an OD600 of 1.0 corresponds to approximately 8 × 10⁸ cells [38,39,40]. The mixture was incubated at 37 °C in a mini tube rotator for 3 min to allow phages to adsorb to the host bacteria. Next, the culture was centrifuged at 17,000 g for 1 min to remove unadsorbed phage particles, and the pellet was washed twice with 1 mL fresh HIB. The mixture was serially diluted to 10⁻⁴ by putting 100 µL of the 10⁻² dilution into 9.9 mL of fresh HIB. To determine the initial phage titer and the time zero, 350 µL of the mixture was aliquoted into two 1.5 mL microtubes. The first microtube was plated immediately in triplicate (100 µL mixture + 100 µL bacterial host in exponential phase) to provide the initial phage titer. The second microtube was centrifuged at 17,000 g for 1 min, and then the supernatant was plated in triplicate (100 µL mixture + 100 µL bacterial host) to determine the time zero. Afterward, 100 µL samples were taken every 5 min for 30 min, centrifuged, and the supernatant was plated. Burst size was calculated using Eq. 1. The latent period was estimated at the midpoint of the exponential curve [41].

$$\text{Burst size }=\frac{{\text{Titer}}_{\text{final}}-{\text{Titer}}_{\text{initial}}}{{\text{Titer}}_{\text{initial}}}$$

(1)

Efficiency of plating assay

The efficiency of plating of phage 66 on strain K-12, compared to strain APEC17, was determined. To achieve this, a phage 66 titration using double agar overlay plaque assay was performed in triplicate on strains K-12 and APEC17, as previously described, and the efficiency of plating was calculated using Eq. 2.

$$\text{Efficiency of plating }=\frac{{\text{Titer}}_{\text{K}-12}}{{\text{Titer}}_{\text{APEC}17}}$$

(2)

Transmission electron microscopy

Phage 66 was prepared for transmission electron microscopy (TEM) as previously described, with a few modifications [42]. Briefly, 8 mL filtered (0.22 µm) phage lysate (10⁸ to 10¹⁰ PFU/mL) was ultracentrifuged for 1 h at 4 °C at 30,000 g. After removing the supernatant to leave 250–500 µL at the bottom, 1.5 mL of 0.1 M ammonium acetate (Thermo Fisher Scientific, Waltham, Massachusetts, USA; 0.1 M, pH 7, cat # AM9070G) was added, and the sample was centrifuged again. This step was repeated once more, and after removing the supernatant, the remaining 50–100 µL was transferred to a clean 1.5 mL microtube and placed at 4 °C until TEM visualization.

Two microliters of 2% uranyl acetate (Electron Microscopy Sciences, Hatfield, Pennsylvania, USA, cat # 22,400–2) was deposited on a glow-discharged copper with carbon film grid (Electron Microscopy Sciences, Hatfield, Pennsylvania, USA, cat # CF200H-Cu- 50) and then mixed with 2 µL of purified phage 66. The liquid was removed after 30 s by softly touching the edge of the grid with blotting paper. The grid was left to dry for at least 24 h at RT before TEM visualization.

Phage 66 morphology was observed using a JEOL JEM- 2100 + transmission electron microscope at 200 kV accelerating voltage (JEOL USA, Peabody, Massachusetts, USA), and pictures were taken with a Gatan OneView camera (AMETEK, Berwyn, Pennsylvania, USA). ImageJ software [43] complemented with Fiji image processing package [44] was used to measure the dimensions of the phages from the means of at least 30 specimens.

Growth kinetics

A subculture of each strain was incubated at 37 °C, 200 RPM until OD > 0.5 (WPA CO 8000 Biowave Cell Density Meter, Biochrom Ltd., Cambridge, United Kingdom). Phage 66 was added at an MOI of 1, or not added to the incubation.

Growth curves were generated using the BioTek EPOCH 2 microplate spectrophotometer (Agilent, Santa Clara, California, USA) with BioTek Gen5 Microplate Reader and Imager Software version 3.12 (Agilent, Santa Clara, California, USA) and 96-well transparent flat-bottom plates (Greiner Bio-One, Monroe, North Carolina, USA, cat # 655,161). The OD was measured at 600 nm every 20 min for 16 h at 37 °C with double-orbital shaking at 200 RPM. The experiment was repeated three times.

Isolation and selection of BIMs

The APEC17 and E. coli K-12 BIMs were isolated using a double agar overlay method, following an adapted protocol [45]. Bacterial subculture (200 µL) in exponential phase (~ 1 × 10⁸ CFU/mL) and 100 µL phage 66 (~ 1 × 10⁹ to 1 × 10¹⁰ PFU/mL) were added to 3 mL soft HIB agar (0.75% agar) kept at 55 °C before being poured uniformly on a thin HIB agar plate. After solidification, plates were incubated O/N at 37 °C. All BIM colonies were collected, inoculated into 8 mL of HIB medium, and incubated for 3.5 h. To confirm resistance, BIMs were subjected to a series of phage 66 amplifications with increasing MOIs (10, 100, and 200), each time using the BIM culture from the previous amplification and incubating in 4 mL fresh HIB medium. If there was growth in the tube from the last amplification, this culture was streaked onto HIB agar and incubated O/N. At this step, if two distinct phenotypes (e.g., mucoid and non-mucoid) were present on the solid medium, both were picked and treated as separate BIMs, and a lower case “a” was added to the identification name given to the non-mucoid BIM. If no growth was observed on the solid medium, the BIM isolation process was repeated from the beginning. An isolated BIM colony was re-streaked two more times, and a spot test using phage 66 was performed to ensure the BIM’s resistance. Glycerol stocks (15%–20% glycerol) of the isolates were made and stored at − 80 °C. Finally, a last spot test was performed using the glycerol stock to confirm the BIM’s resistance to phage 66.

Phage adsorption rate assay

Adsorption rate experiments of phage 66 on strains K-12 and K12-BIM1-M were conducted following a previously described procedure [46], with minor modifications. Briefly, at each time point (5, 10, 15, 20, 30, 40, 50, and 60 min), 100 µL of phage 66 was mixed with 900 µL of bacterial culture at an MOI of 0.01 in a tube under aerobic conditions and incubated at 37 °C with shaking at 200 RPM. The bacterial culture was prepared by centrifuging 2 mL of a culture at OD600 0.6–0.8, followed by resuspension of the pellet in 900 µL of fresh HIB. After incubation, the tube contents were centrifuged to remove adsorbed phages, and the titers of unadsorbed phages in the supernatant were determined using a double agar overlay plaque assay, plated in triplicate. The percentage of phage adsorption at each time point was calculated using Eq. 3.

$$\text{Adsorption rate }(\text{\%}) =\frac{{\text{Titer}}_{\text{initial}}-{\text{Titer}}_{\text{surpernanant}}}{{\text{Titer}}_{\text{initial}}}\times 100$$

(3)

DNA extraction, sequencing, and bioinformatics analysis

Bacterial DNA was extracted using the QIAGEN’s QIAamp PowerFecal Pro DNA kit (Toronto, Ontario, Canada, cat # 51,804) following the manufacturer’s instructions. Approximately half of the bacteria from a streaked agar plate was used, and DNA was extracted into a final volume of 75 µL. The concentration of extracted DNA was measured using PicoGreen (Invitrogen, Waltham, Massachusetts, USA, cat # P7589). DNA samples were stored at − 20 °C until they were sent for sequencing.

The DNA samples of the parental strains were sequenced by Oxford Nanopore. Libraries were prepared with 1.5 µg total gDNA using ligation sequencing gDNA + native barcoding (SQK-LSK_109 with expansion 1–12) following the manufacturer’s instructions. An equal amount of each sample was mixed during the final step of library preparation. The libraries were sequenced on a R9.1 flow cell on a GridIon instrument following the manufacturer’s instructions. Parental and BIM DNA samples were sequenced on a MiSeq apparatus (3 × 300 bp). The Illumina DNA libraries were prepared using the NEBNext UltraII DNA Library prep kit (NEB, Whitby, Ontario, Canada, cat # E7645L). All the sequencing was performed by the Plateforme d’analyse génomique de l’Université Laval.

The raw Illumina sequencing reads of the parental strains (APEC17 and K-12) were filtered with fastp version 0.23.2 [47], while those from Nanopore were filtered using Filtlong version 0.2.1 (https://github.com/rrwick/Filtlong) by keeping the best 90% of reads above 1,000 bp or until only 500 Mbp remained. Both datasets were de novo assembled by Unicycler version 0.5.0 [48], and sequences were annotated using Prokka version 1.14.6 [49]. For BIM analysis, reads were filtered with fastp version 0.23.2 [47], then BRESEQ version 0.37.1 [50] was used to find differences between the sequences of the BIMs and their respective parental strains. The CRISPR-Cas systems were investigated with CRISPR Comparison Toolkit version 1.0.2 [51]. In silico serotyping was performed with ECTyper [52] version 1.0.0. The level of YrfF protein conservation was assessed by the ConfSurf web tool using the default settings [53], except that the maximum number of sequences included in the analysis was changed to 500. Illumina and Nanopore sequencing reads were deposited in the Sequence Read Archive database under the BioProject accession PRJNA1179319, and the assembled genomic sequences of the K-12 and APEC17 strains were also deposited in GenBank under the accession numbers GCA_047038235.1 and GCA_047037185.1, respectively.

Minimum inhibitory concentration

To determine the minimum inhibitory concentration (MIC) of tetracycline, E. coli K-12 wild-type (WT) and K12-BIM4-M subcultures were incubated at 37 °C, 200 RPM to reach an OD corresponding to mid-exponential phase (OD_600nm ~ 0.8). Next, a bacterial lawn was prepared on HIB agar using 1 mL of this subculture, and a tetracycline antibiotic strip (Liofilchem®, Waltham, MA, USA, cat # 22–777–902) was aseptically placed in the center of the plate once the bacterial lawn was dry. The agar plate was incubated O/N at 37 °C, and the MIC was determined following the manufacturer’s guidelines.

Molecular modeling

The WT YrfF protein (identified as a potential key player in phage resistance) and its L643R mutant were studied using MD simulations to gain insights into their structural and functional roles. The initial coordinates for the WT protein were obtained from AlphaFold 2 [54] (ID: AF-P45800-F1). All ionizable residues were assumed to be in their standard protonation states at pH 7. Positioning of Proteins in Membranes 2.0 [55] was used to predict the initial placement of the protein within the membrane. Explicit membrane and solvent systems were built using the CHARMM-GUI v.3.7 [56, 57] and the L643R mutation was introduced using the CHARMM-GUI tools. The bilayers consisted of 200 lipids per leaflet, with a composition similar to the inner membranes of gram-negative bacteria [58,59,60] (Additional file 1). The structures were immersed in neutrally charged rectangular boxes of water containing Na⁺ and Cl⁻ ions at a concentration of 150 mM, resulting in approximately 291,000 atoms per system. The system dimensions measured 116 × 116 × 225 ų.

The MD simulations were performed using NAMD 2.14 and NAMD 2.15alpha2 multicore CUDA [61] with the CHARMM36 all-atom potential energy parameters [62], including bacterial phospholipids [60, 63] and TIP3P water model [64]. Simulations were conducted at 310.15 K under isothermal-isobaric (NPT) ensemble conditions with a 2-fs time step and periodic boundary conditions. Langevin damping with a coefficient of 1 ps⁻¹ maintained constant temperature, while pressure was controlled by a Nosé–Hoover Langevin piston at 1 atm. Bond lengths involving hydrogen atoms were constrained using SETTLE [65] for water molecules and SHAKE [66] for all other molecules. Cutoffs for the short-range electrostatics and the Lennard–Jones interactions were set at 12 Å, with the latter smoothed via a switching function over the range of 10–12 Å. Long-range electrostatic interactions were calculated using the particle mesh Ewald method (PME) [67], with an interpolation order of 6 and a maximum Fast Fourier Transform (FFT) grid spacing of 1.0 Å, at each integration step. Nonbonded pair lists were updated every 10 steps, and coordinates were saved every 1 ns. For both systems (WT and L643R), three independent 300-ns MD trajectories were recorded, retaining the last 100 ns of each trajectory for analysis, for a total of 300 ns for each system. A two-sample t-test was performed to compare the angles of the major periplasmic domains. The angles were measured from 10 structures sampled over the last 100 ns of three independent trajectories per system (YrfF-WT and YrfF-L643R), resulting in a total of 30 angles for each system.

Wordom versions 0.22-rc3.i86 and 0.24. × 86–64 [68] were used to determine residue contacts with the membrane and water. Contacts were defined as any non-hydrogen atom within 4.5 Å of the target residue. PyMOL [69] version 2.5.0 was used for generating images of the systems. Biopython [70] and Scikit-learn [71] were used for extracting and performing linear regression on carbon alpha coordinates to measure the tilt angles of the major periplasmic domain.

Bacterial growth curves, RNA extraction, and sequencing

Bacterial growth curves were performed in duplicate for E. coli K-12 WT and K12-BIM1-M to accurately determine cell concentrations based on OD measurements. For this, 20 mL of subculture, adjusted to an OD_600nm of 0.1 using a concentrated bacterial subculture and fresh HIB medium, was incubated at 37 °C, 200 RPM in a sterile 50 mL Erlenmeyer flask for 4 h. Every 30 min, 0.5 mL was removed for OD measurements during the first five intervals, and 0.25 mL for the last four. Simultaneously, 0.1 mL was taken to perform a serial dilution using HIB medium, and the different dilutions were plated on HIB agar medium using a sterile cell spreader. The inoculated agar plates were incubated O/N at 37 °C. Agar plates with 30 to 300 CFU for each interval used to calculate the bacterial concentration corresponding to its OD measurement. The curves were compared with the CGGC web server [72].

RNA extraction and sequencing (RNA-seq) of E. coli K-12 WT and K12-BIM1-M were carried out. RNA from E. coli K-12 WT and K12-BIM1-M was extracted in triplicate from independent subcultures using QIAGEN’s QIAwave RNA Mini Kit (Toronto, Ontario, Canada, cat # 74,534) following the manufacturer’s instructions, and sent to Genome Quebec Centre of Expertise and Services (Montréal, Québec, Canada) for sequencing on an Illumina Novaseq 6000 apparatus. The sequencing reads were filtered with fastp version 0.23.4 [47] and mapped on the genome sequence of E. coli K-12 generated for this study using bowtie version 2.5.1 [73]. The reads on the genes were counted using featureCounts version 2.0.1 [74]. The differential expression analysis between the genes of the parental strain and those of the BIM strain was performed using DESeq2 version 1.40.1 [75] in R version 4.3.0. Clustering and functional enrichment were carried out using STRING web server version 12 [76]. Illumina sequencing reads were deposited in the Sequence Read Archive database under the BioProject accession PRJNA1179319.

Results

In this study, the relationship between a commercial phage (phage 66) and its bacterial host (E. coli APEC17) was investigated. Phage 66 exhibits an icosahedral head and a long, non-contractile tail, suggesting that it belongs to the class Caudoviricetes (Fig. 1). The head measured 65.00 ± 0.27 nm in diameter (n = 49), and the tail measured 159.0 ± 1.1 nm (n = 30) without the tail fibers. The burst size (number of virions produced) of phage 66 with its bacterial host E. coli APEC17 was 156 ± 7 PFUs per infected cell, and its latent period (time required per round of infection) was 26.0 ± 1.8 min.

Fig. 1

Morphology of phage 66. Transmission electron micrograph of phage 66. Sample was stained with 2% uranyl acetate. 60,000 × magnification, 200 kV accelerating voltage. Scale bar represents 200 nm

Screening of phage 66 against various bacteria revealed that, in addition to its host strain APEC17, this virus can also infect the laboratory strain E. coli K-12 (Additional file 2). A spot test analysis demonstrated that phage 66 has a better lytic capacity on APEC17 (Fig. 2A) than on K-12 (Fig. 2B). The phage 66 titration on K-12 (6.67 × 10⁷ ± 5.49 × 10⁶ PFU/mL) is significantly lower (P < 0.01) than on APEC17 (1.53 × 10⁹ ± 6.89 × 10⁷ PFU/mL). The efficiency of plating of phage 66 on K-12 is 0.040 ± 0.002, indicating that phage 66 forms lysis plaques on K-12 at only 4% of the level observed on APEC17. The growth kinetics of the bacteria in liquid culture in the presence or absence of phage 66 at a MOI of 1 were considerably different (Fig. 2C). The lytic effect of the phage took approximately 150 min to become apparent in the growth of the K-12 strain, whereas the bacterial population of the APEC17 host strain was reduced within minutes. In the absence of the phage, the bacterial growth rates were not significantly different (P > 0.05), suggesting that the observed differences in the presence of the phage were not a result of inherent growth defects in either strain. Another noteworthy observation was that the growth of strain APEC17 with phage 66 rebounded approximately 350 min after the start of infection, suggesting the selection of cells resistant to phage 66. Interestingly, the growth profiles of both APEC17 and K-12 were not significantly different between MOI 0.1 and 1.0 (P > 0.05) (Additional Figures File: Figure S1).

Fig. 2

Spot tests of different dilutions of phage 66 on (A) its host strain E. coli APEC17 and (B) E. coli K-12. C Growth kinetics of strains APEC17 and K-12 with or without phage 66 at an MOI of 1. Error bars represent the standard error of the mean for three replicates

The observed differences in the interaction between phage 66 and the APEC17 and K-12 E. coli strains, particularly in terms of the bacterial resistance to the phage, prompted a more detailed investigation into the protective mechanisms employed by these bacteria against this phage. Isolation of 7 and 8 BIMs for K-12 and APEC strains, respectively, revealed increased mucoidy in some resistant mutants (Additional figures file: Figures S2 and S3), and all the BIMs isolated from the K-12 parental strain exhibited increased mucoidy. In contrast, among the BIMs isolated from APEC17 host strain, only four out of eight BIMs had increased mucoidy. Interestingly, two non-mucoid BIMs (APEC17-BIM1a-NM and APEC17-BIM3a-NM) were co-isolated with two mucoid BIMs (APEC17-BIM1-M and APEC17-BIM3-M) at the final step of BIM isolation (see Materials and Methods section).

The molecular determinants involved in phage resistance were investigated by DNA sequencing and comparative genomics of BIMs and parental strains. DNA sequencing of the APEC17 and K-12 strains revealed two chromosomes of 5,215,601 and 4,654,366 bp. Two plasmids of 138,781 bp and 122,566 bp were also found in the APEC17 strain. In silico serotype prediction suggested that APEC17 would be serotype O78:H4 and K-12 would be O16:H48. The two genomes were each predicted to have an IE-type CRISPR-Cas system with two arrays.

Although each strain possessed two CRISPR-Cas systems, no alterations (such as spacer additions or mutations) were found in the analyzed BIMs. However, numerous mutations were identified in BIMs when compared with the parental strains (Fig. 3, Additional file 3). The mutations are all chromosomal, except for a single mutation in the traG gene (T899 K) found in the 138,781 bp plasmid of APEC17-BIM1-M, APEC17-BIM1a-NM, and APEC17-BIM2-M. While two genes were mutated in both strains (rcsC, yrfF), some were unique to K-12 (rhsC, tmcA, malF, astA) or APEC17 (sitB, sitC, radB, traG, rcsD, ycfZ, lptD, rhlE, yjbG, a gene encoding a putative host specificity protein, and a gene encoding a transcriptional antiterminator). Notably, K12-BIM4-M also had lost a cluster of 12 genes (a total of 8646 bp), some of which provided tetracycline resistance (Additional figures files: Figure S4). This alteration caused increased sensitivity to tetracycline (MIC = 1.5 µg/mL) compared with the parental strain K-12 (MIC ≥ 256 µg/mL), which is considered completely resistant according to the supplier (Liofilchem®). Interestingly, another deletion of 8 genes (a total of 14 607 bp) was identified in a non-mucoid BIM of APEC17, APEC17-BIM4-NM (Additional figures file: Figure S5). This deletion included numerous genes involved in the direct or indirect production of mucoidy, including CA (algC, manC, galF, wcaM, wcaK, wzxC, wcaJ).

Fig. 3

Mutations identified in BIMs isolated from K-12 (A) and APEC17 (B) strains. Non-synonymous, synonymous, and structural mutations (nucleotide additions or deletions) are represented by green, blue, and red bars, respectively. A mutation in the traG gene (T899 K) was also identified in the 138,781 bp plasmid of APEC17-BIM1-M, APEC17-BIM1a-NM, and APEC17-BIM2-M

Of the two genes mutated in both K-12 and APEC, rcsC was mainly found in APEC (4 out of 5 BIMs), while the yrfF gene was mutated twice in both strains. Because yrfF has a propensity to be mutated in both strains and likely reflects a common mechanism of phage resistance, one of the BIMs having a mutation in this gene, K12-BIM1-M, was selected for further investigation. This BIM was chosen because the L643R mutation in YrfF is the only one identified in its genome, making it likely to play a major role in resistance to phage 66.

Firstly, it was interesting to investigate whether K12-BIM1-M's resistance to phage 66 was due to a difference in its adsorption capacity (Additional figures file: Figure S6). During the first 40 min after infection, the adsorption level did not vary significantly (P > 0.05) between the parental K-12 strain (12.890 ± 2.245%) and K12-BIM1-M (18.940 ± 1.853%). However, after 40 min, the number of phages in the supernatant increased in the parental strain, indicating the completion of a lytic cycle, while it remained stable in K12-BIM1-M.

YrfF is a five-pass inner membrane protein known to play a role as a repressor of the regulator of capsule synthesis (Rcs) phosphorelay [77,78,79]. The impact of the L643R mutation on the YrfF protein, as observed in the K12-BIM1-M strain, was inferred from its predicted structure, MD simulations, and transcriptomics. An investigation of the level of conservation of position 643 in relation to homologous sequences available on the UNIREF-90 database revealed that leucine is largely dominant at 95.516%, followed by isoleucine at 2.242%, valine at 1.794%, and serine at 0.448% (Additional file 4). No sequences were detected with an arginine at position 643.

The L643R mutation is located at the interface between the periplasmic domain and the membrane (Fig. 4A). The MD simulations showed that replacing the hydrophobic leucine with the bulkier, hydrophilic arginine led to local protein reorganization in the mutant (Fig. 4B). Typically, R643 remains buried into the protein structure, forming hydrogen bonds with the backbones of neighboring residues L396 and G642 (Fig. 4B, right). The change from leucine to arginine at this position causes the β-sheet regions near the mutation to transition into random coils, increasing water contact for residues that normally interact with the membrane in the WT (Fig. 4B). As a result, the periplasmic domain (residues 358 to 652) in the L643R mutant loses its membrane anchoring, which alters its orientation relative to the membrane plane. In the WT protein, the periplasmic domain has a 20° greater inclination (Fig. 4A) than the L643R variant (p < 0.001), i.e., the WT periplasmic domain is oriented at 40.1 ± 3.4° relative to the bilayer normal, whereas the L643R mutant is at 19.4 ± 4.8° (see also Additional Figures file: Figure S7). Despite this orientation change, the transmembrane domain, composed of five α-helices, remained stable in both the WT and L643R mutant, with only minor displacements observed (Additional file 5).

Fig. 4

Structural analysis of the L643R mutation on the YrfF protein in E. coli K-12. A YrfF-WT (cyan, left) and YrfF-L643R (green, right) in the presence of the membrane. The mutated residue (Leu643 to Arg643) is highlighted in orange spheres in both structures. Red arrows represent the average angle of the periplasmic domain according to the normal of the bilayer (z-axis, black arrow). B Residues in YrfF forming hydrogen bonds with L643 (left) or R643 (right) are shown in orange. Residues in red (G389, V390, V392, D394, K630, I631, and F632) are predicted to lose membrane contact in YrfF-L643R but not in YrfF-WT

Because YrfF is known to be a regulator, transcriptomics was used to verify whether this role was altered. Bacterial growth curves for E. coli K-12 WT and K12-BIM1-M were created prior to the RNA extraction (Additional figures file: Figure S8). The number of CFU and the OD curves were not significantly different (P > 0.01), suggesting that the fitness of the parental strain and the mucoid BIM are quite comparable.

Using RNA-Seq, the differential expression of genes between the parental strain K-12 and K12-BIM1-M was investigated. A total of 476 genes were considered dysregulated in K12-BIM1-M compared with the parental strain (Additional file 6). Of these genes, 325 were down-regulated (log2 FC < − 1, adjusted p-value < 0.01) and 151 were up-regulated (log2 FC > 1, adjusted p-value < 0.01). Proteins with dysregulated genes could be grouped into several clusters (Fig. 5). Proteins with up-regulated genes (Fig. 5A) had an average of 8.15 links formed among them, while proteins with down-regulated genes (Fig. 5B) had an average of 16.2 links formed between them. Therefore, down-regulated genes tended to be twice as connected to each other as up-regulated genes. Verifying the enrichment for certain functional categories shows, among other things, that the 20 genes involved in CA metabolism were up-regulated (Fig. 5A). Among the clustered groupings of down-regulated genes, most were involved in flagellar production (and chemotaxis), xenobiotic transport and metal response, and genes with no known function (Fig. 5B).

Fig. 5

Grouping of proteins according to the STRING database with genes that were (A) up-regulated and (B) down-regulated. Some proteins are colored based on their membership in enriched functional categories. Only functional categories enriched for the main cluster-forming proteins are displayed. The cluster numbers, accompanied by their descriptions, the number of proteins in the network compared with the total number of proteins of this category found in the reference genome, the strength of clustering (log₁₀ observed/expected), and the False Discovery Rate value (FDR, corrected p-value using the Benjamini–Hochberg procedure) are shown

Discussion

The aim of this study was to investigate the ability of phage 66 to infect a strain of APEC and the laboratory strain K-12, and the mechanisms used by these two different strains of E. coli to protect themselves against phage 66. The one-step growth curve revealed that phage 66 has a high burst size (156 ± 7 PFUs per infected cell) and an average latent period (26.0 ± 1.8 min) compared with other APEC-targeting phages, which have burst sizes ranging from as low as 6 to as high as 318 PFUs per infected cell (median = 57, mean = 93 PFUs per infected cell) and latent periods ranging from 10 to 40 min [18, 46, 80]. Phages with a short latent period and large burst size are preferred for phage therapy and biocontrol applications because of their high lytic activity [37, 80, 81].

Although phage 66 can infect both bacterial strains, the infection dynamics vary between the two strains. Phage 66 drastically reduced the bacterial population of strain APEC17, while there was a much longer delay in the reduction of strain K-12. This suggests that phage 66 is better adapted to infect strain APEC17—the bacterium used for its initial isolation by SyntBioLab Inc.—than strain K-12, which might be considered as an alternative host. The efficiency of plating of phage 66 on K-12 is 0.04, indicating that phage 66 forms lysis plaques on K-12 at only 4% of the level observed on APEC17. According to Green et al., this value suggests that the infection capacity of phage 66 is moderate on K-12 compared to its infection on APEC17 [82]. Interestingly, the bacterial population of strain APEC17 increased again after a certain period, unlike that of strain K-12. This observation prompted the investigation of the resistance mechanisms used by the two bacterial strains to protect themselves against phage 66.

Phage 66 was exposed to E. coli strains APEC17 and K-12 to isolate and characterize BIMs. Mucoidy was observed in most of the isolated BIMs, representing a striking phenotypic difference from the parental strains. This mucoidy is known to be caused by an overproduction of CA, an extracellular polysaccharide (EPS) produced by many Enterobacteriaceae such as E. coli and Salmonella spp. [83,84,85,86,87]. In stressful conditions, such as dehydration and osmotic shock, the overproduction of the extracellular matrix creates a protective capsule around the bacterium that favors its survival [25, 83, 88,89,90]. The overproduction of the extracellular matrix can prevent phage adsorption by creating a physical barrier, making it difficult or even impossible for the phage to reach the bacterial receptor to infect the host [25, 85, 90,91,92]. Interestingly, two non-mucoid BIMs (APEC17-BIM1a-NM and APEC17-BIM3a-NM) were co-isolated with two mucoid BIMs (APEC17-BIM1-M and APEC17-BIM3-M) at the final step of BIM isolation. This suggests that these non-mucoid BIMs had equivalent fitness and resistance levels to the mucoid BIMs and that mucoidy is not the only mechanism that provides protection against phage 66; however, the protection strategies used by non-mucoid BIMs are not yet clear and warrant further investigation.

Sequencing of BIM DNA and comparison with sequences from parental strains revealed no alteration of CRISPR-Cas system, such as the addition of spacers. This suggests that CRISPR-Cas systems may be inactive or that the phage has a way of bypassing them [93]. A similar study with other phages is required to confirm this observation. However, several mutations in different genes were observed in the BIMs, with some mutations found only in K-12 BIMs (rhsC, tmcA, malF, astA) and some only in APEC17 BIMs (sitB, sitC, radB, traG, rcsD, ycfZ, lptD, rhlE, yjbG). Realistically, mutations in these genes may reflect the different molecular strategies used by different strains. It is also possible that some mutations, especially those observed only once, were co-selected with the resistance phenotype.

Interestingly, K12-BIM5a-M presents a single nucleotide polymorphism mutation in the tmcA gene, which codes for tRNA(Met) cytidine acetyltransferase TmcA. This protein decreases the DNA-binding ability of H-NS [94], a multifunctional protein that acts as a global transcriptional silencer of genes with high AT content [95, 96] and regulates approximately 5% of all E. coli genes [97]. H-NS represses rcsA expression, which is a gene involved in CA production [98]. Another interesting observation was the loss of a gene cluster in K12-BIM4-M, which results in the loss of a tetracycline resistance gene and decreases the MIC for this antibiotic. It is impossible to determine the impact of the loss of these genes on phage resistance. Interestingly, the non-mucoid BIM APEC17-BIM4-NM also had lost several genes, with the majority involved in mucoidy and CA production [99,100,101,102]. The absence of these genes may contribute to the strain being non-mucoid. For example, in Klebsiella pneumoniae, a frameshift mutation of wcaJ (one of the genes deleted in APEC17-BIM4-NM) has been shown to alter CA production and cause a loss of mucoidy [102].

Two genes were found to be mutated in both strains: rcsC and yrfF. Interestingly, these genes are both involved in the production of CA and, therefore, are realistically involved in the increased mucoidy. Previous studies have shown that a point mutation in the yrfF gene (also known as igaA in Salmonella) leads to an overproduction of capsule material and partial loss of motility in Salmonella mutant [103, 104]. This gene was identified in two large-scale studies and shown to confer protection against various phages in E. coli [91, 105]. Similarly, a truncated rcsC gene in a Salmonella mutant induced mucoidy [106], and rcsC mucoid Salmonella mutants were isolated by resistance to mecillinam [107]. Further supporting this, Majdalani and Gottesman (2005) reported an increase in the expression of capsule synthesis (cps) genes resulting from a rcsC point mutation in E. coli [108]. YrfF negatively regulates RcsC activity, although it is not clear whether this regulation occurs through direct or indirect interaction with RcsC or which specific domain of RcsC was affected [108]. As mentioned in the Results section, yrfF codes for the inner membrane protein YrfF (IgaA homolog) and serves as a negative regulator of the Rcs system, specifically by negatively regulating RcsC [78, 79, 86, 104, 108,109,110]. When RcsF, an outer membrane lipoprotein, detects stress on the cell surface (such as osmotic shock, temperature changes, or pH fluctuations), it signals for YrfF to relieve its inhibition [78]. This action allows the Rcs phosphorelay components—RcsC, RscD, RcsB, and RcsA—to become active and permit the transcription of over 150 genes—including those for biofilm formation and biosynthesis of CA production (cps and wca [111]) [112, 113]—and repress genes involved in flagellar motility [77,78,79, 108].

In this study, all isolated BIMs from both the K-12 and APEC17 strains that had a point mutation in yrfF exhibited a mucoid phenotype, which was not the case for BIMs that had a point mutation in rcsC. Indeed, while a single point mutation in the rcsC gene of K12-BIM2-M led to a mucoid phenotype, some BIMs without increased mucoidy (APEC17-BIM1a-NM and APEC17-BIM3a-NM) also had mutations in rcsC and in other genes. This suggests that another mutation in another gene may compensate for the rcsC mutation, canceling the mucoid phenotype. Because the yrfF gene was more frequently mutated in both strains (K-12 and APEC) and is already known to be important for phage defense, the BIM K12-BIM1-M that had only a single mutation in the yrfF gene (L643R) was used to investigate the impact of the mutation in yrfF in more detail.

Despite the mucoidy produced by K12-BIM1-M compared to the parental strain, no difference in the adsorption level was observed. At this stage, the impact of mucoidy on phage 66's ability to recognize its receptor(s) remains unclear.

Numerous studies indicate that deletion of yrfF, which fully activates the Rcs system, is lethal for bacteria [78, 104, 114]. However, another deletion of either rcsB, rcsC, or rcsD can restore viability [115]. This suggests that the observed yrfF point mutation in K12-BIM1-M likely did not disrupt its regulatory role. However, it may still affect its interactions with other proteins, potentially outside the RcsC pathway. While the mutation in yrfF led to increased CA production, the lack of additional mutations in other genes of the Rcs phosphorelay system suggests that YrfF remains largely functional.

The YrfF-WT protein and its L643R mutant were studied using MD simulations to understand their structural and functional roles. The L643R mutation, located at the interface between the periplasmic domain and the membrane, replaced a nonpolar residue with a positively charged one. This substitution disrupted the secondary structure at the base of the periplasmic domain, inducing a transition from a β-sheet to a random coil conformation. Consequently, residues near the mutation exhibited reduced membrane contact and increased solvent exposure, weakening the membrane anchoring of the YrfF-L643R periplasmic domain and making it more upright. This conformational change may affect the function of YrfF, potentially altering its interactions with other proteins, such as RcsF.

Bacterial growth curves for E. coli K-12 WT and K12-BIM1-M were performed prior to RNA extraction. The number of CFU and the OD curves were not significantly different between the two strains, suggesting that the fitness of the parental strain and the mucoid BIM are quite comparable. This contrasts with the expectation that the overproduction of CA would require significant energy, potentially reducing the fitness of the mucoid BIM compared with the parental strain, as observed by other groups [116, 117]. This discrepancy is likely because the mutations in this study were naturally occurring, and that our isolation protocol selected for BIMs having the best fitness.

The differential expression of K12-BIM1-M genes was investigated in relation to the parental strain. Although the yrfF gene was not significantly dysregulated, we demonstrated that the genes producing CA were well enriched and overexpressed in K12-BIM1-M. As mentioned above, the YrfF protein is known to negatively regulate rcsC, thus minimizing CA production. Given the phenotype, one would expect the rcsC, rcsD, rcsB, and rcsA genes to be overexpressed, because an altered YrfF might fail to repress rcsC, thereby activating the Rcs phosphorelay and leading to increased mucoidy. However, these genes were not significantly dysregulated; therefore, the mucoidy observed in the mutant strain might be a result of a novel interaction involving the altered YrfF protein.

In K12-BIM1-M, the genes yjbE, yjbF, and yjbG were significantly up-regulated, with yjbE showing the highest level of up-regulation. This gene is part of an operon regulated by the Rcs phosphorelay system [118]. Interestingly, there was an enrichment in genes overexpressed in other functional categories (histidine metabolism, fimbrium, oxidase complex, and glycogen metabolic process). The categories for under-expressed genes included iron-sulfur complex, UMP biosynthetic process, flagellar assembly, anaerobic electron transport chain, and xenobiotic transport. This demonstrates the complexity of the interaction network that can exist between genes. Finally, the tet(B) gene and its regulator tetR were overexpressed in K12-BIM1-M; however, the reason why these genes were significantly overexpressed is not clear. It is known that tetR is self-regulating and that its expression level is triggered by an inducer, usually tetracycline [119]. Further analysis is required to determine whether other molecules, such as stress factors, could also induce tet genes transcription.

Conclusion

The objective of this study was to examine phage 66’s capacity to infect an APEC strain and the laboratory strain K-12, as well as the mechanisms these two distinct E. coli strains employ to defend against phage 66. The results demonstrated that phage 66 infects its host strain APEC17 more rapidly than K-12, which could be considered an alternative host. The isolation of BIMs from both strains revealed that most of the resistant bacteria exhibited increased mucoidy. Genomic analysis of these mutant strains identified amino acid changes in several proteins involved in CA production, likely responsible for the increased mucoidy. Further investigation of a mutant strain with a mutation in the YrfF protein (L643R) showed that its interaction with the membrane was altered, and that several genes, including those in the Rcs locus involved in CA production, were dysregulated. It would be interesting to investigate other BIMs with different mutations, such as K12-BIM6-M, which also has a single mutation in the YrfF protein but at a different position (D287N). Several questions remain unanswered. It is not clear, at the molecular level, why phage 66 infects APEC17 more effectively than K-12, despite both strains likely having similar protection mechanisms. Additionally, it would be valuable to explore why none of the BIMs had differences in their CRISPR-Cas systems. Lastly, the involvement of tetracycline resistance genes, whether lost or overexpressed, remains enigmatic and highlights the complexity and interconnectedness of bacterial molecular systems. This study also underscores that, although phage therapy appears promising, many questions remain and warrant further investigation.